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. f 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 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. Κβ•NewlineDelimiter –(cedarcode) style™codešΟn™Kšœύ™…Kšœ Οeœ7™BK™"K™—šΟk ˜ KšœH˜HK˜—š œŸœŸ˜KšŸœ˜#KšŸœ ˜—K˜šŸ˜KšœŸœΟc˜9KšœŸœ #˜AšœŸœ ˜,Kšœ ŸœŸœ#Ÿœ™KKš œŸœŸœ Ÿœ ŸœŸœ ™E—šœ Ÿœ˜'Kš œŸœŸœŸœŸœŸœŸœ™@—šœŸœ˜5Jš œŸœŸœ Ÿœ ŸœŸœŸœ™N—šŸœŸœ Ÿœ˜KšŸœŸœŸœŸœ™ ——K™šœŸœŸœŸœ˜FJšœˆ™ˆšœ ŸœŸœŸœ˜BKšœŸœ˜ šŸœŸœ˜&KšœŸœ5˜>šŸœ Ÿœ˜šœ Ÿœ 5˜TKšœP™P—KšœŸœ ˜KšœŸœ˜šŸœ Ÿ˜˜QK™"—KšŸœ Ÿœ ŸœŸœ˜%K˜KšŸœ˜—K˜—KšŸœ˜KšŸœ˜ —šŸœ˜Kšœ˜K˜Kš ŸœŸœΟb œŸœ‘ œ˜-K˜—K˜—š ŸœŸœŸœŸœŸœŸœŸ˜LK˜KšŸœ˜—Kšœ˜K˜—šœŸœŸœŸœŸœ Ÿœ Ÿœ Ÿœ ŸœŸœ˜’Kšœ9™9KšœR™Ršœ˜KšœŸœ&˜1KšŸœŸœ˜1K˜—KšœŸœ   ˜ K˜Kšœ˜K˜K˜—š œŸœŸœ ŸœŸœ˜pJ™€KšœŸœ˜!K˜šŸœŸœ˜KšŸœŸœ9˜DKšœ˜—šŸœ˜Kšœ ŸœŸœŸœ˜!Kš ŸœŸœŸœŸœŸœ˜8šœ˜KšŸ˜K˜CKšœŸœ˜Kš ŸœŸœŸœŸœŸœ˜8šŸœŸœ˜K˜/KšŸœŸœ+˜JšŸœŸœŸœŸ˜&K˜BKšŸœŸœE˜gKšŸœ˜—K˜>KšŸœ ŸœD˜SKšœ˜—šŸœ˜K˜KšŸœ Ÿœ*˜9K˜—KšŸ˜—Kšœ˜—K˜K˜—š œŸœ Ÿœ˜6KšŸœŸœ˜KšŸ˜Jšœ Ÿœ ŸœŸœ™!Jšœ Ÿœ™KšŸœ ˜KšœŸœ˜$KšœŸœ˜$KšœŸœ˜ Kš œŸœŸœŸœŸœ˜%Kš œŸœŸœŸœŸœ˜%KšœŸœŸœ˜KšœŸœŸœ˜š  œŸœŸœŸœŸœ˜0KšŸ˜KšœŸœ˜#KšœŸœ!˜,KšŸœŸœŸœŸœ˜/KšŸœŸœŸœŸœ˜/KšœŸœ ŸœŸœ ˜*KšŸœ˜—KšŸœŸœŸœ˜KšŸœŸœŸœ˜Kš Ÿœ Ÿœ Ÿœ Ÿœ ŸœŸœ ˜CK˜Jšœ™šŸœŸœ Ÿ˜KšœŸœ˜KšŸœŸœŸœ˜K˜K˜Jšœ™KšŸœŸœŸœ˜KšŸœ˜—Jšœ4™4JšœŸœ ™š ŸœŸœ ŸœŸœ Ÿ˜+KšŸ˜KšœŸœ ˜KšœŸœ ˜KšŸœ˜—K˜K˜šŸœ Ÿ˜KšœŸœ(˜/šŸœŸœŸ˜KšœŸœŸœŸœŸœŸœŸœŸœ˜7—Kš ŸœŸœ ŸœŸœ Ÿœ˜AKšŸœ˜KšŸœ Ÿœ Ÿœ Ÿœ ŸœŸœŸœŸœŸœŸœŸœŸœŸœŸœ˜{KšŸœ˜—KšŸœ˜KšŸœ˜K˜—šœŸœΟdœ’œ’œ’œ’œ’œ’œ’œŸœŸœŸœŸœ˜TKš œ ’œ’œ’œ’œ˜%Kš œ ’œ’œ’œ’œ˜%K˜—K˜KšœŸœ˜1š  œŸœŸœŸœŸœ˜^KšœW™WšŸœ Ÿœ˜(Kšœ Ÿœ$˜2K˜K˜K˜KšŸœ˜K˜—šŸœ˜Kšœ˜Kšœ Ÿœ˜Kšœ˜K˜*K˜'K˜'K˜*Kšœ Ÿœ.˜