GGCirclesImpl.mesa
Author: Eric Bier on June 4, 1985 4:58:33 pm PDT
Last edited by Bier on January 5, 1987 8:24:32 pm PST
Contents: Routines for finding the intersections of various types of circles, lines, and points.
Pier, May 30, 1986 5:24:52 pm PDT
DIRECTORY
GGAngle, GGBasicTypes, GGCircles, GGLines, GGVector, RealFns;
GGCirclesImpl:
CEDAR
PROGRAM
IMPORTS GGAngle, GGLines, GGVector, 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 = GGBasicTypes.Line;
Vector:
TYPE = GGBasicTypes.Vector;
Circles
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 GGVector;
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]] = {
Line: cy -sx - d = 0. (f(x) = 0)
Circle: x2 + y2 = r2. (g(x) = 0)
d, magD, b: REAL;
p, h: Point;
d ← GGLines.SignedLineDistance[circle.origin, line];
magD ← ABS[d];
SELECT magD FROM
> circle.radius => hitCount ← 0;
= circle.radius => {
-- tangent. One intersection point.
hitCount ← 1;
points[1] ← GGVector.Sub[circle.origin, GGVector.Scale[[-line.s, line.c], d]];
};
< circle.radius => {
-- two intersections
hitCount ← 2;
p ← GGVector.Sub[circle.origin, GGVector.Scale[[-line.s, line.c], d]];
b ← RealFns.SqRt[circle.radius*circle.radius - d*d];
h ← GGVector.Scale[[line.c, line.s], b];
points[1] ← GGVector.Sub[p, h];
points[2] ← GGVector.Add[p, h];
};
ENDCASE => ERROR;
};
CircleMeetsEdge:
PUBLIC
PROC [circle: Circle, edge: Edge]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2]] = {
line: Line ← edge.line;
testPoints: ARRAY [1..2] OF Point;
testCount: [0..2];
[testPoints, testCount] ← CircleMeetsLine[circle, line];
hitCount ← 0;
FOR i:
NAT
IN [1..testCount]
DO
IF GGLines.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]] = {
o1ToO2, o1ToO2Hat: Vector;
magO1ToO2: REAL;
IF RatherClose[circle1.origin, circle2.origin]
THEN { -- concentric circles
hitCount ← 0;
RETURN;
};
o1ToO2 ← GGVector.Sub[circle2.origin, circle1.origin];
magO1ToO2 ← GGVector.Magnitude[o1ToO2];
[Artwork node; type 'ArtworkInterpress on' to command tool]
SELECT magO1ToO2 FROM
> circle1.radius + circle2.radius => hitCount ← 0; -- circles far apart
= circle1.radius + circle2.radius => { -- circles just touch as shown
hitCount ← 1;
o1ToO2Hat ← GGVector.Scale[o1ToO2, 1.0/magO1ToO2];
points[1] ← GGVector.Add[circle1.origin, GGVector.Scale[o1ToO2Hat, circle1.radius]];
};
IN (
ABS[circle1.radius - circle2.radius] .. circle1.radius + circle2.radius) => {
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 ← GGVector.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 ← GGVector.Add[circle1.origin, GGVector.Scale[o1ToO2Hat, m]];
}
ELSE {
m ← RealFns.SqRt[a*a - h*h];
p ← GGVector.Sub[circle2.origin, GGVector.Scale[o1ToO2Hat, m]];
};
normal ← [-h*o1ToO2Hat.y, h*o1ToO2Hat.x];
points[1] ← GGVector.Sub[p, normal];
points[2] ← GGVector.Add[p, normal];
};
=
ABS[circle1.radius - circle2.radius] => {
hitCount ← 1;
IF circle1.radius > circle2.radius
THEN {
o1ToO2Hat ← GGVector.Scale[o1ToO2, 1.0/magO1ToO2];
points[1] ← GGVector.Add[circle1.origin, GGVector.Scale[o1ToO2Hat, circle1.radius]];
}
ELSE {
o1ToO2Hat ← GGVector.Scale[o1ToO2, -1.0/magO1ToO2];
points[1] ← GGVector.Add[circle2.origin, GGVector.Scale[o1ToO2Hat, circle2.radius]];
};
};
<
ABS[circle1.radius - circle2.radius] => {
hitCount ← 0;
};
ENDCASE => ERROR;
};
CircleMeetsArc:
PUBLIC
PROC [circle: Circle, arc: Arc]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2]] = {
IF arc.edge # NIL THEN [points, hitCount] ← CircleMeetsEdge[circle, arc.edge]
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount] ← 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;
};
};
ArcMeetsLine:
PUBLIC
PROC [arc: Arc, line: Line]
RETURNS [points:
ARRAY [1..2]
OF Point, hitCount: [0..2]] = {
IF arc.edge #
NIL
THEN {
noHit: BOOL;
[points[1], noHit] ← GGLines.LineMeetsEdge[line, arc.edge];
hitCount ← IF noHit THEN 0 ELSE 1;
}
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount] ← 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]] = {
IF arc.edge #
NIL
THEN {
noHit: BOOL;
[points[1], noHit] ← GGLines.EdgeMeetsEdge[edge, arc.edge];
hitCount ← IF noHit THEN 0 ELSE 1;
}
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount] ← 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]] = {
IF arc1.edge # NIL THEN [points, hitCount] ← ArcMeetsEdge[arc2, arc1.edge]
ELSE {
pts: ARRAY [1..2] OF Point;
hCount: [0..2];
[pts, hCount] ← 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;
};
};
SignedCircleDistance:
PUBLIC
PROC [pt: Point, circle: Circle]
RETURNS [d:
REAL] = {
originToPoint: REAL;
originToPoint ← GGVector.Magnitude[GGVector.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 ← GGVector.Sub[pt, circle.origin];
magOriginToPoint ← GGVector.Magnitude[originToPoint];
IF magOriginToPoint < epsilon
THEN {
projectedPt ← [circle.origin.x + circle.radius, circle.origin.y];
RETURN;
};
projectedPt ← GGVector.Add[circle.origin, GGVector.Scale[originToPoint, circle.radius/magOriginToPoint]];
};
Arcs
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]];
};
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: GGLines.CreateEdge[v0, v2]];
}
ELSE {
vector0 ← GGVector.Sub[v0, arc.circle.origin];
theta0 ← GGVector.AngleFromVector[vector0];
vector1 ← GGVector.Sub[v1, arc.circle.origin];
theta1 ← GGVector.AngleFromVector[vector1];
vector2 ← GGVector.Sub[v2, arc.circle.origin];
theta2 ← GGVector.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 GGAngle.IsInCCWInterval[theta1, theta0, theta2]
THEN {
-- CCW arc
deltaTheta ← GGAngle.CounterClockwiseAngle[theta0, theta2];
arc^ ← [circle: arc.circle, ccw: TRUE, p0: v0, p2: v2, theta0: theta0, deltaTheta: deltaTheta, edge: NIL];
}
ELSE {
deltaTheta ← GGAngle.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 ← GGLines.CreateEmptyEdge[];
GGLines.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 ← GGVector.Sub[pt, arc.circle.origin];
angle: REAL ← GGVector.AngleFromVector[direction];
IF arc.edge # NIL THEN RETURN[GGLines.LinePointOnEdge[pt, arc.edge]];
RETURN[GGAngle.IsInCCWInterval2[angle, arc.theta0, arc.deltaTheta]];
};
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[GGLines.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 GGVector.DistanceSquared[pt, arc.p0] < GGVector.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 ← GGVector.DistanceSquared[pt, arc.p0];
distance2ToPHi ← GGVector.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[GGLines.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[GGLines.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 ← GGVector.Distance[pt, nearEndpoint];
};
};
DistanceSquaredPointToArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [distanceSquared:
REAL] = {
projectedPt: Point;
IF arc.edge # NIL THEN RETURN[GGLines.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];
};
AlmostZero:
PROC [r:
REAL]
RETURNS [
BOOL] = {
epsilon: REAL = 1.0e-5;
RETURN[ABS[r] < epsilon];
};
OnArc:
PUBLIC
PROC [pt: Point, arc: Arc]
RETURNS [
BOOL] = {
d2: REAL;
d2 ← DistanceSquaredPointToArc[pt, arc];
RETURN[AlmostZero[d2]];
};
END.