[转]已知两圆圆心坐标及半径求两圆交点 ，C语言|参数方程求解

在一个二维平面上给定两个圆的圆心横纵坐标、半径共6个参数, 求交点. 这个问题无非是解二元二次方程组. 普通二元二次方程联立消元求解的困难在于, 中间过程里的系数会变得非常复杂, 从而导致容易出错---因为公式毕竟还是要人来推导, 人的出错率比计算机要高得多得多---改用圆的参数方程求解, 可以在显著地减轻这个负担.

现在谈谈这问题的求解过程. 选择圆的参数方程的好处是方程是一次的, 化简方便, 虽然是三角函数方程并且既有正弦也有余弦, 不过到最后可以很方便地求出来.

(下面分析中x^y表示x的y次方)

大家还记得圆的参数方程吧, 圆心 (x0, y0), 半径为 r 的圆的参数方程是:

x = r * cosθ + x0

y = r * sinθ + y0

假设现在两圆参数x1, y1, r1, x2, y2, r2(这些分别表示, 咳, 有谁看不出来它们分别表示什么吗?), 设交点为 (x, y), 代入其中一个圆中的参数方程有

x = r1 * cosθ + x1 且 y = r1 * sinθ + y1

代入另一圆的标准方程, 得到

(r1 * cosθ + x1 - x2)^2 + (r1 * sinθ + y1 - y2)^2 = r2^2

是的, 看起来有关于正余弦二次项, 不过不要惊慌, 展开合并同类项之后, 正好这两项会合并成常数:

左边= (r1 * cosθ)^2 + (r1 * sinθ)^2 +

2 * r1 * (x1 - x2) * cosθ + 2 * r1 * (y1 - y2) * sinθ

= r2^2 - (x1 - x2)^2 - (y1 - y2)^2 =右边

这样就好办了, 把 r1^2 转移到等式右边, 令:

a = 2 * r1 * (x1 - x2)

b = 2 * r1 * (y1 - y2)

c = r2^2 - r1^2 - (x1 - x2)^2 - (y1 - y2)^2

那么方程便成为:

a * cosθ + b * sinθ = c

用 (1 - (cosθ)^2)^(1 / 2) 表示sinθ, 令：

p = a^2 + b^2

q = -2 * a * c

r = c^2 - b^2

便化为一个一元二次方程, 解得:

cosθ = (±(q^2 - 4 * p * r)^(1 / 2) - q) / (2 * p)

然而到此为止还没有结束, 因为刚才将三角方程转化为二次方程时, 等式两边平方了一次, 如果直接这样求对应角的正弦值, 符号总会为正. 为了将纵坐标正确解出, 必须变角. 那么如何变角？方法当然很多, 诱导公式, 或者反过头去把方程变一下, 以正弦作为未知数, 但是这些方法都比较复杂. 在这里可以选择另一种方案, 那就是用验证代替求解: 验证所得的解是不是根, 如果不是, 将纵坐标反号便可以了. 最后注意一下两解的横坐标相同的情况, 这样要先输出正的再输出负的.

下面上代码

#include<math.h> // sqrt fabs

//点

struct point_t {

double x, y;

};

//圆

struct circle_t {

struct point_t center;

double r;

};

//浮点数判同

int double_equals(double const a, double const b)

{

static const double ZERO = 1e-9;

return fabs(a - b) < ZERO;

}

//两点之间距离的平方

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 insect(struct circle_t circles[], struct point_t points[])

{

double d, a, b, c, p, q, r; // a, b, c, p, q, r与上面分析中的量一致

double cos_value[2], sin_value[2]; //交点的在circles[0]上对应的正余弦取值

//余弦值cos_value就是分析中的cosθ

if (double_equals(circles[0].center.x, circles[1].center.x)

&& double_equals(circles[0].center.y, circles[1].center.y)

&& double_equals(circles[0].r, circles[1].r)) {

return -1;

}

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 (double_equals(d, circles[0].r + circles[1].r)

|| double_equals(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]);

points[0].x = circles[0].r * cos_value[0] + circles[0].center.x;

points[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

//在这里验证解是否正确,如果不正确,则将纵坐标符号进行变换

if(!double_equals(distance_sqr(&points[0], &circles[1].center),

circles[1].r * circles[1].r)) {

points[0].y = circles[0].center.y - circles[0].r * sin_value[0];

}

return 1;

}

r = c * c - b * b;

cos_value[0] = (sqrt(q * q - 4.0 * p * r) - q) / p / 2.0;

cos_value[1] = (-sqrt(q * q - 4.0 * p * r) - 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]);

points[0].x = circles[0].r * cos_value[0] + circles[0].center.x;

points[1].x = circles[0].r * cos_value[1] + circles[0].center.x;

points[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

points[1].y = circles[0].r * sin_value[1] + circles[0].center.y;

//验证解是否正确,两个解都需要验证.

if (!double_equals(distance_sqr(&points[0], &circles[1].center),

circles[1].r * circles[1].r)) {

points[0].y = circles[0].center.y - circles[0].r * sin_value[0];

}

if (!double_equals(distance_sqr(&points[1], &circles[1].center),

circles[1].r * circles[1].r)) {

points[1].y = circles[0].center.y - circles[0].r * sin_value[1];

}

//如果求得的两个点坐标相同,则必然其中一个点的纵坐标反号可以求得另一点坐标

if (double_equals(points[0].y, points[1].y)

&& double_equals(points[0].x, points[1].x)) {

if(points[0].y > 0) {

points[1].y = -points[1].y;

} else {

points[0].y = -points[0].y;

}

}

return 2;

}

一个完整的程序:

#include<stdio.h>

#include<math.h>

structpoint_t {

doublex, y;

};

structcircle_t {

structpoint_t center;

doubler;

};

intdouble_equals(doubleconsta,doubleconstb)

{

static const doubleZERO = 1e-9;

returnfabs(a - b) < ZERO;

}

doubledistance_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);

}

doubledistance(struct point_t const* a, struct point_t const* b)

{

returnsqrt(distance_sqr(a, b));

}

intinsect(structcircle_t circles[],structpoint_t points[])

{

doubled, a, b, c, p, q, r;

doublecos_value[2], sin_value[2];

if(double_equals(circles[0].center.x, circles[1].center.x)

&& double_equals(circles[0].center.y, circles[1].center.y)

&& double_equals(circles[0].r, circles[1].r)) {

return -1;

}

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(double_equals(d, circles[0].r + circles[1].r)

|| double_equals(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]);

points[0].x = circles[0].r * cos_value[0] + circles[0].center.x;

points[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

if(!double_equals(distance_sqr(&points[0], &circles[1].center),

circles[1].r * circles[1].r)) {

points[0].y = circles[0].center.y - circles[0].r * sin_value[0];

}

return 1;

}

r = c * c - b * b;

cos_value[0] = (sqrt(q * q - 4.0 * p * r) - q) / p / 2.0;

cos_value[1] = (-sqrt(q * q - 4.0 * p * r) - 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]);

points[0].x = circles[0].r * cos_value[0] + circles[0].center.x;

points[1].x = circles[0].r * cos_value[1] + circles[0].center.x;

points[0].y = circles[0].r * sin_value[0] + circles[0].center.y;

points[1].y = circles[0].r * sin_value[1] + circles[0].center.y;

if(!double_equals(distance_sqr(&points[0], &circles[1].center),

circles[1].r * circles[1].r)) {

points[0].y = circles[0].center.y - circles[0].r * sin_value[0];

}

if(!double_equals(distance_sqr(&points[1], &circles[1].center),

circles[1].r * circles[1].r)) {

points[1].y = circles[0].center.y - circles[0].r * sin_value[1];

}

if(double_equals(points[0].y, points[1].y)

&& double_equals(points[0].x, points[1].x)) {

if(points[0].y > 0) {

points[1].y = -points[1].y;

}else{

points[0].y = -points[0].y;

}

}

return2;

}

intmain()

{

structcircle_t circles[2];

structpoint_t points[2];

while(EOF!= 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)) {

switch(insect(circles, points)) {

case-1:

printf("THE CIRCLES ARE THE SAME\n");

break;

case0:

printf("NO INTERSECTION\n");

break;

case1:

printf("(%.3lf %.3lf)\n", points[0].x, points[0].y);

break;

case2:

printf("(%.3lf %.3lf) (%.3lf %.3lf)\n",

points[0].x, points[0].y,

points[1].x, points[1].y);

}

}

return0;

}