<> <> <> DIRECTORY ImagerBasic USING [Pair], ImagerConic USING [], Real USING [SqRt], RealFns USING [ArcTanDeg, CosDeg, SinDeg]; ImagerConicImpl: CEDAR PROGRAM IMPORTS Real, RealFns EXPORTS ImagerConic ~ BEGIN Pair: TYPE ~ ImagerBasic.Pair; Test: PROC [p0, p1, p2: Pair, r: REAL] RETURNS [list: LIST OF RECORD[p1, p2, p3: Pair]] ~ { Curve: PROC [p1, p2, p3: Pair] ~ { list _ CONS[[p1, p2, p3], list]; }; ToCurves[p0, p1, p2, r, Curve]; }; ArcTest: PROC [p0, p1, p2: Pair] RETURNS [list: LIST OF RECORD[p1, p2: Pair, r: REAL]] ~ { Conic: PROC [p1, p2: Pair, r: REAL] ~ { list _ CONS[[p1, p2, r], list]; }; FromArc[p0, p1, p2, Conic]; }; threshold: REAL _ 0.01; fourThirds: REAL ~ (4.0/3.0); ToCurves: PUBLIC PROC [p0, p1, p2: Pair, r: REAL, curve: PROC [p1, p2, p3: Pair]] ~ { IF NOT (r IN [0.0..1.0]) THEN ERROR; IF ABS[r-0.5]<=threshold THEN { f: REAL ~ fourThirds*r; p01: Pair ~ [p0.x+f*(p1.x-p0.x), p0.y+f*(p1.y-p0.y)]; p21: Pair ~ [p2.x+f*(p1.x-p2.x), p2.y+f*(p1.y-p2.y)]; curve[p01, p21, p2]; } ELSE { m: Pair ~ [(p0.x+p2.x)*0.5, (p0.y+p2.y)*0.5]; p: Pair ~ [m.x+r*(p1.x-m.x), m.y+r*(p1.y-m.y)]; p01: Pair ~ [p0.x+r*(p1.x-p0.x), p0.y+r*(p1.y-p0.y)]; p21: Pair ~ [p2.x+r*(p1.x-p2.x), p2.y+r*(p1.y-p2.y)]; rNew: REAL ~ MIN[1.0/(1.0+Real.SqRt[2.0*(1-r)]), 0.99999]; ToCurves[p0, p01, p, rNew, curve]; ToCurves[p, p21, p2, rNew, curve]; }; }; Add: PROC [a: Pair, b: Pair] RETURNS [Pair] ~ { RETURN[[a.x+b.x, a.y+b.y]] }; Sub: PROC [a: Pair, b: Pair] RETURNS [Pair] ~ { RETURN[[a.x-b.x, a.y-b.y]] }; Sqr: PROC [r: REAL] RETURNS [REAL] ~ { RETURN[r*r] }; DistSqr: PROC [a: Pair, b: Pair] RETURNS [REAL] ~ { RETURN[Sqr[a.x-b.x]+Sqr[a.y-b.y]] }; Dot: PROC [a: Pair, b: Pair] RETURNS [REAL] ~ { RETURN[a.x*b.x+a.y*b.y] }; Cross: PROC [a: Pair, b: Pair] RETURNS [REAL] ~ { RETURN[a.x*b.y-a.y*b.x] }; FromArc: PUBLIC PROC [p0, p1, p2: Pair, conic: PROC [p1, p2: Pair, r: REAL]] ~ { IF p0 = p2 THEN { c: Pair ~ [(p0.x+p1.x)*0.5, (p0.y+p1.y)*0.5]; <> v: Pair ~ [p0.x-c.x, p0.y-c.y]; <> p: Pair ~ [-v.y, v.x]; <> r: REAL ~ Real.SqRt[2.0]-1.0; <> conic[Add[p0,p], Add[c,p], r]; conic[Add[p1,p], p1, r]; conic[Sub[p1,p], Sub[c,p], r]; conic[Sub[p2,p], p2, r]; } ELSE { sgn: REAL ~ IF Cross[Sub[p1,p0], Sub[p2,p0]] >= 0 THEN 1.0 ELSE -1.0; <> distProd: REAL ~ Real.SqRt[DistSqr[p0,p1]*DistSqr[p2,p1]]; IF distProd = 0 THEN { <> conic[p1,p2,0.5]; } ELSE { dot: REAL ~ Dot[Sub[p0,p1], Sub[p2,p1]]; cos: REAL ~ -dot/distProd; <> IF cos>=0.5 THEN { tan: REAL ~ Cross[Sub[p0,p1], Sub[p2,p1]]/dot; m: Pair ~ [(p0.x+p2.x)*0.5, (p0.y+p2.y)*0.5]; p: Pair ~ [tan*(m.y-p0.y)+m.x, tan*(p0.x-m.x)+m.y]; r: REAL ~ cos/(1.0+cos); conic[p, p2, r]; } ELSE { m1: Pair ~ [(p0.x+p1.x)*0.5, (p0.y+p1.y)*0.5]; <> v1: Pair ~ [m1.y-p0.y, p0.x-m1.x]; <> a1: REAL ~ v1.y; b1: REAL ~ -v1.x; c1: REAL ~ a1*m1.x+b1*m1.y; <> m2: Pair ~ [(p0.x+p2.x)*0.5, (p0.y+p2.y)*0.5]; v2: Pair ~ [m2.y-p0.y, p0.x-m2.x]; a2: REAL ~ v2.y; b2: REAL ~ -v2.x; c2: REAL ~ a2*m2.x+b2*m2.y; det: REAL ~ a1*b2-a2*b1; Center: PROC RETURNS [Pair] ~ INLINE { <<*** Possible side effect on p1 in degenerate case>> IF det = 0 OR ABS[det] < ABS[a1*b2]*0.000005 THEN { <> v: Pair _ Add[v1, v2]; m: REAL; IF DistSqr[v1, v2] > v.x*v.x + v.y*v.y THEN v _ Sub[v1, v2]; m _ 1.0E+9/Real.SqRt[v.x*v.x + v.y*v.y]; <> p1 _ [2.0*m*v.x, 2.0*m*v.y]; RETURN [[m*v.x, m*v.y]]; } ELSE RETURN [[(c1*b2-c2*b1)/det, (a1*c2-a2*c1)/det]]; }; center: Pair ~ Center[]; theta0: REAL ~ RealFns.ArcTanDeg[p0.y-center.y, p0.x-center.x]; theta1: REAL _ RealFns.ArcTanDeg[p1.y-center.y, p1.x-center.x]; theta2: REAL _ RealFns.ArcTanDeg[p2.y-center.y, p2.x-center.x]; WHILE theta1 - theta0 < 0 DO theta1 _ theta1 + 360 ENDLOOP; WHILE theta2 - theta0 < 0 DO theta2 _ theta2 + 360 ENDLOOP; IF theta2 < theta1 THEN {theta1 _ theta1 - 360; theta2 _ theta2 - 360}; IF (theta1 - theta0)*(theta2 - theta0) < 0 THEN ERROR ELSE { radius: REAL _ Real.SqRt[DistSqr[center, p0]]; PointOnCircle: PROC [theta: REAL] RETURNS [Pair] ~ { RETURN [[center.x+radius*RealFns.CosDeg[theta], center.y+radius*RealFns.SinDeg[theta]]] }; pMid: Pair _ PointOnCircle[(theta0+theta2)*0.5]; FromArc[p0, PointOnCircle[0.75*theta0+0.25*theta2], pMid, conic]; FromArc[pMid, PointOnCircle[0.25*theta0+0.75*theta2], p2, conic]; }; }; }; }; }; END.