glmとclapackでレイ・トライアングル交錯計算

やあ子どもたち元気にしてたかな。
今日はレイ・トライアングル交錯についての考察だよ。レイトレーシングとかピッキングなどで利用頻度の高い計算の一つだ。今更な話題だけど何度も作り直したりするロジックの一つではないかな。
今回おじさんはglm::lookAt()と、CLAPACKを利用した、barycentric coordinates ベースのアプローチを考えたのでメモしておくぞ。考えたというか実装しただけだけどな。手軽に使えるものは何でも流用しようというわけだね。
GLMとかCLAPACKとかに頼ったので、おかげでストラテジーやコードそのものは至極シンプルなものとなった。それは以下の図（図というか式展開）のようだよ。

そう、レイが三角に交錯したときの状況を、t, λ[3]という4つの未定定数を用いて書き下して作った行列方程式を解くだけだ。よく平面の方程式を考えたりとか、三角形を2ベクトルの線形和で考えたりとか、そういう面倒くさいのは一切なしだ。（いないとは思うがbarycentric coordinatesって何？とかここまで来て今更なこと言ってる人は帰ってくれ。）
さて、これに加えて、今回はオプションサービスがつきます！それは何かというと間便な事前プレフィルターの実装方法も提案しています。それはつまり、glm::lookAt()を使っているという点。そう、つまりレイがちょうどZ軸にくるように全体を座標変換したとき、X軸値の符号が３頂点とも同じ三角、あるいはY軸値の符号が３頂点とも同じ三角については、交錯計算をせずに終了するというギミックを入れている点だよ。
コードは以下のようになった。

int ray_triangle_intersection(
                                 const nurs::ray3& ray,
                                 const glm::vec3& a,
                                 const glm::vec3& b,
                                 const glm::vec3& c,
                                 glm::vec3* pintersection_point
                                 )
{
    // returns
    // 0: hit
    // 1: no hit
    // 2: no hit because of parallel situation
    int result=0;
    const glm::mat4 mat = glm::lookAt( ray.first, ray.first+ray.second, glm::vec3(0,1,0) );
    {
        const glm::vec4 o = mat * glm::vec4( a, 1.0 );
        const glm::vec4 p = mat * glm::vec4( b, 1.0 );
        const glm::vec4 q = mat * glm::vec4( c, 1.0 );
        if(
           (o.x < 0 && p.x< 0 && q.x< 0) ||
           (o.y < 0 && p.y< 0 && q.y < 0) ||
           (o.x > 0 && p.x> 0 && q.x > 0) ||
           (o.y > 0 && p.y> 0 && q.y > 0)
           )
            return 1;
    }
    // 平面との交錯を求めます。
    const glm::vec4 p0 = glm::vec4( ray.first, 1.0 );
    const glm::vec4 v = glm::vec4( ray.second, 1.0 );
    {
        // lapack::dgesv_ 呼び出し準備
        int N=4, NRHS=1, LDA=4, LDB=4;
        const long n=4;
        int IPIV[n], INFO;
        double A[n*n], B[n];
        A[0] = -v.x;
        A[1] = -v.y;
        A[2] = -v.z;
        A[3] = .0;
        A[4] = a.x;
        A[5] = a.y;
        A[6] = a.z;
        A[7] = 1.0;
        A[8] = b.x;
        A[9] = b.y;
        A[10] = b.z;
        A[11] = 1.0;
        A[12] = c.x;
        A[13] = c.y;
        A[14] = c.z;
        A[15] = 1.0;
        B[0] = p0.x;
        B[1] = p0.y;
        B[2] = p0.z;
        B[3] = 1.0;
        {
            // ここが dgesv_ 呼び出し
            ::dgesv_( &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO );
            if( INFO!=0 ){
                return 2;
            }
        }
        const auto& b1 = B[1];
        const auto& b2 = B[2];
        const auto& b3 = B[3];
        if(
           (b1 >= .0 && b1 <= 1.0) &&
           (b2 >= .0 && b2 <= 1.0) &&
           (b3 >= .0 && b3 <= 1.0)
           )
        {
            *pintersection_point = (float)b1*a + (float)b2*b + (float)b3*c;
            result = 0;
        }
        else
        {
            result = 1;
        }
    }
    return result;
}

ちなみに、一本のレイに対して、交錯を計算したい三角が複数ある場合は、上記のglm::lookAt()の計算は一回だけ行えば良いので、このmat4の値だけは外部引数で受け取るようにしておいてもいいかも知れない。おっと、CLAPACKの使い方については、こちらの過去日記の後半部分を参照してくれよ。
ではまたな。チャオ！