diff --git a/lib/fit-curve.js b/lib/fit-curve.js new file mode 100644 index 0000000..ea0ebc0 --- /dev/null +++ b/lib/fit-curve.js @@ -0,0 +1,225 @@ +/* + JavaScript implementation of + Algorithm for Automatically Fitting Digitized Curves + by Philip J. Schneider + "Graphics Gems", Academic Press, 1990 + + The MIT License (MIT) + + Original (C): + https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c + -> Python: + https://github.com/volkerp/fitCurves + -> CoffeeScript/JavaScript + math.js/lodash: + https://github.com/soswow/fit-curves + -> JavaScript (ES6-ish), no dependencies + https://github.com/Sphinxxxx/fit-curve +*/ +(function(root, factory) { + if (typeof define === "function" && define.amd) { + define([], factory); + } else { + if (typeof module === "object" && module.exports) { + module.exports = factory(); + } else { + root.fitCurve = factory(); + } + } +})(this, function() { + function _classCallCheck(instance, Constructor) { + if (!(instance instanceof Constructor)) { + throw new TypeError("Cannot call a class as a function"); + } + } + var maths = function() { + function maths() { + _classCallCheck(this, maths); + } + maths.zeros_Xx2x2 = function zeros_Xx2x2(x) { + var zs = []; + while (x--) { + zs.push([0, 0]); + } + return zs; + }; + maths.mulItems = function mulItems(items, multiplier) { + return [items[0] * multiplier, items[1] * multiplier]; + }; + maths.mulMatrix = function mulMatrix(m1, m2) { + return m1[0] * m2[0] + m1[1] * m2[1]; + }; + maths.subtract = function subtract(arr1, arr2) { + return [arr1[0] - arr2[0], arr1[1] - arr2[1]]; + }; + maths.addArrays = function addArrays(arr1, arr2) { + return [arr1[0] + arr2[0], arr1[1] + arr2[1]]; + }; + maths.addItems = function addItems(items, addition) { + return [items[0] + addition, items[1] + addition]; + }; + maths.sum = function sum(items) { + return items.reduce(function(sum, x) { + return sum + x; + }); + }; + maths.dot = function dot(m1, m2) { + return maths.mulMatrix(m1, m2); + }; + maths.vectorLen = function vectorLen(v) { + var a = v[0], b = v[1]; + return Math.sqrt(a * a + b * b); + }; + maths.divItems = function divItems(items, divisor) { + return [items[0] / divisor, items[1] / divisor]; + }; + maths.squareItems = function squareItems(items) { + var a = items[0], b = items[1]; + return [a * a, b * b]; + }; + maths.normalize = function normalize(v) { + return this.divItems(v, this.vectorLen(v)); + }; + return maths; + }(); + var bezier = function() { + function bezier() { + _classCallCheck(this, bezier); + } + bezier.q = function q(ctrlPoly, t) { + var tx = 1 - t; + var pA = maths.mulItems(ctrlPoly[0], tx * tx * tx), pB = maths.mulItems(ctrlPoly[1], 3 * tx * tx * t), pC = maths.mulItems(ctrlPoly[2], 3 * tx * t * t), pD = maths.mulItems(ctrlPoly[3], t * t * t); + return maths.addArrays(maths.addArrays(pA, pB), maths.addArrays(pC, pD)); + }; + bezier.qprime = function qprime(ctrlPoly, t) { + var tx = 1 - t; + var pA = maths.mulItems(maths.subtract(ctrlPoly[1], ctrlPoly[0]), 3 * tx * tx), pB = maths.mulItems(maths.subtract(ctrlPoly[2], ctrlPoly[1]), 6 * tx * t), pC = maths.mulItems(maths.subtract(ctrlPoly[3], ctrlPoly[2]), 3 * t * t); + return maths.addArrays(maths.addArrays(pA, pB), pC); + }; + bezier.qprimeprime = function qprimeprime(ctrlPoly, t) { + return maths.addArrays(maths.mulItems(maths.addArrays(maths.subtract(ctrlPoly[2], maths.mulItems(ctrlPoly[1], 2)), ctrlPoly[0]), 6 * (1 - t)), maths.mulItems(maths.addArrays(maths.subtract(ctrlPoly[3], maths.mulItems(ctrlPoly[2], 2)), ctrlPoly[1]), 6 * t)); + }; + return bezier; + }(); + function fitCurve(points, maxError) { + var len = points.length, leftTangent = maths.normalize(maths.subtract(points[1], points[0])), rightTangent = maths.normalize(maths.subtract(points[len - 2], points[len - 1])); + return fitCubic(points, leftTangent, rightTangent, maxError); + } + function fitCubic(points, leftTangent, rightTangent, error) { + var MaxIterations = 20; + var bezCurve, u, uPrime, maxError, splitPoint, centerTangent, beziers, dist, i; + if (points.length === 2) { + dist = maths.vectorLen(maths.subtract(points[0], points[1])) / 3; + bezCurve = [points[0], maths.addArrays(points[0], maths.mulItems(leftTangent, dist)), maths.addArrays(points[1], maths.mulItems(rightTangent, dist)), points[1]]; + return [bezCurve]; + } + u = chordLengthParameterize(points); + bezCurve = generateBezier(points, u, leftTangent, rightTangent); + var _computeMaxError = computeMaxError(points, bezCurve, u); + maxError = _computeMaxError[0]; + splitPoint = _computeMaxError[1]; + if (maxError < error) { + return [bezCurve]; + } + if (maxError < error * error) { + for (i = 0;i < MaxIterations;i++) { + uPrime = reparameterize(bezCurve, points, u); + bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent); + var _computeMaxError2 = computeMaxError(points, bezCurve, uPrime); + maxError = _computeMaxError2[0]; + splitPoint = _computeMaxError2[1]; + if (maxError < error) { + return [bezCurve]; + } + u = uPrime; + } + } + beziers = []; + centerTangent = maths.normalize(maths.subtract(points[splitPoint - 1], points[splitPoint + 1])); + beziers = beziers.concat(fitCubic(points.slice(0, splitPoint + 1), leftTangent, centerTangent, error)); + beziers = beziers.concat(fitCubic(points.slice(splitPoint), maths.mulItems(centerTangent, -1), rightTangent, error)); + return beziers; + } + function generateBezier(points, parameters, leftTangent, rightTangent) { + var bezCurve, A, a, C, X, det_C0_C1, det_C0_X, det_X_C1, alpha_l, alpha_r, epsilon, segLength, i, len, tmp, u, ux, firstPoint = points[0], lastPoint = points[points.length - 1]; + bezCurve = [firstPoint, null, null, lastPoint]; + A = maths.zeros_Xx2x2(parameters.length); + for (i = 0, len = parameters.length;i < len;i++) { + u = parameters[i]; + ux = 1 - u; + a = A[i]; + a[0] = maths.mulItems(leftTangent, 3 * u * (ux * ux)); + a[1] = maths.mulItems(rightTangent, 3 * ux * (u * u)); + } + C = [[0, 0], [0, 0]]; + X = [0, 0]; + for (i = 0, len = points.length;i < len;i++) { + u = parameters[i]; + a = A[i]; + C[0][0] += maths.dot(a[0], a[0]); + C[0][1] += maths.dot(a[0], a[1]); + C[1][0] += maths.dot(a[0], a[1]); + C[1][1] += maths.dot(a[1], a[1]); + tmp = maths.subtract(points[i], bezier.q([firstPoint, firstPoint, lastPoint, lastPoint], u)); + X[0] += maths.dot(a[0], tmp); + X[1] += maths.dot(a[1], tmp); + } + det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; + det_C0_X = C[0][0] * X[1] - C[1][0] * X[0]; + det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; + alpha_l = det_C0_C1 === 0 ? 0 : det_X_C1 / det_C0_C1; + alpha_r = det_C0_C1 === 0 ? 0 : det_C0_X / det_C0_C1; + segLength = maths.vectorLen(maths.subtract(firstPoint, lastPoint)); + epsilon = 1E-6 * segLength; + if (alpha_l < epsilon || alpha_r < epsilon) { + bezCurve[1] = maths.addArrays(firstPoint, maths.mulItems(leftTangent, segLength / 3)); + bezCurve[2] = maths.addArrays(lastPoint, maths.mulItems(rightTangent, segLength / 3)); + } else { + bezCurve[1] = maths.addArrays(firstPoint, maths.mulItems(leftTangent, alpha_l)); + bezCurve[2] = maths.addArrays(lastPoint, maths.mulItems(rightTangent, alpha_r)); + } + return bezCurve; + } + function reparameterize(bezier, points, parameters) { + return parameters.map(function(p, i) { + return newtonRaphsonRootFind(bezier, points[i], p); + }); + } + function newtonRaphsonRootFind(bez, point, u) { + var d = maths.subtract(bezier.q(bez, u), point), qprime = bezier.qprime(bez, u), numerator = maths.mulMatrix(d, qprime), denominator = maths.sum(maths.addItems(maths.squareItems(qprime), maths.mulMatrix(d, bezier.qprimeprime(bez, u)))); + if (denominator === 0) { + return u; + } else { + return u - numerator / denominator; + } + } + function chordLengthParameterize(points) { + var u = [], currU, prevU, prevP; + points.forEach(function(p, i) { + currU = i ? prevU + maths.vectorLen(maths.subtract(p, prevP)) : 0; + u.push(currU); + prevU = currU; + prevP = p; + }); + u = u.map(function(x) { + return x / prevU; + }); + return u; + } + function computeMaxError(points, bez, parameters) { + var dist, maxDist, splitPoint, v, i, count, point, u; + maxDist = 0; + splitPoint = points.length / 2; + for (i = 0, count = points.length;i < count;i++) { + point = points[i]; + u = parameters[i]; + v = maths.subtract(bezier.q(bez, u), point); + dist = v[0] * v[0] + v[1] * v[1]; + if (dist > maxDist) { + maxDist = dist; + splitPoint = i; + } + } + return [maxDist, splitPoint]; + } + return fitCurve; +}); diff --git a/lib/fit-curve.min.js b/lib/fit-curve.min.js new file mode 100644 index 0000000..cc0883a --- /dev/null +++ b/lib/fit-curve.min.js @@ -0,0 +1,24 @@ +/* + JavaScript implementation of + Algorithm for Automatically Fitting Digitized Curves + by Philip J. Schneider + "Graphics Gems", Academic Press, 1990 + + The MIT License (MIT) + + Original (C): + https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c + -> Python: + https://github.com/volkerp/fitCurves + -> CoffeeScript/JavaScript + math.js/lodash: + https://github.com/soswow/fit-curves + -> JavaScript (ES6-ish), no dependencies + https://github.com/Sphinxxxx/fit-curve +*/ +(function(u,m){"function"===typeof define&&define.amd?define([],m):"object"===typeof module&&module.exports?module.exports=m():u.fitCurve=m()})(this,function(){function u(a,e){if(!(a instanceof e))throw new TypeError("Cannot call a class as a function");}function m(a,e,f,c){var h,k,g,d,l;if(2===a.length)return c=b.vectorLen(b.subtract(a[0],a[1]))/3,h=[a[0],b.addArrays(a[0],b.mulItems(e,c)),b.addArrays(a[1],b.mulItems(f,c)),a[1]],[h];k=x(a);h=v(a,k,e,f);d=w(a,h,k);g=d[0];d=d[1];if(gl;l++)if(k=y(h,a,k),h=v(a,k,e,f),d=w(a,h,k),g=d[0],d=d[1],gh&&(h=c,k=g); +return[h,k]}var b=function(){function a(){u(this,a)}a.zeros_Xx2x2=function(e){for(var a=[];e--;)a.push([0,0]);return a};a.mulItems=function(e,a){return[e[0]*a,e[1]*a]};a.mulMatrix=function(e,a){return e[0]*a[0]+e[1]*a[1]};a.subtract=function(e,a){return[e[0]-a[0],e[1]-a[1]]};a.addArrays=function(a,b){return[a[0]+b[0],a[1]+b[1]]};a.addItems=function(a,b){return[a[0]+b,a[1]+b]};a.sum=function(a){return a.reduce(function(a,e){return a+e})};a.dot=function(e,b){return a.mulMatrix(e,b)};a.vectorLen=function(a){var b= +a[0];a=a[1];return Math.sqrt(b*b+a*a)};a.divItems=function(a,b){return[a[0]/b,a[1]/b]};a.squareItems=function(a){var b=a[0];a=a[1];return[b*b,a*a]};a.normalize=function(a){return this.divItems(a,this.vectorLen(a))};return a}(),t=function(){function a(){u(this,a)}a.q=function(a,f){var c=1-f,h=b.mulItems(a[0],c*c*c),k=b.mulItems(a[1],3*c*c*f),c=b.mulItems(a[2],3*c*f*f),g=b.mulItems(a[3],f*f*f);return b.addArrays(b.addArrays(h,k),b.addArrays(c,g))};a.qprime=function(a,f){var c=1-f,h=b.mulItems(b.subtract(a[1], +a[0]),3*c*c),c=b.mulItems(b.subtract(a[2],a[1]),6*c*f),k=b.mulItems(b.subtract(a[3],a[2]),3*f*f);return b.addArrays(b.addArrays(h,c),k)};a.qprimeprime=function(a,f){return b.addArrays(b.mulItems(b.addArrays(b.subtract(a[2],b.mulItems(a[1],2)),a[0]),6*(1-f)),b.mulItems(b.addArrays(b.subtract(a[3],b.mulItems(a[2],2)),a[1]),6*f))};return a}();return function(a,e){var f=a.length,c=b.normalize(b.subtract(a[1],a[0])),f=b.normalize(b.subtract(a[f-2],a[f-1]));return m(a,c,f,e)}}); \ No newline at end of file diff --git a/lib/fitCurves.js b/lib/fitCurves.js deleted file mode 100644 index b2d7627..0000000 --- a/lib/fitCurves.js +++ /dev/null @@ -1,182 +0,0 @@ -// Generated by CoffeeScript 1.9.3 -(function() { - " CoffeeScript implementation of\nAlgorithm for Automatically Fitting Digitized Curves\nby Philip J. Schneider\n\"Graphics Gems\", Academic Press, 1990"; - var add, bezier, chain, chordLengthParameterize, computeMaxError, dot, fitCubic, fitCurve, generateBezier, last, lodash, math, multiply, newtonRaphsonRootFind, normalize, reparameterize, subtract, zeros, zip; - - math = (typeof require === "function" ? require('mathjs') : void 0) || (typeof window !== "undefined" && window !== null ? window.mathjs : void 0); - - lodash = (typeof require === "function" ? require('lodash') : void 0) || (typeof window !== "undefined" && window !== null ? window._ : void 0); - - zeros = math.zeros; - - multiply = math.multiply; - - subtract = math.subtract; - - add = math.add; - - chain = math.chain; - - dot = math.dot; - - last = lodash.last; - - zip = lodash.zip; - - bezier = { - q: function(ctrlPoly, t) { - return math.chain(multiply(Math.pow(1.0 - t, 3), ctrlPoly[0])).add(multiply(3 * Math.pow(1.0 - t, 2) * t, ctrlPoly[1])).add(multiply(3 * (1.0 - t) * Math.pow(t, 2), ctrlPoly[2])).add(multiply(Math.pow(t, 3), ctrlPoly[3])).done(); - }, - qprime: function(ctrlPoly, t) { - return math.chain(multiply(3 * Math.pow(1.0 - t, 2), subtract(ctrlPoly[1], ctrlPoly[0]))).add(multiply(6 * (1.0 - t) * t, subtract(ctrlPoly[2], ctrlPoly[1]))).add(multiply(3 * Math.pow(t, 2), subtract(ctrlPoly[3], ctrlPoly[2]))).done(); - }, - qprimeprime: function(ctrlPoly, t) { - return add(multiply(6 * (1.0 - t), add(subtract(ctrlPoly[2], multiply(2, ctrlPoly[1])), ctrlPoly[0])), multiply(6 * t, add(subtract(ctrlPoly[3], multiply(2, ctrlPoly[2])), ctrlPoly[1]))); - } - }; - - fitCurve = function(points, maxError) { - var leftTangent, rightTangent; - leftTangent = normalize(subtract(points[1], points[0])); - rightTangent = normalize(subtract(points[points.length - 2], last(points))); - return fitCubic(points, leftTangent, rightTangent, maxError); - }; - - fitCubic = function(points, leftTangent, rightTangent, error) { - var bezCurve, beziers, centerTangent, dist, i, j, maxError, ref, ref1, splitPoint, u, uPrime; - if (points.length === 2) { - dist = math.norm(subtract(points[0], points[1])) / 3.0; - bezCurve = [points[0], add(points[0], multiply(leftTangent, dist)), add(points[1], multiply(rightTangent, dist)), points[1]]; - return [bezCurve]; - } - u = chordLengthParameterize(points); - bezCurve = generateBezier(points, u, leftTangent, rightTangent); - ref = computeMaxError(points, bezCurve, u), maxError = ref[0], splitPoint = ref[1]; - if (maxError < error) { - return [bezCurve]; - } - if (maxError < Math.pow(error, 2)) { - for (i = j = 0; j < 20; i = ++j) { - uPrime = reparameterize(bezCurve, points, u); - bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent); - ref1 = computeMaxError(points, bezCurve, uPrime), maxError = ref1[0], splitPoint = ref1[1]; - if (maxError < error) { - return [bezCurve]; - } - u = uPrime; - } - } - beziers = []; - centerTangent = normalize(subtract(points[splitPoint - 1], points[splitPoint + 1])); - beziers = beziers.concat(fitCubic(points.slice(0, splitPoint + 1), leftTangent, centerTangent, error)); - beziers = beziers.concat(fitCubic(points.slice(splitPoint), multiply(centerTangent, -1), rightTangent, error)); - return beziers; - }; - - generateBezier = function(points, parameters, leftTangent, rightTangent) { - var A, C, X, alpha_l, alpha_r, bezCurve, det_C0_C1, det_C0_X, det_X_C1, epsilon, i, j, k, len, len1, point, ref, ref1, segLength, tmp, u; - bezCurve = [points[0], null, null, last(points)]; - A = zeros(parameters.length, 2, 2).valueOf(); - for (i = j = 0, len = parameters.length; j < len; i = ++j) { - u = parameters[i]; - A[i][0] = multiply(leftTangent, 3 * Math.pow(1 - u, 2) * u); - A[i][1] = multiply(rightTangent, 3 * (1 - u) * Math.pow(u, 2)); - } - C = zeros(2, 2).valueOf(); - X = zeros(2).valueOf(); - ref = zip(points, parameters); - for (i = k = 0, len1 = ref.length; k < len1; i = ++k) { - ref1 = ref[i], point = ref1[0], u = ref1[1]; - C[0][0] += dot(A[i][0], A[i][0]); - C[0][1] += dot(A[i][0], A[i][1]); - C[1][0] += dot(A[i][0], A[i][1]); - C[1][1] += dot(A[i][1], A[i][1]); - tmp = subtract(point, bezier.q([points[0], points[0], last(points), last(points)], u)); - X[0] += dot(A[i][0], tmp); - X[1] += dot(A[i][1], tmp); - } - det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; - det_C0_X = C[0][0] * X[1] - C[1][0] * X[0]; - det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; - alpha_l = det_C0_C1 === 0 ? 0 : det_X_C1 / det_C0_C1; - alpha_r = det_C0_C1 === 0 ? 0 : det_C0_X / det_C0_C1; - segLength = math.norm(subtract(points[0], last(points))); - epsilon = 1.0e-6 * segLength; - if (alpha_l < epsilon || alpha_r < epsilon) { - bezCurve[1] = add(bezCurve[0], multiply(leftTangent, segLength / 3.0)); - bezCurve[2] = add(bezCurve[3], multiply(rightTangent, segLength / 3.0)); - } else { - bezCurve[1] = add(bezCurve[0], multiply(leftTangent, alpha_l)); - bezCurve[2] = add(bezCurve[3], multiply(rightTangent, alpha_r)); - } - return bezCurve; - }; - - reparameterize = function(bezier, points, parameters) { - var j, len, point, ref, ref1, results, u; - ref = zip(points, parameters); - results = []; - for (j = 0, len = ref.length; j < len; j++) { - ref1 = ref[j], point = ref1[0], u = ref1[1]; - results.push(newtonRaphsonRootFind(bezier, point, u)); - } - return results; - }; - - newtonRaphsonRootFind = function(bez, point, u) { - "Newton's root finding algorithm calculates f(x)=0 by reiterating\nx_n+1 = x_n - f(x_n)/f'(x_n)\nWe are trying to find curve parameter u for some point p that minimizes\nthe distance from that point to the curve. Distance point to curve is d=q(u)-p.\nAt minimum distance the point is perpendicular to the curve.\nWe are solving\nf = q(u)-p * q'(u) = 0\nwith\nf' = q'(u) * q'(u) + q(u)-p * q''(u)\ngives\nu_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)|"; - var d, denominator, numerator, qprime; - d = subtract(bezier.q(bez, u), point); - qprime = bezier.qprime(bez, u); - numerator = math.sum(multiply(d, qprime)); - denominator = math.sum(add(math.dotPow(qprime, 2), multiply(d, bezier.qprimeprime(bez, u)))); - if (denominator === 0) { - return u; - } else { - return u - numerator / denominator; - } - }; - - chordLengthParameterize = function(points) { - var i, j, k, ref, ref1, u; - u = [0]; - for (i = j = 1, ref = points.length; 1 <= ref ? j < ref : j > ref; i = 1 <= ref ? ++j : --j) { - u.push(u[i - 1] + math.norm(subtract(points[i], points[i - 1]))); - } - for (i = k = 0, ref1 = u.length; 0 <= ref1 ? k < ref1 : k > ref1; i = 0 <= ref1 ? ++k : --k) { - u[i] = u[i] / last(u); - } - return u; - }; - - computeMaxError = function(points, bez, parameters) { - var dist, i, j, len, maxDist, point, ref, ref1, splitPoint, u; - maxDist = 0; - splitPoint = points.length / 2; - ref = zip(points, parameters); - for (i = j = 0, len = ref.length; j < len; i = ++j) { - ref1 = ref[i], point = ref1[0], u = ref1[1]; - dist = Math.pow(math.norm(subtract(bezier.q(bez, u), point)), 2); - if (dist > maxDist) { - maxDist = dist; - splitPoint = i; - } - } - return [maxDist, splitPoint]; - }; - - normalize = function(v) { - return math.divide(v, math.norm(v)); - }; - - if (typeof exports !== "undefined" && exports !== null) { - if ((typeof module !== "undefined" && module !== null) && (module.exports != null)) { - module.exports = fitCurve; - } else { - exports.fitCurve = fitCurve; - } - } else if (window) { - window.fitCurve = fitCurve; - } - -}).call(this); diff --git a/src/fit-curve.core.js b/src/fit-curve.core.js new file mode 100644 index 0000000..4430bce --- /dev/null +++ b/src/fit-curve.core.js @@ -0,0 +1,391 @@ +/* + Simplified versions of what we need from math.js + Optimized for our input, which is only numbers and 1x2 arrays (i.e. [x, y] coordinates). +*/ +class maths { + //zeros = logAndRun(math.zeros); + static zeros_Xx2x2(x) { + var zs = []; + while(x--) { zs.push([0,0]); } + return zs + } + + //multiply = logAndRun(math.multiply); + static mulItems(items, multiplier) { + //return items.map(x => x*multiplier); + return [items[0]*multiplier, items[1]*multiplier]; + } + static mulMatrix(m1, m2) { + //https://en.wikipedia.org/wiki/Matrix_multiplication#Matrix_product_.28two_matrices.29 + //Simplified to only handle 1-dimensional matrices (i.e. arrays) of equal length: + // return m1.reduce((sum,x1,i) => sum + (x1*m2[i]), + // 0); + return (m1[0]*m2[0]) + (m1[1]*m2[1]); + } + + //Only used to subract to points (or at least arrays): + // subtract = logAndRun(math.subtract); + static subtract(arr1, arr2) { + //return arr1.map((x1, i) => x1 - arr2[i]); + return [arr1[0]-arr2[0], arr1[1]-arr2[1]]; + } + + //add = logAndRun(math.add); + static addArrays(arr1, arr2) { + //return arr1.map((x1, i) => x1 + arr2[i]); + return [arr1[0]+arr2[0], arr1[1]+arr2[1]]; + } + static addItems(items, addition) { + //return items.map(x => x+addition); + return [items[0]+addition, items[1]+addition]; + } + + //var sum = logAndRun(math.sum); + static sum(items) { + return items.reduce((sum,x) => sum + x); + } + + //chain = math.chain; + + //Only used on two arrays. The dot product is equal to the matrix product in this case: + // dot = logAndRun(math.dot); + static dot(m1, m2) { + return maths.mulMatrix(m1, m2); + } + + //https://en.wikipedia.org/wiki/Norm_(mathematics)#Euclidean_norm + // var norm = logAndRun(math.norm); + static vectorLen(v) { + var a = v[0], b = v[1]; + return Math.sqrt(a*a + b*b); + } + + //math.divide = logAndRun(math.divide); + static divItems(items, divisor) { + //return items.map(x => x/divisor); + return [items[0]/divisor, items[1]/divisor]; + } + + //var dotPow = logAndRun(math.dotPow); + static squareItems(items) { + //return items.map(x => x*x); + var a = items[0], b = items[1]; + return [a*a, b*b]; + } + + static normalize(v) { + return this.divItems(v, this.vectorLen(v)); + } + + //Math.pow = logAndRun(Math.pow); +} + + +class bezier { + //Evaluates cubic bezier at t, return point + static q(ctrlPoly, t) { + var tx = 1.0 - t; + var pA = maths.mulItems( ctrlPoly[0], tx * tx * tx ), + pB = maths.mulItems( ctrlPoly[1], 3 * tx * tx * t ), + pC = maths.mulItems( ctrlPoly[2], 3 * tx * t * t ), + pD = maths.mulItems( ctrlPoly[3], t * t * t ); + return maths.addArrays(maths.addArrays(pA, pB), maths.addArrays(pC, pD)); + } + + //Evaluates cubic bezier first derivative at t, return point + static qprime(ctrlPoly, t) { + var tx = 1.0 - t; + var pA = maths.mulItems( maths.subtract(ctrlPoly[1], ctrlPoly[0]), 3 * tx * tx ), + pB = maths.mulItems( maths.subtract(ctrlPoly[2], ctrlPoly[1]), 6 * tx * t ), + pC = maths.mulItems( maths.subtract(ctrlPoly[3], ctrlPoly[2]), 3 * t * t ); + return maths.addArrays(maths.addArrays(pA, pB), pC); + } + + //Evaluates cubic bezier second derivative at t, return point + static qprimeprime(ctrlPoly, t) { + return maths.addArrays(maths.mulItems( maths.addArrays(maths.subtract(ctrlPoly[2], maths.mulItems(ctrlPoly[1], 2)), ctrlPoly[0]), 6 * (1.0 - t) ), + maths.mulItems( maths.addArrays(maths.subtract(ctrlPoly[3], maths.mulItems(ctrlPoly[2], 2)), ctrlPoly[1]), 6 * t )); + } +} + + +/** + * Fit one or more Bezier curves to a set of points. + * + * @param {Array>} points - Array of digitized points, e.g. [[5,5],[5,50],[110,140],[210,160],[320,110]] + * @param {Number} maxError - Tolerance, squared error between points and fitted curve + * @returns {Array>>} Array of Bezier curves, where each element is [first-point, control-point-1, control-point-2, second-point] and points are [x, y] + */ +function fitCurve(points, maxError) { + var len = points.length, + leftTangent = maths.normalize(maths.subtract(points[1], points[0])), + rightTangent = maths.normalize(maths.subtract(points[len - 2], points[len - 1])); + return fitCubic(points, leftTangent, rightTangent, maxError); +} + +/** + * Fit a Bezier curve to a (sub)set of digitized points. + * Your code should not call this function directly. Use {@link fitCurve} instead. + * + * @param {Array>} points - Array of digitized points, e.g. [[5,5],[5,50],[110,140],[210,160],[320,110]] + * @param {Array} leftTangent - Unit tangent vector at start point + * @param {Array} rightTangent - Unit tangent vector at end point + * @param {Number} error - Tolerance, squared error between points and fitted curve + * @returns {Array>>} Array of Bezier curves, where each element is [first-point, control-point-1, control-point-2, second-point] and points are [x, y] + */ +function fitCubic(points, leftTangent, rightTangent, error) { + const MaxIterations = 20; //Max times to try iterating (to find an acceptable curve) + + var bezCurve, //Control points of fitted Bezier curve + u, //Parameter values for point + uPrime, //Improved parameter values + maxError, //Maximum fitting error + splitPoint, //Point to split point set at if we need more than one curve + centerTangent, //Unit tangent vector at splitPoint + beziers, //Array of fitted Bezier curves if we need more than one curve + dist, i; + + //console.log('fitCubic, ', points.length); + + //Use heuristic if region only has two points in it + if (points.length === 2) { + dist = maths.vectorLen(maths.subtract(points[0], points[1])) / 3.0; + bezCurve = [ + points[0], + maths.addArrays(points[0], maths.mulItems(leftTangent, dist)), + maths.addArrays(points[1], maths.mulItems(rightTangent, dist)), + points[1] + ]; + return [bezCurve]; + } + + //Parameterize points, and attempt to fit curve + u = chordLengthParameterize(points); + bezCurve = generateBezier(points, u, leftTangent, rightTangent); + + //Find max deviation of points to fitted curve + [maxError, splitPoint] = computeMaxError(points, bezCurve, u); + if (maxError < error) { + //console.log('cme ~', maxError, points[splitPoint]); + return [bezCurve]; + } + //If error not too large, try some reparameterization and iteration + if (maxError < (error*error)) { + for (i = 0; i < MaxIterations; i++) { + uPrime = reparameterize(bezCurve, points, u); + bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent); + [maxError, splitPoint] = computeMaxError(points, bezCurve, uPrime); + if (maxError < error) { + //console.log('cme '+i, maxError, points[splitPoint]); + return [bezCurve]; + } + u = uPrime; + } + } + + //Fitting failed -- split at max error point and fit recursively + //console.log('splitting'); + beziers = []; + centerTangent = maths.normalize(maths.subtract(points[splitPoint - 1], points[splitPoint + 1])); + beziers = beziers.concat(fitCubic(points.slice(0, splitPoint + 1), leftTangent, centerTangent, error)); + beziers = beziers.concat(fitCubic(points.slice(splitPoint), maths.mulItems(centerTangent, -1), rightTangent, error)); + return beziers; +}; + +/** + * Use least-squares method to find Bezier control points for region. + * + * @param {Array>} points - Array of digitized points + * @param {Array} parameters - Parameter values for region + * @param {Array} leftTangent - Unit tangent vector at start point + * @param {Array} rightTangent - Unit tangent vector at end point + * @returns {Array>} Approximated Bezier curve: [first-point, control-point-1, control-point-2, second-point] where points are [x, y] + */ +function generateBezier(points, parameters, leftTangent, rightTangent) { + var bezCurve, //Bezier curve ctl pts + A, a, //Precomputed rhs for eqn + C, X, //Matrices C & X + det_C0_C1, det_C0_X, det_X_C1, //Determinants of matrices + alpha_l, alpha_r, //Alpha values, left and right + + epsilon, segLength, + i, len, tmp, u, ux, + firstPoint = points[0], + lastPoint = points[points.length-1]; + + bezCurve = [firstPoint, null, null, lastPoint]; + //console.log('gb', parameters.length); + + //Compute the A's + A = maths.zeros_Xx2x2(parameters.length); + for (i = 0, len = parameters.length; i < len; i++) { + u = parameters[i]; + ux = 1 - u; + a = A[i]; + + a[0] = maths.mulItems(leftTangent, 3 * u * (ux*ux)); + a[1] = maths.mulItems(rightTangent, 3 * ux * (u*u)); + } + + //Create the C and X matrices + C = [[0,0], [0,0]]; + X = [0,0]; + for (i = 0, len = points.length; i < len; i++) { + u = parameters[i]; + a = A[i]; + + C[0][0] += maths.dot(a[0], a[0]); + C[0][1] += maths.dot(a[0], a[1]); + C[1][0] += maths.dot(a[0], a[1]); + C[1][1] += maths.dot(a[1], a[1]); + + tmp = maths.subtract(points[i], bezier.q([firstPoint, firstPoint, lastPoint, lastPoint], u)); + + X[0] += maths.dot(a[0], tmp); + X[1] += maths.dot(a[1], tmp); + } + + //Compute the determinants of C and X + det_C0_C1 = (C[0][0] * C[1][1]) - (C[1][0] * C[0][1]); + det_C0_X = (C[0][0] * X[1] ) - (C[1][0] * X[0] ); + det_X_C1 = (X[0] * C[1][1]) - (X[1] * C[0][1]); + + //Finally, derive alpha values + alpha_l = det_C0_C1 === 0 ? 0 : det_X_C1 / det_C0_C1; + alpha_r = det_C0_C1 === 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. + segLength = maths.vectorLen(maths.subtract(firstPoint, lastPoint)); + epsilon = 1.0e-6 * segLength; + if (alpha_l < epsilon || alpha_r < epsilon) { + //Fall back on standard (probably inaccurate) formula, and subdivide further if needed. + bezCurve[1] = maths.addArrays(firstPoint, maths.mulItems(leftTangent, segLength / 3.0)); + bezCurve[2] = maths.addArrays(lastPoint, maths.mulItems(rightTangent, segLength / 3.0)); + } else { + //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[1] = maths.addArrays(firstPoint, maths.mulItems(leftTangent, alpha_l)); + bezCurve[2] = maths.addArrays(lastPoint, maths.mulItems(rightTangent, alpha_r)); + } + + return bezCurve; +}; + +/** + * Given set of points and their parameterization, try to find a better parameterization. + * + * @param {Array>} bezier - Current fitted curve + * @param {Array>} points - Array of digitized points + * @param {Array} parameters - Current parameter values + * @returns {Array} New parameter values + */ +function reparameterize(bezier, points, parameters) { + /* + var j, len, point, results, u; + results = []; + for (j = 0, len = points.length; j < len; j++) { + point = points[j], u = parameters[j]; + + results.push(newtonRaphsonRootFind(bezier, point, u)); + } + return results; + //*/ + return parameters.map((p, i) => newtonRaphsonRootFind(bezier, points[i], p)); +}; + +/** + * Use Newton-Raphson iteration to find better root. + * + * @param {Array>} bez - Current fitted curve + * @param {Array} point - Digitized point + * @param {Number} u - Parameter value for "P" + * @returns {Number} New u + */ +function newtonRaphsonRootFind(bez, point, u) { + /* + Newton's root finding algorithm calculates f(x)=0 by reiterating + x_n+1 = x_n - f(x_n)/f'(x_n) + We are trying to find curve parameter u for some point p that minimizes + the distance from that point to the curve. Distance point to curve is d=q(u)-p. + At minimum distance the point is perpendicular to the curve. + We are solving + f = q(u)-p * q'(u) = 0 + with + f' = q'(u) * q'(u) + q(u)-p * q''(u) + gives + u_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)| + */ + + var d = maths.subtract(bezier.q(bez, u), point), + qprime = bezier.qprime(bez, u), + numerator = /*sum(*/maths.mulMatrix(d, qprime)/*)*/, + denominator = maths.sum(maths.addItems( maths.squareItems(qprime), maths.mulMatrix(d, bezier.qprimeprime(bez, u)) )); + + if (denominator === 0) { + return u; + } else { + return u - (numerator/denominator); + } +}; + +/** + * Assign parameter values to digitized points using relative distances between points. + * + * @param {Array>} points - Array of digitized points + * @returns {Array} Parameter values + */ +function chordLengthParameterize(points) { + var u = [], currU, prevU, prevP; + + points.forEach((p, i) => { + currU = i ? prevU + maths.vectorLen(maths.subtract(p, prevP)) + : 0; + u.push(currU); + + prevU = currU; + prevP = p; + }) + u = u.map(x => x/prevU); + + return u; +}; + +/** + * Find the maximum squared distance of digitized points to fitted curve. + * + * @param {Array>} points - Array of digitized points + * @param {Array>} bez - Fitted curve + * @param {Array} parameters - Parameterization of points + * @returns {Array} Maximum error (squared) and point of max error + */ +function computeMaxError(points, bez, parameters) { + var dist, //Current error + maxDist, //Maximum error + splitPoint, //Point of maximum error + v, //Vector from point to curve + i, count, point, u; + + maxDist = 0; + splitPoint = points.length / 2; + + for (i = 0, count = points.length; i < count; i++) { + point = points[i]; + u = parameters[i]; + + //len = maths.vectorLen(maths.subtract(bezier.q(bez, u), point)); + //dist = len * len; + v = maths.subtract(bezier.q(bez, u), point); + dist = v[0]*v[0] + v[1]*v[1]; + + if (dist > maxDist) { + maxDist = dist; + splitPoint = i; + } + } + + return [maxDist, splitPoint]; +}; diff --git a/src/fit-curve.umd-wrapper.js b/src/fit-curve.umd-wrapper.js new file mode 100644 index 0000000..f032ca0 --- /dev/null +++ b/src/fit-curve.umd-wrapper.js @@ -0,0 +1,52 @@ +// ==ClosureCompiler== +// @output_file_name fit-curve.min.js +// @compilation_level SIMPLE_OPTIMIZATIONS +// ==/ClosureCompiler== + +/** + * @preserve JavaScript implementation of + * Algorithm for Automatically Fitting Digitized Curves + * by Philip J. Schneider + * "Graphics Gems", Academic Press, 1990 + * + * The MIT License (MIT) + * + * Original (C): + * https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c + * -> Python: + * https://github.com/volkerp/fitCurves + * -> CoffeeScript/JavaScript + math.js/lodash: + * https://github.com/soswow/fit-curves + * -> JavaScript (ES6-ish), no dependencies + * https://github.com/Sphinxxxx/fit-curve + */ + +//UMD pattern from +//https://github.com/umdjs/umd/blob/master/templates/returnExports.js +(function (root, factory) { + if (typeof define === 'function' && define.amd) { + // AMD. Register as an anonymous module. + define([], factory); + } else if (typeof module === 'object' && module.exports) { + // Node. Does not work with strict CommonJS, + // only CommonJS-like environments that support module.exports, like Node. + module.exports = factory(); + } else { + // Browser globals (root is window) + root.fitCurve = factory(); + } +})(this, function () { + + + /* + Insert fit-curve.core.js here. + * Transpile first if necessary, for examle here: + https://babeljs.io/repl/ (Presets: es2015, es2015-loose, stage-2) + + Then minify with Google's Closure Compiler (settings included at the top of this file): + https://closure-compiler.appspot.com/home + */ + + + return fitCurve; +}); diff --git a/src/fitCurves.coffee b/src/fitCurves.coffee deleted file mode 100644 index 3a7d379..0000000 --- a/src/fitCurves.coffee +++ /dev/null @@ -1,220 +0,0 @@ -""" CoffeeScript implementation of - Algorithm for Automatically Fitting Digitized Curves - by Philip J. Schneider - "Graphics Gems", Academic Press, 1990 -""" - -math = require?('mathjs') || window?.mathjs -lodash = require?('lodash') || window?._ - -zeros = math.zeros -multiply = math.multiply -subtract = math.subtract -add = math.add -chain = math.chain -dot = math.dot - -last = lodash.last -zip = lodash.zip - -bezier = - # evaluates cubic bezier at t, return point - q: (ctrlPoly, t) -> - math.chain(multiply((1.0-t)**3, ctrlPoly[0])) - .add(multiply(3*(1.0-t)**2 * t, ctrlPoly[1])) - .add(multiply(3*(1.0-t) * t**2, ctrlPoly[2])) - .add(multiply(t**3, ctrlPoly[3])) - .done() - - # evaluates cubic bezier first derivative at t, return point - qprime: (ctrlPoly, t) -> - math.chain(multiply(3*(1.0-t)**2, subtract(ctrlPoly[1],ctrlPoly[0]))) - .add(multiply(6*(1.0-t) * t, subtract(ctrlPoly[2], ctrlPoly[1]))) - .add(multiply(3*t**2, subtract(ctrlPoly[3],ctrlPoly[2]))) - .done() - - # evaluates cubic bezier second derivative at t, return point - qprimeprime: (ctrlPoly, t) -> - add( - multiply( - 6*(1.0-t), - add( - subtract( - ctrlPoly[2], - multiply(2, ctrlPoly[1]) - ), - ctrlPoly[0] - ) - ), - multiply( - 6*(t), - add( - subtract( - ctrlPoly[3], - multiply(2,ctrlPoly[2]) - ), - ctrlPoly[1] - ) - ) - ) - -# Fit one (ore more) Bezier curves to a set of points -fitCurve = (points, maxError) -> - leftTangent = normalize(subtract(points[1], points[0])) - rightTangent = normalize(subtract(points[points.length - 2], last(points))) - fitCubic(points, leftTangent, rightTangent, maxError) - - -fitCubic = (points, leftTangent, rightTangent, error) -> - # Use heuristic if region only has two points in it - if points.length is 2 - dist = math.norm(subtract(points[0], points[1])) / 3.0 - bezCurve = [ - points[0], - add(points[0], multiply(leftTangent, dist)), - add(points[1], multiply(rightTangent, dist)), - points[1] - ] - return [bezCurve] - - # Parameterize points, and attempt to fit curve - u = chordLengthParameterize(points) - bezCurve = generateBezier(points, u, leftTangent, rightTangent) - # Find max deviation of points to fitted curve - [maxError, splitPoint] = computeMaxError(points, bezCurve, u) - if maxError < error - return [bezCurve] - - # If error not too large, try some reparameterization and iteration - if maxError < error**2 - for i in [0...20] - uPrime = reparameterize(bezCurve, points, u) - bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent) - [maxError, splitPoint] = computeMaxError(points, bezCurve, uPrime) - if maxError < error - return [bezCurve] - u = uPrime - - # Fitting failed -- split at max error point and fit recursively - beziers = [] - centerTangent = normalize subtract points[splitPoint - 1], points[splitPoint + 1] - beziers = beziers.concat fitCubic points[...splitPoint + 1], leftTangent, centerTangent, error - beziers = beziers.concat fitCubic points[splitPoint...], multiply(centerTangent, -1), rightTangent, error - - return beziers - - -generateBezier = (points, parameters, leftTangent, rightTangent) -> - bezCurve = [points[0], null, null, last(points)] - - # compute the A's - A = zeros(parameters.length, 2, 2).valueOf() - - for u, i in parameters - A[i][0] = multiply(leftTangent, 3 * (1 - u)**2 * u) - A[i][1] = multiply(rightTangent, 3 * (1 - u) * u**2) - - # Create the C and X matrices - C = zeros(2, 2).valueOf() - X = zeros(2).valueOf() - - for [point, u], i in zip(points, parameters) - C[0][0] += dot(A[i][0], A[i][0]) - C[0][1] += dot(A[i][0], A[i][1]) - C[1][0] += dot(A[i][0], A[i][1]) - C[1][1] += dot(A[i][1], A[i][1]) - - tmp = subtract(point, bezier.q([points[0], points[0], last(points), last(points)], u)) - - X[0] += dot(A[i][0], tmp) - X[1] += dot(A[i][1], tmp) - - # Compute the determinants of C and X - det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1] - det_C0_X = C[0][0] * X[1] - C[1][0] * X[0] - det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1] - - # Finally, derive alpha values - alpha_l = if det_C0_C1 is 0 then 0 else det_X_C1 / det_C0_C1 - alpha_r = if det_C0_C1 is 0 then 0 else 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. */ - segLength = math.norm(subtract(points[0], last(points))) - epsilon = 1.0e-6 * segLength - if alpha_l < epsilon or alpha_r < epsilon - # fall back on standard (probably inaccurate) formula, and subdivide further if needed. - bezCurve[1] = add(bezCurve[0], multiply(leftTangent, segLength / 3.0)) - bezCurve[2] = add(bezCurve[3], multiply(rightTangent, segLength / 3.0)) - else - # 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[1] = add(bezCurve[0], multiply(leftTangent, alpha_l)) - bezCurve[2] = add(bezCurve[3], multiply(rightTangent, alpha_r)) - - return bezCurve - - -reparameterize = (bezier, points, parameters) -> - (newtonRaphsonRootFind(bezier, point, u) for [point, u] in zip(points, parameters)) - - -newtonRaphsonRootFind = (bez, point, u) -> - """ - Newton's root finding algorithm calculates f(x)=0 by reiterating - x_n+1 = x_n - f(x_n)/f'(x_n) - We are trying to find curve parameter u for some point p that minimizes - the distance from that point to the curve. Distance point to curve is d=q(u)-p. - At minimum distance the point is perpendicular to the curve. - We are solving - f = q(u)-p * q'(u) = 0 - with - f' = q'(u) * q'(u) + q(u)-p * q''(u) - gives - u_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)| - """ - d = subtract(bezier.q(bez, u), point) - qprime = bezier.qprime(bez, u) - numerator = math.sum(multiply(d, qprime)) - denominator = math.sum(add(math.dotPow(qprime, 2), multiply(d, bezier.qprimeprime(bez, u)))) - - if denominator is 0 - u - else - u - numerator / denominator - - -chordLengthParameterize = (points) -> - u = [0] - for i in [1...points.length] - u.push(u[i - 1] + math.norm(subtract(points[i], points[i - 1]))) - - for i in [0...u.length] - u[i] = u[i] / last(u) - - return u - - -computeMaxError = (points, bez, parameters) -> - maxDist = 0 - splitPoint = points.length / 2 - for [point, u], i in zip(points, parameters) - dist = math.norm(subtract(bezier.q(bez, u), point))**2 - if dist > maxDist - maxDist = dist - splitPoint = i - - [maxDist, splitPoint] - -normalize = (v) -> math.divide(v, math.norm(v)) - -if exports? - if module? and module.exports? - module.exports = fitCurve - else - exports.fitCurve = fitCurve -else if window - window.fitCurve = fitCurve diff --git a/src/log-and-run.js b/src/log-and-run.js new file mode 100644 index 0000000..83a3e44 --- /dev/null +++ b/src/log-and-run.js @@ -0,0 +1,42 @@ +/* + Only a utility used for investigating existing code. + Not needed to build fit-curve.js +*/ +function logAndRun(func, log) { + //return func; + + if(!log) { + (this || self)._loggedAndRun = log = {}; + } + + function argType(arg) { + if(Array.isArray(arg)) { + return '['+ arg.map(argType).join(',') +']'; + } + else if(Number.isFinite(arg)) { + return 'num'; + } + else { + return '???' + arg; + } + } + + return function() { + //http://www.codeovertones.com/2011/08/how-to-print-stack-trace-anywhere-in.html + var e = new Error('dummy'); + var stack = e.stack.replace(/^[^\(]+?[\n$]/gm, '') + .replace(/^\s+at\s+/gm, '') + .replace(/^Object.\s*\(/gm, '{anonymous}()@') + .split('\n'); + + var result = func.apply(null, arguments); + + var signature = func.name + ': ' + Array.from(arguments).map(argType).join(', ') + ' ('+ stack[0] +')'; + var context = func.name + ': ' + Array.from(arguments).join(' ') + ' => ' + result; + + log[signature] = log[signature] || []; + log[signature].push(context); + + return result; + } +}