とあるデザイナーのテクニカルノート

Houdiniノート

制御点で描くクロソイド曲線

公開 更新

実装の流れ

2次ベジェ曲線の3つの制御点を利用して連続するクロソイド曲線を描画していく。 ポリラインをベジェ制御点に分割する

ポリラインを3ポイント（2辺）ごとに分割し、

それぞれのガイドカーブにクロソイド曲線を描いていく。

計算

クロソイドカーブの基本についてはこのページを参考に クロソイド曲線について

AとRからクロソイド曲線を計算する。R=1としてクロソイド曲線を描画した後、本来の大きさにスケールする、という流れ。

R = 1 τはIの1n\frac{1}{n}といった一定の値にする。

A=R2τA=R\sqrt{2τ} L=A2RL=\frac{A^2}{R}

1つ目のクロソイドカーブをR＝１のスケールで原点に描画する。

x=A2×2τ(1−12!×5τ2+14!×9τ4+16!×13τ6+⋅⋅⋅)x = \frac{A}{\sqrt{2}}\times2\sqrt{τ}\left( 1-\frac{1}{2!\times5}τ^2+\frac{1}{4!\times9}τ^4+\frac{1}{6!\times13}τ^6+···\right)

y=A2×23ττ(1−τ23!×73+τ45!×113+τ67!×153+⋅⋅⋅)y = \frac{A}{\sqrt{2}}\times\frac{2}{3}τ\sqrt{τ}\left( 1-\frac{τ^2}{\frac{3!\times7}{3}}+\frac{τ^4}{\frac{5!\times11}{3}}+\frac{τ^6}{\frac{7!\times15}{3}}+···\right)

クロソイドカーブの座標X,YはAとτから求まる。τを０から接線角まで刻んで入力して座標を計算する。

つぎに円弧を描画する。円の中心点はクロソイドカーブの終点から接線との直行ベクトルから計算する。円弧の半径は１、角度はthetaとなる。

2つ目のクロソイドカーブを描画し、X方向に反転したものをI（交角）回転させて円弧の最後の座標に移動させる行列をつくる。

クロソイド曲線２つと円弧を合成して、原点から伸びるカーブが出来る。

元のガイドカーブに合わせてスケールする。実際のガイドとの長さを比べて、スケールする値を求める。スケールさせてどちらかのガイドラインの余分な長さをオフセットさせる。

コード

//
// ガイドカーブをもとにクロソイド曲線を生成する
// Run Over: Primitives
// input1: Polyline(Control Points)

// クロソイドと円弧のおおよその割合
float angleRatio = 0.25;

// クロソイド曲線を描く間隔(m)
float resampleLength = 0.1;

// 直線の交点を求める関数
vector CrossPointXY(vector p0; vector p1; vector p2; vector p3)
{
    float a1 = p0.y - p1.y;
    float b1 = p1.x - p0.x;
    float c1 = (p1.x - p0.x) * p0.y - (p1.y - p0.y) * p0.x;
     
    float a2 = p2.y - p3.y;
    float b2 = p3.x - p2.x;
    float c2 = (p3.x - p2.x) * p2.y - (p3.y - p2.y) * p2.x;
     
    float x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
    float y = (a1 * c2 - a2 * c1) / (a1 * b2 - a2 * b1);

    return set(x, y, 0);
}

// クロソイド曲線のX成分を求める関数(tauとAから)
function float clothoid_x(float tau; float A)
{
    return A / sqrt(2) * 2 * sqrt(tau) * (1 - (pow(tau, 2) / 10.0) + (pow(tau, 4) / 216.0) - (pow(tau, 6) / 9360.0));
}

// クロソイド曲線のY成分を求める関数
function float clothoid_y(float tau; float A)
{
    return A / sqrt(2) * (2/3.0) * sqrt(tau) * tau * (1 - (pow(tau, 2) / 14.0) + (pow(tau, 4) / 440.0) - (pow(tau, 6) / 25200.0) );
}

//
// プリミティブごとにクロソイドを描く
//
int npts[] = primpoints(0, @primnum);
if(len(npts) > 2)
{
    //
    // 行列
    //
    vector p0 = point(0, "P", npts[0]);
    vector p1 = point(0, "P", npts[1]);
    vector p2 = point(0, "P", npts[2]);
    
    vector old_p0 = p0;
    vector old_p1 = p1;
    vector old_p2 = p2;
    
    vector tangent = normalize(p1 - p0);
    vector normal = normalize(cross(tangent, p2 - p1));
    vector binormal = cross(normal, tangent);
    matrix worldPlane = maketransform(normal, binormal, p0);
    matrix inversePlane = invert(worldPlane);
    
    p0 *= inversePlane;
    p1 *= inversePlane;
    p2 *= inversePlane;
    
    vector v0 = normalize(p1 - p0);
    vector v1 = normalize(p2 - p1);
    
    // 交角
    float I = acos(dot(v0, v1));
    
    //
    // X座標と接線(tau)からクロソイド曲線を描画する
    //
    float R = 1;

    float tau1 = I * angleRatio / 2;
    float A1 = R * sqrt(tau1*2);
    float L1 = A1 * A1 / R;
    float tau2 = I * angleRatio / 2;
    float A2 = R * sqrt(tau2*2);
    float L2 = A2 * A2 / R;
    float theta = I - (tau1 + tau2);
    
    //
    // クロソイドが描けるか判定
    //
    if(I < 0.001)
    {
        // 角度が浅くほぼ直線の場合はベジェ曲線を描画する
        int prim = addprim(0, "polyline");
    
        int num = 100;
        for(int i = 0; i < num; i++)
        {
            float t = i / float(num-1);
            vector pos = (1-t)*(1-t)*old_p0 + 2*(1-t)*t*old_p1 + t*t*old_p2;
            int pt = addpoint(0, pos);
            addvertex(0, prim, pt);
        }
        
        removeprim(0, @primnum, 1);
    }
    else
    {
        // 接線と同心円の接する座標（L1と同心円の接点）
        float x0 = clothoid_x(tau1, A1);
        float y0 = clothoid_y(tau1, A1);
        
        // クロソイド曲線(L1)
        vector curve_pos0[];
        int num = max(int(L1 / resampleLength), 2);
        for(int i = 0; i < num; i++)
        {
            float t = tau1 / float(num-1) * i;
            
            float x = clothoid_x(t, A1);
            float y = clothoid_y(t, A1);
            
            curve_pos0[i] = set(x, y, 0);
        }
        
        //
        // 円弧
        //
        // 円の中心座標と接線ベクトル
        vector vecTau1 = set(cos(tau1), sin(tau1), 0);
        vector center = set(x0, y0, 0) + cross(set(0,0,1), vecTau1) * R;
        
        vector diff = curve_pos0[-1] - center;
        float startRot = atan2(diff.y, diff.x);
        
        float arcPosLength = R * theta;
        num = max(int(arcPosLength / resampleLength), 2);
        float perAngle = theta / float(num);
        
        vector arcPos[];
        for(int i = 0; i < num+1; i++)
        {
            float angle = startRot + perAngle * i;
            arcPos[i] = set(cos(angle), sin(angle), 0) * R + center;
        }
        
        // クロソイド曲線(L2)
        vector curve_pos1[];
        num = max(int(L2 / resampleLength), 2);
        for(int i = 0; i < num; i++)
        {
            float t = tau2 / float(num-1) * i;
            
            float x = clothoid_x(t, A2);
            float y = clothoid_y(t, A2);
            
            curve_pos1[i] = set(x, y, 0);
        }
        
        // 片方のクロソイド曲線を反転させる行列
        matrix mirrorWorld = maketransform(set(0,0,1), set(0,1,0), set(0,0,0));
        
        scale(mirrorWorld, set(-1, 1, 1));
        translate(mirrorWorld, set(curve_pos1[-1].x, -curve_pos1[-1].y, 0));
        rotate(mirrorWorld, I, set(0,0,1));   // 交角
        
        translate(mirrorWorld, arcPos[-1]);
        
        for(int i = 0; i < num; i++)
            curve_pos1[i] *= mirrorWorld;
        
        // ローカルのガイド
        vector local_p0 = curve_pos0[0];
        vector local_p2 = curve_pos1[0];
        vector local_p1 = CrossPointXY(set(0,0,0), set(1,0,0), curve_pos1[0], curve_pos1[0] + (p1 - p2));
        
        // ローカルのガイドは余白がない、元のガイドは余白を含むのでそこを比較する
        float localRatio = length(local_p1-local_p0)/length(local_p1-local_p2);
        float ratio = length(p1-p0)/length(p1-p2);
        
        // ガイドに合わせてオフセットさせる
        matrix guidWorld = maketransform(set(0,0,1), set(0,1,0), set(0,0,0));
        float scale = 1;
        
        //底辺のほうが長い場合
        if(ratio > localRatio)
        {
            scale = length(p1-p2)/length(local_p1-local_p2);
            float offset = length(p1-p0)-length(local_p1-local_p0)*scale;
            
            scale(guidWorld, set(scale, scale, scale));
            translate(guidWorld, set(offset, 0, 0));
        }
        else
        {
            scale = length(p1-p0)/length(local_p1-local_p0);
            scale(guidWorld, set(scale, scale, scale));
        }
        
        // 端（直線）
        if(ratio > localRatio)
        {
            int prim = addprim(0, "polyline");
            
            int pt = addpoint(0, p0 * worldPlane);
            addvertex(0, prim, pt);
            pt = addpoint(0, curve_pos0[0] * guidWorld * worldPlane);
            addvertex(0, prim, pt);
            setprimgroup(0, "__line", prim, 1);
        }
    
        // クロソイド
        int prim = addprim(0, "polyline");
        for(int i = 0; i < len(curve_pos0); i++)
        {
            int pt = addpoint(0, curve_pos0[i] * guidWorld * worldPlane);
            addvertex(0, prim, pt);
            setprimgroup(0, "__clothoid", prim, 1);
        }
        
        // 円弧
        prim = addprim(0, "polyline");
        for(int i = 0; i < len(arcPos); i++)
        {
            int pt = addpoint(0, arcPos[i] * guidWorld * worldPlane);
            addvertex(0, prim, pt);
            setprimgroup(0, "__arc", prim, 1);
        }
        
        // クロソイド
        prim = addprim(0, "polyline");
        curve_pos1 = reverse(curve_pos1);
        for(int i = 0; i < len(curve_pos1); i++)
        {
            int pt = addpoint(0, curve_pos1[i] * guidWorld * worldPlane);
            addvertex(0, prim, pt);
            setprimgroup(0, "__clothoid", prim, 1);
        }
        
        // 端（直線）
        if(ratio < localRatio)
        {
            prim = addprim(0, "polyline");
            
            int pt = addpoint(0, curve_pos1[-1] * guidWorld * worldPlane);
            addvertex(0, prim, pt);
            pt = addpoint(0, p2 * worldPlane);
            addvertex(0, prim, pt);
            setprimgroup(0, "__line", prim, 1);
        }
        
        removeprim(0, @primnum, 1);
    }
}

Polylineクロソイド曲線線形

← 記事一覧へ