<> <> <> <> <<>> 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.