70

How do I calculate the intersection points of two circles. I would expect there to be either two, one or no intersection points in all cases.

I have the x and y coordinates of the centre-point, and the radius for each circle.

An answer in python would be preferred, but any working algorithm would be acceptable.

1
  • 2
    Another case could result if the center and radius of the two circles are the same.
    – andand
    Commented Jul 28, 2010 at 18:59

6 Answers 6

124

Intersection of two circles

Written by Paul Bourke

The following note describes how to find the intersection point(s) between two circles on a plane, the following notation is used. The aim is to find the two points P3 = (x3, y3) if they exist.

Intersection of 2 circles

First calculate the distance d between the center of the circles. d = ||P1 - P0||.

  • If d > r0 + r1 then there are no solutions, the circles are separate.

  • If d < |r0 - r1| then there are no solutions because one circle is contained within the other.

  • If d = 0 and r0 = r1 then the circles are coincident and there are an infinite number of solutions.

Considering the two triangles P0P2P3 and P1P2P3 we can write

a2 + h2 = r02 and b2 + h2 = r12

Using d = a + b we can solve for a,

a = (r02 - r12 + d2 ) / (2 d)

It can be readily shown that this reduces to r0 when the two circles touch at one point, ie: d = r0 + r1

Solve for h by substituting a into the first equation, h2 = r02 - a2

So

P2 = P0 + a ( P1 - P0 ) / d

And finally, P3 = (x3,y3) in terms of P0 = (x0,y0), P1 = (x1,y1) and P2 = (x2,y2), is

x3 = x2 +- h ( y1 - y0 ) / d

y3 = y2 -+ h ( x1 - x0 ) / d

Source: http://paulbourke.net/geometry/circlesphere/

11
  • 1
    The link is down, can you upload a new one? Also, can anyone explain how a=(r02-r12+d2)/(2*d) was obtained from previous equations? Maybe it's pretty obvious, but I'm just not getting there. Thanks
    – David Neto
    Commented Jan 30, 2014 at 18:24
  • 1
    @DavidNeto I fixed the link. Also, $a$ is obtained by solving for $h^2$ and plugging it into the other equation $h^2 = r_0^2 - a^2$, so now you can solve for $a$ here: $b^2 + r_0^2 - a^2 = r_1^2$ (remember that $d = a + b$)
    – nachocab
    Commented Sep 7, 2014 at 21:00
  • 2
    Great this is a good explanation. Here is the algorithm in python: gist.github.com/xaedes/974535e71009fa8f090e
    – xaedes
    Commented Mar 6, 2016 at 15:17
  • 6
    Using d = a + b we can solve for a. How? This step looks like dark magic to me. Which equation was used to solve for a?
    – Makogan
    Commented Jul 12, 2018 at 0:16
  • 11
    b^2 = r1^2 - h^2 and h^2 = r0^2 - a^2 b^2 = r1^2 - r0^2 + a^2 r0^2 - r1^2 = a^2 - b^2 d^2 + r0^2 - r1^2 = a^2 - b^2 + d^2 Now, use equality (a + b)^2 = d^2 = a^2 + b^2 + 2ab d^2 + r0^2 - r1^2 = a^2 - b^2 + a^2 + b^2 + 2ab d^2 + r0^2 - r1^2 = 2 * a^2 + 2ab From a + b = d you can get that 2 * a^2 + 2ab = 2da if you multiply both sides by 2a Finally, d^2 + r0^2 - r1^2 = 2da
    – Y N
    Commented Aug 5, 2018 at 18:48
16

Why not just use 7 lines of your favorite procedural language (or programmable calculator!) as below.

Assuming you are given P0 coords (x0,y0), P1 coords (x1,y1), r0 and r1 and you want to find P3 coords (x3,y3):

d=sqr((x1-x0)^2 + (y1-y0)^2)
a=(r0^2-r1^2+d^2)/(2*d)
h=sqr(r0^2-a^2)
x2=x0+a*(x1-x0)/d   
y2=y0+a*(y1-y0)/d   
x3=x2+h*(y1-y0)/d       // also x3=x2-h*(y1-y0)/d
y3=y2-h*(x1-x0)/d       // also y3=y2+h*(x1-x0)/d
4
  • What if (r0^2-a^2) < 0?
    – ablaze
    Commented Dec 4, 2017 at 15:58
  • Ah, yes. See Paul Bourke's original answer, you first have to make sure you test to avoid d > r0 + r1 (circles too far apart), d < |r0 - r1| (one circle inside the other) and (d = 0 and r0 = r1) (circles coincident). Outside of these cases, you'll have a solution and the expression you mention will be positive. Commented Dec 5, 2017 at 17:46
  • So it is an absolute difference, I missed that. Thank you!
    – ablaze
    Commented Dec 9, 2017 at 3:32
  • @ablaze: Not really. I should have simply said that if (r0^2-a^2) <= 0 then you have no intersection points between the circles. Commented Sep 1, 2019 at 6:19
15

Here is my C++ implementation based on Paul Bourke's article. It only works if there are two intersections, otherwise it probably returns NaN NAN NAN NAN.

class Point{
    public:
        float x, y;
        Point(float px, float py) {
            x = px;
            y = py;
        }
        Point sub(Point p2) {
            return Point(x - p2.x, y - p2.y);
        }
        Point add(Point p2) {
            return Point(x + p2.x, y + p2.y);
        }
        float distance(Point p2) {
            return sqrt((x - p2.x)*(x - p2.x) + (y - p2.y)*(y - p2.y));
        }
        Point normal() {
            float length = sqrt(x*x + y*y);
            return Point(x/length, y/length);
        }
        Point scale(float s) {
            return Point(x*s, y*s);
        }
};

class Circle {
    public:
        float x, y, r, left;
        Circle(float cx, float cy, float cr) {
            x = cx;
            y = cy;
            r = cr;
            left = x - r;
        }
        pair<Point, Point> intersections(Circle c) {
            Point P0(x, y);
            Point P1(c.x, c.y);
            float d, a, h;
            d = P0.distance(P1);
            a = (r*r - c.r*c.r + d*d)/(2*d);
            h = sqrt(r*r - a*a);
            Point P2 = P1.sub(P0).scale(a/d).add(P0);
            float x3, y3, x4, y4;
            x3 = P2.x + h*(P1.y - P0.y)/d;
            y3 = P2.y - h*(P1.x - P0.x)/d;
            x4 = P2.x - h*(P1.y - P0.y)/d;
            y4 = P2.y + h*(P1.x - P0.x)/d;

            return pair<Point, Point>(Point(x3, y3), Point(x4, y4));
        }

};
3
  • 2
    nope ;) - sounds like a fun exercise to try yourself though ;p
    – Rusty Rob
    Commented May 29, 2014 at 1:54
  • May I ask what's the purpose of Point P2 = P1.sub(P0).scale(a/d).add(P0); this line?
    – padawan
    Commented May 31, 2014 at 15:08
  • it's an arrangement of P2 = P0 + a ( P1 - P0 ) / d
    – Rusty Rob
    Commented May 31, 2014 at 19:59
13

Here's an implementation in Javascript using vectors. The code is well documented, you should be able to follow it. Here's the original source

See live demo here: enter image description here

// Let EPS (epsilon) be a small value
var EPS = 0.0000001;

// Let a point be a pair: (x, y)
function Point(x, y) {
  this.x = x;
  this.y = y;
}

// Define a circle centered at (x,y) with radius r
function Circle(x,y,r) {
  this.x = x;
  this.y = y;
  this.r = r;
}

// Due to double rounding precision the value passed into the Math.acos
// function may be outside its domain of [-1, +1] which would return
// the value NaN which we do not want.
function acossafe(x) {
  if (x >= +1.0) return 0;
  if (x <= -1.0) return Math.PI;
  return Math.acos(x);
}

// Rotates a point about a fixed point at some angle 'a'
function rotatePoint(fp, pt, a) {
  var x = pt.x - fp.x;
  var y = pt.y - fp.y;
  var xRot = x * Math.cos(a) + y * Math.sin(a);
  var yRot = y * Math.cos(a) - x * Math.sin(a);
  return new Point(fp.x+xRot,fp.y+yRot);
}

// Given two circles this method finds the intersection
// point(s) of the two circles (if any exists)
function circleCircleIntersectionPoints(c1, c2) {

  var r, R, d, dx, dy, cx, cy, Cx, Cy;

  if (c1.r < c2.r) {
    r  = c1.r;  R = c2.r;
    cx = c1.x; cy = c1.y;
    Cx = c2.x; Cy = c2.y;
  } else {
    r  = c2.r; R  = c1.r;
    Cx = c1.x; Cy = c1.y;
    cx = c2.x; cy = c2.y;
  }

  // Compute the vector <dx, dy>
  dx = cx - Cx;
  dy = cy - Cy;

  // Find the distance between two points.
  d = Math.sqrt( dx*dx + dy*dy );

  // There are an infinite number of solutions
  // Seems appropriate to also return null
  if (d < EPS && Math.abs(R-r) < EPS) return [];

  // No intersection (circles centered at the 
  // same place with different size)
  else if (d < EPS) return [];

  var x = (dx / d) * R + Cx;
  var y = (dy / d) * R + Cy;
  var P = new Point(x, y);

  // Single intersection (kissing circles)
  if (Math.abs((R+r)-d) < EPS || Math.abs(R-(r+d)) < EPS) return [P];

  // No intersection. Either the small circle contained within 
  // big circle or circles are simply disjoint.
  if ( (d+r) < R || (R+r < d) ) return [];

  var C = new Point(Cx, Cy);
  var angle = acossafe((r*r-d*d-R*R)/(-2.0*d*R));
  var pt1 = rotatePoint(C, P, +angle);
  var pt2 = rotatePoint(C, P, -angle);
  return [pt1, pt2];

}
1
  • 2
    Demo link is broken
    – brainjam
    Commented Oct 29, 2019 at 17:33
0

Try this;

    def ri(cr1,cr2,cp1,cp2):
        int1=[]
        int2=[]
        ori=0
        if cp1[0]<cp2[0] and cp1[1]!=cp2[1]:
            p1=cp1
            p2=cp2
            r1=cr1
            r2=cr2
            if cp1[1]<cp2[1]:
                ori+=1
            elif cp1[1]>cp2[1]:
                ori+=2        
        elif cp1[0]>cp2[0] and cp1[1]!=cp2[1]:
            p1=cp2
            p2=cp1
            r1=cr2
            r2=cr1
            if p1[1]<p2[1]:
                ori+=1
            elif p1[1]>p2[1]:
                ori+=2
        elif cp1[0]==cp2[0]:
            ori+=4
            if cp1[1]>cp2[1]:
                p1=cp1
                p2=cp2
                r1=cr1
                r2=cr2
            elif cp1[1]<cp2[1]:
                p1=cp2
                p2=cp1
                r1=cr2
                r2=cr1
        elif cp1[1]==cp2[1]:
            ori+=3
            if cp1[0]>cp2[0]:
                p1=cp2
                p2=cp1
                r1=cr2
                r2=cr1
            elif cp1[0]<cp2[0]:
                p1=cp1
                p2=cp2
                r1=cr1
                r2=cr2
        if ori==1:#+
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            thta=math.degrees(math.acos(A/r1))
            rs=p2[1]-p1[1]
            rn=p2[0]-p1[0]
            gd=rs/rn
            yint=p1[1]-((gd)*p1[0])
            dty=calc_dist(p1,[0,yint])

            aa=p1[1]-yint
            bb=math.degrees(math.asin(aa/dty))
            d=90-bb
            e=180-d-thta
            g=(dty/math.sin(math.radians(e)))*math.sin(math.radians(thta))
            f=(g/math.sin(math.radians(thta)))*math.sin(math.radians(d))
            oty=yint+g
            h=f+r1
            i=90-e
            j=180-90-i
            l=math.sin(math.radians(i))*h
            k=math.cos(math.radians(i))*h
            iy2=oty-l
            ix2=k
            int2.append(ix2)
            int2.append(iy2)

            m=90+bb
            n=180-m-thta
            p=(dty/math.sin(math.radians(n)))*math.sin(math.radians(m))
            o=(p/math.sin(math.radians(m)))*math.sin(math.radians(thta))
            q=p+r1
            r=90-n
            s=math.sin(math.radians(r))*q
            t=math.cos(math.radians(r))*q
            otty=yint-o
            iy1=otty+s
            ix1=t
            int1.append(ix1)
            int1.append(iy1)
        elif ori==2:#-
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            thta=math.degrees(math.acos(A/r1))
            rs=p2[1]-p1[1]
            rn=p2[0]-p1[0]
            gd=rs/rn
            yint=p1[1]-((gd)*p1[0])
            dty=calc_dist(p1,[0,yint])

            aa=yint-p1[1]
            bb=math.degrees(math.asin(aa/dty))
            c=180-90-bb
            d=180-c-thta
            e=180-90-d
            f=math.tan(math.radians(e))*p1[0]
            g=math.sqrt(p1[0]**2+f**2)
            h=g+r1
            i=180-90-e
            j=math.sin(math.radians(e))*h
            jj=math.cos(math.radians(i))*h
            k=math.cos(math.radians(e))*h
            kk=math.sin(math.radians(i))*h
            l=90-bb
            m=90-e
            tt=l+m+thta
            n=(dty/math.sin(math.radians(m)))*math.sin(math.radians(thta))
            nn=(g/math.sin(math.radians(l)))*math.sin(math.radians(thta))
            oty=yint-n
            iy1=oty+j
            ix1=k
            int1.append(ix1)
            int1.append(iy1)

            o=bb+90
            p=180-o-thta
            q=90-p
            r=180-90-q
            s=(dty/math.sin(math.radians(p)))*math.sin(math.radians(o))
            t=(s/math.sin(math.radians(o)))*math.sin(math.radians(thta))
            u=s+r1
            v=math.sin(math.radians(r))*u
            vv=math.cos(math.radians(q))*u
            w=math.cos(math.radians(r))*u
            ww=math.sin(math.radians(q))*u
            ix2=v
            otty=yint+t
            iy2=otty-w
            int2.append(ix2)
            int2.append(iy2)

        elif ori==3:#y
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            b=math.sqrt(r1**2-A**2)
            int1.append(p1[0]+A)
            int1.append(p1[1]+b)
            int2.append(p1[0]+A)
            int2.append(p1[1]-b)
        elif ori==4:#x
            D=calc_dist(p1,p2)
            tr=r1+r2
            el=tr-D
            a=r1-el
            b=r2-el
            A=a+(el/2)
            B=b+(el/2)
            b=math.sqrt(r1**2-A**2)
            int1.append(p1[0]+b)
            int1.append(p1[1]-A)
            int2.append(p1[0]-b)
            int2.append(p1[1]-A)
        return [int1,int2]
    def calc_dist(p1,p2):
        return math.sqrt((p2[0] - p1[0]) ** 2 +
                         (p2[1] - p1[1]) ** 2)
2
  • It uses the Cartesian axis to calculate the intersect points using trigonometry. Just calculate the intercept point from a pair of coordinates and any origin location. then just work out the angles until you find the position of each. Commented May 30, 2018 at 15:02
  • Does that mean any Long/Lat in decimal degrees need to be converted to Cartesian coordinates? And the results are in Cartesian coordinates that then need to be re-converted into Long/Lat?
    – simpson
    Commented Sep 6, 2023 at 21:57
-2

A quick way to determine the intersection points P1 and P2 is to take the vector w between the center points A and B of the circles.

When now imagining a rectangle that is spanned by two vectors a and b between point A and P1, we can say

P1 = A + a + b
P2 = A + a - b

The only remaining question is, how long the vectors a and b are:

We know that |a|^2 + |b|^2 = r_A^2 and (|w| - |a|)^2 + |b|^2 = r_B^2, setting them equal and solving for |a| yields

|a| = (r_A^2 - r_b^2 + |w|^2) / (2|w|)
|b| = |b| = +- sqrt(r_A^2 - |a|^2)

Now the lengths of the vectors can be used to construct the actual vectors a and b by using the normalized vector w.

When implementing this idea, the solution is pretty straightforward

w = {
    x: B.x - A.x,
    y: B.y - A.y
}
    
d = hypot(w.x, w.y)

if (d <= A.r + B.r && abs(B.r - A.r) <= d) {

    w.x/= d;
    w.y/= d;

    a = (A.r * A.r - B.r * B.r + d * d) / (2 * d);
    b = Math.sqrt(A.r * A.r - a * a);

    P1 = {
        x: A.x + a * w.x - b * w.y,
        y: A.y + a * w.y + b * w.x
    }

    P2 = {
        x: A.x + a * w.x + b * w.y,
        y: A.y + a * w.y - b * w.x
    }

} else {
    P1 = P2 = null
}
0

Not the answer you're looking for? Browse other questions tagged or ask your own question.