3次元での三角形の外心を求めるプログラム

C++ math


やあ子どもたち。今年の冬は寒いけどみんな元気でやっているか。
２次元での三角形の外心を求める計算はよく紹介されているけれども、今日は２次元ではなく3次元空間における三角形の外心を求めるプログラムを作成したので紹介するぞ。
以下ソース冒頭のコメントにも書いたけど3次元空間における三角形の外心を計算する問題は、最終的には3平面の交錯点を求める問題に帰着されるんだね。なので3x3行列の行列方程式になる。とここまでくればもうあとはLAPACKさんで解決だろ。LAPACKさんのー？、はい、dgesv_() 関数だよな。
じゃ今回の関数、CalcCircumcenter()のソースを見ていこうか。

/*
 　・外心の計算ストラテジーの考案
 　　・３次元空間内の三角形の外心を計算する方法は、
 　　　p0-p1(法線n0)と、p1-p2（法線n1）の、
　　　それぞれ垂直二等分面と、
 　　　p0,p1,p2の３点を通る、平面（法線n2）で、
 　　　これら３平面の交点ということで計算してやる。
　　　得られる方程式は、
　　　 p dot n0 = m0 dot n0
 　　　p dot n1 = m1 dot n1
 　　　p dot n2 = p1 dot n2
 　　　という連立方程式を p(x, y, z) について解けばよく。
 即ち
　　　 [ n0x n0y n0z ] [x] = [ n0 dot m0 ]
　　　 [ n1x n1y n1z ] [y] = [ n1 dot m1 ]
 　　　[ n2x n2y n2z ] [z] = [ n2 dot p1 ]
 　　・これを解くために、LAPACK::dgesv_　を使って、解を一発で求める作戦を考えています。
 */
bool CalcCircumcenter( const glm::vec3& p0,
                      const glm::vec3& p1,
                      const glm::vec3& p2, float* radii, glm::vec3* centre )
{
    // 3点の外心とその半径を求めます。
    // 計算が正常に完了した場合は 0を、
    // 3点が直線状に並び外心が計算されない場合は!0を、それぞれ返します。
    const float thres = 1e-12;
    glm::vec3 n0 = p1 - p0;
    glm::vec3 n1 = p2 - p1;
    glm::vec3 m0 = .5f*(p0+p1);
    glm::vec3 m1 = .5f*(p1+p2);
    n0 = glm::normalize( n0 );
    n1 = glm::normalize( n1 );
    glm::vec3 n2 = glm::cross( n0, n1 );
    {
        float norm = glm::length( n2 );
        if( norm*norm < thres ){
            return !0; // 3点が絶望的に直線状に並んでいた場合。
        }
    }
    n2 = glm::normalize( n2 );
    {
        // lapack::dgesv_ 呼び出し準備
        long N=3, NRHS=1, LDA=3, LDB=3;
        const long n=3;
        long IPIV[n], INFO;
        double A[n*n], B[n];
        A[0] = n0.x;
        A[1] = n1.x;
        A[2] = n2.x;
        A[3] = n0.y;
        A[4] = n1.y;
        A[5] = n2.y;
        A[6] = n0.z;
        A[7] = n1.z;
        A[8] = n2.z;
        B[0] = glm::dot( n0, m0 );
        B[1] = glm::dot( n1, m1 );
        B[2] = glm::dot( n2, p1 );
        {
            // ここが dgesv_ 呼び出し
            ::dgesv_( &N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO );
            if( INFO!=0 ){
                return !0; // 計算は失敗しました。
            }
        }
        *centre = glm::vec3( B[0], B[1], B[2] );
    }
    // 三角形の外接円半径の計算
    *radii = glm::distance( p1, *centre );
    return 0;
}

おっといい忘れたけどglmベースだったね。GLMはなければ取ってくるか、あるいはそこは別の数値ライブラリで置き換えてくれ。
さて上記コードはMacOSX なら、AccelerateFramework をインポートした状態でそのまま使えるし（今はvecLibではもうないらしいね。内包してるだけなのかも知れんが）、WindowsでもCLAPACKさえ使える環境なら、そのままコンパイルするはずだよ。ちなみに、後者の場合は、こんな注意事項もあった気がするから気をつけてくれ。もう昔の話になってればいいのだけどな。

わざわざLAPACKを引っ張りだしてみたけど使えるものは使っとけだ。これの計算スピードはね、速いよ。もう一瞬でおわるから。だからこれで十分実用なんだよねえ。どうかなアリストテレスもびっくりじゃないかこの手抜き加減。
子どもたちの中でもっといい方法があるよって知ってる人がいたら一応教えてくれよな。じゃ今日はこの辺で。