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 { 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. NImagerConicImpl.mesa Michael Plass, February 16, 1984 10:04:11 am PST Edited by Doug Wyatt, November 22, 1983 11:32 am center of the circle vector from center to p0 v rotated by +90 degrees conic ratio for a 90 degree arc positive if arc is counterclockwise p1 coincides with p0 or p2 - shame, shame cos of 1/2 the angle subtended by the arc midpoint of segment from p0 to p1 direction vector of perpendicular bisector a1*x + b1*y = c1 is the equation of the perpendicular bisector *** Possible side effect on p1 in degenerate case The determinant is zero, or the calculation has lost most of its significant digits; in this case the points are practically collinear, but p1 is not between the other two. Pick as the center a point way out in the boonies, but in the right direction. Not so big that later calls on FromArc will overflow. Ê’˜J™J™0J™0J˜šÏk ˜ Jšœ œ˜Jšœ œ˜Jšœœ˜Jšœœ˜*—J˜Jšœ ˜Jšœ˜Jšœ ˜Jšœ˜J˜šœœ˜J˜—šÏnœœœœœœœ˜[šžœœ˜"Jšœœ˜ Jšœ˜—Jšœ˜Jšœ˜J˜—šžœœœœœœœ˜Zšžœœœ˜'Jšœœ˜Jšœ˜—Jšžœ˜Jšœ˜J˜—Jšœ œ˜Jšœ œ ˜š žœœœœ œ˜UJš œœœ œœ˜$šœœœ˜Jšœœ˜J˜5J˜5Jšœ˜Jšœ˜—šœ˜J˜-J˜/J˜5J˜5Jšœœœ*˜:Jšœ"˜"Jšœ"˜"Jšœ˜—Jšœ˜J˜—Jšžœœœ œ˜NJšžœœœ œ˜NJš žœœœœœœ˜6Jš žœœœœœ˜YJš žœœœœœ˜KJš žœœœœœ˜Mš žœœœœœ˜Pšœ œ˜˜-J™—˜J™—˜J™—šœœ˜J™—Jšœ˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜—šœ˜š œœœ$œœ˜EJ™#—Jšœ œ,˜:šœœ˜Jšœ)™)Jšœ˜Jšœ˜—šœ˜Jšœœ˜(šœœ˜J™)—šœ œ˜Jšœœ%˜.J˜-Jšœ3˜3Jšœœ˜Jšœ˜Jšœ˜—šœ˜˜.J™!—˜"J™*—Jšœœ˜Jšœœ ˜šœœ˜J™>—J˜.J˜"Jšœœ˜Jšœœ ˜Jšœœ˜Jšœœ˜šžœœœ œ˜&J™1š œ œœœœ˜3J™üJ˜Jšœœ˜Jšœ%œ˜<šœ(˜(Jšœ5™5—Jšœ˜Jšœ˜Jšœ˜—Jšœœ*˜5Jšœ˜—Jšœ˜Jšœœ3˜?Jšœœ3˜?Jšœœ3˜?Jšœœœ˜;Jšœœœ˜;Jšœœ0˜GJšœ)œ˜5šœ˜Jšœœ"˜.šž œœ œœ ˜4JšœQ˜WJšœ˜—Jšœ0˜0JšœA˜AJšœA˜AJšœ˜—Jšœ˜—Jšœ˜—Jšœ˜—Jšœ˜J˜—Jšœ˜—…—"