<> <> <> <> <> DIRECTORY Angles2d, GGBasicTypes, GGCircles, Lines2d, GGUtility, Vectors2d, RealFns; 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 = GGBasicTypes.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] = { <> 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; <> 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; <> <> 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] = { <> <> 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; 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) => { <> << [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; }; }; 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; }; }; 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] = { <> 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]]; }; <<>> <> 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: 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] = { <> 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]]; }; NearestEndpoint: PUBLIC PROC [pt: Point, arc: Arc] RETURNS [endpoint: Point] = { <> 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] = { <> 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] = { <> 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.