GGCubic2Impl.mesa
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. Hopefully, these routines will eventually find their way back into Cubic2, and CubicPaths.
Copyright Ó 1988 by Xerox Corporation. All rights reserved.
Bier, May 11, 1989 11:04:19 pm PDT
DIRECTORY
Cubic2, CubicPaths, GGCubic2, GGSegmentTypes, Polynomial, Real, Vector2;
GGCubic2Impl: CEDAR PROGRAM
IMPORTS Cubic2, Polynomial, Vector2
EXPORTS GGCubic2 =
BEGIN
Bezier: TYPE = Cubic2.Bezier; -- RECORD[b0,b1,b2,b3: VEC]
Coeffs: TYPE = Cubic2.Coeffs; -- RECORD[c0,c1,c2,c3: Imager.VEC];
Path: TYPE = CubicPaths.Path; -- REF PathRec
PathRec: TYPE = RECORD[bounds: Imager.Rectangle, cubics: REF PathSequence];
PathSequence: TYPE = RECORD[bezier: SEQUENCE length: NAT OF Bezier];
PointProc: TYPE = CubicPaths.PointProc;
PointProc: TYPE = PROC [p: VEC] RETURNS [stop: BOOLEAN ← FALSE];
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
ShortRealRootRec: TYPE = RECORD [nRoots: NAT, realRoot: ARRAY [0..5) OF REAL];
VEC: TYPE = Vector2.VEC;
VEC: TYPE ~ RECORD [x, y: REAL];
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.
subdivide: PROC[bezier: Cubic2.Bezier] RETURNS [quit: BOOLEAN] = {
quit ← FALSE;
IF Cubic2.Flat[bezier, epsilon] THEN {
len: REAL ← Vector2.Length[Vector2.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 ← Vector2.Add[Vector2.Mul[bezier.b0, 1-alpha], Vector2.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] = {
Finds the closest point of path to pt using sub-division.
THIS ROUTINE IS VERY SLOW. Use it only to help in debugging ClosestPointAnalytic.
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];
};
CheapRealRootsInInterval: PUBLIC PROC [poly: Polynomial.Ref, lo, hi: REAL] RETURNS [roots: ShortRealRootRec] = {
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.
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;
};
PairSequence: TYPE ~ GGSegmentTypes.PairSequence;
AsPolyline: PUBLIC PROC [bezier: Bezier, tolerance: REAL] RETURNS [polyline: PairSequence] = {
Approximate bezier with a polyline that deviates by no more than tolerance from bezier.
IF Cubic2.Flat[bezier, tolerance] THEN {
polyline ← NEW[GGSegmentTypes.PairSequenceRep[2]];
polyline.length ← 2;
polyline[0] ← bezier.b0;
polyline[1] ← bezier.b3;
RETURN;
}
ELSE {
bezier1, bezier2: Bezier;
totalLength: NAT;
poly1, poly2: PairSequence;
[bezier1, bezier2] ← Cubic2.Split[bezier];
poly1 ← AsPolyline[bezier1, tolerance];
poly2 ← AsPolyline[bezier2, tolerance];
totalLength ← poly1.length+poly2.length-1;
polyline ← NEW[GGSegmentTypes.PairSequenceRep[totalLength]];
polyline.length ← totalLength;
FOR i: NAT IN [0..poly1.length) DO
polyline[i] ← poly1[i];
ENDLOOP;
FOR i: NAT IN [1..poly2.length) DO -- skip the first vertex of poly2
polyline[i-1+poly1.length] ← poly2[i];
ENDLOOP;
};
};
END.