三角形の交差判定

programming 3D c++

3Dプログラミングにおいて三角形（ポリゴン）の交差判定は重要である。このルーチンの品質・速度は全体のパフォーマンスに大きく影響する。
三角形の交差判定でよく使われるのが以下の手法だと思う。

このうち「Arenbergの手法」は行列を用いて、三角形を変形してやることで判定を行なう。
「Arenbergの手法」では逆行列を用いてるが、特別に逆行列が必要なわけではない。
移動して回転してシアーしてスケーリングすれば良い。
以下の型が定義されていたとする。

typedef double real;

class vector3;//

typedef vector3 vertices[4];

enum{
    VERTEX_NORMAL,
    VERTEX_BINORMAL_U,
    VERTEX_BINORMAL_V
};

typedef vector3 ray[2];
enum{
    RAY_ORIGIN,
    RAY_DIRECTION
};

typedef vector3 intersection[2];
enum{
    INTERSECTION_NORMAL,
    INTERSECTION_TUV
};

このうえで、移動して回転してシアーしてスケーリングをうまくまとめまくると、三角形の交差判定は以下のようになる。引数は左から交差情報を格納する配列、頂点座標の配列、光線情報、最遠値である。

bool intersect(intersection is, const vertices vtx,const ray r, real dist){
    static const real EPSILON = 1e-10;//お好みで
    real t, u, v;
    
    vector3 edge1(vtx[1]-vtx[0]); 
    vector3 edge2(vtx[2]-vtx[0]); 
    vector3 r_org(r[RAY_ORIGIN]-vtx[0]);//
    
    vector3 normal(cross(edge1, edge2));
    
    real nd =  dot(r_org, normal);
    real nv = -dot(r[RAY_DIRECTION], normal);
        
    if(nv < EPSILON)return false;
        
    t = nd/nv;
    
    if( t < 0 || dist < t )return false;
    
    vector3 pos(r[RAY_DIRECTION]*t + r_org);
    
    real el = dot(edge1,edge2);
    //real e1 = dot(edge1,edge1);
    //real e2 = dot(edge2,edge2);   
    
    vector3 nx = edge1*dot(edge2,edge2)-edge2*(el);              //cross(edge2,normal);//x
    vector3 ny = edge1*el              -edge2*(dot(edge1,edge1));//cross(edge1,normal);//y
    
    u = dot(pos,nx)/dot(edge1,nx); //
    v = dot(pos,ny)/dot(edge2,ny); // 
    
    if( (EPSILON>u)||(EPSILON>v)||(u+v>1-EPSILON) )return false;
    
    
    is[INTERSECTION_NORMAL] = normal;   
    is[INTERSECTION_TUV]    = vector3(t,u,v);
    
    return true;
    
}

どうだろう？もし内積・外積が速いのなら複雑な処理がない分速くなるだろう。
しかし上記のコードは内積・外積が特別速くない限りは「Tomas Mollerの手法」に速度で勝てない。「教科書の手法」にわずかに勝る程度だ。
ただし、上記のコードはおもしろい特徴がある。それは途中に計算する頂点値から得られる４つのベクトルが事前に計算しておくことが可能であるということだ。つまり事前に頂点座標からベクトルを計算しておき、呼びだし側で頂点座標の代わりに計算したベクトルを関数に与えてやることで高速化が期待できる。

vector3 vtx[4]={...}

//省略

   vector3 edge1(vtx[1]-vtx[0]); 
   vector3 edge2(vtx[2]-vtx[0]); 
vector3 normal(cross(edge1,edge2));
   vector3 nx(cross(edge2,normal));//x
   vector3 ny(cross(edge1,normal));//y
vector3 binormal_u(nx/dot(edge1,nx));
vector3 binormal_v(ny/dot(edge2,ny));

//vtx[0]=vtx[0];
vtx[1]=normal;//
vtx[2]=binormal_u;
vtx[3]=binormal_v;

//省略

if( intersect_ex(insect,vtx,r,100000000.0) ){...}//呼び出し

とすると交差判定はいかのコードになる。

bool intersect_ex(intersection is, const vertices vtx, const ray r, real dist){
    static const real EPSILON = 1e-10;
    
    real t, u, v;
    
    vector3 r_org(r[RAY_ORIGIN]-vtx[0]);
    
    real nd =  dot(r_org,            vtx[VERTEX_NORMAL]);
    real nv = -dot(r[RAY_DIRECTION], vtx[VERTEX_NORMAL]);
        
    if( nv < EPSILON )return false;
        
    t = nd/nv;
    
    if( t < 0 || dist < t )return false;
    
    vector3 pos(r[RAY_DIRECTION]*t + r_org);
    
    u = dot(pos,vtx[VERTEX_BINORMAL_U]);
    v = dot(pos,vtx[VERTEX_BINORMAL_V]);
    
    if( (EPSILON > u)||(EPSILON > v)||(u+v > 1-EPSILON) )return false;
    
    
    is[INTERSECTION_NORMAL] = vtx[VERTEX_NORMAL]);   
    is[INTERSECTION_TUV]    = vector3(t,u,v);
    
    return true;
    
}

これならベクトルに関して4つの内積といくつかの加減算ですむ。
速度的には「Tomas Mollerの手法」の1.3倍くらいの速度となる。たぶんこれより速い三角形交差ルーチンはないと思う。
弱点はポリゴンごとにデータを作っておかなければならないため、インデックス頂点式に比べ、メモリ上のデータがかさんでしまうことだろうか。
精度に関してだが厳密に調査したわけではないのでわからない（EPSILONと前計算の精度次第）。
上記のコードを用いて、ポリゴンを描写したものを以下に示す。

以下は比較用。多少明度をあげて見やすくした。



実はもともとは「Arenbergの手法」の詳細を知らずに求めていったのだが、結局「Arenbergの手法」とほぼ同一のものになってしまった。
「Arenbergの手法」において精度が気になるのなら、必要なベクトルを求める際、逆行列を使わず、上記のコードを用いて逐次的に求めていったほうがよいだろう。その際できるだけ高精度で計算すればよい。
また既存のコードへの変更はメッシュのデータ構造に対してのみなので、レイトレースシステム全体に影響する光線に手を入れるPlucker coordinateの手法を使うぐらいだったら、こちらを選んだ方が良いと思われる。

追記：
いや、そもそも「Arenbergの手法」の手法と呼んでいいのだろうか？わからない。
・・・ここなんか見ても理論としてはそのまんまですし・・・
どうやら「Barycentric Coordinatesを用いた手法」というのがより一般的な呼称の気がする。
Barycentric Coordinatesとは日本語で「重心座標系」のことだ。要するに重心比をわかりやすいよう変形した座標系のこと。つうか上のコード中にで出したu,vっつうのは重心座標系の値そのものなわけで、それを交差判定に直接利用すれば一石二鳥なわけで。
くわしくはこの記事の後半
あらかじめ頂点ベクトルから別のベクトルを求めるという作業は、重心座標を求めていたに過ぎないわけだ。
どうやらこの重心座標を求める計算は以下のようにしたほうが簡単に求められるようだ。

vector3 e1 = cross(vtx[2],vtx[0]);
vector3 e2 = cross(vtx[0],vtx[1]);
vector3 normal( cross( (vtx[1]-vtx[0]),(vtx[2]-vtx[0]) ) );//cross(vtx[1],vtx[2])+e1+e2
real  D = dot(e2,vtx[2]);

vector3 m1 = e1/D;
vector3 m2 = e2/D;

//vtx[0]=vtx[0];
vtx[1]=normal;
vtx[2]=m1;
vtx[3]=m2;

*1:だれが考案したものなのだろう？知っている方教えてください。