三角プリミティブを対角線で分割する

//
// 三角形ポリゴンの長いエッジから隣り合うプリミティブを探す
// Run OVer: Primitive
//
int prim_vtxs[] = primvertices(0, @primnum);
int pts[];

// 三角形なら
if (len(prim_vtxs) == 3)
{
    float max_length = 0;
    int longest_edge;
    foreach (int vtx; prim_vtxs)
    {
        int edge = vertexhedge(0, vtx);
        int src_pt = hedge_srcpoint(0, edge);
        int dst_pt = hedge_dstpoint(0, edge);
        vector src_pos = point(0, "P", src_pt);
        vector dst_pos = point(0, "P", dst_pt);
        float length = distance(src_pos, dst_pos);
        if (length > max_length)
        {
            longest_edge = edge;
            max_length = length;
            pts[0] = src_pt;
            pts[1] = dst_pt;
        }
    }  
    
    // プライマリエッジを取得する
    // 隣接した面のエッジと同じかどうかを比較するにはプライマリエッジを取得する必要あり
    int edge = hedge_nextequiv(0, longest_edge);
    int prim = hedge_prim(0, edge);
    
    if(prim != @primnum)
        i@__pare_prim = prim;
    else
        i@__pare_prim = -1;
        
    // エッジのポイントインデックス
    i[]@__common_pts = pts;
    // 三角形のポイントインデックス
    int primpoints[] = primpoints(0, @primnum);
    i[]@__primpoints = primpoints;
    
    // 
    for(int i = 0; i < len(pts); i++)
        removevalue(primpoints, pts[i]);
    
    i@__corner_pt = primpoints[0];
}

//
// 線分の交差判定
// Run Over: Primitive
//
int IsIntersectLines(vector A; vector B; vector C; vector D; export vector Cross)
{
    // 直線にタッチしているポイントも考慮するための許容誤差
    float epsilon = 0.0001;
    
    int result = -1;    

    float denominator = (B.z - A.z) * (D.x - C.x) - (B.x - A.x) * (D.z - C.z);
    
    if(denominator != 0)
    {
        float s = ((B.x - A.x) * (C.z - A.z) - (B.z - A.z) * (C.x - A.x)) / denominator;
        if(s >= 0 - epsilon && s <= 1.0 + epsilon)
        {
            float r = ((A.x - C.x) * (D.z - C.z) - (A.z - C.z) * (D.x - C.x)) / denominator;
            
            if(r >= 0 - epsilon && r <= 1.0 + epsilon)
            {
                result = 1;
                Cross = A + (B - A) * r;
            }
        }
    }
    
    return result;
}

// ペアのプリミティブがある場合
if(i@__pare_prim != -1)
{
    // ペアのプリミティブから見たペアの番号
    int pare_prim = prim(0, "__pare_prim", i@__pare_prim);
    
    // 互いがペアの場合
    if(pare_prim == @primnum && @primnum > i@__pare_prim)
    {
        // 色を塗ってみる
        //@Cd = set(1,0,0);
        //setprimattrib(0, "Cd", i@pare_prim, set(1, 0.5, 0));
        
        // 中間点の計算
        int cornerPoint = prim(0, "__corner_pt", i@__pare_prim);
        vector p0 = point(0, "P", i@__corner_pt);
        vector p1 = point(0, "P", i[]@__common_pts[0]);
        vector p2 = point(0, "P", i[]@__common_pts[1]);
        vector p3 = point(0, "P", cornerPoint);
        
        vector cross = (p1 + p2) / 2;   // とりあえず中間点
        
        // 2次元的に交差ポイントを計算する
        IsIntersectLines(p0, p3, p1, p2, cross);
        
        // 対角線との一番近いポイントを探す
        vector ab = normalize(p2 - p1);
        float d = abs(dot(p1-cross, ab));
        cross = p1 + ab * d; // 交点座標
        
        int prim = addprim(0, "poly");
        int center_pt = addpoint(0, cross);
        addvertex(0, prim, i@__corner_pt);
        addvertex(0, prim, i[]@__common_pts[0]);
        addvertex(0, prim, center_pt);
        
        prim = addprim(0, "poly");
        addvertex(0, prim, i@__corner_pt);
        addvertex(0, prim, center_pt);
        addvertex(0, prim, i[]@__common_pts[1]);
        
        // 
        int common_pts[] =  prim(0, "__common_pts", i@__pare_prim);
        
        prim = addprim(0, "poly");
        addvertex(0, prim, cornerPoint);
        addvertex(0, prim, common_pts[0]);
        addvertex(0, prim, center_pt);
        
        prim = addprim(0, "poly");
        addvertex(0, prim, cornerPoint);
        addvertex(0, prim, center_pt);
        addvertex(0, prim, common_pts[1]);

        // 元のプリミティブを削除
        removeprim(0, @primnum, 1);
        removeprim(0, i@__pare_prim, 1);
    }
}