DIRECTORY
Cubic2, CubicPaths, CubicPathsExtras, CubicSplines, Imager, ImagerPath, ImagerTransformation, Polynomial, Real, RealFns, Vector2;

CubicPathsImpl: CEDAR PROGRAM
IMPORTS Cubic2, ImagerTransformation, Polynomial, RealFns, Vector2
EXPORTS CubicPaths, CubicPathsExtras
~ BEGIN OPEN CubicPaths;

Bezier: TYPE = Cubic2.Bezier;
Coeffs: TYPE = Cubic2.Coeffs;
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
VEC: TYPE = Imager.VEC;
X: NAT = CubicSplines.X;
Y: NAT = CubicSplines.Y;
PathFromCubic: PUBLIC PROC [coeffs: CubicSplines.CoeffsSequence, newBounds: BOOL _ TRUE] RETURNS [Path] = {
path: Path _ NEW[PathRec];
path.cubics _ NEW[PathSequence[coeffs.length]];
FOR i: NAT IN [0..coeffs.length) DO
cfs: Cubic2.Coeffs _ [
c3: [coeffs[i].t3[X], coeffs[i].t3[Y]],
c2: [coeffs[i].t2[X], coeffs[i].t2[Y]],
c1: [coeffs[i].t1[X], coeffs[i].t1[Y]],
c0: [coeffs[i].t0[X], coeffs[i].t0[Y]]];
path.cubics[i] _ Cubic2.CoeffsToBezier[cfs];
ENDLOOP;
IF newBounds THEN UpdateBounds[path];
RETURN [path];
};

EnumeratePath: PUBLIC PROC [path: Path, moveTo: ImagerPath.MoveToProc, curveTo: ImagerPath.CurveToProc] = {
moveTo[path.cubics[0].b0];
FOR i: NAT IN [0..path.cubics.length) DO
curveTo[path.cubics[i].b1, path.cubics[i].b2, path.cubics[i].b3];
ENDLOOP;
};

WalkPath: PUBLIC PROC [path: Path, tol: REAL, proc: PointProc] = {
subdivide: PROC [bezier: Cubic2.Bezier] RETURNS [quit: BOOLEAN]= {
quit _ FALSE;
IF Cubic2.Flat[bezier, tol] THEN {
len: REAL _ Vector2.Length[Vector2.Sub[bezier.b0, bezier.b3]];
IF len > 0 THEN {
step: REAL _ tol/len;
alpha: REAL _ step;
pt: VEC;
UNTIL alpha > 1 DO
pt _ Vector2.Add[Vector2.Mul[bezier.b3,alpha],Vector2.Mul[bezier.b0,1-alpha]];
IF proc[pt] THEN {quit _ TRUE; EXIT};
alpha _ alpha+step;
ENDLOOP;
}
ELSE quit _ proc[bezier.b3];
RETURN [quit]}
ELSE {
b1, b2: Cubic2.Bezier;
[b1,b2] _ Cubic2.Split[bezier];
IF NOT subdivide[b1] THEN [] _ subdivide[b2];
};
};
IF NOT proc[path.cubics[0].b0] THEN FOR i: NAT IN [0..path.cubics.length) DO
[] _ subdivide[path.cubics[i]];
ENDLOOP;
};
 
PointOnPath: PUBLIC PROC [pt: VEC, path: Path, tol: REAL] RETURNS [BOOLEAN] = {
tolSq: REAL _ tol*tol;
found: BOOLEAN _ FALSE;
hit: PointProc = {
IF Vector2.Square[Vector2.Sub[p, pt]] <= tolSq THEN found _ TRUE;
RETURN [found];
};
IF pt.x < path.bounds.x OR pt.y < path.bounds.y OR pt.x>(path.bounds.x+path.bounds.w)
OR pt.y > (path.bounds.y+path.bounds.h) THEN RETURN [FALSE]
ELSE WalkPath[path, tol, hit];
RETURN [found];
}; 

TranslatePath: PUBLIC PROC [path: Path, offset: VEC] = {
FOR i: NAT IN [0..path.cubics.length) DO
path.cubics[i].b0 _ Vector2.Add[offset, path.cubics[i].b0];
path.cubics[i].b1 _ Vector2.Add[offset, path.cubics[i].b1];
path.cubics[i].b2 _ Vector2.Add[offset, path.cubics[i].b2];
path.cubics[i].b3 _ Vector2.Add[offset, path.cubics[i].b3];
ENDLOOP;
path.bounds.x _ path.bounds.x+offset.x;
path.bounds.y _ path.bounds.y+offset.y;
};

TransformPath: PUBLIC PROC [path: Path, tranformation: ImagerTransformation.Transformation, newBounds: BOOL _ TRUE] = {
FOR i: NAT IN [0..path.cubics.length) DO
path.cubics[i].b0 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b0];
path.cubics[i].b1 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b1];
path.cubics[i].b2 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b2];
path.cubics[i].b3 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b3];
ENDLOOP;
IF newBounds THEN UpdateBounds[path];
}; 

<<UpdateBounds: PUBLIC PROC [path: Path] = {
minX, minY: REAL _ 10E6;
maxX, maxY: REAL _ -10E6;
pointProc: PointProc = {
IF p.x > maxX THEN maxX _ p.x;
IF p.y > maxY THEN maxY _ p.y;
IF p.x < minX THEN minX _ p.x;
IF p.y < minY THEN minY _ p.y;
};
WalkPath[path, 1.0, pointProc];
path.bounds _ [x: minX, y: minY, w: maxX-minX, h: maxY-minY];
};
>>
Evaluate: PROC [a3, a2, a1, a0, b3, b2, b1, b0: REAL, u: REAL] RETURNS [pt: VEC] = {
pt.x _ (((a3)*u + a2)*u + a1)*u + a0;
pt.y _ (((b3)*u + b2)*u + b1)*u + b0;
};

EvaluateX: PROC [coeffs: Cubic2.Coeffs, t: REAL] RETURNS [x: REAL] = {
x _ (((coeffs.c3.x)*t + coeffs.c2.x)*t + coeffs.c1.x)*t + coeffs.c0.x;
};

EvaluateY: PROC [coeffs: Cubic2.Coeffs, t: REAL] RETURNS [y: REAL] = {
y _ (((coeffs.c3.y)*t + coeffs.c2.y)*t + coeffs.c1.y)*t + coeffs.c0.y;
};

UpdateBounds: PUBLIC PROC [path: Path] = {
minX, minY: REAL _ 10E6;
maxX, maxY: REAL _ -10E6;
loEndX, hiEndX, loEndY, hiEndY: REAL;
coeffs: Cubic2.Coeffs;
tArray: ARRAY [1..2] OF REAL;
ptX, ptY: REAL;
rootCount: NAT;
FOR i: NAT IN [0..path.cubics.length) DO
coeffs _ Cubic2.BezierToCoeffs[path.cubics[i]];
[tArray, rootCount] _ QuadraticFormula[3.0*coeffs.c3.x, 2.0*coeffs.c2.x, coeffs.c1.x];
FOR j: NAT IN [1..rootCount] DO
IF tArray[j] >= 0.0 AND tArray[j] <= 1.0 THEN {
ptX _ EvaluateX[coeffs, tArray[j]];
minX _ MIN[minX, ptX];
maxX _ MAX[maxX, ptX];
};
ENDLOOP;
[tArray, rootCount] _ QuadraticFormula[3.0*coeffs.c3.y, 2.0*coeffs.c2.y, coeffs.c1.y];
FOR j: NAT IN [1..rootCount] DO
IF tArray[j] >= 0.0 AND tArray[j] <= 1.0 THEN {
ptY _ EvaluateY[coeffs, tArray[j]];
minY _ MIN[minY, ptY];
maxY _ MAX[maxY, ptY];
};
ENDLOOP;
loEndX _ path.cubics[i].b0.x;
loEndY _ path.cubics[i].b0.y;
hiEndX _ path.cubics[i].b3.x;
hiEndY _ path.cubics[i].b3.y;
minX _ MIN[minX, loEndX, hiEndX];
minY _ MIN[minY, loEndY, hiEndY];
maxX _ MAX[maxX, loEndX, hiEndX];
maxY _ MAX[maxY, loEndY, hiEndY];
ENDLOOP;
path.bounds _ [x: minX, y: minY, w: maxX-minX, h: maxY-minY];
};

ClosestPoint: PUBLIC PROC [pt: VEC, path: Path] RETURNS [VEC] = {
closest: VEC _ path.cubics[0].b0;
minD: REAL _ 10E6;
test: PointProc = {
thisD: REAL _ Vector2.Square[Vector2.Sub[p, pt]];
IF thisD < minD THEN {closest _ p; minD _ thisD};
};
WalkPath[path, 1.0, test];
RETURN [closest];
};

ClosestPointAnalytic: PUBLIC PROC [pt: VEC, path: Path, tolerance: REAL _ 9999.0] RETURNS [closest: VEC, success: BOOL] = {
bezier: Bezier;
c: Coeffs;
bestDist2, thisDist2: REAL;
thisPt: VEC;
thisSuccess: BOOL;

bestDist2 _ Real.LargestNumber;
success _ FALSE;

FOR i: NAT IN [0..path.cubics.length) DO
bezier _ path.cubics[i];
c _ Cubic2.BezierToCoeffs[bezier];
[thisPt, thisDist2, thisSuccess] _ ClosestPointCubic[c.c3.x, c.c2.x, c.c1.x, c.c0.x, c.c3.y, c.c2.y, c.c1.y, c.c0.y, pt];
IF NOT thisSuccess THEN ERROR;
IF thisDist2 < bestDist2 THEN {
bestDist2 _ thisDist2;
closest _ thisPt;
success _ TRUE;
};
ENDLOOP;
};

ClosestPointCubic: PROC [a3, a2, a1, a0, b3, b2, b1, b0: REAL, q: VEC] RETURNS [closest: VEC, bestD2: REAL, success: BOOL] = {
c5, c4, c3, c2, c1, c0: REAL;
poly: Polynomial.Ref;
roots: ShortRealRootRec;
thisD2: REAL;
thisPt: VEC;
c5 _ 3*(a3*a3 + b3*b3);
c4 _ 5*(a2*a3 + b2*b3);
c3 _ 4*(a1*a3 + b1*b3) + 2*(a2*a2 + b2*b2);
c2 _ 3*(a1*a2 + b1*b2) + 3*(a3*(a0 - q.x) + b3*(b0 - q.y));
c1 _ 2*(a2*(a0 - q.x) + b2*(b0 - q.y)) + (a1*a1 + b1*b1);
c0 _ (a1*(a0 - q.x) + b1*(b0 - q.y));
poly _ Polynomial.Quintic[[c0, c1, c2, c3, c4, c5]];
roots _ CheapRealRootsInInterval[poly, 0.0, 1.0];
bestD2 _ Real.LargestNumber;
success _ FALSE;
FOR i: NAT IN [0..roots.nRoots) DO
IF roots.realRoot[i] > 0.0 AND roots.realRoot[i] <= 1.0 THEN {
thisPt _ Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, roots.realRoot[i]];
thisD2 _ Vector2.Square[Vector2.Sub[thisPt, q]];
IF thisD2 < bestD2 THEN {
bestD2 _ thisD2;
closest _ thisPt;
success _ TRUE;
};
};
ENDLOOP;
FOR u: REAL _ 0, u+1.0 UNTIL u > 1.0 DO
thisPt _ Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, u];
thisD2 _ Vector2.Square[Vector2.Sub[thisPt, q]];
IF thisD2 < bestD2 THEN {
bestD2 _ thisD2;
closest _ thisPt;
success _ TRUE;
};
ENDLOOP;
};

CopyPath: PUBLIC PROC [path: Path] RETURNS [Path] = {
new: Path _ NEW[PathRec];
new.cubics _ NEW[PathSequence[path.cubics.length]];
FOR i: NAT IN [0..path.cubics.length) DO
new.cubics[i] _ path.cubics[i];
ENDLOOP;
new.bounds _ path.bounds;
RETURN [new];
};

ReversePath: PUBLIC PROC [path: Path] = {
last: NAT _ path.cubics.length-1;
FOR i: NAT IN [0..path.cubics.length/2) DO
temp: Cubic2.Bezier _ path.cubics[i];
path.cubics[i] _ path.cubics[last-i];
path.cubics[last-i] _ temp;
ENDLOOP;
FOR i: NAT IN [0..path.cubics.length) DO
temp: VEC _ path.cubics[i].b0;
path.cubics[i].b0 _ path.cubics[i].b3;
path.cubics[i].b3 _ temp;
temp _ path.cubics[i].b1;
path.cubics[i].b1 _ path.cubics[i].b2;
path.cubics[i].b2 _ temp;
ENDLOOP;
};
QuadraticFormula: PUBLIC PROC [a, b, c: REAL] RETURNS [roots: ARRAY [1..2] OF REAL, rootCount: NAT] = {
discriminant, temp: REAL;
IF a = 0 THEN {
[roots[1], rootCount] _ LinearFormula [b, c];
RETURN;
};
IF c = 0 THEN {
roots[1] _ 0.0;
[roots[2], rootCount] _ LinearFormula[a, b];
rootCount _ rootCount + 1;
IF roots[1] > roots[2] THEN { -- swap them
temp _ roots[1];
roots[1] _ roots[2];
roots[2] _ temp;
};
RETURN};
discriminant _ b*b - 4*a*c; -- 3 mult, 1 add
SELECT discriminant FROM
=0 => {rootCount _ 1;
roots[1] _ -b/(2*a);
RETURN};	-- 1 mult, 1 div (4 mult, 1 div, 1 add total)
<0 => {rootCount _ 0;
RETURN}; -- (3 mult, 1 add total)
>0 => {
sqrtRadical: REAL _ RealFns.SqRt[discriminant];
term: REAL;
rootCount _ 2;
term _ IF b < 0 THEN sqrtRadical - b ELSE -sqrtRadical - b;
roots[1] _ term/(2.0*a);
roots[2] _ (2.0*c)/term;
IF roots[1] > roots[2] THEN { -- swap them
temp _ roots[1];
roots[1] _ roots[2];
roots[2] _ temp;
};
RETURN}; -- 1 SqRt, 2 mult, 2 div, 2 add (1 SqRt, 5 mult, 2 div, 3 add total)
ENDCASE => ERROR;
}; -- end of QuadraticFormula

LinearFormula: PROC [a, b: REAL] RETURNS [root: REAL, rootCount: NAT] = {
IF a = 0 THEN {
IF b = 0 THEN {root _ 0.0; rootCount _ 0; RETURN}
ELSE {root _ 0.0; rootCount _ 0; RETURN}	-- inconsistent case
}
ELSE {
IF b = 0  THEN {root _ 0.0; rootCount _ 1; RETURN}
ELSE {root _ -b/a; rootCount _ 1; RETURN};
};
};

CheapRealRootsInInterval: PUBLIC PROC [poly: Polynomial.Ref, lo, hi: REAL] RETURNS [roots: ShortRealRootRec] = {
d: NAT _ Polynomial.Degree[poly];
roots.nRoots _ 0;
IF d<=1 THEN {
IF d=1 THEN {roots.nRoots _ 1; roots.realRoot[0] _ -poly[0]/poly[1]}
}
ELSE {
savedCoeff: ARRAY [0..5] OF REAL;
FOR i: NAT IN [0..d] DO savedCoeff[i] _ poly[i] ENDLOOP;
Polynomial.Differentiate[poly];
BEGIN
extrema: ShortRealRootRec _ CheapRealRootsInInterval[poly, lo, hi];
x: REAL;
FOR i: NAT IN [0..d] DO poly[i] _ savedCoeff[i] ENDLOOP;
IF extrema.nRoots>0 THEN {
x _ RootBetween[poly, lo, extrema.realRoot[0]];
IF x<= extrema.realRoot[0] THEN {roots.nRoots _ 1; roots.realRoot[0] _ x};
FOR i: NAT IN [0..extrema.nRoots-1) DO
x _ RootBetween[poly, extrema.realRoot[i], extrema.realRoot[i+1]];
IF x <= extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1};
ENDLOOP;
x _ RootBetween[poly, extrema.realRoot[extrema.nRoots-1], hi];
IF x <= hi THEN {roots.realRoot[roots.nRoots] _ x; roots.nRoots _ roots.nRoots + 1}
}
ELSE {
x _ RootBetween[poly, lo, hi];
IF x <= hi THEN {roots.realRoot[0] _ x; roots.nRoots _ 1}
};
END
};
};

RootBetween: PROC [poly: Polynomial.Ref, x0, x1: REAL]
RETURNS [x: REAL] =
BEGIN
OPEN Polynomial;
y0: REAL _ Polynomial.Eval[poly,x0];
y1: REAL _ Polynomial.Eval[poly,x1];
xx: REAL;
xx0: REAL _ IF y0<y1 THEN x0 ELSE x1;
xx1: REAL _ IF y0<y1 THEN x1 ELSE x0;
yy0: REAL _ MIN[y0,y1];
yy1: REAL _ MAX[y0,y1];
NewtonStep: PROC [x: REAL] RETURNS [newx:REAL] =
BEGIN
y: REAL _ Polynomial.Eval[poly, x];
deriv: REAL _ Polynomial.EvalDeriv[poly, x];
IF y<0 THEN {IF y>yy0 THEN {xx0 _ x; yy0 _ y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 _ x; yy1 _ y}};
newx _ IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) THEN RETURN [x1+100];
xx _ x _ (x0+x1)/2;
WHILE x IN [x0..x1] DO
newx:REAL _ NewtonStep[x];
IF x=newx THEN RETURN;
x _ NewtonStep[newx];
xx _ NewtonStep[xx];
IF xx=x THEN EXIT;
ENDLOOP;
IF xx0 IN [x0..x1] AND xx1 IN [x0..x1] THEN
BEGIN
x0 _ MIN[xx0,xx1];
x1 _ MAX[xx0,xx1];
END;
y0 _ Polynomial.Eval[poly,x0];
y1 _ Polynomial.Eval[poly,x1];
THROUGH [0..500) DO
y: REAL _ Polynomial.Eval[poly, x _ (x0+x1)/2];
IF x=x0 OR x=x1 THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF (y > 0 AND y0 < 0) OR (y < 0 AND y0 > 0) THEN {x1 _ x; y1 _ y}
ELSE {x0 _ x; y0 _ y};
IF (y0 > 0 AND y1 > 0) OR (y0 < 0 AND y1 < 0) OR (y0=0 OR y1=0) THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;

GetParam: PUBLIC PROC [bezier: Cubic2.Bezier, pt: VEC] RETURNS [u: REAL] = {
p0, p1, p2, p3: VEC;
a,b,c,d: REAL;
ignoreX, ignoreY: BOOL _ FALSE;
polyX, polyY: Polynomial.Ref;
epsilon: REAL =  0.01;
rootsX, rootsY: ShortRealRootRec;
p0 _ bezier.b0;
p1 _ bezier.b1;
p2 _ bezier.b2;
p3 _ bezier.b3;
d _ p0.x - pt.x;
c _ 3*(p1.x - p0.x);
b _ 3*p0.x - 6*p1.x + 3*p2.x;
a _ -p0.x + 3*(p1.x - p2.x) + p3.x;
IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreX _ TRUE
ELSE {
polyX _ Polynomial.Cubic[[d,c,b,a]];
rootsX _ CheapRealRootsInInterval[polyX, 0.0, 1.0];
};
d _ p0.y - pt.y;
c _ 3*(p1.y - p0.y);
b _ 3*p0.y - 6*p1.y + 3*p2.y;
a _ -p0.y + 3*(p1.y - p2.y) + p3.y;
IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreY _ TRUE
ELSE {
polyY _ Polynomial.Cubic[[d,c,b,a]];
rootsY _ CheapRealRootsInInterval[polyY, 0.0, 1.0];
};
IF ignoreX THEN {
IF ignoreY THEN ERROR;
RETURN [rootsY.realRoot[0]];
}
ELSE IF ignoreY THEN {
RETURN [rootsX.realRoot[0]];
}
ELSE {
FOR i: NAT IN [0..rootsX.nRoots) DO
FOR j: NAT IN [0..rootsY.nRoots) DO
IF ABS[rootsX.realRoot[i] - rootsY.realRoot[j]] < epsilon THEN
RETURN [(rootsX.realRoot[i] + rootsY.realRoot[j]) / 2.0];
ENDLOOP;
ENDLOOP;
ERROR;
};
};

CubicMeetsLine: PUBLIC PROC [bezier: Bezier, a, b, c: REAL] RETURNS [points: ARRAY [0..2] OF VEC, hitCount: [0..3], tangency: ARRAY [0..2] OF BOOL] = {
f: Coeffs;
poly: Polynomial.Ref;
roots: ShortRealRootRec;
lastT: REAL;
success: BOOL;
epsilon: REAL = 0.0001;

hitCount _ 0;
tangency[0] _ tangency[1] _ tangency[2] _ FALSE;
f _ Cubic2.BezierToCoeffs[bezier];
poly _ Polynomial.Cubic[[a*f.c0.x + b*f.c0.y + c, a*f.c1.x + b*f.c1.y, a*f.c2.x + b*f.c2.y, a*f.c3.x + b*f.c3.y]];
roots _ CheapRealRootsInInterval[poly, 0.0, 1.0];

lastT _ 2.0; -- a value outside the range [0..1];
success _ FALSE;
FOR i: NAT IN [0..roots.nRoots) DO
IF roots.realRoot[i] > 0.0 AND roots.realRoot[i] <= 1.0 THEN {
IF ABS[roots.realRoot[i]-lastT] > epsilon THEN {
lastT _ roots.realRoot[i];
points[hitCount] _ Evaluate[f.c3.x, f.c2.x, f.c1.x, f.c0.x, f.c3.y, f.c2.y, f.c1.y, f.c0.y, lastT];
hitCount _ hitCount + 1;
}
ELSE {
IF lastT # 2.0 THEN tangency[hitCount-1] _ TRUE;
};
};
ENDLOOP;
};

END.
CubicPathsImpl.mesa
Copyright Σ 1985, 1987 by Xerox Corporation.  All rights reserved.
Written by:  Maureen Stone August 5, 1985 8:49:21 pm PDT
Kurlander August 15, 1986 10:08:11 am PDT
Russ Atkinson (RRA) February 2, 1987 10:44:17 pm PST
Pier, August 5, 1987 3:50:54 pm PDT
Bier, February 6, 1989 4:44:55 pm PST


path.cubics[0].b0;
Find vertical points (find dx/dt and solve for t)
Find horizontal points (find dy/dt and solve for t)
Consider Endpoints
CubicPathsExtras all below plus made UpdateBounds Public
Find the fifth degree equation in the parameter u and solve using an iterative root finder.
The spline equation is first represented as:
x(u) = a3*u3 + a2*u2 + a1*u + a0
y(u) = b3*u3 + b2*u2 + b1*u + b0
for 0 < u < 1.
a3, a2, a1, a0, b3, b2, b1, b0: REAL;
p0, p1, p2, p3: VEC;
p0 _ bezier.b0; p1 _ bezier.b1; p2 _ bezier.b2; p3 _ bezier.b3;
a3 _ -p0.x + 3*p1.x - 3*p2.x + p3.x;
a2 _ 3*p0.x - 6*p1.x + 3*p2.x;
a1 _ -3*p0.x + 3*p1.x;
a0 _ p0.x;
b3 _ -p0.y + 3*p1.y - 3*p2.y + p3.y;
b2 _ 3*p0.y - 6*p1.y + 3*p2.y;
b1 _ -3*p0.y + 3*p1.y;
b0 _ p0.y;
Given a cubic parametric curve
x(u) = a3*u3 + a2*u2 + a1*u + a0
y(u) = b3*u3 + b2*u2 + b1*u + b0,
for 0 < u < 1.
we wish to find the closest point on this curve to q.
First, we find the extrema of (x - q.x)2 + (y - q.y)2 by setting
dx/du(x - q.x) + dy/du(y - q.y) = 0.
This gives us up to 5 real roots.  We eliminate those for which u < 0 or u > 1.  We evaluate the spline at the remain values of u and compute (x - q.x)2 + (y - q.y)2.  We also compute the distance from q to the spline endpoints at u = 0 and u = 1.  We choose the best point.
reverses it in place

The solution to the equation "ax2+bx+c=0".  If a=0 this is just a linear equation.  If
c = 0, one root is zero and we solve a linear equation.  Otherwise, we use the
quadratic formula in either the form [-b+-(b2-4ac)1/2]/2a or (2c)/[-b-+(b2-4ac)1/2]
depending on the sign of b.  We will arrange the roots so that roots[1] < roots[2].
The solution to the equation "ax + b = 0".
Uses a modified Newton-Raphelson iteration to find the real roots of a poly.  It finds only those roots in the interval [lo..hi].  If the interval is relatively small (e.g., [0..1]) this routine is significantly faster than Polynomial.CheapRootsInInterval.
debug  xi: ARRAY [0..50) OF REAL;
debug  i: NAT _ 0;
first try Newton's method
debug  xi[i] _ xx; i _ i+1;
If it oscillated or went out of range, try bisection
debug  SIGNAL Bisection;
Assumes that pt is a point on bezier.  Finds a parameter, u, corresponding to that point.
ShortRealRootRec: TYPE = RECORD [nRoots: NAT, realRoot: ARRAY [0..5) OF REAL];
Take advantage of the roots of both polynomials to improve accuracy.
The spline equation is represented as:
x(u) = a3*t3 + a2*t2 + a1*t + a0
y(u) = b3*t3 + b2*t2 + b1*t + b0
for 0 < u < 1.
The line is represented as:
ax + by + c = 0.
Thus, the points of intersection are at the t values that satisfy:
a*(a3*t3 + a2*t2 + a1*t + a0) + b*(b3*t3 + b2*t2 + b1*t + b0) + c = 0.
or equivalently:
(a*a3 + b*b3)*t3 + (a*a2 + b*b2)*t2 + (a*a1 + b*b1)*t + (a*a0 + b*b0 + c) = 0.
Κœ˜
codešœΟk
œ™Kš
œB™BKš
œ8™8K™)K™4K™#K™%K™—š
	˜	Jš
œ
˜
K˜—šΟnœœ
˜Jšœ;˜BKšœ˜$Kšœ
œ˜K˜Kšœœ˜Kšœœ˜Jšœœ˜5Kšœœ
œ
˜Jšœœ˜Jšœœ˜—K™šž
œœ
œ2œœœ˜kKšœ
œ
˜Kšœœ˜/šœœ
œ˜#˜K˜'K˜'K˜'K˜(—K˜,Kšœ
˜—Kšœœ˜%Kšœ˜K˜K˜—šž
œœQ˜kJ˜šœœ
œ˜(Kš
œA˜AKšœ
˜—K˜K˜—šžœœœ˜Bšœœœœ˜BJšœœ
˜
šœœ˜"Kšœœ5˜>˜Kšœœ˜Kšœœ˜Kšœœ
˜šœ˜Kš
œN˜NKšœ
œ	œœ˜%K˜Kšœ
˜—K˜
—Kšœ˜Jšœ˜—šœ˜Jš
œ˜J˜Jšœ
œœ˜-J˜—J˜—šœ
œœœ
œ˜LKš
œ˜Kšœ
˜—Kš
œ˜Kš
œ
˜
—šžœœœœœœ˜OKšœœ˜Kšœœœ
˜˜Kšœ-œ	œ
˜AKšœ	˜K˜—šœœœ#˜UKšœ&œ
œœ
˜;—Kšœ˜Kšœ	˜Kš
œ˜K˜—šž
œœœ˜8šœœ
œ˜(K˜;K˜;K˜;K˜;Kšœ
˜—K˜'K˜'K˜K˜—šž
œœMœœ˜wšœœ
œ˜(Kš
œU˜UKš
œU˜UKš
œU˜UKš
œU˜UKšœ
˜—Kšœœ˜%Kš
œ˜K˜—šžœœ˜,Kšœœ˜Kšœœ	˜š
œ˜Kšœœ˜Kšœœ˜Kšœœ˜Kšœœ˜Kš
œ˜—K˜K˜=K˜Kš
ž˜—šžœœΟd
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œœœœœ˜TJš	œŸ
œŸ
œŸ
œŸ
œ
˜%Jš	œŸ
œŸ
œŸ
œŸ
œ
˜%J˜J˜—š
ž	œœœœœ˜FJš
œF˜FJ˜J˜—š
ž	œœœœœ˜FJš
œF˜FJ˜J˜—šžœœ
œ˜*Kš
œ™Kšœœ˜Kšœœ	˜Kšœ œ
˜%Kš
œ˜Kšœœœ
œ
˜Kšœ
œ
˜Kšœœ
˜šœœ
œ˜(Kš
œ/˜/Kš
Οb1™1Kš
œV˜Všœœ
œ˜šœœœ˜/Kš
œ#˜#Kšœœ˜Kšœœ˜K˜—Kšœ
˜—Kš
 3™3Kš
œV˜Všœœ
œ˜šœœœ˜/Kš
œ#˜#Kšœœ˜Kšœœ˜K˜—Kšœ
˜—Kš
 ™Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ˜Kšœœ˜!Kšœœ˜!Kšœœ˜!Kšœœ˜!Kšœ
˜—K˜=Kš
œ˜—K˜K™8š
žœœœœœ˜AKšœ	œ˜!Kšœœ˜š
œ˜Kšœœ&˜1Kšœœ˜1K˜—K˜Kšœ˜K˜K˜—šžœœ
œœœœœœ˜{K™[K™,KšœŸ
œΟu
œŸ
œ‘
œŸ
œŸ
™ KšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ
™ KšœΟm
œ’
œ™Kš
œ˜K˜
Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œœ
™%Kšœ
Ÿ
œŸ
œŸ
œŸ
œœ
™Kšœœ
˜Kšœœ
˜Kšœ
œ
˜K˜Kš
œ˜Kšœ
œ
˜K˜šœœ
œ˜(Kš
œ˜Kš	œ
Ÿ
œŸ
œŸ
œŸ
œ
™?Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ™$Kš	œ
Ÿ
œŸ
œŸ
œŸ
œ™Kšœ
Ÿ
œŸ
œŸ
œ™Kšœ
Ÿ
œŸ
Ϫ
Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ™$Kš	œ
Ÿ
œŸ
œŸ
œŸ
œ™Kšœ
Ÿ
œŸ
œŸ
œ™Kšœ
Ÿ
œŸ
Ϫ
Kš
œ"˜"Kš
œy˜yKšœ
œ
œ
œ
˜šœœ˜Kš
œ˜Kš
œ˜Kšœ
œ
˜K˜—Kšœ
˜—K˜K˜—šžœœ"œœœœ
œœ˜~K™KšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ
™ KšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ™!Kšœ’
œ’
œ™Kš
œ5™5Kšœ'‘
œ‘
œ™@Kš’
œ’
œ
’
œ’
œ™$Kšœ—
‘
œ‘
œm™’Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œŸ
œœ
˜Kš
œ˜Kš
œ˜Kšœœ
˜
Kšœœ
˜Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ˜Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ˜Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œ˜+Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œ	˜;Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ
Ÿ
œŸ
œŸ
œŸ
œ˜9Kšœ
Ÿ
œŸ
œŸ
œŸ
œŸ
œ	˜%Kš
œ4˜4Kš
œ1˜1Kš
œ˜Kšœ
œ
˜šœœ
œ˜"šœœœ˜>KšœŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œ˜EKš
œ0˜0šœœ˜Kš
œ˜Kš
œ˜Kšœ
œ
˜K˜—K˜—Kšœ
˜—šœœœ	˜'KšœŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œŸ
œ˜5Kš
œ0˜0šœœ˜Kš
œ˜Kš
œ˜Kšœ
œ
˜K˜—Kšœ
˜—K˜K˜—šžœœœ˜5Kšœœ
˜Kšœ
œ#˜3šœœ
œ˜(Kš
œ˜Kšœ
˜—K˜Kšœ˜
K˜K˜—šžœœ˜)Kš
œ™Kšœœ˜!šœœ
œ˜*K˜%Kš
œ%˜%Kš
œ˜Kšœ
˜—šœœ
œ˜(Kšœœ˜K˜&K˜K˜K˜&K˜Kšœ
˜—K˜—K™
procšžœœ
œœœ	œœ
œ
œ˜gIlead1šœ ‘
ϱ
‘
œ‘œ‘
œ‘œU™ΝLšœœ
˜šœœ˜Lš
œ-˜-Lšœ
˜L˜—šœœ˜Lš
œ˜Lš
œ,˜,Lš
œ˜šœœΟc˜*L˜L˜L˜L˜—Lšœ˜—Lšœ£˜,Lšœ˜š
œ˜Lš
œ˜Lšœ£-˜6—š
œ˜Lšœ£˜!—š
œ˜Lšœ
œ˜/Lšœœ
˜Lš
œ˜Lšœœœœ˜;Lš
œ˜L˜šœœ£˜*L˜L˜L˜L˜—Lšœ£D˜M—Lšœœ
˜Lšœ£˜L˜—šž
œœœœœ
œ˜IL™*šœœ˜Lšœœœ
˜1Lšœœ£˜=L˜
—šœ˜Lšœœœ
˜2Lšœœ˜*L˜—Lš
œ˜—š
žœœ
œ œœ˜pJ™€Jšœœ˜!J˜šœœ˜Jšœœ9˜DJš
œ
˜
—šœ˜Jšœœœ
œ
˜!Jš
œœ
œœœ
˜8š
œ˜Jš
˜Jš
œC˜CJšœœ
˜Jš
œœ
œœœ
˜8šœœ˜Jš
œ/˜/Jšœœ+˜Jšœœ
œ˜&Jš
œB˜BJšœœE˜gJšœ
˜—Jš
œ>˜>Jšœ	œD˜SJš
œ
˜
—šœ˜Jš
œ˜Jšœ	œ*˜9J˜—Jš
˜—Jš
œ˜—J˜J˜—šžœœ œ
˜6Jšœœ˜Jš
˜Jšœœ	œ
œ
™!Jšœ
œ™Jšœ˜Jšœœ˜$Jšœœ˜$Jšœœ
˜	Jš	œœœœœ˜%Jš	œœœœœ˜%Jšœœœ˜Jšœœœ˜š
ž
œœœœœ˜0Jš
˜Jšœœ˜#Jšœœ!˜,Jšœœœœ˜/Jšœœœœ˜/Jšœœ	œœ˜*Jšœ
˜—Jšœœ
œ˜Jšœœ
œ˜Jšœ	œ	œ	œ	œ
œ
˜CJ˜Jš
œ™šœœ
˜Jšœœ˜Jšœœ
œ
˜J˜J˜Jš
œ™Jšœœ
œ
˜Jšœ
˜—Jš
œ4™4Jšœœ™š	œœ
œœ
˜+Jš
˜Jšœœ
˜Jšœœ
˜Jšœ
˜—Jš
œ˜Jš
œ˜šœ
˜Jšœœ(˜/šœœ˜Jšœ
œ
œœœ
œœ
œ˜7—Jš
œœ	œœ	œ˜AJšœ˜Jšœ	œ	œ	œ	œœœœ
œœœ
œœ
œ˜{Jšœ
˜—Jšœ
˜Jšœ
˜—J˜šžœœ
œœœœ˜LKš
œY™YKšœœ
˜Kšœ	œ
˜Kšœœœ
˜Kš
œ˜Kšœ	œ	˜Kš
œ!˜!Jš
œœœ
œœœ
œ™NK˜K˜K˜K˜Kš
œ˜K˜K˜Kš
œ#˜#Kšœ
œœ
œœ
œœ
œœ˜fšœ˜Kš
œ$˜$Kš
œ3˜3K˜—Kš
œ˜K˜K˜Kš
œ#˜#Kšœ
œœ
œœ
œœ
œœ˜fšœ˜Kš
œ$˜$Kš
œ3˜3K˜—šœ	œ˜Kšœ	œ
œ
˜Kšœ˜K˜
—šœ
œ	œ˜Kšœ˜K˜
—šœ˜K™Dšœœ
œ˜#šœœ
œ˜#šœ
œ4˜>Kšœ3˜9—Kšœ
˜—Kšœ
˜—Kšœ
˜K˜—K˜K˜—šžœœ
œœœ
œœ
œœœ
œ˜—
K™&KšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ
™ KšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ
™ Kšœ’
œ’
œ™K™K™K™BKšœŸ
œ‘
œŸ
œ‘
œŸ
œŸ
œŸ
œ‘
œŸ
œ‘
œŸ
œŸ
œ
™FK™KšœŸœŸœ‘
ŸœŸœŸ
œ‘
œŸ
œŸ
œ
Ÿ
œŸ
œ
™NK˜
Kš
œ˜Kš
œ˜Kšœœ
˜Kšœ	œ
˜Kšœ	œ
˜K˜Kš
œ
˜
Kšœ*œ
˜0Kš
œ"˜"Kš
œr˜rKš
œ1˜1K˜Kšœ
£$˜1Kšœ
œ
˜šœœ
œ˜"šœœœ˜>šœ
œ$œ˜0Kš
œ˜Kš
œc˜cKš
œ˜K˜
—šœ˜Kšœ
œœ
˜0K˜—K˜—Kšœ
˜—L˜K˜—Kšœ
˜—…—82[ά