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 PointProc: TYPE = CubicPaths.PointProc; ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec; VEC: TYPE = Vector2.VEC; WalkPath: PUBLIC PROC [path: Path, epsilon: REAL, proc: PointProc] = { 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 alpha: REAL _ alphaStep; pt: VEC; UNTIL alpha > 1 DO pt _ Vector2.Add[Vector2.Mul[bezier.b0, 1-alpha], Vector2.Mul[bezier.b3, alpha]]; 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] = { 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] = { 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 y0yy0 THEN {xx0 _ x; yy0 _ y}}; IF y>0 THEN {IF y 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; 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] = { 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. `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 PathRec: TYPE = RECORD[bounds: Imager.Rectangle, cubics: REF PathSequence]; PathSequence: TYPE = RECORD[bezier: SEQUENCE length: NAT OF Bezier]; PointProc: TYPE = PROC [p: VEC] RETURNS [stop: BOOLEAN _ FALSE]; ShortRealRootRec: TYPE = RECORD [nRoots: NAT, realRoot: ARRAY [0..5) OF REAL]; VEC: TYPE ~ RECORD [x, y: REAL]; Walk along the curve in steps of size epsilon or less, calling proc. Where possible, approximate the curve with straight line segments. Choose this value so that each step (= alphaStep * len) will be of size epsilon. Go alpha of the way from b0 to b3. Finds the closest point of path to pt using sub-division. THIS ROUTINE IS VERY SLOW. Use it only to help in debugging ClosestPointAnalytic. 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; Approximate bezier with a polyline that deviates by no more than tolerance from bezier. Κ ˜IcodešΟn™šœύ™…K™šžœ žœ˜šœ žœŸ5˜TKšœP™P—Kšœžœ ˜Kšœžœ˜šžœ ž˜šœQ˜QK™"—Kšžœ žœ žœžœ˜%Kšœ˜Kšžœ˜—K˜—Kšžœ˜Jšžœ˜ —šžœ˜Jšœ˜J˜Jš žœžœΟb œžœ  œ˜-J˜—J˜—š žœžœžœžœžœžœž˜LKšœ˜Kšžœ˜—Kšœ˜K˜—šœžœžœžœžœ žœ žœ žœ žœžœ˜’Kšœ9™9KšœR™Ršœ˜Kšœžœ&˜1Kšžœžœ˜1K˜—Kšœžœ Ÿ ˜ Kšœ˜Kšœ˜K˜K˜—š œžœžœ žœžœ˜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˜—šœžœΟdœ‘œ‘œ‘œ‘œ‘œ‘œ‘œžœžœžœžœ˜TJš œ ‘œ‘œ‘œ‘œ˜%Jš œ ‘œ‘œ‘œ‘œ˜%J˜—J˜Jšœžœ˜1š  œžœžœžœžœ˜^KšœW™Wšžœ žœ˜(Kšœ žœ$˜2Kšœ˜Kšœ˜Kšœ˜Kšžœ˜K˜—šžœ˜Kšœ˜Kšœ žœ˜Kšœ˜Kšœ*˜*Kšœ'˜'Kšœ'˜'Kšœ*˜*Kšœ žœ.˜