# Circles Intersection

Posted at 2012-01-10 12:30:40 | Updated at 2020-11-30 20:06:53

We are going to calculate points of intersections for 2 circels, given their centers and radiuses on a 2D plane. Say, we will implement the following C function
`int intersect(struct circle_t const circles[], struct point_t intersections[]);`
where `point_t` and `circle_t` are
`struct point_t {    double x;    double y;};struct circle_t {    point_t center;    double r;};`
The function `intersect` takes 2 circles as parameter, returns the number of points of intersection, and the details of the points will be put in parameter `intersection`. Since 2 circles may have up to 2 points of intersection, we may call that function in following way:
`#include <stdio.h>int main(void){    struct circle_t circles; // 2 circles    struct point_t points; // and up to 2 points    /* here we read the information via stdin     * in this format     *   (x0, y0, r0)     *   (x1, y1, r1)     */    scanf("%lf%lf%lf%lf%lf%lf",          &circles.center.x, &circles.center.y, &circles.r,          &circles.center.x, &circles.center.y, &circles.r);    /* 2 circles are the same */    /* NOTE: since x, y and r are all of type double,     *       use == to compare them may cause problem.     */    if (circles.center.x == circles.center.x     && circles.center.y == circles.center.y     && circles.r == circles.r)    {       puts("The circles are the same.");       return 0;    }    /* we call //intersect// here.     * both //circles// and //points// are an array of 2.     * if only one point of intersection, it will be stored in //points//.     */    switch (intersect(circles, points)) {        case 0:            puts("No intersection.");            break;        case 1:            printf("(%.3lf %.3lf)\n", points.x, points.y);            break;        case 2:            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",                   points.x, points.y,                   points.x, points.y);    }    return 0;}`
First, let's consider a simple situation: there is no intersection
`#include <math.h> // we need //sqrt// and //fabs//double distance_sqr(struct point_t const* a, struct point_t const* b){    return (a->x - b->x) * (a->x - b->x) + (a->y - b->y) * (a->y - b->y);}double distance(struct point_t const* a, struct point_t const* b){    return sqrt(distance_sqr(a, b));}int intersect(struct circle_t const circles[], struct point_t intersections[]){    double d = distance(&circles.center, &circles.center);    if (d > circles.r + circles.r     || d < fabs(circles.r - circles.r))    {        return 0;    }    /* more codes here */}`
To solute the rest of this problem, besides the usual way of representing circles like
• (x - x0)2 + (y - y0)2 = r2, where (x0, y0) is the center and r is the radius
we will also use parametric equation for that
• x = r * cosθ + x0
• y = r * sinθ + y0

Now, let (xa, ya) be one point of intersections, then for circle0 (that is (x0, y0) and r0), there is some θa that meets
• xa = r0 * cosθa + x0
• ya = r0 * sinθa + y0
then substitute (xa, ya) in the equation of circle1
(r0 * cosθa + x0 - x1)2 + (r0 * sinθa + y0 - y1)2 = r12

Don't panic though the equation is full of horrible calls to trigonometrical function. We do some tricks on the left and will get
//LEFT// = ((r0 * cosθa) + (x0 - x1))2 + ((r0 * sinθa) + (y0 - y1))2
= (r0 * cosθa)2 + (x0 - x1)2 + 2 * (r0 * cosθa) * (x0 - x1) +
(r0 * sinθa)2 + (y0 - y1)2 + 2 * (r0 * sinθa) * (y0 - y1)
= (r0 * cosθa)2 + (r0 * sinθa)2 + (x0 - x1)2 + (y0 - y1)2 +
2 * r0 * (x0 - x1) * cosθa + 2 * r0 * (y0 - y1) * sinθa
= r02 + (x0 - x1)2 + (y0 - y1)2 +
2 * r0 * (x0 - x1) * cosθa + 2 * r0 * (y0 - y1) * sinθa
= r12 = //RIGHT//

Then let
• a = 2 * r0 * (x0 - x1)
• b = 2 * r0 * (y0 - x1)
• c = r12 - (r02 + (x0 - x1)2 + (y0 - y1)2)
the equation will be transformed to
a * cosθa + b * sinθa = c
and then represent sinθa as ±(1 - (cosθa)2)1/2
b * (±(1 - (cosθa)2))1/2 = c - a * cosθa
square both sides (then the annoying ± is gone)
b2 * (1 - (cosθa)2) = c2 + a2 * (cosθa)2 - 2 * c * a * cosθa
=> b2 - b2 * (cosθa)2 = c2 + a2 * (cosθa)2 - 2 * c * a * cosθa
=> (a2 + b2) * (cosθa)2 - (2 * a * c) * cosθa + (c2 - b2) = 0

See, we are about to solute a quadratic equation of cosθa, which won't be very hard
cosθ = (±(q2 - 4 * p * s)1/2 - q) / (2 * p)
where
• p = a2 + b2
• q = -2 * a * c
• s = c2 - b2

You may notice the sign ±, so how may we determine which it should really be? We have 2 choices here, one is to calculate for more details, until we find whether it should be positive or negative, the other is, we just //verify// the solution, by using the magical computing power of the computer. In the latter way, what we will do is substitute the following 4 solutions
• (((q2 - 4 * p * s)1/2 - q) / (2 * p), (1 - (cosθa)2)1/2)
• ((-(q2 - 4 * p * s)1/2 - q) / (2 * p), (1 - (cosθa)2)1/2)
• (((q2 - 4 * p * s)1/2 - q) / (2 * p), -(1 - (cosθa)2)1/2)
• ((-(q2 - 4 * p * s)1/2 - q) / (2 * p), -(1 - (cosθa)2)1/2)
in the origin equation, to see which will meet.
By the way, if there is only one point of intersection, (q2 - 4 * p * s)1/2 will be 0 and the sign is not important.

Here is the complete code
`#include <stdio.h>#include <math.h>struct point_t {    double x;    double y;};struct circle_t {    struct point_t center;    double r;};double distance_sqr(struct point_t const* a, struct point_t const* b){    return (a->x - b->x) * (a->x - b->x) + (a->y - b->y) * (a->y - b->y);}double distance(struct point_t const* a, struct point_t const* b){    return sqrt(distance_sqr(a, b));}int intersect(struct circle_t const circles[], struct point_t intersections[]){    double a, b, c, p, q, s; // we have talked about those a, b, c, p, q, s    double cos_value, sin_value; // values for cosÎ¸a sinÎ¸a    double d // distance between centers of the two circles             = distance(&circles.center, &circles.center);    if (d > circles.r + circles.r        || d < fabs(circles.r - circles.r))    {        return 0;    }    a = 2.0 * circles.r * (circles.center.x - circles.center.x);    b = 2.0 * circles.r * (circles.center.y - circles.center.y);    c = circles.r * circles.r - circles.r * circles.r        - distance_sqr(&circles.center, &circles.center);    p = a * a + b * b;    q = -2.0 * a * c;    // there is only one point of intersection    if (d == circles.r + circles.r     || d == fabs(circles.r - circles.r))    {        cos_value = -q / p / 2.0;        sin_value = sqrt(1 - cos_value * cos_value);        intersections.x = circles.r * cos_value + circles.center.x;        intersections.y = circles.r * sin_value + circles.center.y;        // verify the solution here, if it is not, negate sinÎ¸a        if (distance_sqr(&intersections, &circles.center)            != circles.r * circles.r)        {            intersections.y = circles.center.y                                 - circles.r * sin_value;        }        return 1;    }    s = c * c - b * b;    cos_value = (sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;    cos_value = (-sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;    sin_value = sqrt(1 - cos_value * cos_value);    sin_value = sqrt(1 - cos_value * cos_value);    intersections.x = circles.r * cos_value + circles.center.x;    intersections.x = circles.r * cos_value + circles.center.x;    intersections.y = circles.r * sin_value + circles.center.y;    intersections.y = circles.r * sin_value + circles.center.y;    // verify and correct both solutions, still by negating sinÎ¸a    if (distance_sqr(&intersections, &circles.center)        != circles.r * circles.r)    {        intersections.y = circles.center.y - circles.r * sin_value;    }    if (distance_sqr(&intersections, &circles.center)        != circles.r * circles.r)    {        intersections.y = circles.center.y - circles.r * sin_value;    }    /* if we got 2 same solutions, pick one and negate its Y coord     * it is because cosÎ¸a is used to represent X coord of the solutions,     * so when we got same solutions the problem must be on the Y coord.     */    if (intersections.y == intersections.y     && intersections.x == intersections.x)    {        intersections.y = -intersections.y;    }    return 2;}int main(void){    struct circle_t circles;    struct point_t points;    scanf("%lf%lf%lf%lf%lf%lf",          &circles.center.x, &circles.center.y, &circles.r,          &circles.center.x, &circles.center.y, &circles.r);    if (circles.center.x == circles.center.x     && circles.center.y == circles.center.y     && circles.r == circles.r)    {       puts("The circles are the same.");       return 0;    }    switch (intersect(circles, points)) {        case 0:            puts("No intersection.");            break;        case 1:            printf("(%.3lf %.3lf)\n", points.x, points.y);            break;        case 2:            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",                   points.x, points.y,                   points.x, points.y);    }    return 0;}`

Post tags: Your comment will be licensed under CC-NC-ND 3.0 Back to Bit Focus NijiPress - Copyright (C) Neuron Teckid @ Bit Focus About this site