Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!elroy.jpl.nasa.gov!ncar!gatech!udel!rochester!pt.cs.cmu.edu!OREO.GRAPHICS.CS.CMU.EDU!gleicher
From: gleicher@CS.CMU.EDU (Michael Gleicher)
Newsgroups: comp.graphics
Subject: Re: HELP - intersection of 2 circles
Message-ID: <GLEICHER.90Nov17122218@OREO.GRAPHICS.CS.CMU.EDU>
Date: 17 Nov 90 17:22:18 GMT
References: <GLEICHER.90Nov16004901@OREO.GRAPHICS.CS.CMU.EDU>
Sender: gleicher@OREO.GRAPHICS.CS.CMU.EDU
Distribution: comp
Organization: School of Computer Science, Carnegie Mellon University
Lines: 146
In-reply-to: gleicher@CS.CMU.EDU's message of 16 Nov 90 05:49:01 GMT


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

------------------