!!!交差判定 レイトレース処理での一番大事な部分、交差判定について記述します。 ここではレイ（視点位置と視線ベクトルを持つ）とポリゴン（三角形）との交差判定になりますね。 !!単純な交差判定 三次元空間上のポリゴンをX/Y/Z軸を圧縮する形で２次元に投影してしまいます。 これは、X-Y平面への投影・X-Z平面への投影・Z-Y平面への投影の３つの投影があります。一番確実なのは（誤差を少なくするのは）それぞれの面に投影した場合の面積を計算して、一番面積の大きい面に投影するとするといいです。 {{ref_image touei_20040821.png}} 上図の場合は、X-Z平面に投影しています（三角形の頂点座標のうち、X/Z成分のみを取り出します）。 また、レイの方向ベクトルと面の法線ベクトルにより「直線と面の交点位置」を求めます（これは３次元空間での処理）。このときに交点が求まります。が、三角形内に内包されているかは分かりません。これを、X-Z平面に三角形を投影している場合は交点座標のX/Z成分のみを取り出すとよいです。 さて、後は２次元空間の処理になります。 {{ref_image inner_poly_20040821.png}} 空間としては、横方向プラスは右向き・縦方向プラスは上向きで三角形の頂点の並びが逆時計周り、としていますが下記方法だと特に意識する必要はありません。 ポリゴン頂点座標を0-1-2、交点位置をPとします。 P-0、0-1の順にベクトルを作成（上図の灰色の矢印）してこれの外積を求めます。 P-0のベクトル = (ax, ay) 0-1のベクトル = (bx, by) 外積 = ax * by - ay * bx この外積が>0の場合はカウンタを+1します。それ以外は-1します。 同様に、P-1/1-2のベクトル・P-2/2-0のベクトルも外積を求めてカウンタを操作します。 この結果、カウンタが-3または+3の場合では三角形内に交点Pが内包されている、と結論を下すことができます。 !!Tomas Mollerの交差判定 「Practical Analysis of Optimized Ray-Triangle Intersection 」の論文の交差判定です。この場合は３次元のポリゴンを２次元に投影する必要はありません。 http://www.ce.chalmers.se/old/staff/tomasm/raytri/ また、研究されつくしているようないわゆる「枯れた」手法ですのでよく使われます。 これは、レイの始点位置・レイの方向ベクトル（単位ベクトル済み）と三角形の３頂点が与えられたときの視点からの距離Tと交点情報U/Vを求めます。 Origがレイの開始位置、dirがレイの方向ベクトル（単位ベクトルに正規化済み）、 v0/v1/v2が三角形の頂点座標とし、結果をpRetT/pRetU/pRetVで示されるポインタに 返す関数は以下のようになります。 typedef struct { float x, y, z; } VECTOR3; bool TriangleIntersect(VECTOR3 Orig, VECTOR3 dir, VECTOR3 v0, VECTOR3 v1, VECTOR3 v2, float *pRetT, float *pRetU, float *pRetV) { VECTOR3 e1, e2, pvec, tvec, qvec; float det; float t, u, v; float inv_det; Vec3Sub(&e1, &v1, &v0); vec3Sub(&e2, &v2, &v0); Vec3Cross(&pvec, &dir, &e2); det = Vec3Dot(&e1, &pvec); if (det > (1e-3)) { Vec3Sub(&tvec, &Orig, &v0); u = Vec3Dot(&tvec, &pvec); if (u < 0.0f || u > det) return false; Vec3Cross(&qvec, &tvec, &e1); v = Vec3Dot(&dir, &qvec); if (v < 0.0 || u + v > det) return false; } else if (det < -(1e-3)) { Vec3Sub(&tvec, &Orig, &v0); u = Vec3Dot(&tvec, &pvec); if (u > 0.0 || u < det) return false; Vec3Cross(&qvec, &tvec, &e1); v = Vec3Dot(&dir, &qvec); if (v > 0.0 || u + v < det) return false; } else { return false; } inv_det = 1.0f / det; t = Vec3Dot(&e2, &qvec); t *= inv_det; u *= inv_det; v *= inv_det; if(pRetT) *pRetT = t; if(pRetU) *pRetU = u; if(pRetV) *pRetV = v; return true; //hit!! } ちなみに、各関数は3次元のベクトル操作関数です。 [[ベクトル操作関数|よく使うベクトル・行列計算関数_3DCG]]は 自前で用意しておくのがいいかと思います（よく使います）。 Vec3Sub(&c, &a, &b); ... c = a - b を計算 e = Vec3Dot(&a, &b); ... aとbの内積を求め、結果はe（float） Vec3Cross(&c, &a, &b); ... aとbの外積を求め、結果はc（VECTOR3） 結果、始点Origからの距離Tで交わります。また、交点情報U/Vが求まります。 U/Vは、法線やテクスチャのUV（同じUVでややこしいですが別物です）を求める際に使用します。 交点座標は以下の感じで求まります。これは、レイの方向ベクトルが単位ベクトル（長さ1.0）なんで単純ですね。 px = dir.x * T + Orig.x; py = dir.y * T + Orig.y; pz = dir.z * T + Orig.z; 法線を求める場合は、頂点ごとに「頂点法線(VECTOR3 n1/n2/n3とする)」が定義されているとして以下の感じです。 float fDat = 1.0f - U - V; n1 = n1 * fDat; n2 = n2 * U; n3 = n3 * V; n = n1 + n2 + n3; ここでのnが交点位置での法線ベクトルになります（正規化されたものです）。 !!Arenbergの交差判定 http://www.acm.org/tog/resources/RTNews/html/rtnews5b.html#art3 上記の方法よりも、より計算処理を省略することができます。 ただし、精度としては少し荒く感じます。 この手法では三角形をレイに交差させる前に、三角形自体を(0, 0, 0) - (0, 1, 0) - (1, 0, 0)の形に変換し、同時にレイの位置・方向ベクトルも変換（アフィン変換）します。 {{ref_image Arenberg_20040822.png}} 計算処理が少なくて済むというのが一番のメリットです。 たとえば、レイの開始位置を(ox, oy, oz)、レイの方向ベクトルを(dx, dy, dz)とし、 三角形を(x0, y0, z0)-(x1, y1, z1)-(x2, y2, z2)とします。 これより、まずポリゴンの法線ベクトルを求めて正規化します。 float ax, ay, az, bx, by, bz, cx, cy, cz; float nx, ny, nz; double dnx, dny, dnz, da; ax = x0 - x2; ay = y0 - y2; az = z0 - z2; bx = x1 - x2; by = y1 - y2; bz = z1 - z2; dnx = ay * bz - az * by; dny = az * bx - ax * bz; dnz = ax * by - ay * bx; da = sqrt(dnx * dnx + dny * dny + dnz * dnz); if(da <= 1e-6) return 0; da = 1.0f / da; dnx *= da; dny *= da; dnz *= da; nx = (float)dnx; ny = (float)dny; nz = (float)dnz; 次に、レイを三角形(0, 0, 0) - (0, 1, 0) - (1, 0, 0)の空間に変換するための 変換行列を求めます。 [(x0 - x2) (x1 - x2) (nx - x2)] Q = [(y0 - y2) (y1 - y2) (ny - y2)] [(z0 - z2) (z1 - z2) (nz - z2)] の「逆行列(Q'とします）」が求める変換行列になります。 MATRIX mat1, mat2; MatrixIdentity(&mat1); //単位行列にする mat1.m[0][0] = x0 - x2; mat1.m[0][1] = x1 - x2; mat1.m[0][2] = nx - x2; mat1.m[1][0] = y0 - y2; mat1.m[1][1] = y1 - y2; mat1.m[1][2] = ny - y2; mat1.m[2][0] = z0 - z2; mat1.m[2][1] = z1 - z2; mat1.m[2][2] = nz - z2; MatrixInverse(&mat2, &mat1); //逆行列を求める ax = mat2.m[0][0]; bx = mat2.m[0][1]; nx = mat2.m[0][2]; ay = mat2.m[1][0]; by = mat2.m[1][1]; ny = mat2.m[1][2]; az = mat2.m[2][0]; bz = mat2.m[2][1]; nz = mat2.m[2][2]; cx = x2; cy = y2; cz = z2; 以上が前処理です。 ここで求まった(ax, ay, az)、(bx, by, bz)、(nx, ny, nz)、(cx, cy, cz) が処理に必要な情報となります（高速に処理する場合は上記のfloat * 3 * 4をバッファに蓄えておく）。三角形の頂点座標は不要ですので破棄してOKです。 さて、それでは交差判定を行います。 レイ開始位置(ox, oy, oz)、レイ方向ベクトル(dx, dy, dz)が渡されたときに このレイ情報をアフィン変換する必要があります。 P' = (Q') * (P - C) で座標位置PをP'に変換します。 Q'は前処理で計算した行列、Cは同じく前処理で計算した(cx, cy, cz)です。 fx = ox - cx; fy = oy - cy; fz = oz - cz; ox2 = ax * fx + bx * fy + nx * fz; oy2 = ay * fx + by * fy + ny * fz; oz2 = az * fx + bz * fy + nz * fz; この(ox2, oy2, oz2)がレイの開始位置の変換後の座標になります。 ベクトルは以下の式で変換します。 V' = (Q') * V Vが方向ベクトルになりV'に変換します。 dx2 = ax * dx + bx * dy + nx * dz; dy2 = ay * dx + by * dy + ny * dz; dz2 = az * dx + bz * dy + nz * dz; この(dx2, dy2, dz2)がレイの方向ベクトルになります。 注意点としてここで、「単位ベクトルに正規化しないでください」。 これは、すでに三角形(0, 0, 0) - (0, 1, 0) - (1, 0, 0)の空間にレイを変換しているからです。この状態で変換前の距離Tを算出するのでそのままにしておく必要があるわけです。 また、方向ベクトル(dx2, dy2, dz2)のdx2が0.0の場合は 三角形と交わらないことになりますので「交差なし」でスキップできます。 if(dz2 >= -(1e-6) && dz2 <= (1e-6)) return 0; 距離Tを求めます。 t = -oz2 / dz2; if(t <= (1e-5)) return 0; 交点情報のU/Vを求めます。 u = ox2 + t * dx2; v = oy2 + t * dy2; if(u < 0.0f || v < 0.0f || (u + v) > 1.0f) return 0; t/u/vは、「Tomas Mollerの交差判定」のときと同じ結果です。 tがレイの始点からの距離でu/vが交点情報になります。 同様にU/Vより交点位置と法線ベクトルなどを計算することができます。 逆行列を求めてから後はアフィン変換とT/U/Vの算出だけ、ということで 非常に低コストですね。 ただし、精度は多少粗が見えます（微妙に一定周期のノイズが乗っています）。 しかしリアルタイムでレイトレを行うなどの場合には大変有効なアルゴリズムであると 思います（オフラインレンダラでは品質優先の場合は厳しいかな）。 ちなみにこの方法を使ったリアルタイムレイトレとして「SaarCOR http://www.saarcor.de/」があります。これはZバッファベースのOpenGL/DirectX の手法ではなく、レイトレーシングを行う方法でハードウェア（FPGA）にて 交差判定を含むレイシューティングを行っています。 サンプルを見る限りは比較的大きなポリゴンが多いですので違和感がないですが、 細かいオフラインレンダラのシーン（極小ポリゴンが密集しているシーンなど）は 品質として厳しくなるかもしれないですね。 （後、仮数部16bitの24ビット浮動小数点を使っているとすると結構厳しいかなと。 最低でも32bit精度がないと私の場合は品質として見るに耐えないと感じました） と話題がそれましたが、私としては「Tomas Moller」の交差判定手法をお勧めします。