OpenCLでいくぜ！ (4)

　ようやくレイトレースの話。
　視線と球の交差判定なんですが、これは視線の式と球の式を連立させた方程式を解くってことなんですな。

　で、視線をどう数式にするかというと
　3次元空間上の座標（x, y, z）と、X軸、Y軸、Z軸に平行なそれぞれの移動量dx, dy, dzによって、視線が表現できるわけですわ。



　（dx, dy, dz）を視線方向ベクトルといいます。
　で、この視線方向ベクトルを、方向はそのままで長さが1になるように調整したもの（単位ベクトルっていいます）を使い、視点である座標（x, y, z）から、視線方向単位ベクトルを何倍すれば球に当たるか考えるわけです。
　もちろん当たらない場合もある。



　これを倍率をtとして、球の表面との交差点（X, Y, Z）との関係式で表すと以下のようになるわけですよ。



　一方、球の表面上の座標（X,Y,Z）を表現する式は、球の半径をr、球の中心座標を（cx, cy, cz）とすると以下のようになるわけで



　視線が球と交差するなら交差点（X, Y, Z）はどちらの式でも共通であるわけだから視線の式を球の式に代入できるんですな。これを実際にやったのが以下の式



　であるなら未知の変数はtだけとなり、tについて方程式を組み直すと以下のような2次方程式となるわけですよ。



　で、高校で習ったとおり（中学だっけ？）2次方程式は解の公式を使い解く事ができるわけですよ。



　ここで判別式が0か正の値ならtは実数として存在し、すなわち視線ベクトルと球は接触することになるわけです。



　この演算を以下の4つのパラメータ

sphere_cente：球の中心座標
sphere_radius：球の半径
eye_position：視点座標
eye_direction：視線単位ベクトル

　を引数で受け取り、交差するならtの値、しないなら-1を返すようにした関数がcalc_tという関数なわけです。
　注意）座標、ベクトルに関してはX,Y,Zの3要素があるんですが、これを別々の変数にはせずX,Y,Z順のfloata型配列としとります。例えば球の中心座標のX座標はsphere_cente[0]でY座標はsphere_cente[1]となる。

static float calc_t(
const float sphere_center[3], // 球の中心座標
const float sphere_radius, // 球の半径
const float eye_position[3], // 視点座標
const float eye_direction[3]) // 視線単位ベクトル
{
float dx = eye_position[0] - sphere_center[0];
float dy = eye_position[1] - sphere_center[1];
float dz = eye_position[2] - sphere_center[2];
float A = 1.0; // eye_directionは単位ベクトルなので1になる
float B = eye_direction[0] * dx
+ eye_direction[1] * dy
+ eye_direction[2] * dz;
float C = dx * dx
+ dy * dy
+ dz * dz
- sphere_radius * sphere_radius;
// 返り値を-1で初期化
float t = -1;
// 解の公式の値が実数をもつか判別
float D = B * B - A * C;
if (D >= 0) {
// 2つの解を計算
float t0 = (-B - sqrt(D)) / A;
float t1 = (-B + sqrt(D)) / A;
// 小さい方を選択
if (t0 > t1) {
t0 = t1;
}
// 視点、または視点より前にあればtを設定する
if (t0 > 0) {
t = t0;
}
}
return t;
}

　まずは判定式で実数解をもつ（接触する）事を確認し、接触するなら解の公式でtを求めればよいわけです。



　上記の式からわかるように解tは2つあるんですが、値の小さい方が0より大きいなら球は、視点の視線方向先にあり、tはその交差点までの視線方向単位ベクトルの倍率を示す事になるわけです。



　ちなみに値の大きい方のtは、球の後ろ側の交差点までの大きさを示す。



　２つのtの値が一致する場合は、視線は球の接線上となるわけです。



　また、小さい方のtの値が0より小さいなら、球は視線より後ろ、または視点は球の中となるし



　視点が球の内部かそうでないかは、大きい方のtの値が負の値かどうかで判断できるわけです。



　今回の場合は、視点より前の球だけ描画対象とするとし、小さい方のtの値が0以下の場合はすべて-1になるようにしました。

　rendering関数では視点はZ軸上、スクリーンはXY平面に重ねてZ軸が中心を通るようしています。
　この状態でスクリーンの全画素についてcalc_t関数呼び出しをおこない、返り値が-1以外のところは球であるとし、その画素は赤になるようにすればレイトレーシングによる描画が完成するってわけです。


視点と、スクリーンのどの画素をチェックするかで視線方向が決定するわけだ

static void rendering(CGContextRef context)
{
int width = CGBitmapContextGetWidth(context);
int height = CGBitmapContextGetHeight(context);
unsigned char* baseAddress = (unsigned char*) CGBitmapContextGetData(context);
size_t bytesPerRow = CGBitmapContextGetBytesPerRow(context);
float eye_position[] = {0.0, 0.0, -width};
float sphere_radius = width / 2;
float sphere_center[3] = {0.0, 0.0, sphere_radius * 1.1};
for (int v = 0; v < height; v++) {
for (int h = 0; h < width; h++) {
unsigned char* p = baseAddress + (bytesPerRow * v) + h * 4;
float x = h - (width / 2);
float y = v - (height / 2);
renderingPixel(p, x, y, 0.0, eye_position, sphere_center, sphere_radius);
}
}
}

　renderingPixel関数では視点とスクリーン画素の座標を元に視線ベクトル（単位ベクトル）を計算してcalc_t関数を呼び出し画素の色を設定しています。

static void renderingPixel(
unsigned char* p, // 画素のバッファ
float x, // 画素の3次元空間上のx座標
float y, // 画素の3次元空間上のy座標
float z, // 画素の3次元空間上のz座標
const float eye_position[3], // 視点座標
const float sphere_center[3], // 球中心座標
const float sphere_radius) // 球の半径
{
// 視線方向ベクトルの決定
float eye_direction[3]; // 視線方向ベクトル
eye_direction[0] = x - eye_position[0];
eye_direction[1] = y - eye_position[1];
eye_direction[2] = z - eye_position[2];
// 視線方向ベクトルを単位ベクトルにする。
float l = eye_direction[0] * eye_direction[0]
+ eye_direction[1] * eye_direction[1]
+ eye_direction[2] * eye_direction[2];
l = sqrt(l);
eye_direction[0] /= l;
eye_direction[1] /= l;
eye_direction[2] /= l;

float t = calc_t(sphere_center, sphere_radius, eye_position, eye_direction);
if (t >= 0) {
// 交差するので赤色にする
*p = 255;
}
}

　calc_t関数に渡すための視線方向の単位ベクトル（単位ベクトルは長さが1のベクトル）は以下のように計算。
　視点からスクリーン画素の座標へ向かう方向を視線ベクトルとし

// 視線方向ベクトルの決定
float eye_direction[3]; // 視線方向ベクトル
eye_direction[0] = x - eye_position[0];
eye_direction[1] = y - eye_position[1];
eye_direction[2] = z - eye_position[2];

　これを単位ベクトルに変換するためにベクトルの大きさでX,Y,Z成分をベクトルの大きさで割る。
　ベクトルの大きさは三平方の定理を利用して求める。

// 視線方向ベクトルを単位ベクトルにする。
float l = eye_direction[0] * eye_direction[0]
+ eye_direction[1] * eye_direction[1]
+ eye_direction[2] * eye_direction[2];
l = sqrt(l);
eye_direction[0] /= l;
eye_direction[1] /= l;
eye_direction[2] /= l;

　これで描画される画面はこんな感じ。



　残念感でいっぱいです。

　ま

if (t >= 0) {
// 交差するので赤色にする
*p = 255;
}
}

　なんてふうに赤一色で描いてるんだからあたりまえ。
　この赤色の濃淡を球体の陰影にあわせて決定してやることで、もうちょっとリアルな球体になる。

if (t >= 0) {

// 法線ベクトルを計算
float normal[3]; // 法線ベクトル
normal[0] = eye_direction[0] * t + eye_position[0] - sphere_center[0];
normal[1] = eye_direction[1] * t + eye_position[1] - sphere_center[1];
normal[2] = eye_direction[2] * t + eye_position[2] - sphere_center[2];
// 法線ベクトルを単位ベクトルにする。
normal[0] /= sphere_radius;
normal[1] /= sphere_radius;
normal[2] /= sphere_radius;
// 光の方向ベクトルを設定
// 光の方向ベクトルを設定
float light[] = {0.0, -1.0, 0.0};
float l = sqrt(light[0] * light[0] + light[1] * light[1] + light[2] * light[2]);
light[0] /= l;
light[1] /= l;
light[2] /= l;
// 球表面の明るさを決定
float bright = -normal[0] * light[0] + -normal[1] * light[1] + -normal[2] * light[2];
if (bright < 0)
bright = 0;
else if (bright > 1.0)
bright = 1.0;
*p = 255 * bright;
}
}

　これで以下のような結果が得られる。


　人類、火星に立つ！

　というか、球体の陰影はついてんだけど宇宙空間に置かれた球体のようになっているわけだ…
　ここらへんの対策や、なぜ上の計算で球体の陰影が計算できるのか等は、まて次回！

　それにしても、この1000x1000のレイトレース画像を一瞬で演算しちゃうCore i7って…
　その昔、Z80 4MHzのマシンでHu-BASICを使い数時間かけて…
　しつこいですか、そうですか。

------------
サンプルプロジェクト：opencl02.zip