mirror of
https://github.com/scratchfoundation/paper.js.git
synced 2025-01-05 20:32:00 -05:00
Replace slow simpson() method with insanely fast Gauss-Legendre Numerical Integration by Jim Armstrong which was further optimised.
This commit is contained in:
Jürg Lehni
parent
c4ad95b0ac
commit
29e57cc521
2 changed files with 81 additions and 233 deletions
|
|
@ -125,22 +125,21 @@ var Curve = this.Curve = Base.extend({
|
||||||
|| this._path.closed && curves[curves.length - 1]) || null;
|
|| this._path.closed && curves[curves.length - 1]) || null;
|
||||||
},
|
},
|
||||||
|
|
||||||
// Calculates arclength of a cubic using adaptive simpson integration.
|
getLength: function() {
|
||||||
getLength: function(goal) {
|
|
||||||
var z0 = this._segment1._point,
|
var z0 = this._segment1._point,
|
||||||
z1 = this._segment2._point,
|
z1 = this._segment2._point,
|
||||||
c0 = z0.add(this._segment1._handleOut),
|
c0 = z0.add(this._segment1._handleOut),
|
||||||
c1 = z1.add(this._segment2._handleIn);
|
c1 = z1.add(this._segment2._handleIn);
|
||||||
// TODO: Check for straight lines and handle separately.
|
// TODO: Check for straight lines and handle separately.
|
||||||
|
|
||||||
// Calculate the coefficients of a Bezier derivative, divided by 3.
|
// Calculate the coefficients of a Bezier derivative.
|
||||||
var ax = 3 * (c0.x - c1.x) - z0.x + z1.x,
|
var ax = 9 * (c0.x - c1.x) + 3 * (z1.x - z0.x),
|
||||||
bx = 2 * (z0.x + c1.x) - 4 * c0.x,
|
bx = 6 * (z0.x + c1.x) - 12 * c0.x,
|
||||||
cx = c0.x - z0.x,
|
cx = 3 * (c0.x - z0.x),
|
||||||
|
|
||||||
ay = 3 * (c0.y - c1.y) - z0.y + z1.y,
|
ay = 9 * (c0.y - c1.y) + 3 * (z1.y - z0.y),
|
||||||
by = 2 * (z0.y + c1.y) - 4 * c0.y,
|
by = 6 * (z0.y + c1.y) - 12 * c0.y,
|
||||||
cy = c0.y - z0.y;
|
cy = 3 * (c0.y - z0.y);
|
||||||
|
|
||||||
function ds(t) {
|
function ds(t) {
|
||||||
// Calculate quadratic equations of derivatives for x and y
|
// 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);
|
return Math.sqrt(dx * dx + dy * dy);
|
||||||
}
|
}
|
||||||
|
|
||||||
var integral = MathUtils.simpson(ds, 0.0, 1.0, MathUtils.EPSILON, 1.0);
|
return MathUtils.gauss(ds, 0.0, 1.0, 8);
|
||||||
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;
|
|
||||||
},
|
},
|
||||||
|
|
||||||
|
/**
|
||||||
|
* 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() {
|
clone: function() {
|
||||||
return new Curve(this._segment1, this._segment2);
|
return new Curve(this._segment1, this._segment2);
|
||||||
},
|
},
|
||||||
|
|
|
||||||
|
|
@ -1,220 +1,47 @@
|
||||||
var MathUtils = {
|
var MathUtils = new function() {
|
||||||
// TODO: 64 bit: 53 digits, 32 bit: 24 digits
|
var maxDepth = 53; // MANT_DIGITS: 64 bit: 53 digits, 32 bit: 24 digits
|
||||||
// TODO: Naming? What is MANT standing for?
|
|
||||||
MANT_DIGITS: 53,
|
|
||||||
EPSILON: Math.pow(2, -52),
|
|
||||||
|
|
||||||
// Compute a numerical approximation to an integral via adaptive Simpson's
|
// Gauss-Legendre Numerical Integration, a Gauss.as port from Singularity:
|
||||||
// Rule This routine ignores underflow.
|
//
|
||||||
|
// 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.
|
var abscissa = [
|
||||||
// f: Pointer to function to be integrated.
|
-0.5773502692, 0.5773502692,
|
||||||
// a, b: Lower, upper limits of integration.
|
-0.7745966692, 0.7745966692, 0,
|
||||||
// accuracy: Desired relative accuracy of integral,
|
-0.8611363116, 0.8611363116, -0.3399810436, 0.3399810436,
|
||||||
// try to make |error| <= accuracy*abs(integral).
|
-0.9061798459, 0.9061798459, -0.5384693101, 0.5384693101, 0.0000000000,
|
||||||
// dxmax: Maximum limit on the width of a subinterval
|
-0.9324695142, 0.9324695142, -0.6612093865, 0.6612093865, -0.2386191861, 0.2386191861,
|
||||||
// For periodic functions, dxmax should be
|
-0.9491079123, 0.9491079123, -0.7415311856, 0.7415311856, -0.4058451514, 0.4058451514, 0.0000000000,
|
||||||
// set to the period or smaller to prevent
|
-0.9602898565, 0.9602898565, -0.7966664774, 0.7966664774, -0.5255324099, 0.5255324099, -0.1834346425, 0.1834346425
|
||||||
// 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;
|
|
||||||
|
|
||||||
// Have estimate est of integral on (alpha, alpha+da).
|
var weight = [
|
||||||
// Bisect and compute estimates on left and right half intervals.
|
1, 1,
|
||||||
// integral is the best value for the integral.
|
0.5555555556, 0.5555555556, 0.8888888888,
|
||||||
for (;;) {
|
0.3478548451, 0.3478548451, 0.6521451549, 0.6521451549,
|
||||||
var dx = 0.5 * da,
|
0.2369268851, 0.2369268851, 0.4786286705, 0.4786286705, 0.5688888888,
|
||||||
arg = alpha + 0.5 * dx;
|
0.1713244924, 0.1713244924, 0.3607615730, 0.3607615730, 0.4679139346, 0.4679139346,
|
||||||
fv1 = f(arg);
|
0.1294849662, 0.1294849662, 0.2797053915, 0.2797053915, 0.3818300505, 0.3818300505, 0.4179591837,
|
||||||
fv3 = f(arg + dx);
|
0.1012285363, 0.1012285363, 0.2223810345, 0.2223810345, 0.3137066459, 0.3137066459, 0.3626837834, 0.3626837834
|
||||||
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;
|
|
||||||
|
|
||||||
if (index >= table.length)
|
return {
|
||||||
success = false;
|
EPSILON: Math.pow(2, -maxDepth + 1),
|
||||||
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];
|
|
||||||
}
|
|
||||||
|
|
||||||
} else {
|
gauss: function(f, a, b, n) {
|
||||||
// Raise level and store information for processing right-half
|
n = Math.min(Math.max(n, 2), 8);
|
||||||
// interval.
|
|
||||||
p = table[++index] = {
|
var l = n == 2 ? 0 : n * (n - 1) / 2 - 1,
|
||||||
left: false,
|
sum = 0,
|
||||||
f1t: fv2,
|
mul = 0.5 * (b - a),
|
||||||
f2t: fv3,
|
ab2 = mul + a;
|
||||||
f3t: fv4,
|
for(var i = 0; i < n; i++)
|
||||||
dat: dx,
|
sum += f(ab2 + mul * abscissa[l + i]) * weight[l + i];
|
||||||
estr: estr
|
|
||||||
};
|
return mul * sum;
|
||||||
da = dx;
|
|
||||||
est = estl;
|
|
||||||
fv4 = fv2;
|
|
||||||
fv2 = fv1;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
},
|
|
||||||
|
|
||||||
// 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 };
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue