About
RSS

Bit Focus


已知两圆圆心坐标及半径求两圆交点

Posted at 2009-11-15 14:56:09 | Updated at 2018-05-12 05:27:03

在一个二维平面上给定两个圆的圆心横纵坐标、半径共 6 个参数, 求交点. 即实现下面的 C 函数

int intersect(struct circle_t const circles[], struct point_t intersections[]);

其中 point_tcircle_t 的定义分别是

Code Snippet 0-0

struct point_t {
    double x;
    double y;
};

struct circle_t {
    point_t center;
    double r;
};

函数 intersect 的输入参数为两个圆, 返回交点个数, 另外, 交点详细信息将被存入另一个参数数组 intersections 中. 由于两个圆之多有 2 个交点, 因此该函数可以如下形式调用:

Code Snippet 0-1

#include <stdio.h>

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

    /* 从 stdin 输入圆参数
     * 按照如下格式
     * $(x_0, y_0, r_0)$
     * $(x_1, y_1, r_1)$
     */
    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);

    // 如果两个圆相同
    // *注意:* 由于 x, y, 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;
}

用程序实现求交点方法可以有多种, 比较极端的甚至可以用到二分逼近, 反正计算机的计算能力跟笔算是两个世界; 不过本文还是老实一点, 仍试着来解方程. 在求解过程中, 除了用到圆的标准方程

$ (x - x_0)^2 + (y - y_0)^2 = r_0^2

(其中 $(x_0, y_0)$ 是其中一个圆的圆心, $r_0$ 是其半径; 当然对另一圆也是如此, 此处省略之) 圆的参数方程也将会被用到

$ x = r_0 * cos(A) + x_0

$ y = r_0 * sin(A) + y_0

现在, 设其中其中至少有一个交点 (没有交点的情况可以简单利用圆心距大于半径之和来判定), 换言之存在某个 $A_s$ 使得存在对应的点 $(x_s, y_s)$ 于两个圆的方程都成立, 即

$ x_s = r_0 * cos(A_s) + x_0

$ y_s = r_0 * sin(A_s) + y_0

$ (x_s - x_1)^2 + (y_s - y_1)^2 = r_1^2

那么, 将上面两个等式代入第三个式子, 就得到

$ (r_0 * cos(A_s) + x_0 - x_1)^2 + (r_0 * sin(A_s) + y_0 - y_1)^2 = r_1^2

虽然这里满是恐怖的三角函数, 以及它们的平方, 但别惊慌, 通过一些变换后可令其简化很多

$ ((r_0 * cos(A_s)) + (x_0 - x_1))^2 + ((r_0 * sin(A_s)) + (y_0 - y_1))^2

$ = (r_0 * cos(A_s))^2 + (x_0 - x_1)^2 + 2 * (r_0 * cos(A_s)) * (x_0 - x_1) +

$ (r_0 * sin(A_s))^2 + (y_0 - y_1)^2 + 2 * (r_0 * sin(A_s)) * (y_0 - y_1)

$ = (r_0 * cos(A_s))^2 + (r_0 * sin(A_s))^2 + (x_0 - x_1)^2 + (y_0 - y_1)^2 +

$ 2 * r_0 * (x_0 - x_1) * cos(A_s) + 2 * r_0 * (y_0 - y_1) * sin(A_s)

$ = r_0^2 + (x_0 - x_1)^2 + (y_0 - y_1)^2 +

$ 2 * r_0 * (x_0 - x_1) * cos(A_s) + 2 * r_0 * (y_0 - y_1) * sin(A_s)

$ = r_1^2

现在, 令

$ a = 2 * r_0 * (x_0 - x_1)

$ b = 2 * r_0 * (y_0 - x_1)

$ c = r_1^2 - (r_0^2 + (x_0 - x_1)^2 + (y_0 - y_1)^2)

则方程被转换成了

$ a * cos(A_s) + b * sin(A_s) = c

然后代入 $sin(A_s) = +-sqrt(1 - (cos(A_s))^2)$

$ b * (+-sqrt(1 - (cos(A_s))^2)) = c - a * cos(A_s)

两边同时平方 (± 就消失了) 得

$ b^2 * (1 - (cos(A_s))^2) = c^2 + a^2 * (cos(A_s))^2 - 2 * c * a * cos(A_s)

$ b^2 - b^2 * (cos(A_s))^2 = c^2 + a^2 * (cos(A_s))^2 - 2 * c * a * cos(A_s)

$ (a^2 + b^2) * (cos(A_s))^2 - (2 * a * c) * cos(A_s) + (c^2 - b^2) = 0

这样就变换为了 $cos(A_s)$ 的一元二次方程; 再令

$ p = a^2 + b^2

$ q = -2 * a * c

$ s = c^2 - b^2

很容易将方程的解表示为

$ cos(A_s) = (+-sqrt(q^2 - 4 * p * s) - q) / (2 * p)

对于计算机求解的另一个好处在于, 以上使用了许多 (a, b, c, p, q, s) 中间量来表示方程, 在写代码时, 同样可以同名临时变量替换之. (另吐槽一下数学里基本上不怎么认真对待变量和函数命名, 无论什么解答过程通通是 f, g, h, 函数形参通通是 x, y, z, 太容易混淆了 :-)

然而, 这个解集中又出现了讨厌的 ±, 应该如何断定此时它到底是正还是负呢? 最直接的当然是列更多的方程, 找到更准确的解. 不过, 现在到底是在计算机上实现解方程, 而有胜过笔纸几个数量级的超强计算能力, 继续列方程求解是完全没必要的. 此处采用暴力计算验证的方式, 即, 将 $cos(A_s)$ 分别取正负号的情况, 以及每个余弦值对应两个正弦值总共 4 组解代入原来的圆方程, 验证一下哪些是正确的即可. 而如果圆本身只有一个交点, 那么 $sqrt(q^2 - 4 * p * s)$ 必然将会是 0, 此时验证两组解即可.

下面是完整代码

Code Snippet 0-2

#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; // 如同在上面讨论到的中间量 a, b, c, p, q, s
    double cos_value[2], sin_value[2]; // 存 $cos(A_s)$ 和 $sin(A_s)$ 的值
    // 圆心距离
    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;
    }

    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;

    // 只有一个交点
    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;

        // 验证解, 如果验证发现有误, 则将 $sin(A_s)$ 取反
        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;
    cos_value[1] = (-sqrt(q * q - 4.0 * p * s) - q) / p / 2;
    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[0].y = circles[0].r * sin_value[0] + circles[0].center.y;
    intersections[1].x = circles[0].r * cos_value[1] + circles[0].center.x;
    intersections[1].y = circles[0].r * sin_value[1] + circles[0].center.y;

    // 同样, 验证失败时, 将 $sin(A_s)$ 取反
    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 (intersections[0].y == intersections[1].y
        && intersections[0].x == intersections[1].x)
    {
        intersections[1].y = circles[0].center.y - circles[0].r * sin_value[1];
    }
    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:   Algorithm  C  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