如何计算两个圆的交点。我期望在所有情况下都存在两个、一个或者零个交点。
我有每个圆心点的 x 和 y 坐标以及每个圆的半径。
首选 Python 的答案，但任何可行的算法都可以接受。
6个回答
123
两个圆的交点
作者：Paul Bourke
以下说明如何在平面上找到两个圆的交点，使用以下符号表示。目标是找到两个点P3 = (x3, y3)，如果它们存在的话。首先计算圆心之间的距离d = ||P1 - P0||。
如果d > r0 + r1，则没有解，圆是分开的。
如果d < |r0 - r1|，则没有解，因为一个圆包含在另一个圆内。
如果d = 0且r0 = r1，则圆重合，并且有无限多个解。
考虑三角形P0P2P3和P1P2P3，我们可以写出a² + h² = r0²和b² + h² = r1²。
使用d = a + b，我们可以解出a，即a = (r0² - r1² + d²) / (2d)。
通过将a代入第一个方程式中，我们可以解出h² = r0² - a²。
所以P2 = P0 + a(P1 - P0) / d。
最后，P3 = (x3, y3)用P0 = (x0, y0)，P1 = (x1, y1)和P2 = (x2, y2)表示为：
x3 = x2 ± h(y1 - y0) / d
y3 = y2 ± h(x1 - x0) / d
来源：http://paulbourke.net/geometry/circlesphere/。这篇文章与圆和球体的几何学相关，介绍了一些公式和算法。
- Jacob
10
1链接已失效，请上传新的链接。此外，有人可以解释一下如何从先前的公式得出a=(r02-r12+d2)/(2*d)吗？也许这很明显，但我还没有理解。谢谢！ - David Neto
1@DavidNeto 我修复了链接。此外，$a$ 是通过解 $h^2$ 并将其插入到另一个方程 $h^2 = r_0^2 - a^2$ 中获得的，因此现在您可以在这里解决 $a$：$b^2 + r_0^2 - a^2 = r_1^2$（请记住 $d = a + b$）。 - nachocab
2很好，这是一个很好的解释。这是Python中的算法：
https://gist.github.com/xaedes/974535e71009fa8f090e - xaedes
6我们可以使用 d = a + b 来解出 a。如何操作呢？这一步看起来对我来说就像黑魔法一样。用哪个方程式来解出 a 呢？ - Makogan
9
b^2 = r1^2 - h^2 和 h^2 = r0^2 - a^2b^2 = r1^2 - r0^2 + a^2r0^2 - r1^2 = a^2 - b^2
d^2 + r0^2 - r1^2 = a^2 - b^2 + d^2现在，使用等式 (a + b)^2 = d^2 = a^2 + b^2 + 2abd^2 + r0^2 - r1^2 = a^2 - b^2 + a^2 + b^2 + 2ab
d^2 + r0^2 - r1^2 = 2 * a^2 + 2ab从 a + b = d 可得，如果两边乘以 2a，则有 2 * a^2 + 2ab = 2da最终结果为，
d^2 + r0^2 - r1^2 = 2da - Y N显示剩余5条评论
17
为什么不直接使用你最喜欢的过程式语言（或可编程计算器！）如下所示，只需使用7行。
假设给定P0坐标（x0，y0）、P1坐标（x1，y1）、r0和r1，并且要找到P3坐标（x3，y3）：
d=sqr((x1-x0)^2 + (y1-y0)^2)
a=(r0^2-r1^2+d^2)/(2*d)
h=sqr(r0^2-a^2)
x2=x0+a*(x1-x0)/d
y2=y0+a*(y1-y0)/d
x3=x2+h*(y1-y0)/d // also x3=x2-h*(y1-y0)/d
y3=y2-h*(x1-x0)/d // also y3=y2+h*(x1-x0)/d
- CuriousMarc
4
如果(r0^2-a^2) < 0会怎样？ - ablaze
啊，是的。请参考Paul Bourke的原始答案，首先要确保您进行测试以避免d > r0 + r1（圆太远），d < |r0 - r1|（一个圆在另一个圆内）和（d = 0且r0 = r1）（圆重合）。除了这些情况之外，您将有一个解决方案，并且您提到的表达式将为正数。 - CuriousMarc
所以这是一个绝对的差异，我错过了。谢谢！ - ablaze
@ablaze：不是很对。我应该简单地说，如果(r0^2-a^2) <= 0，则两个圆之间没有交点。 - CuriousMarc
15
这是基于Paul Bourke的文章的C++实现。 它只在有两个交点时起作用，否则可能会返回NaN NAN NAN NAN。
class Point{
public:
float x, y;
Point(float px, float py) {
x = px;
y = py;
}
Point sub(Point p2) {
return Point(x - p2.x, y - p2.y);
}
Point add(Point p2) {
return Point(x + p2.x, y + p2.y);
}
float distance(Point p2) {
return sqrt((x - p2.x)*(x - p2.x) + (y - p2.y)*(y - p2.y));
}
Point normal() {
float length = sqrt(x*x + y*y);
return Point(x/length, y/length);
}
Point scale(float s) {
return Point(x*s, y*s);
}
};
class Circle {
public:
float x, y, r, left;
Circle(float cx, float cy, float cr) {
x = cx;
y = cy;
r = cr;
left = x - r;
}
pair<Point, Point> intersections(Circle c) {
Point P0(x, y);
Point P1(c.x, c.y);
float d, a, h;
d = P0.distance(P1);
a = (r*r - c.r*c.r + d*d)/(2*d);
h = sqrt(r*r - a*a);
Point P2 = P1.sub(P0).scale(a/d).add(P0);
float x3, y3, x4, y4;
x3 = P2.x + h*(P1.y - P0.y)/d;
y3 = P2.y - h*(P1.x - P0.x)/d;
x4 = P2.x - h*(P1.y - P0.y)/d;
y4 = P2.y + h*(P1.x - P0.x)/d;
return pair<Point, Point>(Point(x3, y3), Point(x4, y4));
}
};
- Rusty Rob
3
2不行 ;) - 但是这听起来像是一个有趣的练习，你可以自己尝试一下 ;p - Rusty Rob
我可以问一下
Point P2 = P1.sub(P0).scale(a/d).add(P0); 这行代码的目的是什么？ - padawan这是一个关于编程的内容：P2 = P0 + a ( P1 - P0 ) / d 的排列。 - Rusty Rob
13
这是一个使用向量的Javascript实现。代码有详细的注释，您应该能够理解它。这是原始来源
查看演示：这里:
// Let EPS (epsilon) be a small value
var EPS = 0.0000001;
// Let a point be a pair: (x, y)
function Point(x, y) {
this.x = x;
this.y = y;
}
// Define a circle centered at (x,y) with radius r
function Circle(x,y,r) {
this.x = x;
this.y = y;
this.r = r;
}
// Due to double rounding precision the value passed into the Math.acos
// function may be outside its domain of [-1, +1] which would return
// the value NaN which we do not want.
function acossafe(x) {
if (x >= +1.0) return 0;
if (x <= -1.0) return Math.PI;
return Math.acos(x);
}
// Rotates a point about a fixed point at some angle 'a'
function rotatePoint(fp, pt, a) {
var x = pt.x - fp.x;
var y = pt.y - fp.y;
var xRot = x * Math.cos(a) + y * Math.sin(a);
var yRot = y * Math.cos(a) - x * Math.sin(a);
return new Point(fp.x+xRot,fp.y+yRot);
}
// Given two circles this method finds the intersection
// point(s) of the two circles (if any exists)
function circleCircleIntersectionPoints(c1, c2) {
var r, R, d, dx, dy, cx, cy, Cx, Cy;
if (c1.r < c2.r) {
r = c1.r; R = c2.r;
cx = c1.x; cy = c1.y;
Cx = c2.x; Cy = c2.y;
} else {
r = c2.r; R = c1.r;
Cx = c1.x; Cy = c1.y;
cx = c2.x; cy = c2.y;
}
// Compute the vector <dx, dy>
dx = cx - Cx;
dy = cy - Cy;
// Find the distance between two points.
d = Math.sqrt( dx*dx + dy*dy );
// There are an infinite number of solutions
// Seems appropriate to also return null
if (d < EPS && Math.abs(R-r) < EPS) return [];
// No intersection (circles centered at the
// same place with different size)
else if (d < EPS) return [];
var x = (dx / d) * R + Cx;
var y = (dy / d) * R + Cy;
var P = new Point(x, y);
// Single intersection (kissing circles)
if (Math.abs((R+r)-d) < EPS || Math.abs(R-(r+d)) < EPS) return [P];
// No intersection. Either the small circle contained within
// big circle or circles are simply disjoint.
if ( (d+r) < R || (R+r < d) ) return [];
var C = new Point(Cx, Cy);
var angle = acossafe((r*r-d*d-R*R)/(-2.0*d*R));
var pt1 = rotatePoint(C, P, +angle);
var pt2 = rotatePoint(C, P, -angle);
return [pt1, pt2];
}
- will.fiset
1
2演示链接已损坏。 - brainjam
0
试试这个；
def ri(cr1,cr2,cp1,cp2):
int1=[]
int2=[]
ori=0
if cp1[0]<cp2[0] and cp1[1]!=cp2[1]:
p1=cp1
p2=cp2
r1=cr1
r2=cr2
if cp1[1]<cp2[1]:
ori+=1
elif cp1[1]>cp2[1]:
ori+=2
elif cp1[0]>cp2[0] and cp1[1]!=cp2[1]:
p1=cp2
p2=cp1
r1=cr2
r2=cr1
if p1[1]<p2[1]:
ori+=1
elif p1[1]>p2[1]:
ori+=2
elif cp1[0]==cp2[0]:
ori+=4
if cp1[1]>cp2[1]:
p1=cp1
p2=cp2
r1=cr1
r2=cr2
elif cp1[1]<cp2[1]:
p1=cp2
p2=cp1
r1=cr2
r2=cr1
elif cp1[1]==cp2[1]:
ori+=3
if cp1[0]>cp2[0]:
p1=cp2
p2=cp1
r1=cr2
r2=cr1
elif cp1[0]<cp2[0]:
p1=cp1
p2=cp2
r1=cr1
r2=cr2
if ori==1:#+
D=calc_dist(p1,p2)
tr=r1+r2
el=tr-D
a=r1-el
b=r2-el
A=a+(el/2)
B=b+(el/2)
thta=math.degrees(math.acos(A/r1))
rs=p2[1]-p1[1]
rn=p2[0]-p1[0]
gd=rs/rn
yint=p1[1]-((gd)*p1[0])
dty=calc_dist(p1,[0,yint])
aa=p1[1]-yint
bb=math.degrees(math.asin(aa/dty))
d=90-bb
e=180-d-thta
g=(dty/math.sin(math.radians(e)))*math.sin(math.radians(thta))
f=(g/math.sin(math.radians(thta)))*math.sin(math.radians(d))
oty=yint+g
h=f+r1
i=90-e
j=180-90-i
l=math.sin(math.radians(i))*h
k=math.cos(math.radians(i))*h
iy2=oty-l
ix2=k
int2.append(ix2)
int2.append(iy2)
m=90+bb
n=180-m-thta
p=(dty/math.sin(math.radians(n)))*math.sin(math.radians(m))
o=(p/math.sin(math.radians(m)))*math.sin(math.radians(thta))
q=p+r1
r=90-n
s=math.sin(math.radians(r))*q
t=math.cos(math.radians(r))*q
otty=yint-o
iy1=otty+s
ix1=t
int1.append(ix1)
int1.append(iy1)
elif ori==2:#-
D=calc_dist(p1,p2)
tr=r1+r2
el=tr-D
a=r1-el
b=r2-el
A=a+(el/2)
B=b+(el/2)
thta=math.degrees(math.acos(A/r1))
rs=p2[1]-p1[1]
rn=p2[0]-p1[0]
gd=rs/rn
yint=p1[1]-((gd)*p1[0])
dty=calc_dist(p1,[0,yint])
aa=yint-p1[1]
bb=math.degrees(math.asin(aa/dty))
c=180-90-bb
d=180-c-thta
e=180-90-d
f=math.tan(math.radians(e))*p1[0]
g=math.sqrt(p1[0]**2+f**2)
h=g+r1
i=180-90-e
j=math.sin(math.radians(e))*h
jj=math.cos(math.radians(i))*h
k=math.cos(math.radians(e))*h
kk=math.sin(math.radians(i))*h
l=90-bb
m=90-e
tt=l+m+thta
n=(dty/math.sin(math.radians(m)))*math.sin(math.radians(thta))
nn=(g/math.sin(math.radians(l)))*math.sin(math.radians(thta))
oty=yint-n
iy1=oty+j
ix1=k
int1.append(ix1)
int1.append(iy1)
o=bb+90
p=180-o-thta
q=90-p
r=180-90-q
s=(dty/math.sin(math.radians(p)))*math.sin(math.radians(o))
t=(s/math.sin(math.radians(o)))*math.sin(math.radians(thta))
u=s+r1
v=math.sin(math.radians(r))*u
vv=math.cos(math.radians(q))*u
w=math.cos(math.radians(r))*u
ww=math.sin(math.radians(q))*u
ix2=v
otty=yint+t
iy2=otty-w
int2.append(ix2)
int2.append(iy2)
elif ori==3:#y
D=calc_dist(p1,p2)
tr=r1+r2
el=tr-D
a=r1-el
b=r2-el
A=a+(el/2)
B=b+(el/2)
b=math.sqrt(r1**2-A**2)
int1.append(p1[0]+A)
int1.append(p1[1]+b)
int2.append(p1[0]+A)
int2.append(p1[1]-b)
elif ori==4:#x
D=calc_dist(p1,p2)
tr=r1+r2
el=tr-D
a=r1-el
b=r2-el
A=a+(el/2)
B=b+(el/2)
b=math.sqrt(r1**2-A**2)
int1.append(p1[0]+b)
int1.append(p1[1]-A)
int2.append(p1[0]-b)
int2.append(p1[1]-A)
return [int1,int2]
def calc_dist(p1,p2):
return math.sqrt((p2[0] - p1[0]) ** 2 +
(p2[1] - p1[1]) ** 2)
- kulprit.001
2
它使用笛卡尔坐标系，通过三角函数计算交点。仅从一对坐标和任意原点位置计算拦截点，然后计算角度，直到找到每个位置的位置。 - kulprit.001
这是否意味着所有以十进制度表示的经纬度都需要转换为笛卡尔坐标？而且结果将以笛卡尔坐标表示，然后需要再次转换为经纬度坐标？ - undefined
-2
确定交点P1和P2的快速方法是取圆心A和B之间的向量w。
现在想象一个由点A和P1之间的两个向量a和b所围成的矩形，我们可以说
唯一剩下的问题是，向量a和b的长度有多长：
我们知道 |a|^2 + |b|^2 = r_A^2 和 (|w| - |a|)^2 + |b|^2 = r_B^2，将它们相等并解出 |a| 的值。
现在可以使用向量的长度来通过使用归一化向量w来构建实际的向量a和b。
在实施这个想法时，解决方案非常直接。
现在想象一个由点A和P1之间的两个向量a和b所围成的矩形，我们可以说
P1 = A + a + b
P2 = A + a - b
唯一剩下的问题是，向量a和b的长度有多长：
我们知道 |a|^2 + |b|^2 = r_A^2 和 (|w| - |a|)^2 + |b|^2 = r_B^2，将它们相等并解出 |a| 的值。
|a| = (r_A^2 - r_b^2 + |w|^2) / (2|w|)
|b| = |b| = +- sqrt(r_A^2 - |a|^2)
现在可以使用向量的长度来通过使用归一化向量w来构建实际的向量a和b。
在实施这个想法时，解决方案非常直接。
w = {
x: B.x - A.x,
y: B.y - A.y
}
d = hypot(w.x, w.y)
if (d <= A.r + B.r && abs(B.r - A.r) <= d) {
w.x/= d;
w.y/= d;
a = (A.r * A.r - B.r * B.r + d * d) / (2 * d);
b = Math.sqrt(A.r * A.r - a * a);
P1 = {
x: A.x + a * w.x - b * w.y,
y: A.y + a * w.y + b * w.x
}
P2 = {
x: A.x + a * w.x + b * w.y,
y: A.y + a * w.y - b * w.x
}
} else {
P1 = P2 = null
}
- Robert Eisele
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接