About

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

在一个二维平面上给定两个圆的圆心横纵坐标、半径共 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


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
#include

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

Tags: Algorithm C Geometry

Loading comments


. Back to Tobary book
Tobary book - Copyright (C) ZhePlus @ Bit Focus
About