曲率と曲率半径を計算する

//
// ポリラインの角ポイントの曲率を計算する
// Run Over: Primitives
//
float step = 5.0;   // 曲率を計算する幅（片側）

int pts[] = primpoints(0, @primnum);

float curveLength = primintrinsic(0, "measuredperimeter", @primnum);

// ポリラインが連環しているか始点のつながりで確認する
int npt = neighbourcount(0, pts[0]);

for(int i = 0; i < len(pts); i++)
{
    vector pos = point(0, "P", pts[i]);
    int index;
    vector uv;
    xyzdist(0, pos, index, uv);
    
    // 端はクランプする
    float prevU = clamp(uv.x - (step / curveLength), 0, 1);
    float nextU = clamp(uv.x + (step / curveLength), 0, 1);
    
    vector p0 = primuv(0, "P", index, set(prevU, 0));
    vector p1 = point(0, "P", pts[i]);
    vector p2 = primuv(0, "P", index, set(nextU, 0));
    
    // XZ平面で計算する場合はこのチェックを外す
    //p0.y = 0;
    //p1.y = 0;
    //p2.y = 0;
    
    // 弧長
    float c = length(p2 - p0);
    // p1から一番近いP0-P2上の点
    vector p = p0 + dot((p1 - p0), normalize(p2 - p0)) * normalize(p2 - p0);
    // 矢高
    float h = length(p - p1);
    
    // ほとんど矢向が０に近い場合は無限半径（＝０）、またはほぼ直線とみなす
    if(h < 0.0001)
    {
        setpointattrib(0, "__radius", pts[i], 0.0);
        setpointattrib(0, "__theta", pts[i], PI);
        setpointattrib(0, "__curvature", pts[i], 0.0);
    }
    else
    {
        // 半径
        float radius = (c * c / 4 + h * h) / (2 * h);

        float theta = 2 * asin(c / (2 * radius));
        
        if(theta > radians(0.0))
        {
            setpointattrib(0, "__radius", pts[i], radius);
            setpointattrib(0, "__theta", pts[i], theta);
            setpointattrib(0, "__curvature", pts[i], 1.0/radius);
        }
    }
}

//
// ポリラインの曲率を視覚化する
// Run Over: Primitives
//
float minRadius = chf("radius"); // 上限の最小曲率半径(m)（赤く表示される値）

float maxCurvature = 1.0/minRadius;    // 上限の曲率
float step = 1.0;   // 曲率を計算する幅（片側）

float curveLength = primintrinsic(0, "measuredperimeter", @primnum);

int pts[] = primpoints(0, @primnum);
for(int i = 0; i < len(pts); i++)
{
    vector pos = point(0, "P", pts[i]);
    int index;
    vector uv;
    xyzdist(0, pos, index, uv);
    
    // 端はクランプする
    float prevU = clamp(uv.x - (step / curveLength), 0, 1);
    float nextU = clamp(uv.x + (step / curveLength), 0, 1);
    
    vector p0 = primuv(0, "P", index, set(prevU, 0));
    vector p1 = point(0, "P", pts[i]);
    vector p2 = primuv(0, "P", index, set(nextU, 0));
    
    // 弧長
    float c = length(p2 - p0);
    // p1から一番近いP0-P2上の点
    vector p = p0 + dot((p1 - p0), normalize(p2 - p0)) * normalize(p2 - p0);
    // 矢高
    float h = length(p - p1);
    
    // ほとんど矢向が０に近い場合は無限半径（＝０）とする
    if(h < 0.0001)
    {
        setpointattrib(0, "__radius", pts[i], 0.0);     // 曲率半径
        setpointattrib(0, "__theta", pts[i], PI);       // 中心角
        setpointattrib(0, "__curvature", pts[i], 0.0);  // 曲率
        
        vector color =  hsvtorgb(0.5, 1.0, 1.0);
        setpointattrib(0, "Cd", pts[i], color);
    }
    else
    {
        // 半径
        float radius = (c * c / 4 + h * h) / (2 * h);
        float curvature = 1.0/radius;
        
        // 中心角
        float theta = 2 * asin(c / (2 * radius));
        
        if(theta > radians(0.0))
        {
            setpointattrib(0, "__radius", pts[i], radius);         // 曲率半径
            setpointattrib(0, "__theta", pts[i], theta);           // 中心角
            setpointattrib(0, "__curvature", pts[i], curvature);   // 曲率
        
            float hue = clamp(curvature/maxCurvature, 0, 1.0) * 0.5 * -1 + 0.5;
            vector color =  hsvtorgb(hue, 1.0, 1.0);
            
            setpointattrib(0, "Cd", pts[i], color);
        }
    }
}