GGCirclesImpl.mesa
Author: Eric Bier on June 4, 1985 4:58:33 pm PDT
Last edited by Bier on July 16, 1987 1:58:44 pm PDT
Contents: Routines for finding the intersections of various types of circles, lines, and points.
Pier, September 18, 1990 11:24 am PDT
Bier, May 8, 1989 2:14:49 pm PDT
DIRECTORY
Angles2d, GGBasicTypes, GGCircles, GGCoreTypes, GGUtility, Lines2d, RealFns, Vectors2d;
GGCirclesImpl:
CEDAR
PROGRAM
IMPORTS Angles2d, Lines2d, Vectors2d, RealFns
EXPORTS GGCircles = BEGIN
Arc: TYPE = REF ArcObj;
ArcObj: TYPE = GGBasicTypes.ArcObj;
Circle: TYPE = REF CircleObj;
CircleObj: TYPE = GGBasicTypes.CircleObj;
Point: TYPE = GGBasicTypes.Point;
Edge: TYPE = GGBasicTypes.Edge;
Line: TYPE = GGCoreTypes.Line;
Vector:
TYPE = GGBasicTypes.Vector;
CreateEmptyCircle:
PUBLIC
PROC []
RETURNS [circle: Circle] = {
circle ←
NEW[CircleObj ← [
origin: [0.0, 0.0],
radius: 0.0
]];
};
CopyCircle:
PUBLIC
PROC [from: Circle, to: Circle] = {
to.origin ← from.origin;
to.radius ← from.radius;
};
FillCircleFromPointAndRadius:
PUBLIC
PROC [pt: Point, radius:
REAL, circle: Circle] = {
circle.origin ← pt;
circle.radius ← radius;
};
CircleFromPointAndRadius:
PUBLIC
PROC [pt: Point, radius:
REAL]
RETURNS [circle: Circle] = {
circle ← CreateEmptyCircle[];
FillCircleFromPointAndRadius[pt, radius, circle];
};
CircleFrom3Points:
PUBLIC
PROC [p0, p1, p2: Point]
RETURNS [circle: Circle, linear:
BOOL] = {
circle ← CreateEmptyCircle[];
linear ← FillCircleFrom3Points[p0, p1, p2, circle];
};
FillCircleFrom3Points:
PUBLIC
PROC [p0, p1, p2: Point, circle: Circle]
RETURNS [linear:
BOOL] = {
Find the perpendicular bisectors of [p0, p1] and [p1, p2]. Their intersection is the center of circle. The distance of p0 or p1 or p2 from the center is the radius.
OPEN Vectors2d;
IF p0 = p2
THEN {
circle.origin ← [(p0.x+p1.x)*0.5, (p0.y+p1.y)*0.5]; -- center of the circle
circle.radius ← Distance[circle.origin, p0]; -- vector from center to p0
linear ← FALSE;
}
ELSE {
m1: Point ← [(p0.x+p1.x)*0.5, (p0.y+p1.y)*0.5]; -- midpoint of segment from p0 to p1
v1: Point ← [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 first perpendicular bisector
m2: Point ~ [(p0.x+p2.x)*0.5, (p0.y+p2.y)*0.5]; -- midpoint of segment from p1 to p2
v2: Point ~ [m2.y-p0.y, p0.x-m2.x];
a2: REAL ~ v2.y;
b2: REAL ~ -v2.x;
c2:
REAL ~ a2*m2.x+b2*m2.y;
a2*x + b2*y = c2 is the equation of the second perpendicular bisector
These lines are now two simultaneous linear equations. The determinant is:
det: REAL ~ a1*b2-a2*b1;
IF det = 0
OR
ABS[det] <
ABS[a1*b2]*0.000005
THEN {
circle.origin ← [0.0, 0.0];
circle.radius ← 12.0; -- magic number for debugging
linear ← TRUE;
}
ELSE {
circle.origin ← [(c1*b2-c2*b1)/det, (a1*c2-a2*c1)/det]; -- solve the linear system
circle.radius ← Distance[circle.origin, p0];
linear ← FALSE;
};
};
};
RatherClose:
PROC [p1, p2: Point]
RETURNS [
BOOL] = {
epsilon: REAL = 1.0e-5;
RETURN[ABS[p1.x - p2.x] < epsilon AND ABS[p1.y - p2.y] < epsilon];
};
CircleMeetsLine:
PUBLIC
PROC [circle: Circle, line: Line]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL ←
FALSE] = {
Line: cy -sx - d = 0. (f(x) = 0)
Circle: x2 + y2 = r2. (g(x) = 0)
d, magD, b: REAL;
p, h: Point;
epsilon: REAL = GGUtility.epsilonInPoints;
d ← Lines2d.SignedLineDistance[circle.origin, line];
magD ← ABS[d];
SELECT magD - circle.radius FROM
> epsilon => hitCount ← 0;
IN [-epsilon..epsilon] => {
-- tangent. One intersection point.
hitCount ← 1;
points[1] ← Vectors2d.Sub[circle.origin, Vectors2d.Scale[[-line.s, line.c], d]];
tangent ← TRUE;
};
< -epsilon => {
-- two intersections
hitCount ← 2;
p ← Vectors2d.Sub[circle.origin, Vectors2d.Scale[[-line.s, line.c], d]];
b ← RealFns.SqRt[circle.radius*circle.radius - d*d];
h ← Vectors2d.Scale[[line.c, line.s], b];
points[1] ← Vectors2d.Sub[p, h];
points[2] ← Vectors2d.Add[p, h];
};
ENDCASE => ERROR;
};
CircleMeetsEdge:
PUBLIC
PROC [circle: Circle, edge: Edge]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL] = {
line: Line ← edge.line;
testPoints: ARRAY [1..2] OF Point;
testCount: [0..2];
[testPoints, testCount, tangent] ← CircleMeetsLine[circle, line];
hitCount ← 0;
FOR i:
NAT
IN [1..testCount]
DO
IF Lines2d.OnEdge[testPoints[i], edge]
THEN {
hitCount ← hitCount + 1;
points[hitCount] ← testPoints[i];
};
ENDLOOP;
};
CircleMeetsCircle:
PUBLIC
PROC [circle1, circle2: Circle]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL ←
FALSE] = {
o1ToO2, o1ToO2Hat: Vector;
epsilon: REAL = GGUtility.epsilonInPoints;
magO1ToO2, outerTangent, innerTangent: REAL;
IF RatherClose[circle1.origin, circle2.origin]
THEN {
-- concentric circles
hitCount ← 0;
points[1] ← points[2] ← [0.0, 0.0]; -- avoids compiler warnings
RETURN;
};
o1ToO2 ← Vectors2d.Sub[circle2.origin, circle1.origin];
magO1ToO2 ← Vectors2d.Magnitude[o1ToO2];
outerTangent ← circle1.radius + circle2.radius;
innerTangent ← ABS[circle1.radius - circle2.radius];
[Artwork node; type 'ArtworkInterpress on' to command tool]
SELECT magO1ToO2 FROM
> outerTangent+epsilon => hitCount ← 0; -- circles far apart
IN [outerTangent-epsilon .. outerTangent+epsilon] => {
-- circles just touch as shown
hitCount ← 1;
tangent ← TRUE;
o1ToO2Hat ← Vectors2d.Scale[o1ToO2, 1.0/magO1ToO2];
points[1] ← Vectors2d.Add[circle1.origin, Vectors2d.Scale[o1ToO2Hat, circle1.radius]];
};
IN (innerTangent+epsilon .. outerTangent-epsilon) => {
The two circles overlap. We expect two roots. In the picture below, point C is one of the circle intersection points. Segment AB is the segment O1O2 in the picture above. We wish to find the altitude h of the triangle. b = r1 and a = r2. From Heron's formula for the area of a triangle, area K = sqrt(s(s-a)(s-b)(s-c)), where s = (a+b+ c)/2. We also know that the area K = 0.5ch. So h = 2K/c. This gives us the y coordinate of the intersection points. The x coordinate is the point where the altitude hits the base, which we get from the Pythagorean Theorem.
[Artwork node; type 'ArtworkInterpress on' to command tool]
p: Point;
normal: Vector;
s, h, m: REAL;
b: REAL ← circle1.radius;
a: REAL ← circle2.radius;
hitCount ← 2;
o1ToO2Hat ← Vectors2d.Scale[o1ToO2, 1.0/magO1ToO2];
s ← 0.5*(b + a + magO1ToO2);
h ← 2.0*RealFns.SqRt[s*(s-magO1ToO2)*(s-b)*(s-a)]/magO1ToO2;
IF b > a
THEN {
m ← RealFns.SqRt[b*b - h*h];
p ← Vectors2d.Add[circle1.origin, Vectors2d.Scale[o1ToO2Hat, m]];
}
ELSE {
m ← RealFns.SqRt[a*a - h*h];
p ← Vectors2d.Sub[circle2.origin, Vectors2d.Scale[o1ToO2Hat, m]];
};
normal ← [-h*o1ToO2Hat.y, h*o1ToO2Hat.x];
points[1] ← Vectors2d.Sub[p, normal];
points[2] ← Vectors2d.Add[p, normal];
};
IN [innerTangent-epsilon .. innerTangent+epsilon] => {
hitCount ← 1;
tangent ← TRUE;
IF circle1.radius > circle2.radius
THEN {
o1ToO2Hat ← Vectors2d.Scale[o1ToO2, 1.0/magO1ToO2];
points[1] ← Vectors2d.Add[circle1.origin, Vectors2d.Scale[o1ToO2Hat, circle1.radius]];
}
ELSE {
o1ToO2Hat ← Vectors2d.Scale[o1ToO2, -1.0/magO1ToO2];
points[1] ← Vectors2d.Add[circle2.origin, Vectors2d.Scale[o1ToO2Hat, circle2.radius]];
};
};
< innerTangent-epsilon => {
hitCount ← 0;
};
ENDCASE => ERROR;
};
CircleMeetsArc:
PUBLIC
PROC [circle: Circle, arc: Arc]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL] = {
IF arc.edge # NIL THEN [points, hitCount, tangent] ← CircleMeetsEdge[circle, arc.edge]
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount, tangent] ← CircleMeetsCircle[circle, arc.circle];
hitCount ← 0;
FOR i:
NAT
IN [1..hCount]
DO
IF OnArc[pts[i], arc]
THEN {
hitCount ← hitCount + 1;
points[hitCount] ← pts[i];
};
ENDLOOP;
};
};
SignedCircleDistance:
PUBLIC
PROC [pt: Point, circle: Circle]
RETURNS [d:
REAL] = {
originToPoint: REAL;
originToPoint ← Vectors2d.Magnitude[Vectors2d.Sub[pt, circle.origin]];
d ← originToPoint - circle.radius;
};
CircleDistance:
PUBLIC
PROC [pt: Point, circle: Circle]
RETURNS [d:
REAL] = {
d ← ABS[SignedCircleDistance[pt, circle]];
};
PointProjectedOntoCircle:
PUBLIC
PROC [pt: Point, circle: Circle]
RETURNS [projectedPt: Point] = {
We drop a normal from the point onto the circle and find where it hits.
epsilon: REAL = 1.0e-5;
originToPoint: Vector;
magOriginToPoint: REAL;
originToPoint ← Vectors2d.Sub[pt, circle.origin];
magOriginToPoint ← Vectors2d.Magnitude[originToPoint];
IF magOriginToPoint < epsilon
THEN {
projectedPt ← [circle.origin.x + circle.radius, circle.origin.y];
RETURN;
};
projectedPt ← Vectors2d.Add[circle.origin, Vectors2d.Scale[originToPoint, circle.radius/magOriginToPoint]];
};
PointIsInCircle:
PUBLIC
PROC [pt: Point, circle: Circle]
RETURNS [
BOOL] = {
deltaX, deltaY: REAL;
deltaX ← pt.x - circle.origin.x;
deltaY ← pt.y - circle.origin.y;
RETURN[deltaX*deltaX+deltaY*deltaY < circle.radius*circle.radius];
};
CreateArc:
PUBLIC
PROC [v0, v1, v2: Point]
RETURNS [arc: Arc] = {
arc ← CreateEmptyArc[];
FillArc[v0, v1, v2, arc];
};
CreateEmptyArc:
PUBLIC
PROC []
RETURNS [arc: Arc] = {
arc ← NEW[ArcObj ← [circle: CreateEmptyCircle[], ccw: FALSE, p0: [0.0,0.0], p2: [0.0,0.0], theta0: 0.0, deltaTheta: 0.0]];
};
ReverseArc:
PUBLIC
PROC [arc: Arc] = {
Invariants: arc.deltaTheta is always a positive, counterclockwise angle. All arcs can be thought of as counterclockwise if you start at the proper end. arc.theta0 always correpondes to arc.p0. arc.ccw says whether the arc is really counterclockwise or not.
arc.ccw ← NOT arc.ccw;
IF arc.edge # NIL THEN arc.edge.startIsFirst ← NOT arc.edge.startIsFirst; -- This should be a call to Lines2d.ReverseEdge, but there is no Lines2d.ReverseEdge yet.
};
FillArc:
PUBLIC
PROC [v0, v1, v2: Point, arc: Arc] = {
theta0, theta1, theta2, deltaTheta: REAL;
vector0, vector1, vector2: Vector;
linear: BOOL;
linear ← FillCircleFrom3Points[v0, v1, v2, arc.circle];
IF linear
THEN {
arc^ ← [circle: arc.circle, ccw: TRUE, p0: v0, p2: v2, theta0: 0.0, deltaTheta: 0.0, edge: Lines2d.CreateEdge[v0, v2]];
}
ELSE {
vector0 ← Vectors2d.Sub[v0, arc.circle.origin];
theta0 ← Vectors2d.AngleFromVector[vector0];
vector1 ← Vectors2d.Sub[v1, arc.circle.origin];
theta1 ← Vectors2d.AngleFromVector[vector1];
vector2 ← Vectors2d.Sub[v2, arc.circle.origin];
theta2 ← Vectors2d.AngleFromVector[vector2];
IF theta0 = theta2
THEN {
deltaTheta ← 360.0;
arc^ ← [circle: arc.circle, ccw: TRUE, p0: v0, p2: v2, theta0: theta0, deltaTheta: deltaTheta, edge: NIL];
}
ELSE {
IF Angles2d.IsInCCWInterval[theta1, theta0, theta2]
THEN {
-- CCW arc
deltaTheta ← Angles2d.CounterClockwiseAngle[theta0, theta2];
arc^ ← [circle: arc.circle, ccw: TRUE, p0: v0, p2: v2, theta0: theta0, deltaTheta: deltaTheta, edge: NIL];
}
ELSE {
deltaTheta ← Angles2d.CounterClockwiseAngle[theta2, theta0];
arc^ ← [circle: arc.circle, ccw: FALSE, p0: v2, p2: v0, theta0: theta2, deltaTheta: deltaTheta, edge: NIL];
};
};
};
};
CopyArc:
PUBLIC
PROC [from: Arc, to: Arc] = {
CopyCircle[from.circle, to.circle];
to.ccw ← from.ccw;
to.p0 ← from.p0;
to.p2 ← from.p2;
to.theta0 ← from.theta0;
to.deltaTheta ← from.deltaTheta;
IF from.edge #
NIL
THEN {
to.edge ← Lines2d.CreateEmptyEdge[];
Lines2d.CopyEdge[from: from.edge, to: to.edge];
};
};
CirclePointOnArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [
BOOL] = {
Assumes pt is on arc.circle. Is it on arc?
direction: Vector ← Vectors2d.Sub[pt, arc.circle.origin];
angle: REAL ← Vectors2d.AngleFromVector[direction];
IF arc.edge # NIL THEN RETURN[Lines2d.LinePointOnEdge[pt, arc.edge]];
RETURN[Angles2d.IsInCCWInterval2[angle, arc.theta0, arc.deltaTheta]];
};
ArcMeetsLine:
PUBLIC
PROC [arc: Arc, line: Line]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL ←
FALSE] = {
IF arc.edge #
NIL
THEN {
noHit: BOOL;
[points[1], noHit] ← Lines2d.LineMeetsEdge[line, arc.edge];
hitCount ← IF noHit THEN 0 ELSE 1;
}
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount, tangent] ← CircleMeetsLine[arc.circle, line];
hitCount ← 0;
FOR i:
NAT
IN [1..hCount]
DO
IF OnArc[pts[i], arc]
THEN {
hitCount ← hitCount + 1;
points[hitCount] ← pts[i];
};
ENDLOOP;
};
};
ArcMeetsEdge:
PUBLIC
PROC [arc: Arc, edge: Edge]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL ←
FALSE] = {
IF arc.edge #
NIL
THEN {
noHit: BOOL;
[points[1], noHit] ← Lines2d.EdgeMeetsEdge[edge, arc.edge];
hitCount ← IF noHit THEN 0 ELSE 1;
}
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount, tangent] ← CircleMeetsEdge[arc.circle, edge];
hitCount ← 0;
FOR i:
NAT
IN [1..hCount]
DO
IF OnArc[pts[i], arc]
THEN {
hitCount ← hitCount + 1;
points[hitCount] ← pts[i];
};
ENDLOOP;
};
};
ArcMeetsArc:
PUBLIC
PROC [arc1, arc2: Arc]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2], tangent:
BOOL] = {
IF arc1.edge # NIL THEN [points, hitCount, tangent] ← ArcMeetsEdge[arc2, arc1.edge]
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount, tangent] ← CircleMeetsArc[arc1.circle, arc2];
hitCount ← 0;
FOR i:
NAT
IN [1..hCount]
DO
IF OnArc[pts[i], arc1]
THEN {
hitCount ← hitCount + 1;
points[hitCount] ← pts[i];
};
ENDLOOP;
};
};
NearestEndpoint:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [endpoint: Point] = {
Look for a winner with the infinity norm first. If that fails, use the 2 norm.
IF arc.edge # NIL THEN RETURN[Lines2d.NearestEndpoint[pt, arc.edge]];
BEGIN
IF
ABS[pt.x-arc.p0.x] <=
ABS[pt.x-arc.p2.x]
THEN {
IF ABS[pt.y-arc.p0.y] <= ABS[pt.y-arc.p2.y] THEN RETURN[arc.p0]
ELSE GOTO DoMath
}
ELSE
IF ABS[pt.y-arc.p0.y] > ABS[pt.y-arc.p2.y] THEN RETURN[arc.p2]
ELSE GOTO DoMath;
EXITS
DoMath =>
IF Vectors2d.DistanceSquared[pt, arc.p0] < Vectors2d.DistanceSquared[pt, arc.p2]
THEN endpoint ← arc.p0
ELSE endpoint ← arc.p2;
END;
};
DistanceSquaredToNearestEndpoint:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [distanceSquared:
REAL] = {
distance2ToPLo, distance2ToPHi: REAL;
distance2ToPLo ← Vectors2d.DistanceSquared[pt, arc.p0];
distance2ToPHi ← Vectors2d.DistanceSquared[pt, arc.p2];
RETURN[MIN[distance2ToPLo, distance2ToPHi]];
};
NearestPointOnArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [onArc: Point] = {
The nearest point on arc.circle to pt, if the nearest point is on arc. Otherwise, the nearest endpoint.
projectedPt: Point;
IF arc.edge # NIL THEN RETURN[Lines2d.NearestPointOnEdge[pt, arc.edge]];
projectedPt ← PointProjectedOntoCircle[pt, arc.circle];
IF CirclePointOnArc[projectedPt, arc] THEN onArc ← projectedPt
ELSE onArc ← NearestEndpoint[pt, arc];
};
DistancePointToArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [distance:
REAL] = {
Distance[pt, NearestPointOnArc[pt, arc]].
projectedPt, nearEndpoint: Point;
IF arc.edge # NIL THEN RETURN[Lines2d.DistancePointToEdge[pt, arc.edge]];
projectedPt ← PointProjectedOntoCircle[pt, arc.circle];
IF CirclePointOnArc[projectedPt, arc] THEN distance ← ABS[CircleDistance[pt, arc.circle]]
ELSE {
nearEndpoint ← NearestEndpoint[pt, arc];
distance ← Vectors2d.Distance[pt, nearEndpoint];
};
};
DistanceSquaredPointToArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [distanceSquared:
REAL] = {
projectedPt: Point;
IF arc.edge # NIL THEN RETURN[Lines2d.DistanceSquaredPointToEdge[pt, arc.edge]];
projectedPt ← PointProjectedOntoCircle[pt, arc.circle];
IF CirclePointOnArc[projectedPt, arc]
THEN {distanceSquared ← CircleDistance[pt, arc.circle];
distanceSquared ← distanceSquared*distanceSquared}
ELSE distanceSquared ← DistanceSquaredToNearestEndpoint[pt, arc];
};
OnArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [
BOOL] = {
d2: REAL;
epsilon: REAL = 1.0e-5;
d2 ← DistanceSquaredPointToArc[pt, arc];
RETURN[d2 < epsilon];
};
END.