レイと三角形の衝突判定

ray-triangle.png

1997年に Möller-Trumboreさんという人が「Fast, Minimum Storage Ray/Triangle Intersection」という論文を発表しました[1]。
その論文にはレイと三角形の衝突判定の効率的な方法が示されています。
そしてその手法が現在でも使用されています。
どんな方法かというと
レイと三角形の交点が、レイの方程式と、三角形の重心座標の両方で表されることに着目して、
2つの方程式を解くことにより、交点を求める、という。。ざっくり説明するとそんなかんじです。
これから詳しくみていきます。
レイと三角形の交点をPとします。
三角形の頂点はA,B,Cとします。
Pは三角形の重心座標で表現できます。

(1)
\begin{equation} P = wA+uB+vC \end{equation}

この(w,u,v)が重心座標です。
u+v+w=1です。
なので、w=1-u-vです。

(2)
\begin{equation} P = (1-u-v)A+uB+vC \end{equation}

この式を変形すると

(3)
\begin{equation} P=A-uA-vA+uB+vC=A+u(B-A)+v(C-A) \end{equation}

(B-A)と(C-A)は三角形ABCを構成する辺のベクトルです。
交点Pはレイの方程式でも表現できます。

(4)
\begin{equation} P=O+tD \end{equation}

Oはレイの原点
tはPがどのぐらい原点Oから離れているかを示すパラメータ
Dはレイの方向ベクトルです。
重心座標で示されたPの式と、レイの方程式で表された‘Pの式を合わせてみます

(5)
\begin{align} O+tD = A+u(B-A)+v(C-A)\\ O-A = -tD+u(B-A)+v(C-A) \end{align}

未知の定数は'(t,u,v)で右側にまとめました。
これをベクトルと行列形式で表現すると

(6)
\begin{align} \begin{bmatrix} -D&(B-A)&(C-A)\\ \end{bmatrix} \begin{bmatrix} t\\ u\\ v \end{bmatrix} = O-A \end{align}

既知の項(O-A),B-A,C-A,Dはベクトルです。
一方、未知の項 t,u,vはスカラー値です。

もし三角形に回転などの変形を施したとしたら、
交点Pはデカルト座標系(x,y,z)としての値は変わりますが、
三角形ABCの重心座標としてのPの値は、どんなに三角形に変形を施したとしても変わりません。
Mollerさんのアルゴリズムはこの利点を活かしています。

O-Aという項は原点がOからAへ移動させるという意味にもとれます。
つまり三角形の頂点Aが重心座標の原点です。

クラメルの公式を利用して解く 

この式はクラメルの公式により、こんな風に解けます

(7)
\begin{bmatrix} t\\ u\\ v \end{bmatrix} = \frac{1}{| -D B-A C-A|} \begin{bmatrix} | (O-A) (B-A) (C-A)|\\ |-D (O-A) (C-A)|\\ |-D (B-A) (O-A)| \end{bmatrix}

行列式の値を求める

1x3の行列式の求め方は

(8)
\begin{align} |ABC| = -(A\times C)\cdot B = -(C \times B)\cdot A \end{align}

3x3行列の行列式の値は三重積です。
なので

(9)
\begin{bmatrix} t\\ u\\ v \end{bmatrix} = \frac{1}{|D \times (B-A)) \cdot (C-A)|} \begin{bmatrix} ((O-A) \times (B-A)) \cdot (C-A)\\ (D \times (O-A)) \cdot (C-A)\\ ((O-A) \times (B-A)) \cdot D \end{bmatrix}

なのでこういう風になります

(10)
\begin{align} t=((O-A) \times (B-A)) \cdot (C-A)/d \newline u=(D \times (O-A)) \cdot (C-A)/d \newline v = ((O-A) \times (B-A)) \cdot D/d \end{align}

$d=(D \times (C-A)) \cdot (B-A)$

ソースコード

vec3 hit(vec3 orig,vec3 dir,vec3 a,vec3 b,vec3 c){
   float eps= 0.0000001;
   vec3 ab=b-a;
   vec3 ac=c-a;
 
   vec3 n=cross(dir,ac);
 
   float det=dot(ab,n);
 
   if(det<=eps){ return vec3(0.3,0.3,0.3);}
 
   vec3 ao=orig-a;
   float u=dot(ao,n)/det;
   if(u<0.0 || u>1.0){ return vec3(0.3,0.3,0.3);}
 
   vec3 e=cross(ao,ab);
   float v=dot(dir,e)/det;
    if(v<0.0||u+v>1.0){ return vec3(0.3,0.3,0.3);}
 
    float t= dot(ac,e)/det;
    return vec3(u,v,t);
}

レイと三角形の衝突判定を利用したラスター化デモ


フラグメントシェーダでレイと三角形の衝突判定をして、
当たったら描画、することにより三角形を描画しました。
フルソースコード(Shadertoy)


raytrace-rasterization raytracing

サポートサイト Wikidot.com raytrace-rasterizationraytracing