GGCubic2Impl.mesa
Copyright © 1986 by Xerox Corporation. All rights reserved.
Last edited by Bier on December 1, 1986 4:48:43 pm PST
Contents: The routines in Cubic2 and CubicPaths are inadequate for Gargoyle purposes. The Closest Point routine is too slow and doesn't properly document its accuracy. There is also a need for routines that are slow but very accurate so we can test how good our fast routines are (maybe using rational arithmetic?). Hopefully, these routines will eventually find their way back into Cubic2, and CubicPaths.
DIRECTORY
Cubic2, CubicPaths, GGCubic2, Polynomial, Real, Vector2;
GGCubic2Impl: CEDAR PROGRAM
IMPORTS Cubic2, Polynomial, Vector2
EXPORTS GGCubic2 =
BEGIN
Bezier: TYPE = Cubic2.Bezier;
Path: TYPE = CubicPaths.Path;
PointProc: TYPE = CubicPaths.PointProc;
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
VEC: TYPE = Vector2.VEC;
[Artwork node; type 'ArtworkInterpress on' to command tool]
Flat: PUBLIC PROC [bezier: Bezier, epsilon: REAL] RETURNS [BOOL] = {
Trivial Rejection 1: If b1 or b2 are more than epsilon from the bounding box of b0 and b3, then we assume the curve is not flat.
Trivial Rejection 2: If b0 is within epsilon of b3 then assume the curve is flat. This relies on the assumption that the curve is not a very high frequency one in this region.
OPEN Vector2;
dx,dy: REAL;
d1,d2,d,bl,bh: VEC;
oh: VEC=[epsilon, epsilon];
Trivial Rejection 1
bh ← Add[UpperRight[bezier.b0,bezier.b3],oh];
bl ← Sub[LowerLeft[bezier.b0,bezier.b3],oh];
IF ~In[bezier.b1,bl,bh] OR ~In[bezier.b2,bl,bh] THEN RETURN[FALSE];
Trivial Rejection 2
d ← Sub[bezier.b3,bezier.b0];
dx ← ABS[d.x]; dy ← ABS[d.y];
IF dx+dy < epsilon THEN RETURN[TRUE];
Compute Pseudo-Altitudes, A and B
d1 ← Sub[bezier.b1,bezier.b0];
d2 ← Sub[bezier.b2,bezier.b0];
IF dy < dx THEN {
Compute Vertical Altitudes as shown in the figure.
dydx: REAL ← d.y/d.x;
RETURN[ABS[d2.y-d2.x*dydx]<epsilon AND ABS[d1.y-d1.x*dydx]<epsilon] }
ELSE {
Compute Horizontal Altitudes (not shown).
dxdy: REAL ← d.x/d.y;
RETURN[ABS[d2.x-d2.y*dxdy]<epsilon AND ABS[d1.x-d1.y*dxdy]<epsilon] };
};
LowerLeft: PROC[a: VEC, b: VEC] RETURNS[VEC] = INLINE { RETURN[[MIN[a.x,b.x],MIN[a.y,b.y]]] };
UpperRight: PROC[a: VEC, b: VEC] RETURNS[VEC] = INLINE { RETURN[[MAX[a.x,b.x],MAX[a.y,b.y]]] };
In: PROC[a: VEC, b: VEC, c: VEC] RETURNS[BOOLEAN] = INLINE { RETURN[a.x IN[b.x..c.x] AND a.y IN[b.y..c.y]] };
WalkPath: PUBLIC PROC[path: Path, epsilon: REAL, proc: PointProc] = {
Walk along the curve in steps of size epsilon or less, calling proc. Where possible, approximate the curve with straight line segments.
OPEN Vector2;
subdivide: PROC[bezier: Cubic2.Bezier] RETURNS [quit: BOOLEAN] = {
quit ← FALSE;
IF Flat[bezier, epsilon] THEN {
len: REAL ← Length[Sub[bezier.b0, bezier.b3]];
IF len > 0 THEN {
alphaStep: REAL ← epsilon/len; -- what fraction of the curve to advance at each step
Choose this value so that each step (= alphaStep * len) will be of size epsilon.
alpha: REAL ← alphaStep;
pt: VEC;
UNTIL alpha > 1 DO
pt ← Add[Mul[bezier.b0, 1-alpha], Mul[bezier.b3, alpha]];
Go alpha of the way from b0 to b3.
IF proc[pt] THEN {quit ← TRUE; EXIT};
alpha ← alpha+alphaStep;
ENDLOOP;
}
ELSE quit ← proc[bezier.b3];
RETURN[quit]}
ELSE {
b1, b2: 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;
};
ClosestPointSubdivide: PUBLIC PROC[pt: VEC, path: Path, epsilon: REAL, tolerance: REAL ← 9999.0] RETURNS [closest: VEC, success: BOOL ← TRUE] = {
This is very inefficient because it uses WalkPath. A new version should be written that prunes the search.
test: PointProc = {
thisD: REAL ← Vector2.Square[Vector2.Sub[p, pt]];
IF thisD < minD THEN {closest ← p; minD ← thisD};
};
minD: REAL ← 10E6; -- "infinity"
closest ← path.cubics[0].b0;
WalkPath[path, epsilon, test];
};
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: Polynomial.ShortRealRootRec;
ShortRealRootRec: TYPE = RECORD [nRoots: NAT, realRoot: ARRAY [0..5) OF REAL];
p0 ← bezier.b0;
p1 ← bezier.b1;
p2 ← bezier.b2;
p3 ← bezier.b3;
d ← p0.x - pt.x;
c ← -3*p0.x + 3*p1.x;
b ← 3*p0.x - 6*p1.x + 3*p2.x;
a ← -p0.x + 3*p1.x - 3*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 ← Polynomial.CheapRealRoots[polyX];
};
d ← p0.y - pt.y;
c ← -3*p0.y + 3*p1.y;
b ← 3*p0.y - 6*p1.y + 3*p2.y;
a ← -p0.y + 3*p1.y - 3*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 ← Polynomial.CheapRealRoots[polyY];
};
IF ignoreX THEN {
IF ignoreY THEN ERROR;
FOR j: NAT IN [0..rootsY.nRoots) DO
IF rootsY.realRoot[j] > 0.0 AND rootsY.realRoot[j] < 1.0 THEN RETURN [rootsY.realRoot[j]];
ENDLOOP;
}
ELSE IF ignoreY THEN {
FOR i: NAT IN [0..rootsX.nRoots) DO
IF rootsX.realRoot[i] > 0.0 AND rootsX.realRoot[i] < 1.0 THEN RETURN [rootsX.realRoot[i]];
ENDLOOP;
}
ELSE {
FOR i: NAT IN [0..rootsX.nRoots) DO
FOR j: NAT IN [0..rootsY.nRoots) DO
IF rootsX.realRoot[i] > 0.0
AND rootsX.realRoot[i] < 1.0 AND ABS[rootsX.realRoot[i] - rootsY.realRoot[j]] < epsilon THEN
RETURN [(rootsX.realRoot[i] + rootsY.realRoot[j]) / 2.0];
ENDLOOP;
ENDLOOP;
ERROR;
};
};
AlphaSplit: PUBLIC PROC [bezier: Bezier, alpha: REAL] RETURNS[Bezier,Bezier] = {
Fraction: PROC[u,v: VEC, alpha: REAL, beta: REAL] RETURNS[VEC] = {
RETURN[[u.x*beta+v.x*alpha, u.y*beta+v.y*alpha]] };
a,b,c,ab,bc,p: VEC;
oneMinusAlpha: REAL;
oneMinusAlpha ← 1.0-alpha;
a ← Fraction[bezier.b0, bezier.b1, alpha, oneMinusAlpha];
b ← Fraction[bezier.b1, bezier.b2, alpha, oneMinusAlpha];
c ← Fraction[bezier.b2, bezier.b3, alpha, oneMinusAlpha];
ab ← Fraction[a, b, alpha, oneMinusAlpha];
bc ← Fraction[b, c, alpha, oneMinusAlpha];
p ← Fraction[ab, bc, alpha, oneMinusAlpha];
RETURN[[bezier.b0, a, ab, p],[p, bc, c, bezier.b3]];
};
ClosestPointAnalytic: PUBLIC PROC[pt: VEC, path: Path, tolerance: REAL ← 9999.0] RETURNS [closest: VEC, success: BOOL] = {
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 d u d 1.
bezier: Bezier;
a3, a2, a1, a0, b3, b2, b1, b0: REAL;
p0, p1, p2, p3: VEC;
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];
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;
[thisPt, thisDist2, thisSuccess] ← ClosestPointCubic[a3, a2, a1, a0, b3, b2, b1, b0, 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] = {
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 d u d 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
x/u(x - q.x) + y/u(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.
c5, c4, c3, c2, c1, c0: REAL;
poly: Polynomial.Ref;
roots: Polynomial.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;
};
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
debug xi: ARRAY [0..50) OF REAL;
debug i: NAT ← 0;
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;
first try Newton's method
WHILE x IN [x0..x1] DO
newx:REAL ← NewtonStep[x];
IF x=newx THEN RETURN;
x ← NewtonStep[newx];
xx ← NewtonStep[xx];
debug xi[i] ← xx; i ← i+1;
IF xx=x THEN EXIT;
ENDLOOP;
If it oscillated or went out of range, try bisection
debug SIGNAL Bisection;
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;
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;
};
END.