Circle-circle intersection points

Question

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.

Answer 1

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 P₃ = (x₃, y₃) if they exist.

First calculate the distance d between the center of the circles. d = ||P₁ - P₀||.

• If d > r₀ + r₁ then there are no solutions, the circles are separate.

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

• If d = 0 and r₀ = r₁ then the circles are coincident and there are an infinite number of solutions.

Considering the two triangles P₀P₂P₃ and P₁P₂P₃ we can write

a² + h² = r₀² and b² + h² = r₁²

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

a = (r₀² - r₁² + d² ) / (2 d)

It can be readily shown that this reduces to r₀ when the two circles touch at one point, ie: d = r₀ + r₁ Solve for h by substituting a into the first equation, h² = r₀² - a²

So

P₂ = P₀ + a ( P₁ - P₀ ) / d

And finally, P₃ = (x₃,y₃) in terms of P₀ = (x₀,y₀), P₁ = (x₁,y₁) and P₂ = (x₂,y₂), is

x₃ = x₂ +- h ( y₁ - y₀ ) / d

y₃ = y₂ -+ h ( x₁ - x₀ ) / d

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

Answer 2

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));         }  };