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, 1992 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 [a
3, a
2, a
1, a
0, b
3, b
2, b
1, b
0:
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.