About
RSS

Bit Focus


Circles Intersection

Posted at 2012-01-10 12:30:40 | Updated at 2018-10-22 20:49:43

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]; // 2 circles
    struct point_t points[2]; // 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[0].center.x, &circles[0].center.y, &circles[0].r,
          &circles[1].center.x, &circles[1].center.y, &circles[1].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[0].center.x == circles[1].center.x
     && circles[0].center.y == circles[1].center.y
     && circles[0].r == circles[1].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[0]//.
     */
    switch (intersect(circles, points)) {
        case 0:
            puts("No intersection.");
            break;
        case 1:
            printf("(%.3lf %.3lf)\n", points[0].x, points[0].y);
            break;
        case 2:
            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",
                   points[0].x, points[0].y,
                   points[1].x, points[1].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[0].center, &circles[1].center);
    if (d > circles[0].r + circles[1].r
     || d < fabs(circles[0].r - circles[1].r))
    {
        return 0;
    }
    /* more codes here */
}
To solute the rest of this problem, besides the usual way of representing circles like
we will also use parametric equation for that

Now, let (xa, ya) be one point of intersections, then for circle0 (that is (x0, y0) and r0), there is some θa that meets
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
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

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
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[2], sin_value[2]; // values for cosθa sinθa
    double d // distance between centers of the two circles
             = distance(&circles[0].center, &circles[1].center);

    if (d > circles[0].r + circles[1].r
        || d < fabs(circles[0].r - circles[1].r))
    {
        return 0;
    }

    a = 2.0 * circles[0].r * (circles[0].center.x - circles[1].center.x);
    b = 2.0 * circles[0].r * (circles[0].center.y - circles[1].center.y);
    c = circles[1].r * circles[1].r - circles[0].r * circles[0].r
        - distance_sqr(&circles[0].center, &circles[1].center);
    p = a * a + b * b;
    q = -2.0 * a * c;

    // there is only one point of intersection
    if (d == circles[0].r + circles[1].r
     || d == fabs(circles[0].r - circles[1].r))
    {
        cos_value[0] = -q / p / 2.0;
        sin_value[0] = sqrt(1 - cos_value[0] * cos_value[0]);

        intersections[0].x = circles[0].r * cos_value[0] + circles[0].center.x;
        intersections[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

        // verify the solution here, if it is not, negate sinθa
        if (distance_sqr(&intersections[0], &circles[1].center)
            != circles[1].r * circles[1].r)
        {
            intersections[0].y = circles[0].center.y
                                 - circles[0].r * sin_value[0];
        }
        return 1;
    }

    s = c * c - b * b;
    cos_value[0] = (sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;
    cos_value[1] = (-sqrt(q * q - 4.0 * p * s) - q) / p / 2.0;
    sin_value[0] = sqrt(1 - cos_value[0] * cos_value[0]);
    sin_value[1] = sqrt(1 - cos_value[1] * cos_value[1]);

    intersections[0].x = circles[0].r * cos_value[0] + circles[0].center.x;
    intersections[1].x = circles[0].r * cos_value[1] + circles[0].center.x;
    intersections[0].y = circles[0].r * sin_value[0] + circles[0].center.y;
    intersections[1].y = circles[0].r * sin_value[1] + circles[0].center.y;

    // verify and correct both solutions, still by negating sinθa
    if (distance_sqr(&intersections[0], &circles[1].center)
        != circles[1].r * circles[1].r)
    {
        intersections[0].y = circles[0].center.y - circles[0].r * sin_value[0];
    }
    if (distance_sqr(&intersections[1], &circles[1].center)
        != circles[1].r * circles[1].r)
    {
        intersections[1].y = circles[0].center.y - circles[0].r * sin_value[1];
    }

    /* 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[0].y == intersections[1].y
     && intersections[0].x == intersections[1].x)
    {
        intersections[1].y = -intersections[1].y;
    }
    return 2;
}

int main(void)
{
    struct circle_t circles[2];
    struct point_t points[2];

    scanf("%lf%lf%lf%lf%lf%lf",
          &circles[0].center.x, &circles[0].center.y, &circles[0].r,
          &circles[1].center.x, &circles[1].center.y, &circles[1].r);

    if (circles[0].center.x == circles[1].center.x
     && circles[0].center.y == circles[1].center.y
     && circles[0].r == circles[1].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[0].x, points[0].y);
            break;
        case 2:
            printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",
                   points[0].x, points[0].y,
                   points[1].x, points[1].y);
    }
    return 0;
}

Post tags:   C  Algorithm  Geometry

Leave a comment:




Creative Commons License 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