FitCurve.c を JavaScript に移植

はじめに

Graphics Gems 1 にある AN ALGORITHM FOR AUTOMATICALLY FITTING DIGITIZED CURVES という
アルゴリズムをC言語で書いた FitCurve.c とソースがあります。
それを2024年12月22日に JavaScript に移植したときのログのような記事です。
移植前を手元に用意して本記事をご覧ください。

全般的な処理

ソースが書かれたのが1990年で且つC言語なので、全体的に少しだけ直す必要があります。

・return (a) の括弧を外しました。
・メソッド名は camelCase に変更した。e.g. GenerateBezier -> generateBezier
・変数名は大文字で始まろうともそのままにしました。
・変数のスコープを狭くしても良いものはスコープを狭くしました。
・malloc(sizeof(T) * n) は new Array(n) で置き換えました。
　空配列で初期化しようとも考えたが new Array(n) だとメモリの再割り当ても無いだろうと判断しました。
・ベクトルの演算は Vector2 のメソッドで代替しました。

それでは関数ごとに見ていきたいと思います。

FitCurve

・FitCurve -> fitCurves と camelCaseにしただけではなく、Curve -> Curves と複数形に変更しました。
・FitCurve のメソッド内メソッドとしてほとんどの関数を作成しました。
　curves というベジェ曲線配列を作成。ベジェ曲線が確定した時点でこの配列に格納していきます。

/**
 * Fit a Bezier curve to a set of digitized points 
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} error User-defined error squared
 * @returns {Array<Bezier2>} Array of cubic bezier curve
 */
static fitCurves(d, error) {
    const curves = [];      // Array of cubic bezier curve
    const nPts = d.length;  // Number of digitized points

    // Unit tangent vectors at endpoints
    const tHat1 = computeLeftTangent(d, 0);
    const tHat2 = computeRightTangent(d, nPts - 1);
    fitCubic(d, 0, nPts - 1, tHat1, tHat2, error);

    return curves;
}

FitCubic

・FitCubic の uPrime は for文のスコープとしました。

/**
 * Fit a Bezier curve to a (sub)set of digitized points
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} first Index of first point in region
 * @param {number} last Index of last point in region
 * @param {Vector2} tHat1 Unit tangent vector at first point
 * @param {Vector2} tHat2 Unit tangent vector at last point
 * @param {number} error User-defined error squared
 */
function fitCubic(d, first, last, tHat1, tHat2, error) {
    let bezCurve,   // Control points of fitted Bezier curve
        u,  //  Parameter values for point
        maxError,   // Maximum fitting error
        splitPoint; // Point to split point set at
    const nPts = last - first + 1;  // Number of points in subset
    const maxIterations = 4;    // Max times to try iterating
    const iterationError = error * 4.0;	// Error below which you try iterating - fixed issue 23            

    // Use heuristic if region only has two points in it
    if(nPts === 2) {
        const dist = Vector2.distance(d[last], d[first]) / 3.0;
        const bezCurve = new Array(4);
        bezCurve[0] = d[first];
        bezCurve[3] = d[last];
        bezCurve[1] = Vector2.add(bezCurve[0], Vector2.scale(tHat1, dist));
        bezCurve[2] = Vector2.add(bezCurve[3], Vector2.scale(tHat2, dist));
        curves.push(bezCurve);
        return;
    }            

    // Parameterize points, and attempt to fit curve
    u = chordLengthParameterize(d, first, last);
    bezCurve = generateBezier(d, first, last, u, tHat1, tHat2);

    // Find max deviation of points to fitted curve */
    [maxError, splitPoint] = computeMaxError(d, first, last, bezCurve, u);
    if(maxError < error) {
        curves.push(bezCurve);
        return;
    }

    // If error not too large, try some reparameterization and iteration
    if(maxError < iterationError) {
        for(let i = 0; i < maxIterations; i += 1) {
            const uPrime = reparameterize(d, first, last, u, bezCurve); // Improved parameter values
            bezCurve = generateBezier(d, first, last, uPrime, tHat1, tHat2);
            [maxError, splitPoint] = computeMaxError(d, first, last, bezCurve, uPrime);
            if(maxError < error) {
                curves.push(bezCurve);
                return;
            }
            u = uPrime;
        }
    }

    // Fitting failed -- split at max error point and fit recursively
    let tHatCenter = computeCenterTangent(d, splitPoint); // Unit tangent vector at splitPoint
    fitCubic(d, first, splitPoint, tHat1, tHatCenter, error);
    tHatCenter = Vector2.negate(tHatCenter);
    fitCubic(d, splitPoint, last, tHatCenter, tHat2, error);
}

GenerateBezier

・GenerateBezier の V2Scale(&tHat1, alpha_l) で tHat1 はアドレスを渡して値が変更されているが
　そもそも GenerateBezier 呼び出しの時点で tHat1 は参照渡しではないため、
　JSでは Vector2.scale(tHat1, dist) のように実装しました。本関数呼び出しでtHat1は変わりません。tHat2も同様です。
・GenerateBezier で使用される B0, B1, B2, B3 メソッドを Scalar.bernsteinBasis(3, r, u) で置換しました。

/**
 * Use least-squares method to find Bezier control points for region.
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} first Index of first point in region
 * @param {number} last Index of last point in region
 * @param {Array<number>} uPrime Parameter values for region
 * @param {Vector2} tHat1 Unit tangent vector at first point
 * @param {Vector2} tHat2 Unit tangent vector at last point
 * @returns {Array<Vector2>} Bezier control points for region
 */
function generateBezier(d, first, last, uPrime, tHat1, tHat2) {
    const nPts = last - first + 1;
    const A = new Array(nPts);
    const bezCurve = new Array(4);            

    // Compute the A's
    for(let i = 0; i < nPts; i += 1) {
        const v1 = Vector2.scale(tHat1, Scalar.bernsteinBasis(3, 1, uPrime[i]));
        const v2 = Vector2.scale(tHat2, Scalar.bernsteinBasis(3, 2, uPrime[i]));
        A[i] = [v1, v2];
    }

    // Create the C and X matrices
    const C = [ [0, 0], [0, 0] ];
    const X = [0, 0];

    for (let i = 0; i < nPts; i++) {
        C[0][0] += Vector2.dot(A[i][0], A[i][0]);
        C[0][1] += Vector2.dot(A[i][0], A[i][1]);
        // C[1][0] += Vector2.dot(A[i][0], A[i][1]);	
        C[1][0] = C[0][1];
        C[1][1] += Vector2.dot(A[i][1], A[i][1]);

        const tmp = Vector2.subtract(
            d[first + i],
            Vector2.add(
                Vector2.scale(d[first], Scalar.bernsteinBasis(3, 0, uPrime[i])),
                Vector2.add(
                    Vector2.scale(d[first], Scalar.bernsteinBasis(3, 1, uPrime[i])),
                    Vector2.add(
                        Vector2.scale(d[last], Scalar.bernsteinBasis(3, 2, uPrime[i])),
                        Vector2.scale(d[last], Scalar.bernsteinBasis(3, 3, uPrime[i]))
                    )
                )
            )
        );

        X[0] += Vector2.dot(A[i][0], tmp);
        X[1] += Vector2.dot(A[i][1], tmp);
    }

    // Compute the determinants of C and X
    const det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
    const det_C0_X  = C[0][0] * X[1]    - C[1][0] * X[0];
    const det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];

    // Finally, derive alpha values
    const alpha_l = (det_C0_C1 === 0) ? 0.0 : det_X_C1 / det_C0_C1;
    const alpha_r = (det_C0_C1 === 0) ? 0.0 : det_C0_X / det_C0_C1;

    // If alpha negative, use the Wu/Barsky heuristic (see text)
    // (if alpha is 0, you get coincident control points that lead to
    // divide by zero in any subsequent NewtonRaphsonRootFind() call.
    const segLength = Vector2.distance(d[last], d[first]);
    const epsilon = 1.0e-6 * segLength;
    if (alpha_l < epsilon || alpha_r < epsilon) {
        // fall back on standard (probably inaccurate) formula, and subdivide further if needed. */
        const dist = segLength / 3.0;
        bezCurve[0] = Vector2.copy(d[first]);
        bezCurve[3] = Vector2.copy(d[last]);
        bezCurve[1] = Vector2.add(bezCurve[0], Vector2.scale(tHat1, dist));
        bezCurve[2] = Vector2.add(bezCurve[3], Vector2.scale(tHat2, dist));
        return bezCurve;
    }

    // First and last control points of the Bezier curve are
    // positioned exactly at the first and last data points
    // Control points 1 and 2 are positioned an alpha distance out
    // on the tangent vectors, left and right, respectively
    bezCurve[0] = Vector2.copy(d[first]);
    bezCurve[3] = Vector2.copy(d[last]);
    bezCurve[1] = Vector2.add(bezCurve[0], Vector2.scale(tHat1, alpha_l));
    bezCurve[2] = Vector2.add(bezCurve[3], Vector2.scale(tHat2, alpha_r));
    return bezCurve;
}

Reparameterize

reparameterize で始点、終点におけるパラメータは再計算する必要がないので u[0] = 0, uPrime[nPts - 1] = 1 としました。

/**
 * Given set of points and their parameterization, try to find a better parameterization.
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} first Index of first point in region
 * @param {number} last Index of last point in region
 * @param {Array<number>} u Current parameter values
 * @param {Array<Vector2>} bezCurve Current fitted curve
 * @returns {Array<number>} Parameterization of points
 */
function reparameterize(d, first, last, u, bezCurve) {
    const nPts = last - first + 1;
    const uPrime = new Array(nPts); // New parameter values
    
    // No need to calculate the values <200b><200b>corresponding to the start and end points
    uPrime[0] = 0;          // The value corresponding to the start point is 0
    uPrime[nPts - 1] = 1;   // The value corresponding to the end point is 1

    // Calculate the value corresponding to a point excluding the start and end points
    for(let i = first + 1; i <= last - 1; i += 1) {
        uPrime[i - first] = newtonRaphsonRootFind(bezCurve, d[i], u[i - first]);
    }
    return uPrime;
}

NewtonRaphsonRootFind

/**
 * Use Newton-Raphson iteration to find better root.
 * @param {Array<Vector2>} Q Current fitted curve
 * @param {Vector2} P Digitized point
 * @param {number} u Parameter value for "P"
 * @returns {number} Improved parameter value for "P"
 */
function newtonRaphsonRootFind(Q, P, u) {
    const Q1 = new Array(3);    // Q'
    const Q2 = new Array(2);    // Q''

    // Compute Q(u)
    const Q_u = bezierII(3, Q, u);  // u evaluated at Q

    // Generate control vertices for Q'
    for(let i = 0; i <= 2; i += 1) {
        Q1[i] = Vector2.scale(Vector2.subtract(Q[i + 1], Q[i]), 3.0);
    }

    // Generate control vertices for Q''
    for(let i = 0; i <= 1; i += 1) {
        Q2[i] = Vector2.scale(Vector2.subtract(Q1[i + 1], Q1[i]), 2.0);
    }

    // Compute Q'(u) and Q''(u)
    const Q1_u = bezierII(2, Q1, u);    // u evaluated at Q'
    const Q2_u = bezierII(1, Q2, u);    // u evaluated at Q''

    // Compute f(u)/f'(u)
    const numerator = (Q_u[0] - P[0]) * (Q1_u[0]) + (Q_u[1] - P[1]) * (Q1_u[1]);
    const denominator = (Q1_u[0]) * (Q1_u[0]) + (Q1_u[1]) * (Q1_u[1]) +
                        (Q_u[0] - P[0]) * (Q2_u[0]) + (Q_u[1] - P[1]) * (Q2_u[1]);
    if(denominator === 0.0) return u;

    // u = u - f(u)/f'(u)
    const uPrime = u - numerator / denominator; // Improved u
    return uPrime;
}

BezierII

/**
 * Evaluate a Bezier curve at a particular parameter value
 * @param {number} degree The degree of the bezier curve
 * @param {Array<Vector2>} V Array of control points
 * @param {number} t Parametric value to find point for
 * @returns {Vector2} Point on curve at parameter t
 */
function bezierII(degree, V, t) {
    // Copy array
    const Vtemp = new Array(degree + 1);  // Local copy of control points
    for(let i = 0; i <= degree; i += 1) {
        Vtemp[i] = Vector2.copy(V[i]);
    }
    // Triangle computation
    for(let i = 1; i <= degree; i += 1) {	
        for(let j = 0; j <= degree - i; j += 1) {
            Vector2.lerp(Vtemp[j], Vtemp[j + 1], t, Vtemp[j]);
        }
    }
    const Q = Vtemp[0]; // Point on curve at parameter t
    return Q;
}

ComputeLeftTangent

・ComputeLeftTangent の第二引数の名称を end から start に変更しました。

/**
 * Approximate unit tangents at endpoints and "left" of digitized curve
 * @param {Array<Vector2>} d Digitized points
 * @param {number} start Index to "left" end of region
 * @returns {Vector2} Unit tangents at endpoints
 */
function computeLeftTangent(d, start) {
    if(start + 1 > d.length - 1) { throw 'Invalid Parameters. (computeLeftTangent)'; }
    const tHat1 = Vector2.subtract(d[start + 1], d[start]);
    Vector2.normalize(tHat1, tHat1);
    return tHat1;
}

ComputeRightTangent

/**
 * Approximate unit tangents at endpoints and "right" of digitized curve
 * @param {Array<Vector2>} d Digitized points
 * @param {number} end Index to "right" end of region
 * @returns {Vector2} Unit tangents at endpoints
 */
function computeRightTangent(d, end) {
    if(end - 1 < 0) { throw 'Invalid Parameters. (computeRightTangent)'; }
    const tHat2 = Vector2.subtract(d[end - 1], d[end]);
    Vector2.normalize(tHat2, tHat2);
    return tHat2;
}

ComputeCenterTangent

・ComputeCenterTangent の tHatCenter は / 2 する必要は無いので / 2 の演算を省略しました。

/**
 * Approximate unit tangents at endpoints and "center" of digitized curve
 * @param {Array<Vector2>} d Digitized points
 * @param {number} center Index to "center" end of region
 * @returns {Vector2} Unit tangents at endpoints
 */
function computeCenterTangent(d, center) {
    if(center - 1 < 0 || center + 1 > d.length - 1) { throw 'Invalid Parameters. (computeCenterTangent)'; }
    const V1 = Vector2.subtract(d[center - 1], d[center]);
    const V2 = Vector2.subtract(d[center], d[center + 1]);
    const tHatCenter = Vector2.add(V1, V2);
    Vector2.normalize(tHatCenter, tHatCenter);
    return tHatCenter;
}

ChordLengthParameterize

/**
 * Assign parameter values to digitized points using relative distances between points.
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} first Indices defining region
 * @param {number} last Indices defining region
 * @returns {Array<number>} Parameterization of points
 */
function chordLengthParameterize(d, first, last) {
    if(first < 0 || last > d.length - 1) { throw 'Invalid Parameters. (chordLengthParameterize)'; }
    const u = new Array(last - first + 1);  // Parameterization
    u[0] = 0.0;
    for(let i = first + 1; i <= last; i += 1) {
        u[i - first] = u[i - first - 1] + Vector2.distance(d[i], d[i - 1]);
    }
    for(let i = first + 1; i <= last; i += 1) {
        u[i - first] = u[i - first] / u[last - first];
    }
    return u;
}

ComputeMaxError

・maxDist と splitPoint を返すよう変更しました。

/**
 * Find the maximum squared distance of digitized points to fitted curve.
 * @param {Array<Vector2>} d Array of digitized points
 * @param {number} first Indices defining region
 * @param {number} last Indices defining region
 * @param {Array<Vector2>} bezCurve Fitted Bezier curve
 * @param {number} u Parameterization of points
 */
function computeMaxError(d, first, last, bezCurve, u) {
    let splitPoint = parseInt((last - first + 1) / 2);
    let maxDist = 0.0;  // Maximum error
    for(let i = first + 1; i < last; i += 1) {
        const P = bezierII(3, bezCurve, u[i - first]);    // Point on curve
        const v = Vector2.subtract(P, d[i]);    // Vector from point to curve
        const dist = Vector2.squaredLength(v);  // Current error
        if(dist >= maxDist) {
            maxDist = dist;
            splitPoint = i;
        }
    }
    return [ maxDist, splitPoint, ];
}