在一个二维平面上给定两个圆的圆心横纵坐标、半径共 6 个参数, 求交点. 即实现下面的 C 函数
int intersect(struct circle_t const circles[], struct point_t intersections[]);其中 point_t 与 circle_t 的定义分别是
struct point_t {
double x;
double y;
};
struct circle_t {
point_t center;
double r;
};函数 intersect 的输入参数为两个圆, 返回交点个数, 另外, 交点详细信息将被存入另一个参数数组 intersections 中. 由于两个圆之多有 2 个交点, 因此该函数可以如下形式调用:
#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, 此时验证两组解即可.
下面是完整代码
#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;
}