问题描述:
给两个相交的圆,第一个圆的圆心为(x1,y1),半径为r1,第二个圆的圆心为(x2,y2),半径为r2,求两个圆的交点。
问题分析:
《训练指南》上求两圆交点的模板用了atan2,acos等库函数,精度损失比较严重。
下面介绍一种精度损失较小的做法:
原文地址
首先回顾一下圆的两种表示方法:
圆的标准方程:(x−x0)2+(y−y0)2=r2
圆的参数方程:{x=x0+r⋅cosθy=y0+r⋅sinθ
将第一个圆的参数方程带入第二个圆的标准方程:
(x1+r1cosθ−x2)2+(y1+r1sinθ−y2)2=r22
展开后得到:
2r1(x1−x2)cosθ+2r1(y1−y2)sinθ=r22−r21−(x1−x2)2−(y1−y2)2
令:
a=2r1(x1−x2)
b=2r1(y1−y2)
c=r22−r21−(x1−x2)2−(y1−y2)2
原式变为:
acosθ+bsinθ=c
令cosθ=x,sinθ=1−x2−−−−−√,关于sinθ的正负后面再判断。
代入方程,得到,ax+b1−x2−−−−−√=c
移项再两边平方,(ax−c)2=b2(1−x2)
整理得,(a2+b2)x2−2acx+c2−b2=0
下面就是解一元二次方程了。
将sinθ和cosθ代回到第一个圆的参数方程,能得到交点的坐标。
如果该点不在第二个圆上,说明sinθ是个负数,还需要对这个交点稍作调整。
还有一种特殊情况就是,如果已经确定有两个不同的交点,但解出来的cosθ值只有一个。
说明对应的sinθ值必然一正一负。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
typedef long double LD;
const LD eps = 1e-10;
int dcmp(LD x) {
if(fabs(x) < eps) return 0;
return x < 0 ? -1 : 1;
}
LD sqr(LD x) { return x * x; }
struct Point
{
LD x, y;
Point(LD x = 0, LD y = 0):x(x), y(y) {}
void read() { cin >> x >> y; }
};
Point operator - (const Point& A, const Point& B) {
return Point(A.x - B.x, A.y - B.y);
}
bool operator == (const Point& A, const Point& B) {
return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.x) == 0;
}
LD Dot(const Point& A, const Point& B) {
return A.x * B.x + A.y * B.y;
}
LD Length(const Point& A) { return sqrt(Dot(A, A)); }
struct Circle
{
Point c;
LD r;
Circle() {}
Circle(Point c, LD r):c(c), r(r) {}
};
int getCircleIntersection(Circle C1, Circle C2) {
LD &r1 = C1.r, &r2 = C2.r;
LD &x1 = C1.c.x, &x2 = C2.c.x, &y1 = C1.c.y, &y2 = C2.c.y;
LD d = Length(C1.c - C2.c);
if(dcmp(fabs(r1-r2) - d) > 0) return -1;
if(dcmp(r1 + r2 - d) < 0) return 0;
LD d2 = Dot(C1.c - C2.c, C1.c - C2.c);
LD a = r1*(x1-x2)*2, b = r1*(y1-y2)*2, c = r2*r2-r1*r1-d*d;
LD p = a*a+b*b, q = -a*c*2, r = c*c-b*b;
LD cosa, sina, cosb, sinb;
//One Intersection
if(dcmp(d - (r1 + r2)) == 0 || dcmp(d - fabs(r1 - r2)) == 0) {
cosa = -q / p / 2;
sina = sqrt(1 - sqr(cosa));
Point p(x1 + C1.r * cosa, y1 + C1.r * sina);
if(!OnCircle(p, C2)) p.y = y1 - C1.r * sina;
inter.push_back(p);
return 1;
}
//Two Intersections
LD delta = sqrt(q * q - p * r * 4);
cosa = (delta - q) / p / 2;
cosb = (-delta - q) / p / 2;
sina = sqrt(1 - sqr(cosa));
sinb = sqrt(1 - sqr(cosb));
Point p1(x1 + C1.r * cosa, y1 + C1.r * sina);
Point p2(x1 + C1.r * cosb, y1 + C1.r * sinb);
if(!OnCircle(p1, C2)) p1.y = y1 - C1.r * sina;
if(!OnCircle(p2, C2)) p2.y = y1 - C1.r * sinb;
if(p1 == p2) p1.y = y1 - C1.r * sina;
inter.push_back(p1);
inter.push_back(p2);
return 2;
}
给两个相交的圆,第一个圆的圆心为(x1,y1),半径为r1,第二个圆的圆心为(x2,y2),半径为r2,求两个圆的交点。
问题分析:
《训练指南》上求两圆交点的模板用了atan2,acos等库函数,精度损失比较严重。
下面介绍一种精度损失较小的做法:
原文地址
首先回顾一下圆的两种表示方法:
圆的标准方程:(x−x0)2+(y−y0)2=r2
圆的参数方程:{x=x0+r⋅cosθy=y0+r⋅sinθ
将第一个圆的参数方程带入第二个圆的标准方程:
(x1+r1cosθ−x2)2+(y1+r1sinθ−y2)2=r22
展开后得到:
2r1(x1−x2)cosθ+2r1(y1−y2)sinθ=r22−r21−(x1−x2)2−(y1−y2)2
令:
a=2r1(x1−x2)
b=2r1(y1−y2)
c=r22−r21−(x1−x2)2−(y1−y2)2
原式变为:
acosθ+bsinθ=c
令cosθ=x,sinθ=1−x2−−−−−√,关于sinθ的正负后面再判断。
代入方程,得到,ax+b1−x2−−−−−√=c
移项再两边平方,(ax−c)2=b2(1−x2)
整理得,(a2+b2)x2−2acx+c2−b2=0
下面就是解一元二次方程了。
将sinθ和cosθ代回到第一个圆的参数方程,能得到交点的坐标。
如果该点不在第二个圆上,说明sinθ是个负数,还需要对这个交点稍作调整。
还有一种特殊情况就是,如果已经确定有两个不同的交点,但解出来的cosθ值只有一个。
说明对应的sinθ值必然一正一负。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
typedef long double LD;
const LD eps = 1e-10;
int dcmp(LD x) {
if(fabs(x) < eps) return 0;
return x < 0 ? -1 : 1;
}
LD sqr(LD x) { return x * x; }
struct Point
{
LD x, y;
Point(LD x = 0, LD y = 0):x(x), y(y) {}
void read() { cin >> x >> y; }
};
Point operator - (const Point& A, const Point& B) {
return Point(A.x - B.x, A.y - B.y);
}
bool operator == (const Point& A, const Point& B) {
return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.x) == 0;
}
LD Dot(const Point& A, const Point& B) {
return A.x * B.x + A.y * B.y;
}
LD Length(const Point& A) { return sqrt(Dot(A, A)); }
struct Circle
{
Point c;
LD r;
Circle() {}
Circle(Point c, LD r):c(c), r(r) {}
};
int getCircleIntersection(Circle C1, Circle C2) {
LD &r1 = C1.r, &r2 = C2.r;
LD &x1 = C1.c.x, &x2 = C2.c.x, &y1 = C1.c.y, &y2 = C2.c.y;
LD d = Length(C1.c - C2.c);
if(dcmp(fabs(r1-r2) - d) > 0) return -1;
if(dcmp(r1 + r2 - d) < 0) return 0;
LD d2 = Dot(C1.c - C2.c, C1.c - C2.c);
LD a = r1*(x1-x2)*2, b = r1*(y1-y2)*2, c = r2*r2-r1*r1-d*d;
LD p = a*a+b*b, q = -a*c*2, r = c*c-b*b;
LD cosa, sina, cosb, sinb;
//One Intersection
if(dcmp(d - (r1 + r2)) == 0 || dcmp(d - fabs(r1 - r2)) == 0) {
cosa = -q / p / 2;
sina = sqrt(1 - sqr(cosa));
Point p(x1 + C1.r * cosa, y1 + C1.r * sina);
if(!OnCircle(p, C2)) p.y = y1 - C1.r * sina;
inter.push_back(p);
return 1;
}
//Two Intersections
LD delta = sqrt(q * q - p * r * 4);
cosa = (delta - q) / p / 2;
cosb = (-delta - q) / p / 2;
sina = sqrt(1 - sqr(cosa));
sinb = sqrt(1 - sqr(cosb));
Point p1(x1 + C1.r * cosa, y1 + C1.r * sina);
Point p2(x1 + C1.r * cosb, y1 + C1.r * sinb);
if(!OnCircle(p1, C2)) p1.y = y1 - C1.r * sina;
if(!OnCircle(p2, C2)) p2.y = y1 - C1.r * sinb;
if(p1 == p2) p1.y = y1 - C1.r * sina;
inter.push_back(p1);
inter.push_back(p2);
return 2;
}