<> <> <> <> <> <> <> <<>> DIRECTORY Cubic2, CubicPaths, CubicPathsExtras, CubicSplines, Imager, ImagerPath, ImagerTransformation, Polynomial, Real, RealFns, Vector2; CubicPathsImpl: CEDAR PROGRAM IMPORTS Cubic2, ImagerTransformation, Polynomial, RealFns, Vector2 EXPORTS CubicPaths, CubicPathsExtras ~ BEGIN OPEN CubicPaths; Bezier: TYPE = Cubic2.Bezier; Coeffs: TYPE = Cubic2.Coeffs; ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec; VEC: TYPE = Imager.VEC; X: NAT = CubicSplines.X; Y: NAT = CubicSplines.Y; <<>> PathFromCubic: PUBLIC PROC [coeffs: CubicSplines.CoeffsSequence, newBounds: BOOL _ TRUE] RETURNS [Path] = { path: Path _ NEW[PathRec]; path.cubics _ NEW[PathSequence[coeffs.length]]; FOR i: NAT IN [0..coeffs.length) DO cfs: Cubic2.Coeffs _ [ c3: [coeffs[i].t3[X], coeffs[i].t3[Y]], c2: [coeffs[i].t2[X], coeffs[i].t2[Y]], c1: [coeffs[i].t1[X], coeffs[i].t1[Y]], c0: [coeffs[i].t0[X], coeffs[i].t0[Y]]]; path.cubics[i] _ Cubic2.CoeffsToBezier[cfs]; ENDLOOP; IF newBounds THEN UpdateBounds[path]; RETURN [path]; }; EnumeratePath: PUBLIC PROC [path: Path, moveTo: ImagerPath.MoveToProc, curveTo: ImagerPath.CurveToProc] = { moveTo[path.cubics[0].b0]; FOR i: NAT IN [0..path.cubics.length) DO curveTo[path.cubics[i].b1, path.cubics[i].b2, path.cubics[i].b3]; ENDLOOP; }; WalkPath: PUBLIC PROC [path: Path, tol: REAL, proc: PointProc] = { subdivide: PROC [bezier: Cubic2.Bezier] RETURNS [quit: BOOLEAN]= { quit _ FALSE; IF Cubic2.Flat[bezier, tol] THEN { len: REAL _ Vector2.Length[Vector2.Sub[bezier.b0, bezier.b3]]; IF len > 0 THEN { step: REAL _ tol/len; alpha: REAL _ step; pt: VEC; UNTIL alpha > 1 DO pt _ Vector2.Add[Vector2.Mul[bezier.b3,alpha],Vector2.Mul[bezier.b0,1-alpha]]; IF proc[pt] THEN {quit _ TRUE; EXIT}; alpha _ alpha+step; ENDLOOP; } ELSE quit _ proc[bezier.b3]; RETURN [quit]} ELSE { b1, b2: Cubic2.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; }; PointOnPath: PUBLIC PROC [pt: VEC, path: Path, tol: REAL] RETURNS [BOOLEAN] = { tolSq: REAL _ tol*tol; found: BOOLEAN _ FALSE; hit: PointProc = { IF Vector2.Square[Vector2.Sub[p, pt]] <= tolSq THEN found _ TRUE; RETURN [found]; }; IF pt.x < path.bounds.x OR pt.y < path.bounds.y OR pt.x>(path.bounds.x+path.bounds.w) OR pt.y > (path.bounds.y+path.bounds.h) THEN RETURN [FALSE] ELSE WalkPath[path, tol, hit]; RETURN [found]; }; TranslatePath: PUBLIC PROC [path: Path, offset: VEC] = { FOR i: NAT IN [0..path.cubics.length) DO path.cubics[i].b0 _ Vector2.Add[offset, path.cubics[i].b0]; path.cubics[i].b1 _ Vector2.Add[offset, path.cubics[i].b1]; path.cubics[i].b2 _ Vector2.Add[offset, path.cubics[i].b2]; path.cubics[i].b3 _ Vector2.Add[offset, path.cubics[i].b3]; ENDLOOP; path.bounds.x _ path.bounds.x+offset.x; path.bounds.y _ path.bounds.y+offset.y; }; TransformPath: PUBLIC PROC [path: Path, tranformation: ImagerTransformation.Transformation, newBounds: BOOL _ TRUE] = { FOR i: NAT IN [0..path.cubics.length) DO path.cubics[i].b0 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b0]; path.cubics[i].b1 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b1]; path.cubics[i].b2 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b2]; path.cubics[i].b3 _ ImagerTransformation.Transform[tranformation, path.cubics[i].b3]; ENDLOOP; IF newBounds THEN UpdateBounds[path]; }; < maxX THEN maxX _ p.x; IF p.y > maxY THEN maxY _ p.y; IF p.x < minX THEN minX _ p.x; IF p.y < minY THEN minY _ p.y; }; WalkPath[path, 1.0, pointProc]; path.bounds _ [x: minX, y: minY, w: maxX-minX, h: maxY-minY]; }; >> 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; }; EvaluateX: PROC [coeffs: Cubic2.Coeffs, t: REAL] RETURNS [x: REAL] = { x _ (((coeffs.c3.x)*t + coeffs.c2.x)*t + coeffs.c1.x)*t + coeffs.c0.x; }; EvaluateY: PROC [coeffs: Cubic2.Coeffs, t: REAL] RETURNS [y: REAL] = { y _ (((coeffs.c3.y)*t + coeffs.c2.y)*t + coeffs.c1.y)*t + coeffs.c0.y; }; UpdateBounds: PUBLIC PROC [path: Path] = { <> minX, minY: REAL _ 10E6; maxX, maxY: REAL _ -10E6; loEndX, hiEndX, loEndY, hiEndY: REAL; coeffs: Cubic2.Coeffs; tArray: ARRAY [1..2] OF REAL; ptX, ptY: REAL; rootCount: NAT; FOR i: NAT IN [0..path.cubics.length) DO coeffs _ Cubic2.BezierToCoeffs[path.cubics[i]]; <> [tArray, rootCount] _ QuadraticFormula[3.0*coeffs.c3.x, 2.0*coeffs.c2.x, coeffs.c1.x]; FOR j: NAT IN [1..rootCount] DO IF tArray[j] >= 0.0 AND tArray[j] <= 1.0 THEN { ptX _ EvaluateX[coeffs, tArray[j]]; minX _ MIN[minX, ptX]; maxX _ MAX[maxX, ptX]; }; ENDLOOP; <> [tArray, rootCount] _ QuadraticFormula[3.0*coeffs.c3.y, 2.0*coeffs.c2.y, coeffs.c1.y]; FOR j: NAT IN [1..rootCount] DO IF tArray[j] >= 0.0 AND tArray[j] <= 1.0 THEN { ptY _ EvaluateY[coeffs, tArray[j]]; minY _ MIN[minY, ptY]; maxY _ MAX[maxY, ptY]; }; ENDLOOP; <> loEndX _ path.cubics[i].b0.x; loEndY _ path.cubics[i].b0.y; hiEndX _ path.cubics[i].b3.x; hiEndY _ path.cubics[i].b3.y; minX _ MIN[minX, loEndX, hiEndX]; minY _ MIN[minY, loEndY, hiEndY]; maxX _ MAX[maxX, loEndX, hiEndX]; maxY _ MAX[maxY, loEndY, hiEndY]; ENDLOOP; path.bounds _ [x: minX, y: minY, w: maxX-minX, h: maxY-minY]; }; <> ClosestPoint: PUBLIC PROC [pt: VEC, path: Path] RETURNS [VEC] = { closest: VEC _ path.cubics[0].b0; minD: REAL _ 10E6; test: PointProc = { thisD: REAL _ Vector2.Square[Vector2.Sub[p, pt]]; IF thisD < minD THEN {closest _ p; minD _ thisD}; }; WalkPath[path, 1.0, test]; RETURN [closest]; }; ClosestPointAnalytic: PUBLIC PROC [pt: VEC, path: Path, tolerance: REAL _ 9999.0] RETURNS [closest: VEC, success: BOOL] = { <> <> <> <> <> bezier: Bezier; c: Coeffs; <> <> 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]; <> <> <> <> <> <> <> <> <> c _ Cubic2.BezierToCoeffs[bezier]; [thisPt, thisDist2, thisSuccess] _ ClosestPointCubic[c.c3.x, c.c2.x, c.c1.x, c.c0.x, c.c3.y, c.c2.y, c.c1.y, c.c0.y, 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: 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; }; CopyPath: PUBLIC PROC [path: Path] RETURNS [Path] = { new: Path _ NEW[PathRec]; new.cubics _ NEW[PathSequence[path.cubics.length]]; FOR i: NAT IN [0..path.cubics.length) DO new.cubics[i] _ path.cubics[i]; ENDLOOP; new.bounds _ path.bounds; RETURN [new]; }; ReversePath: PUBLIC PROC [path: Path] = { <> last: NAT _ path.cubics.length-1; FOR i: NAT IN [0..path.cubics.length/2) DO temp: Cubic2.Bezier _ path.cubics[i]; path.cubics[i] _ path.cubics[last-i]; path.cubics[last-i] _ temp; ENDLOOP; FOR i: NAT IN [0..path.cubics.length) DO temp: VEC _ path.cubics[i].b0; path.cubics[i].b0 _ path.cubics[i].b3; path.cubics[i].b3 _ temp; temp _ path.cubics[i].b1; path.cubics[i].b1 _ path.cubics[i].b2; path.cubics[i].b2 _ temp; ENDLOOP; }; <<>> QuadraticFormula: PUBLIC PROC [a, b, c: REAL] RETURNS [roots: ARRAY [1..2] OF REAL, rootCount: NAT] = { <> discriminant, temp: REAL; IF a = 0 THEN { [roots[1], rootCount] _ LinearFormula [b, c]; RETURN; }; IF c = 0 THEN { roots[1] _ 0.0; [roots[2], rootCount] _ LinearFormula[a, b]; rootCount _ rootCount + 1; IF roots[1] > roots[2] THEN { -- swap them temp _ roots[1]; roots[1] _ roots[2]; roots[2] _ temp; }; RETURN}; discriminant _ b*b - 4*a*c; -- 3 mult, 1 add SELECT discriminant FROM =0 => {rootCount _ 1; roots[1] _ -b/(2*a); RETURN}; -- 1 mult, 1 div (4 mult, 1 div, 1 add total) <0 => {rootCount _ 0; RETURN}; -- (3 mult, 1 add total) >0 => { sqrtRadical: REAL _ RealFns.SqRt[discriminant]; term: REAL; rootCount _ 2; term _ IF b < 0 THEN sqrtRadical - b ELSE -sqrtRadical - b; roots[1] _ term/(2.0*a); roots[2] _ (2.0*c)/term; IF roots[1] > roots[2] THEN { -- swap them temp _ roots[1]; roots[1] _ roots[2]; roots[2] _ temp; }; RETURN}; -- 1 SqRt, 2 mult, 2 div, 2 add (1 SqRt, 5 mult, 2 div, 3 add total) ENDCASE => ERROR; }; -- end of QuadraticFormula LinearFormula: PROC [a, b: REAL] RETURNS [root: REAL, rootCount: NAT] = { <> IF a = 0 THEN { IF b = 0 THEN {root _ 0.0; rootCount _ 0; RETURN} ELSE {root _ 0.0; rootCount _ 0; RETURN} -- inconsistent case } ELSE { IF b = 0 THEN {root _ 0.0; rootCount _ 1; RETURN} ELSE {root _ -b/a; rootCount _ 1; RETURN}; }; }; 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; 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: ShortRealRootRec; <> p0 _ bezier.b0; p1 _ bezier.b1; p2 _ bezier.b2; p3 _ bezier.b3; d _ p0.x - pt.x; c _ 3*(p1.x - p0.x); b _ 3*p0.x - 6*p1.x + 3*p2.x; a _ -p0.x + 3*(p1.x - 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 _ CheapRealRootsInInterval[polyX, 0.0, 1.0]; }; d _ p0.y - pt.y; c _ 3*(p1.y - p0.y); b _ 3*p0.y - 6*p1.y + 3*p2.y; a _ -p0.y + 3*(p1.y - 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 _ CheapRealRootsInInterval[polyY, 0.0, 1.0]; }; IF ignoreX THEN { IF ignoreY THEN ERROR; RETURN [rootsY.realRoot[0]]; } ELSE IF ignoreY THEN { RETURN [rootsX.realRoot[0]]; } ELSE { <> FOR i: NAT IN [0..rootsX.nRoots) DO FOR j: NAT IN [0..rootsY.nRoots) DO IF ABS[rootsX.realRoot[i] - rootsY.realRoot[j]] < epsilon THEN RETURN [(rootsX.realRoot[i] + rootsY.realRoot[j]) / 2.0]; ENDLOOP; ENDLOOP; ERROR; }; }; CubicMeetsLine: PUBLIC PROC [bezier: Bezier, a, b, c: REAL] RETURNS [points: ARRAY [0..2] OF VEC, hitCount: [0..3], tangency: ARRAY [0..2] OF BOOL] = { <> <> <> <> <> <> <> <> <> <<(a*a3 + b*b3)*t3 + (a*a2 + b*b2)*t2 + (a*a1 + b*b1)*t + (a*a0 + b*b0 + c) = 0.>> f: Coeffs; poly: Polynomial.Ref; roots: ShortRealRootRec; lastT: REAL; success: BOOL; epsilon: REAL = 0.0001; hitCount _ 0; tangency[0] _ tangency[1] _ tangency[2] _ FALSE; f _ Cubic2.BezierToCoeffs[bezier]; poly _ Polynomial.Cubic[[a*f.c0.x + b*f.c0.y + c, a*f.c1.x + b*f.c1.y, a*f.c2.x + b*f.c2.y, a*f.c3.x + b*f.c3.y]]; roots _ CheapRealRootsInInterval[poly, 0.0, 1.0]; lastT _ 2.0; -- a value outside the range [0..1]; success _ FALSE; FOR i: NAT IN [0..roots.nRoots) DO IF roots.realRoot[i] > 0.0 AND roots.realRoot[i] <= 1.0 THEN { IF ABS[roots.realRoot[i]-lastT] > epsilon THEN { lastT _ roots.realRoot[i]; points[hitCount] _ Evaluate[f.c3.x, f.c2.x, f.c1.x, f.c0.x, f.c3.y, f.c2.y, f.c1.y, f.c0.y, lastT]; hitCount _ hitCount + 1; } ELSE { IF lastT # 2.0 THEN tangency[hitCount-1] _ TRUE; }; }; ENDLOOP; }; END.