classic fastest method --> Tomas Moller's method

以前にもポリゴンとレイとの交差について考えた*1。いろいろ試してみたけど、最近ではいわゆる[classic fastest method]と呼ばれている手法を使っている。classicというから、歴史は古いのだろう。この方法は数学的にとてもすっきり導出できる。インデックス式だし、重心座標も求めることが出来るし,~~Tomas Mollerの方法より速い結局結局余計なことはするなということなんだろう。~~。これを推し進めると結局Tomas Mollerのコードになる。（はじめから気づけよ！＞俺。多分速度が以前おそかったのは両面判定にする段階であまりよろしくないことをしてたのかもしれないなあ。とりあえずこの記事はTomas Mollerの導出として残しておきます。thx>骸骨さん）

まず以下の図のような三角形を考える。それぞれの頂点をp0,p1,p2とする。またp0p1をe1,p0p2をe2とする。

このときの三角形上の任意の点をPとするとPは以下の式で表せる。

（式１）
ただし

一方、光線を光線の初期位置ベクトルorgと方向ベクトルdirで表すと、光線上の点は距離（スカラー値）tを用いて以下のように表せる。
（式２）

以上より（式１）と（式２）の右辺を等号すると方程式が出来る。
（式３）
この式を少し整理してやると以下のようになる。
（式４）

ここで（式４）は以下のような形で表されている。
（式５）
（式５）のようにある（３次）ベクトルが他の３つの（３次）ベクトルのスカラー倍の和で表されているとき、その３つのスカラー値はクラメールの公式を用いることによって解くことが出来る。以下に公式を表す。ただしdet(v1,v2,v3)は３つのベクトルを縦ベクトルとしその並びを行列として行列式を解くものとする。
（式６）

３角形の交差判定はクラメールの公式を式４に適用しu,v,tを求め、その値が三角形の交差条件となる値か調べることで実現できる。
その際

式６はそれぞれの式の形がよく似ていて、また右辺の分子と分母もよく似てていることが分かる。これをうまく利用して無駄な計算は省く。
割り算は非常に重たい計算なので出来るだけ避ける。

以上の２点を抑えることで実際のコードが高速になる。
以下はそれを踏まえてC++を用いて書いたコード。

class vector3;//ベクタークラス
class ray;//光線クラス
//上記クラスは適当に実装すべし

//ヒットしたときにはその情報を格納
struct hit_info{
	float distance;
	
	vector3 position;//交点
	vector3 normal;//法線
};

//-----
//３角形交差判定ルーチン
//true :三角形にヒットした
//false:三角形にヒットしなかった
//-----
bool test_polygon(
	const vector3 & p0, const vector3& p1, const vector3& p2,	//３つの頂点
	hit_info * info,//出力・・・ヒットしたときはその情報を格納
	const ray& r, float distance,//r:光線,
    //distance:それまでに最も近いポリゴン面までの距離
    //この値よりも大きな距離で交差したら無視する。
){
	static const float EPSILON = value::epsilon()*1000;
	float u, v ,t;
	
	//vA
	//-e1 = p0-p1
	float A = p0.x - p1.x;
	float B = p0.y - p1.y;
	float C = p0.z - p1.z;
	
	//vB
	//-e2 = p0-p2 
	float D = p0.x - p2.x;
	float E = p0.y - p2.y;
	float F = p0.z - p2.z;
	
	//vC
	//dir
	float G = r.direction().x;
	float H = r.direction().y;
	float I = r.direction().z;
	
	//vD
	//p0-org
	float J = (*p_p0).x - r.origin().x;
	float K = (*p_p0).y - r.origin().y;
	float L = (*p_p0).z - r.origin().z;
	
	float P = E*I-H*F;
	float Q = G*F-D*I;
	float R = D*H-E*G;
	
	float iM = (A*P + B*Q + C*R);//分母

	float S,T,U;
	if(EPSILON < iM){//マイナスだと表面
		u = (J*P + K*Q + L*R);
		if(u < 0 || iM < u)return false;
		
		//--------------

		S = B*L-K*C;//BLKC
		T = J*C-A*L;//JCAL
		U = A*K-J*B;//AKJB

		//--------------
        
		v = (G*S + H*T + I*U);
		if(v < 0 || iM  < u+v)return false;

		t = -(D*S + E*T + F*U);
		if(t < 0)return false;
	}else if(iM < -EPSILON){//裏面　裏面だとヒットしないのなら（片面のみなら）ここでfalseを返せばよい。
		u = (J*P + K*Q + L*R);
		if(u > 0 || iM > u)return false;
		
		//--------------

		S = B*L-K*C;//BLKC
		T = J*C-A*L;//JCAL
		U = A*K-J*B;//AKJB

		//--------------
        
		v = (G*S + H*T + I*U);
		if(v > 0 || iM  > u+v)return false;

		t = -(D*S + E*T + F*U);
		if(t < 0)return false;
	}else{
		//分母が0のとき、これは
		return false;
	}

	iM = float(1.0)/iM;//ここでようやく割り算

	t *= iM;
	if(distance < t)return false;
	
	info->distance = t;
	info->position = r.origin()+t*r.direction();//場所
	if(iM<0){//表面
		info->normal = normalize(cross(vecto3(-A,-B,-C),vector3(-D,-E,-F)));//
	}else{	//裏面なら　裏面だとヒットしないのなら（片面のみなら）このルーチンいらない
		info->normal = normalize(cross(vecto3( A, B, C),vector3( D, E, F)));//
	}


	return true;
}

さらにdot,crossの演算でまとめると

bool test_polygon(
	const vector3 & p0, const vector3& p1, const vector3& p2,	//３つの頂点
	hit_info * info,//出力・・・ヒットしたときはその情報を格納
	const ray& r, float distance,//r:光線,
    //distance:それまでに最も近いポリゴン面までの距離
    //この値よりも大きな距離で交差したら無視する。
){
	static const float EPSILON = value::epsilon()*1000;
	float u, v ,t;
	
	//vA
	//-e1 = p0-p1
	vector3 vA(p0-p1);
	
	//vB
	//-e2 = p0-p2 
	vector3 vB(p0-p2);
	
	//vC
	//dir
	const vector3& vC = r.direction();
		
	//vD
	vector3 vD(cross(vC,vB));	
	
	float iM = dot(vA,vD);//分母
	
	if(EPSILON < iM){//マイナスだと表面	
		vector3 vE(p0-r.origin());
		
		u = dot(vD,vE);
		if(u < 0 || iM < u)return false;
		
		vector3 vF(cross(vA,vE));
		v = dot(vC,vF);
		if(v < 0 || iM  < u+v)return false;

		t = -dot(vB,vF);
		if(t < 0)return false;
	}else if(iM < -EPSILON){//裏面　裏面だとヒットしないのなら（片面のみなら）ここでfalseを返せばよい。
		vector3 vE(p0-r.origin());
		
		u = dot(vD,vE);
		if(u > 0 || iM > u)return false;
		
		vector3 vF(cross(vA,vE));
		v = dot(vC,vF);
		if(v > 0 || iM  > u+v)return false;

		t = -dot(vB,vF);
		if(t < 0)return false;
	}else{
		//分母が0のとき、これは
		return false;
	}

	iM = float(1.0)/iM;//ここでようやく割り算

	t *= iM;
	if(distance < t)return false;
	
	info->distance = t;
	info->position = r.origin()+t*r.direction();//場所
	if(iM<0){//表面
		info->normal = -normalize(cross(vA,vB));//
	}else{	//裏面なら　裏面だとヒットしないのなら（片面のみなら）このルーチンいらない
		info->normal =  normalize(cross(vA,vB));//
	}

	return true;
}

Tomas Mollerの手法とほぼ同様になる。

ポリゴンメッシュへの交差判定は三角形単体のみでいくら高速であろうとも、その上のアクセル構造（UG,Octree,BSP-Tree...）が高速でなければ、意味がない。１つの三角形交差判定のスピードの違いというのは、アクセル構造トラバースのスピードからすれば誤差のようなものだろう。でもやっぱりこだわりたいところではある。

参考図書：