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