gleicher@CS.CMU.EDU (Michael Gleicher) (11/16/90)
Help! I can't seem to find a good way to do this. I have 2 circles in the plane (x1,y1,r1) and (x2,y2,r2), I know they aren't the same. I would like to know how many intersection points they have (0,1 or 2) and where these intersection points are. I'd like an analytic solution (closed form). I believe that one should exist, but i can't seem to work out the geometry. Any help would be GREATLY appreciated (a reference, a sketch of a proof, an equation, some C code, ...). Thanks, Mike
stam@dgp.toronto.edu (Jos Stam) (11/17/90)
In article <GLEICHER.90Nov16004901@OREO.GRAPHICS.CS.CMU.EDU> gleicher@CS.CMU.EDU (Michael Gleicher) writes: > >Help! I can't seem to find a good way to do this. > >I have 2 circles in the plane (x1,y1,r1) and (x2,y2,r2), I know they aren't >the same. >I would like to know how many intersection points they have (0,1 or 2) and >where these intersection points are. I'd like an analytic solution (closed >form). I believe that one should exist, but i can't seem to work out the >geometry. Any help would be GREATLY appreciated (a reference, a sketch of a >proof, an equation, some C code, ...). > >Thanks, > Mike Your problem can be stated more generally as: find the intersections (common roots) of two algebraic curves: F(x,y) = 0 G(x,y) = 0 If F is of degree n and G is of degree m, then there are nm solutions by Bezout's theorem (counting complex roots and roots "at infinity"). Of course assuming F and G have no common components... There are several methods to solve the above system, one way is to calculate the resultant, which can be given by the determinant of the Sylvester matrix of F and G. Kajiya has applied this stuff to ray-trace parametric patches: J.T. Kajiya, "Ray Tracing Parametric Patches", SIGGRAPH 82 Proceedings, 245-259. In the case of your circles, the resultant will be a polynomial of degree 4 in y (or x). Solving the polynomial will give you y_1 and y_2 (you are only interested in real solutions I assume...). Plugging these values back in the circle equations will give you the corresponding x_1 and x_2. As I'm lazy I will not calculate the resultant for you. good luck. Jos (there is probably a more clever way to do things for your special case...)
ph@miro.Berkeley.EDU (Paul Heckbert) (11/17/90)
In article gleicher@CS.CMU.EDU (Michael Gleicher) writes: >I have 2 circles in the plane (x1,y1,r1) and (x2,y2,r2). >I would like to [find their intersection points]. Hi Mike. A simple approach is to first rotate and translate to a new coordinate system (u,v) where circle 1 is centered at (0,0) and circle 2 is centered at (u2,0), where u2 is the distance between circle centers: u2^2=(x2-x1)^2+(y2-y1)^2. Now we solve u^2+v^2=r1^2, (u-u2)^2+v^2=r2^2. Subtracting these two equations, the u^2 and v^2 terms drop and we're left with a linear equation for u: u = u2/2 + (r1^2-r2^2)/(2*u2). v can then be found as +-sqrt(r1^2-u^2) and you'll get 0, 1, or 2 intersection points when r1^2-u^2 is <0, =0, or >0, respectively. It should be possible to do the transformations without trig. A second, nearly equivalent way to set up the problem is to note that you know the three side lengths of the triangle defined by the two centers and one of the intersection points, so you could use the law of cosines to find the angles and do a little trig to compute the answer. A third approach is to crank through the algebra explicitly from the original equations: (x-x1)^2+(y-y1)^2=r1^2, (x-x2)^2+(y-y2)^2=r2^2. Solve for y in one equation, substitute into the other, expand the squares, move the remaining sqrt term to one side of the equation, yielding sqrt(quadratic in x)=(linear in x), so if you square both sides you'll get a quadratic equation which you can solve for x. Paul Heckbert, Computer Science Dept. 570 Evans Hall, UC Berkeley INTERNET: ph@miro.berkeley.edu Berkeley, CA 94720 UUCP: ucbvax!miro.berkeley.edu!ph ps: we don't need no stinkin' resultants.
gleicher@CS.CMU.EDU (Michael Gleicher) (11/18/90)
Thanks to everyone who responded.
I sorta feel silly though. Immediately after posting I found an extremely
simple solution. Thats what I get for posting in frustration. My office was
covered with picture of intersecting circles, I just needed to look at the
diagrams a little differently. I'm sorry if I wasted people's time.
The following morning I recieved a few responses, one of which was the very
same solution, some others were simple variants.
In case you're curious Here's my solution. It's nice since it uses no trig,
and didn't require me to solve any non-linear equations.
Given and x1,y1,r1,x2,y2,r2
Determine whether or not the circles intersect by comparing the distance
between the centers and the sum and differences of the radii. The distance
must fall between these two quantities.
Assuming they do intersect:
let c and d be the centers * mark the intersections, and p be the point where
the segment connecting the intersection points meet.
*
r1/|\r2 let V be the vector c->d
/ | \ let d be |V|
c--p--d let a be the distance from 1 to p,
\ | / the distance 2 to p is d-p
\|/ let b be the distance from an intersection to p
*
so p^2 = r1^2 - a^2 = r2^2 - (d-a)^2,
a little algrebra leaves:
a = (r1^2 - r2^2 + d^2) / (2*d)
then b can be computed:
b = sqrt(r1^2-a^2)
So,
p = c + a/d V,
and the intersections are:
i1 = p + b/d N
and i2 = p - b/d N (where N is the normal vector, perpedicular to V)
or, in code (with a test driver for an SGI Iris):
(BTW: this is C++, though there's nothing magic being used)
#include "stdio.h"
#include "math.h"
typedef float Real;
inline Real ABS(const Real x) { return( ((x) >= 0) ? (x) : -(x)); }
const Real CI_TOLL = .03; /* don't put intersections closer than this */
int circInters(Real x1, Real y1, Real r1,
Real x2, Real y2, Real r2,
Real& p1x, Real& p1y,
Real& p2x, Real& p2y)
{
/* the vector between the centers */
Real vx = x2-x1;
Real vy = y2-y1;
/* its length */
Real d = hypot(vx,vy);
if ((ABS(r1-r2) > d) ||
( (r1+r2) < d) ) return 0;
/* a unit vector along v */
Real uvx = vx / d;
Real uvy = vy / d;
/* compute the distance along v where the line between the points is */
Real a = (r1*r1 - r2*r2 + d*d) / (2 * d);
if (a==0) return 0; /* infinite # of intersections */
/* compute the point A */
Real ax = x1 + uvx*a;
Real ay = y1 + uvy*a;
/* compute the distance from A to the intersections */
Real ad = sqrt(r1*r1 - a*a);
if (ad < CI_TOLL) {
p1x = p2x = ax;
p1y = p2y = ay;
return 1;
}
/* the unit normal vector */
Real nx = uvy;
Real ny = -uvx;
/* the points of intersection */
p1x = ax + ad*nx;
p1y = ay + ad*ny;
p2x = ax - ad*nx;
p2y = ay - ad*ny;
return 2;
}
#ifdef TESTDRIVE
/* a test driver (for an iris) */
#include "gl.h"
main()
{
Real x1,y1,r1,x2,y2,r2,p1x,p1y,p2x,p2y;
int n;
keepaspect(1,1);
foreground();
winopen("circles");
gconfig();
ortho2(-5,5,-5,5);
color(WHITE);
clear();
while(x1 != -999) {
printf("x1 y1 r1 x2 y2 z2 ?\n");
scanf("%f %f %f %f %f %f",&x1,&y1,&r1,&x2,&y2,&r2);
n = circInters(x1,y1,r1,x2,y2,r2,p1x,p1y,p2x,p2y);
printf("%d intersections (%5.2f %5.2f) (%5.2f %5.2f)\n",
n,p1x,p1y,p2x,p2y);
color(WHITE);
clear();
color(BLUE);
linewidth(3);
circ(x1,y1,r1);
circ(x2,y2,r2);
color(RED);
if (n>0) sboxf(p1x-.1,p1y-.1,p1x+.1,p1y+.1);
if (n>1) sboxf(p2x-.1,p2y-.1,p2x+.1,p2y+.1);
}
}
#endif
------------------nujy@vax5.cit.cornell.edu (11/18/90)
The number of intersections is easy. Given the circles (x1,y1) radius r1
and (x2,y2) radius r2 (where r2 > r1), then the distance between
the centers is d=sqrt((x2-x1)^2+(y2-y1)^2) and:
abs(d-r2) > r1 --> 0 intersections
abs(d-r2) = r1 --> 1 intersection (2 degenerate intersections)
abs(d-r2) < r1 --> 2 intersections
abs(d-r2) is the distance from the center of the smaller circle to the
nearest point on the larger circle.
Finding the actual intersections is a bit harder. The points must lie
on a line (Lp) perpendicular to the line connecting the centers of the
circles (Lc).
/ <-- C2
___ | + : center of a circle
C1 --> _-- -*_ * : intersection
/ |! ! : perpendicular line (Lp)
| +-|!--|-------+
_ |!_/ ^--- Lc That's not too confusing, is it?
--___-*
| Fig. 1
Furthermore these points are the same distance from Lc, call it h. This
proof is left as an exercise for the reader (LEFR). Now all we need to
find is the location of the intersection of Lp and Lc (x',y'), and h.
Figure 2 shows the geometry of one case. We need to solve for p1.
*_
r1 /! --__ r2
/ !h --__
+-------------+ Fig. 2
|p1|<---p2--->|
|<-----d----->|
r1^2 - p1^2 = h^2 (1) Pythagorean Thm.
r2^2 - p2^2 = h^2 (2) Ditto
p1^2 - p2^2 = r1^2 - r2^2 (3) from 1 & 2
p1 - p2 = (r1^2 - r2^2) / d (4) from 3 and p1+p2=d
p1 = (d^2 + r1^2 - r2^2) / 2d (5) add d to 4, isolate x1
Using the parametric equation of a line, and letting t=p1/d:
x' = x1 + (x2-x1)*p1 / d,
y' = y1 + (y2-y1)*p1 / d.
The intersections are at:
(x'+x",y'+y"), (x'-x",y'-y"),
where
x" = h*(y2-y1) / d,
y" = h*(x1-x2) / d,
is LEFR.
The other possible configuration is:
*_
!--__ r2
h! r1--__
---+-------+
|p |<--d-->|
However, the above solution also applies here. Isn't that nice.
-Chris
Chris Schoeneman | "I was neat, clean, shaved, and sober,
nujy@vax5.ccs.cornell.edu | and I didn't care who knew it."
| - Raymond Chandler