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 -> fitCurves と camelCaseにしただけではなく、Curve -> Curves と複数形に変更しました。 ・FitCurve のメソッド内メソッドとしてほとんどの関数を作成しました。
curves というベジェ曲線配列を作成。ベジェ曲線が確定した時点でこの配列に格納していきます。
static fitCurves(d, error) { const curves = []; const nPts = d.length; const tHat1 = computeLeftTangent(d, 0); const tHat2 = computeRightTangent(d, nPts - 1); fitCubic(d, 0, nPts - 1, tHat1, tHat2, error); return curves; }
・FitCubic の uPrime は for文のスコープとしました。
function fitCubic(d, first, last, tHat1, tHat2, error) { let bezCurve, u, maxError, splitPoint; const nPts = last - first + 1; const maxIterations = 4; const iterationError = error * 4.0; 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; } u = chordLengthParameterize(d, first, last); bezCurve = generateBezier(d, first, last, u, tHat1, tHat2); [maxError, splitPoint] = computeMaxError(d, first, last, bezCurve, u); if(maxError < error) { curves.push(bezCurve); return; } if(maxError < iterationError) { for(let i = 0; i < maxIterations; i += 1) { const uPrime = reparameterize(d, first, last, u, bezCurve); 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; } } let tHatCenter = computeCenterTangent(d, splitPoint); fitCubic(d, first, splitPoint, tHat1, tHatCenter, error); tHatCenter = Vector2.negate(tHatCenter); fitCubic(d, splitPoint, last, tHatCenter, tHat2, error); }
・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) で置換しました。
function generateBezier(d, first, last, uPrime, tHat1, tHat2) { const nPts = last - first + 1; const A = new Array(nPts); const bezCurve = new Array(4); 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]; } 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] = 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); } 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]; 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; const segLength = Vector2.distance(d[last], d[first]); const epsilon = 1.0e-6 * segLength; if (alpha_l < epsilon || alpha_r < epsilon) { 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; } 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 で始点、終点におけるパラメータは再計算する必要がないので u[0] = 0, uPrime[nPts - 1] = 1 としました。
function reparameterize(d, first, last, u, bezCurve) { const nPts = last - first + 1; const uPrime = new Array(nPts); uPrime[0] = 0; uPrime[nPts - 1] = 1; for(let i = first + 1; i <= last - 1; i += 1) { uPrime[i - first] = newtonRaphsonRootFind(bezCurve, d[i], u[i - first]); } return uPrime; }
function newtonRaphsonRootFind(Q, P, u) { const Q1 = new Array(3); const Q2 = new Array(2); const Q_u = bezierII(3, Q, u); for(let i = 0; i <= 2; i += 1) { Q1[i] = Vector2.scale(Vector2.subtract(Q[i + 1], Q[i]), 3.0); } for(let i = 0; i <= 1; i += 1) { Q2[i] = Vector2.scale(Vector2.subtract(Q1[i + 1], Q1[i]), 2.0); } const Q1_u = bezierII(2, Q1, u); const Q2_u = bezierII(1, Q2, 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; const uPrime = u - numerator / denominator; return uPrime; }
function bezierII(degree, V, t) { const Vtemp = new Array(degree + 1); for(let i = 0; i <= degree; i += 1) { Vtemp[i] = Vector2.copy(V[i]); } 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]; return Q; }
・ComputeLeftTangent の第二引数の名称を end から start に変更しました。
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; }
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 の tHatCenter は / 2 する必要は無いので / 2 の演算を省略しました。
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; }
function chordLengthParameterize(d, first, last) { if(first < 0 || last > d.length - 1) { throw 'Invalid Parameters. (chordLengthParameterize)'; } const u = new Array(last - first + 1); 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; }
・maxDist と splitPoint を返すよう変更しました。
function computeMaxError(d, first, last, bezCurve, u) { let splitPoint = parseInt((last - first + 1) / 2); let maxDist = 0.0; for(let i = first + 1; i < last; i += 1) { const P = bezierII(3, bezCurve, u[i - first]); const v = Vector2.subtract(P, d[i]); const dist = Vector2.squaredLength(v); if(dist >= maxDist) { maxDist = dist; splitPoint = i; } } return [ maxDist, splitPoint, ]; }