ImagerConicImpl.mesa
Michael Plass, February 16, 1984 10:04:11 am PST
Edited by Doug Wyatt, November 22, 1983 11:32 am
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];
center of the circle
v: Pair ~ [p0.x-c.x, p0.y-c.y];
vector from center to p0
p: Pair ~ [-v.y, v.x];
v rotated by +90 degrees
r:
REAL ~ Real.SqRt[2.0]-1.0;
conic ratio for a 90 degree arc
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;
positive if arc is counterclockwise
distProd: REAL ~ Real.SqRt[DistSqr[p0,p1]*DistSqr[p2,p1]];
IF distProd = 0
THEN {
p1 coincides with p0 or p2 - shame, shame
conic[p1,p2,0.5];
}
ELSE {
dot: REAL ~ Dot[Sub[p0,p1], Sub[p2,p1]];
cos:
REAL ~ -dot/distProd;
cos of 1/2 the angle subtended by the arc
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];
midpoint of segment from p0 to p1
v1: Pair ~ [m1.y-p0.y, p0.x-m1.x];
direction vector of perpendicular bisector
a1: REAL ~ v1.y;
b1: REAL ~ -v1.x;
c1:
REAL ~ a1*m1.x+b1*m1.y;
a1*x + b1*y = c1 is the equation of the perpendicular bisector
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 {
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.
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];
Not so big that later calls on FromArc will overflow.
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.