From 29e57cc52169a3cb110aa514f1472940c046b626 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrg=20Lehni?= Date: Sun, 6 Mar 2011 23:24:33 +0000 Subject: [PATCH] Replace slow simpson() method with insanely fast Gauss-Legendre Numerical Integration by Jim Armstrong which was further optimised. --- src/path/Curve.js | 63 +++++++---- src/util/MathUtils.js | 251 +++++++----------------------------------- 2 files changed, 81 insertions(+), 233 deletions(-) diff --git a/src/path/Curve.js b/src/path/Curve.js index 4be5a6a1..8a4e0fd1 100644 --- a/src/path/Curve.js +++ b/src/path/Curve.js @@ -125,22 +125,21 @@ var Curve = this.Curve = Base.extend({ || this._path.closed && curves[curves.length - 1]) || null; }, - // Calculates arclength of a cubic using adaptive simpson integration. - getLength: function(goal) { + getLength: function() { var z0 = this._segment1._point, z1 = this._segment2._point, c0 = z0.add(this._segment1._handleOut), c1 = z1.add(this._segment2._handleIn); // TODO: Check for straight lines and handle separately. - // Calculate the coefficients of a Bezier derivative, divided by 3. - var ax = 3 * (c0.x - c1.x) - z0.x + z1.x, - bx = 2 * (z0.x + c1.x) - 4 * c0.x, - cx = c0.x - z0.x, + // Calculate the coefficients of a Bezier derivative. + var ax = 9 * (c0.x - c1.x) + 3 * (z1.x - z0.x), + bx = 6 * (z0.x + c1.x) - 12 * c0.x, + cx = 3 * (c0.x - z0.x), - ay = 3 * (c0.y - c1.y) - z0.y + z1.y, - by = 2 * (z0.y + c1.y) - 4 * c0.y, - cy = c0.y - z0.y; + ay = 9 * (c0.y - c1.y) + 3 * (z1.y - z0.y), + by = 6 * (z0.y + c1.y) - 12 * c0.y, + cy = 3 * (c0.y - z0.y); function ds(t) { // Calculate quadratic equations of derivatives for x and y @@ -149,20 +148,42 @@ var Curve = this.Curve = Base.extend({ return Math.sqrt(dx * dx + dy * dy); } - var integral = MathUtils.simpson(ds, 0.0, 1.0, MathUtils.EPSILON, 1.0); - if (integral == null) - throw new Error('Nesting capacity exceeded in Path#getLenght()'); - // Multiply by 3 again, as derivative was divided by 3 - var length = 3 * integral; - if (goal == undefined || goal < 0 || goal >= length) - return length; - var result = MathUtils.unsimpson(goal, ds, 0, goal / integral, - 100 * MathUtils.EPSILON, integral, Math.sqrt(MathUtils.EPSILON), 1); - if (!result) - throw new Error('Nesting capacity exceeded in computing arctime'); - return -result.b; + return MathUtils.gauss(ds, 0.0, 1.0, 8); }, + /** + * Checks if this curve is linear, meaning it does not define any curve + * handle. + + * @return {@true if the curve is linear} + */ + isLinear: function() { + return this._segment1._handleOut.isZero() + && this._segment2._handleIn.isZero(); + }, + + // TODO: getParameter(length) + // TODO: getParameter(point, precision) + // TODO: getLocation + // TODO: getIntersections + // TODO: adjustThroughPoint + + transform: function(matrix) { + return new Curve( + matrix.transform(this._segment1._point), + matrix.transform(this._segment1._handleOut), + matrix.transform(this._segment2._handleIn), + matrix.transform(this._segment2._point)); + }, + + reverse: function() { + return new Curve(this._segment2.reverse(), this._segment1.reverse()); + }, + + // TODO: divide + // TODO: split + // TODO: getPartLength(fromParameter, toParameter) + clone: function() { return new Curve(this._segment1, this._segment2); }, diff --git a/src/util/MathUtils.js b/src/util/MathUtils.js index e2721232..57277cf2 100644 --- a/src/util/MathUtils.js +++ b/src/util/MathUtils.js @@ -1,220 +1,47 @@ -var MathUtils = { - // TODO: 64 bit: 53 digits, 32 bit: 24 digits - // TODO: Naming? What is MANT standing for? - MANT_DIGITS: 53, - EPSILON: Math.pow(2, -52), +var MathUtils = new function() { + var maxDepth = 53; // MANT_DIGITS: 64 bit: 53 digits, 32 bit: 24 digits - // Compute a numerical approximation to an integral via adaptive Simpson's - // Rule This routine ignores underflow. + // Gauss-Legendre Numerical Integration, a Gauss.as port from Singularity: + // + // Copyright (c) 2006-2007, Jim Armstrong (www.algorithmist.net) + // All Rights Reserved. + // + // Simplified and further optimised by Juerg Lehni. - // Returns approximate value of the integral if successful. - // f: Pointer to function to be integrated. - // a, b: Lower, upper limits of integration. - // accuracy: Desired relative accuracy of integral, - // try to make |error| <= accuracy*abs(integral). - // dxmax: Maximum limit on the width of a subinterval - // For periodic functions, dxmax should be - // set to the period or smaller to prevent - // premature convergence of Simpson's rule. - simpson: function(f, a, b, accuracy, dxmax) { - var table = new Array(MathUtils.MANT_DIGITS); - var p = table[0] = { - left: true, - psum: 0 - }; - var index = 0, - success = true, - alpha = a, - da = b - a, - fv0 = f(alpha), - fv1, - fv2 = f(alpha + 0.5 * da), - fv3, - fv4 = f(alpha + da), - fv5, - wt = da * (1 / 6), - est = wt * (fv0 + 4 * fv2 + fv4), - area = est; + var abscissa = [ + -0.5773502692, 0.5773502692, + -0.7745966692, 0.7745966692, 0, + -0.8611363116, 0.8611363116, -0.3399810436, 0.3399810436, + -0.9061798459, 0.9061798459, -0.5384693101, 0.5384693101, 0.0000000000, + -0.9324695142, 0.9324695142, -0.6612093865, 0.6612093865, -0.2386191861, 0.2386191861, + -0.9491079123, 0.9491079123, -0.7415311856, 0.7415311856, -0.4058451514, 0.4058451514, 0.0000000000, + -0.9602898565, 0.9602898565, -0.7966664774, 0.7966664774, -0.5255324099, 0.5255324099, -0.1834346425, 0.1834346425 + ]; - // Have estimate est of integral on (alpha, alpha+da). - // Bisect and compute estimates on left and right half intervals. - // integral is the best value for the integral. - for (;;) { - var dx = 0.5 * da, - arg = alpha + 0.5 * dx; - fv1 = f(arg); - fv3 = f(arg + dx); - var wt = dx * (1 / 6), - estl = wt * (fv0 + 4 * fv1 + fv2), - estr = wt * (fv2 + 4 * fv3 + fv4), - integral = estl + estr, - diff = est - integral; - area -= diff; + var weight = [ + 1, 1, + 0.5555555556, 0.5555555556, 0.8888888888, + 0.3478548451, 0.3478548451, 0.6521451549, 0.6521451549, + 0.2369268851, 0.2369268851, 0.4786286705, 0.4786286705, 0.5688888888, + 0.1713244924, 0.1713244924, 0.3607615730, 0.3607615730, 0.4679139346, 0.4679139346, + 0.1294849662, 0.1294849662, 0.2797053915, 0.2797053915, 0.3818300505, 0.3818300505, 0.4179591837, + 0.1012285363, 0.1012285363, 0.2223810345, 0.2223810345, 0.3137066459, 0.3137066459, 0.3626837834, 0.3626837834 + ]; - if (index >= table.length) - success = false; - if (!success || (Math.abs(diff) <= accuracy * Math.abs(area) - && da <= dxmax)) { - // Accept approximate integral. - // If it was a right interval, add results to finish at this - // level. If it was a left interval, process right interval. - for (;;) { - if (!p.left) { // process right-half interval - alpha += da; - p.left = true; - p.psum = integral; - fv0 = p.f1t; - fv2 = p.f2t; - fv4 = p.f3t; - da = p.dat; - est = p.estr; - break; - } - integral += p.psum; - if (--index <= 0) - return success ? integral : null; - p = table[index]; - } + return { + EPSILON: Math.pow(2, -maxDepth + 1), - } else { - // Raise level and store information for processing right-half - // interval. - p = table[++index] = { - left: false, - f1t: fv2, - f2t: fv3, - f3t: fv4, - dat: dx, - estr: estr - }; - da = dx; - est = estl; - fv4 = fv2; - fv2 = fv1; - } + gauss: function(f, a, b, n) { + n = Math.min(Math.max(n, 2), 8); + + var l = n == 2 ? 0 : n * (n - 1) / 2 - 1, + sum = 0, + mul = 0.5 * (b - a), + ab2 = mul + a; + for(var i = 0; i < n; i++) + sum += f(ab2 + mul * abscissa[l + i]) * weight[l + i]; + + return mul * sum; } - }, - - // Use adaptive Simpson integration to determine the upper limit of - // integration required to make the definite integral of a continuous - // non-negative function close to a user specified sum. - // This routine ignores underflow. - // integral: Given value for the integral. - // f: Pointer to function to be integrated. - // a, b: Lower, upper limits of integration (a <= b). - // The value of b provided on entry is used - // as an initial guess; somewhat faster if the - // given value is an underestimation. - // accuracy: Desired relative accuracy of b. - // Try to make |integral-area| <= accuracy*integral. - // area: Computed integral of f(x) on [a,b]. - // dxmin: Lower limit on sampling width. - // dxmax: Maximum limit on the width of a subinterval - // For periodic functions, dxmax should be - // set to the period or smaller to prevent - // premature convergence of Simpson's rule. - unsimpson: function(integral, f, a, b, accuracy, area, dxmin, dxmax) { - var table = new Array(MathUtils.MANT_DIGITS); - var p = table[0] = { - psum: 0 - }; - var index = 0, - alpha = a, - parea = 0, - pdiff = 0, - fv0, fv1, fv2, fv3, fv4, fv5; - for (;;) { - p.left = true; - var da = b - alpha; - fv0 = f(alpha); - fv2 = f(alpha + 0.5 * da); - fv4 = f(alpha + da); - var wt = da * (1 / 6), - est = area = wt * (fv0 + 4 * fv2 + fv4); - // Have estimate est of integral on (alpha, alpha+da). - // Bisect and compute estimates on left and right half intervals. - // Sum is better value for integral. - var cont = true; - while(cont) { - var dx = 0.5 * da, - arg = alpha + 0.5 * dx; - fv1 = f(arg); - fv3 = f(arg + dx); - var wt = dx * (1 / 6), - estl = wt * (fv0 + 4 * fv1 + fv2), - estr = wt * (fv2 + 4 * fv3 + fv4), - sum = estl + estr, - diff = est - sum; - area = parea + sum; - var b2 = alpha + da; - if (Math.abs(Math.abs(integral - area) - Math.abs(pdiff)) - + Math.abs(diff) <= fv4 * accuracy * (b2 - a)) { - return { b: b2, area: area }; - } - if (Math.abs(integral - area) > Math.abs(pdiff + diff)) { - if (integral <= area) { - index = 0; - p = table[0]; - p.left = true; - p.psum = parea; - } else { - if ((Math.abs(diff) <= fv4 * accuracy * da - || dx <= dxmin) && da <= dxmax) { - // Accept approximate integral sum. - // If it was a right interval, add results to finish - // at this level. If it was a left interval, process - // right interval. - pdiff += diff; - for (;;) { - if (!p.left) { // process right-half interval - parea += sum; - alpha += da; - p.left = true; - p.psum = sum; - fv0 = p.f1t; - fv2 = p.f2t; - fv4 = p.f3t; - da = p.dat; - est = p.estr; - break; - } - sum += p.psum; - parea -= p.psum; - if (--index <= 0) { - index = 0; - p = table[0]; - p.psum = parea = sum; - alpha += da; - b += b - a; - cont = false; - break; - } else { - p = table[index]; - } - } - continue; - } - } - } - if (index >= table.length) - return null; - // Raise level and store information for processing right-half - // interval. - da = dx; - est = estl; - p = table[++index] = { - left: false, - f1t: fv2, - f2t: fv3, - f3t: fv4, - dat: dx, - estr: estr, - psum: 0 - }; - fv4 = fv2; - fv2 = fv1; - } - } - return { b: b, area: area }; } };