<> <> <> <> <<>> DIRECTORY Cubic2, CubicPaths, GGCubic2, Polynomial, Real, Vector2; GGCubic2Impl: CEDAR PROGRAM IMPORTS Cubic2, Polynomial, Vector2 EXPORTS GGCubic2 = BEGIN Bezier: TYPE = Cubic2.Bezier; Path: TYPE = CubicPaths.Path; PointProc: TYPE = CubicPaths.PointProc; ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec; VEC: TYPE = Vector2.VEC; << [Artwork node; type 'ArtworkInterpress on' to command tool] >> Flat: PUBLIC PROC [bezier: Bezier, epsilon: REAL] RETURNS [BOOL] = { <> <> OPEN Vector2; dx,dy: REAL; d1,d2,d,bl,bh: VEC; oh: VEC=[epsilon, epsilon]; <> bh _ Add[UpperRight[bezier.b0,bezier.b3],oh]; bl _ Sub[LowerLeft[bezier.b0,bezier.b3],oh]; IF ~In[bezier.b1,bl,bh] OR ~In[bezier.b2,bl,bh] THEN RETURN[FALSE]; <> d _ Sub[bezier.b3,bezier.b0]; dx _ ABS[d.x]; dy _ ABS[d.y]; IF dx+dy < epsilon THEN RETURN[TRUE]; <> d1 _ Sub[bezier.b1,bezier.b0]; d2 _ Sub[bezier.b2,bezier.b0]; IF dy < dx THEN { <> dydx: REAL _ d.y/d.x; RETURN[ABS[d2.y-d2.x*dydx]> dxdy: REAL _ d.x/d.y; RETURN[ABS[d2.x-d2.y*dxdy]> OPEN Vector2; subdivide: PROC[bezier: Cubic2.Bezier] RETURNS [quit: BOOLEAN] = { quit _ FALSE; IF Flat[bezier, epsilon] THEN { len: REAL _ Length[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 _ Add[Mul[bezier.b0, 1-alpha], 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]; }; GetParam: PUBLIC PROC [bezier: Cubic2.Bezier, pt: VEC] RETURNS [u: REAL] = { p0, p1, p2, p3: VEC; a,b,c,d: REAL; ignoreX, ignoreY: BOOL _ FALSE; polyX, polyY: Polynomial.Ref; epsilon: REAL = 0.01; rootsX, rootsY: Polynomial.ShortRealRootRec; <> p0 _ bezier.b0; p1 _ bezier.b1; p2 _ bezier.b2; p3 _ bezier.b3; d _ p0.x - pt.x; c _ -3*p0.x + 3*p1.x; b _ 3*p0.x - 6*p1.x + 3*p2.x; a _ -p0.x + 3*p1.x - 3*p2.x + p3.x; IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreX _ TRUE ELSE { polyX _ Polynomial.Cubic[[d,c,b,a]]; rootsX _ Polynomial.CheapRealRoots[polyX]; }; d _ p0.y - pt.y; c _ -3*p0.y + 3*p1.y; b _ 3*p0.y - 6*p1.y + 3*p2.y; a _ -p0.y + 3*p1.y - 3*p2.y + p3.y; IF ABS[d] < epsilon AND ABS[c] < epsilon AND ABS[b] < epsilon AND ABS[a] < epsilon THEN ignoreY _ TRUE ELSE { polyY _ Polynomial.Cubic[[d,c,b,a]]; rootsY _ Polynomial.CheapRealRoots[polyY]; }; IF ignoreX THEN { IF ignoreY THEN ERROR; FOR j: NAT IN [0..rootsY.nRoots) DO IF rootsY.realRoot[j] > 0.0 AND rootsY.realRoot[j] < 1.0 THEN RETURN [rootsY.realRoot[j]]; ENDLOOP; } ELSE IF ignoreY THEN { FOR i: NAT IN [0..rootsX.nRoots) DO IF rootsX.realRoot[i] > 0.0 AND rootsX.realRoot[i] < 1.0 THEN RETURN [rootsX.realRoot[i]]; ENDLOOP; } ELSE { FOR i: NAT IN [0..rootsX.nRoots) DO FOR j: NAT IN [0..rootsY.nRoots) DO IF rootsX.realRoot[i] > 0.0 AND rootsX.realRoot[i] < 1.0 AND ABS[rootsX.realRoot[i] - rootsY.realRoot[j]] < epsilon THEN RETURN [(rootsX.realRoot[i] + rootsY.realRoot[j]) / 2.0]; ENDLOOP; ENDLOOP; ERROR; }; }; AlphaSplit: PUBLIC PROC [bezier: Bezier, alpha: REAL] RETURNS[Bezier,Bezier] = { Fraction: PROC[u,v: VEC, alpha: REAL, beta: REAL] RETURNS[VEC] = { RETURN[[u.x*beta+v.x*alpha, u.y*beta+v.y*alpha]] }; a,b,c,ab,bc,p: VEC; oneMinusAlpha: REAL; oneMinusAlpha _ 1.0-alpha; a _ Fraction[bezier.b0, bezier.b1, alpha, oneMinusAlpha]; b _ Fraction[bezier.b1, bezier.b2, alpha, oneMinusAlpha]; c _ Fraction[bezier.b2, bezier.b3, alpha, oneMinusAlpha]; ab _ Fraction[a, b, alpha, oneMinusAlpha]; bc _ Fraction[b, c, alpha, oneMinusAlpha]; p _ Fraction[ab, bc, alpha, oneMinusAlpha]; RETURN[[bezier.b0, a, ab, p],[p, bc, c, bezier.b3]]; }; ClosestPointAnalytic: PUBLIC PROC[pt: VEC, path: Path, tolerance: REAL _ 9999.0] RETURNS [closest: VEC, success: BOOL] = { <> <> <> <> <> bezier: Bezier; a3, a2, a1, a0, b3, b2, b1, b0: REAL; p0, p1, p2, p3: VEC; bestDist2, thisDist2: REAL; thisPt: VEC; thisSuccess: BOOL; bestDist2 _ Real.LargestNumber; success _ FALSE; FOR i: NAT IN [0..path.cubics.length) DO bezier _ path.cubics[i]; p0 _ bezier.b0; p1 _ bezier.b1; p2 _ bezier.b2; p3 _ bezier.b3; a3 _ -p0.x + 3*p1.x - 3*p2.x + p3.x; a2 _ 3*p0.x - 6*p1.x + 3*p2.x; a1 _ -3*p0.x + 3*p1.x; a0 _ p0.x; b3 _ -p0.y + 3*p1.y - 3*p2.y + p3.y; b2 _ 3*p0.y - 6*p1.y + 3*p2.y; b1 _ -3*p0.y + 3*p1.y; b0 _ p0.y; [thisPt, thisDist2, thisSuccess] _ ClosestPointCubic[a3, a2, a1, a0, b3, b2, b1, b0, pt]; IF NOT thisSuccess THEN ERROR; IF thisDist2 < bestDist2 THEN { bestDist2 _ thisDist2; closest _ thisPt; success _ TRUE; }; ENDLOOP; }; ClosestPointCubic: PROC[a3, a2, a1, a0, b3, b2, b1, b0: REAL, q: VEC] RETURNS [closest: VEC, bestD2: REAL, success: BOOL] = { <> <> <> <> <> <> <<>> < 1. We evaluate the spline at the remain values of u and compute (x - q.x)2 + (y - q.y)2. We also compute the distance from q to the spline endpoints at u = 0 and u = 1. We choose the best point.>> c5, c4, c3, c2, c1, c0: REAL; poly: Polynomial.Ref; roots: Polynomial.ShortRealRootRec; thisD2: REAL; thisPt: VEC; c5 _ 3*(a3*a3 + b3*b3); c4 _ 5*(a2*a3 + b2*b3); c3 _ 4*(a1*a3 + b1*b3) + 2*(a2*a2 + b2*b2); c2 _ 3*(a1*a2 + b1*b2) + 3*(a3*(a0 - q.x) + b3*(b0 - q.y)); c1 _ 2*(a2*(a0 - q.x) + b2*(b0 - q.y)) + (a1*a1 + b1*b1); c0 _ (a1*(a0 - q.x) + b1*(b0 - q.y)); poly _ Polynomial.Quintic[[c0, c1, c2, c3, c4, c5]]; roots _ CheapRealRootsInInterval[poly, 0.0, 1.0]; bestD2 _ Real.LargestNumber; success _ FALSE; FOR i: NAT IN [0..roots.nRoots) DO IF roots.realRoot[i] > 0.0 AND roots.realRoot[i] <= 1.0 THEN { thisPt _ Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, roots.realRoot[i]]; thisD2 _ Vector2.Square[Vector2.Sub[thisPt, q]]; IF thisD2 < bestD2 THEN { bestD2 _ thisD2; closest _ thisPt; success _ TRUE; }; }; ENDLOOP; FOR u: REAL _ 0, u+1.0 UNTIL u > 1.0 DO thisPt _ Evaluate[a3, a2, a1, a0, b3, b2, b1, b0, u]; thisD2 _ Vector2.Square[Vector2.Sub[thisPt, q]]; IF thisD2 < bestD2 THEN { bestD2 _ thisD2; closest _ thisPt; success _ TRUE; }; ENDLOOP; }; 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; }; END.