[转]已知两圆圆心坐标及半径求两圆交点 (C语言|参数方程求解)
2011-08-19 19:16
543 查看
在一个二维平面上给定两个圆的圆心横纵坐标、半径共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>
struct point_t
{
double x, y;
};
struct
circle_t
{
struct point_t
center;
double r;
};
int
double_equals(double consta, 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;
double cos_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;
}
}
return 2;
}
int
main()
{
struct circle_t
circles[2];
struct point_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;
case 0:
printf("NO
INTERSECTION\n");
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^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>
struct point_t
{
double x, y;
};
struct
circle_t
{
struct point_t
center;
double r;
};
int
double_equals(double consta, 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;
double cos_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;
}
}
return 2;
}
int
main()
{
struct circle_t
circles[2];
struct point_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;
case 0:
printf("NO
INTERSECTION\n");
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;
}
相关文章推荐
- 已知两圆圆心坐标及半径求两圆交点 (C语言|参数方程求解)
- 已知两圆圆心坐标及半径求两圆交点 (C语言|参数方程求解)
- 已知两圆圆心坐标及半径求两圆交点 (C语言|参数方程求解)
- 已知两圆圆心坐标及半径求两圆交点 (C语言|参数方程求解)
- 已知两圆圆心坐标及半径求两圆交点
- [转]已知两圆圆心坐标及半径求两圆…
- c语言:有4个圆塔,已知圆心和半径,塔以外无建筑物。输入任一点坐标,求该点的建筑高度
- c语言:有4个圆塔,已知圆心和半径,塔以外无建筑物。输入任一点坐标,求该点的建筑高度
- 求内切圆的圆心和半径(已知三个点的坐标)
- WPF arcgis中已知圆心坐标和半径,求圆上一点的坐标(当前点的x坐标和圆心相等)
- 已知圆心,半径,角度,求圆上的点坐标
- 已知圆弧上两点坐标和半径求圆心坐标的算法(C++)
- 利用mfc动态画圆(已知圆心坐标,半径,以及每次转过的角度)
- 已知圆心,半径,角度,求圆上的点坐标。
- 已知圆心,半径,角度,求圆上的点坐标
- 已知圆心,半径,角度,求圆上的点坐标
- 已知圆心,半径,角度,求圆上的点坐标
- 已知两点坐标和半径求圆心坐标程序C++
- 已知圆弧上两点坐标及圆半径,计算圆心坐标
- 高德已知圆心,半径,计算圆弧坐标