寻找三个球之间的交点

Question

寻找三个球之间的交点

math

15

我正在寻找一种算法来找到3个球体之间的公共交点。

如果没有完整的算法，详细描述数学方法也会非常有帮助。

到目前为止，这是我找到的唯一有用的资源： http://mathforum.org/library/drmath/view/63138.html

但是那里描述的方法都不够详细，无法编写算法。

我更喜欢第二篇帖子中描述的纯代数方法，但任何方法都可以。

- Adam

1

你能否确认你指的是球面而不是实体，然后再加一些内容，使其不仅仅是数学问题，而且与编程有关？ - Pete Kirkham

好的，我想要一个C++算法来完成这个任务，但是我首先需要理解它背后的数学原理。至于你问题的另一部分，是的，只涉及到球体的表面。 - Adam

只需在谷歌上搜索“三边测量”即可：http://en.wikipedia.org/wiki/Trilateration - thalm

6个回答

9

可能比构建3D圆更容易，因为主要是在线和平面上操作：

对于每一对球体，通过减去球体方程（每个都是X^2+Y^2+Z^2+aX+bY+cZ+d=0的形式），得到包含它们交点圆的平面方程。然后你将得到三个平面P12 P23 P31。

这些平面有一个公共的线L，垂直于Q平面，并通过三个球体的中心。你要寻找的两个点在这条线上。点的中间是L和Q的交点H。

要实现这个：

- 计算P12 P23 P32的方程（球体方程的差） - 计算Q的方程（解线性系统或计算叉积） - 计算这四个平面的交点H的坐标。（解线性系统） - 从其方程中获取Q的法向量U（归一化向量） - 计算H和解X之间的距离t：t^2=R1^2-HC1^2，其中（C1，R1）是第一个球体的中心和半径。 - 解为H+tU和H-tU alt text

Cabri 3D构造展示了各种平面和线L。

- Eric Bainville

三个平面不一定有公共直线。相交的两个平面会形成一个平面或一条直线。相交的三个平面会形成一个平面、一条直线或一个点。 - ldog

1

在这种情况下，由于平面的方程不是独立的，因此交点可能是一个平面或一条直线（当三个中心共线时可能在无穷远处）。在一般情况下，交点将是一条直线。 - Eric Bainville

7

更新

这个答案的Python实现和使用示例可以在此GitHub仓库中找到。

使用此方法，解析解实际上非常好，可以告诉您何时存在解决方案，何时不存在解决方案（也可能只有一个解决方案）。没有理由使用牛顿法。

在我的测试中，我认为这比下面给出的三边测量技术更容易理解和简单。但是，两种技术都可以给出正确的答案。

原始答案

考虑两个球体的交点。为了可视化它，请考虑连接两个球体中心的3D线段N。考虑这个横截面

_{（来源：googlepages.com）}

其中红线是法线为N的平面的横截面。由于对称性，您可以从任意角度旋转此横截面，并且红线段的长度不会改变。这意味着两个球体的交点形成的结果曲线是一个圆，并且必须位于法线为N的平面中。

话虽如此，让我们开始寻找交点。首先，我们想要描述两个球体相交的结果圆。您不能使用1个方程式来描述3D中的圆，因为3D中的圆本质上是3D中的曲线，而您无法通过1个方程式来描述3D中的曲线。

考虑图片
_{（来源：googlepages.com）}

让P成为蓝色和红色线的交点。让h是从点P向上沿着红色线的线段长度。让两个球体的距离表示为d。让x是从小圆心到P的距离。然后我们必须有

x^2 +h^2 = r1^2
(d-x)^2 +h^2 = r2^2
==> h = sqrt(r1^2 - 1/d^2*(r1^2-r2^2+d^2)^2)

即你可以解决h，这是交点圆的半径。您可以从x沿着连接2个圆心的线N找到圆的中心点C。

然后，您可以将圆完全描述为（X、C、U、V都是向量）

X = C + (h * cos t) U + (h * sin t) V for t in [0,2*PI)

其中U和V是垂直的向量，位于法向量N所在的平面上。

最后一部分最容易了。只需找到此圆与最终球体的交点即可。这只需要将圆的参数方程中的x、y、z代入最后一个方程式中，并解出t即可。

编辑 ---

您将得到的方程实际上非常复杂，您将有许多正弦和余弦等于某个值。要解决这个问题，可以有两种方法：

使用等式

e^(it) = cos t + i sin t

将余弦和正弦表示为指数形式，然后将所有e^(it)项分组，您应该会得到一个e^(it)的二次方程，可以使用二次公式解出e^(it)，然后解出t。这将给您精确的解。这种方法实际上会告诉您是否存在一个解，两个解或存在一个解，具体取决于二次方法中有多少个点是实数。
使用牛顿法解出t，这种方法并不精确，但其计算方法更容易理解，并且对于这种情况非常有效。

- ldog

嗯，你确定吗？我记得在大学里，我以为这是正确的，但我的教授标记为错误了...也许我记错了，让我仔细考虑一下。 - ldog

2

只有当你将一个点定义为圆时才是正确的，不过这可能会让教授们感到不满 ;) - Alice Purcell

我之前使用的主机（Google Pages）可能已经取消了我的托管。 - ldog

Teensy注意，我认为您在h的方程中错过了一个sqrt（2）项。例如，对于r1 = r2 = d = 1的平凡情况，您的公式给出h = sqrt（1-1/1 *（1-1 + 1）^ 2）= 0。它应该是h = sqrt（r1 ^ 2-（（r1 ^ 2-r2 ^ 2 + d ^ 2）/ 2d）^ 2）。（括号地狱请原谅），上述平凡情况的答案是sqrt（3）/ 2。 - tanantish

这种方法太复杂了，建议搜索三边定位或查看我的答案... - thalm

显示剩余4条评论

5

基本上，您需要分三步完成此操作。假设您有三个球，S1、S2和S3。

C12是由S1和S2的交集创建的圆。
C23是由S2和S3的交集创建的圆。
P1、P2是C12和C13的交点。

这里唯一真正困难的部分是球体相交，幸运的是Mathworld已经解决了这个问题。事实上，Mathworld还提供了圆形相交的解决方案。

根据这些信息，您应该能够创建一个算法。

- ReaperUnreal

3

在搜索网络后，这是最先出现的结果之一，我在这里发布了我在研究几个小时后找到的最干净和最简单的解决方案：

三边测量

这个wiki网站包含了一个快速易懂的向量方法的完整描述，因此可以轻松编写代码。

- thalm

这是一个非常不错的看待问题的方式。我认为它并不比我的提议方法更简单，只是将“复杂性”转移到了代数运算中。需要指出的一件事是它大量使用了球体关于每个轴的旋转对称性，这使得这种方法无法扩展到任何没有此属性的曲线或表面。不过还是很好的发现。 - ldog

1

这里是对Eric上面发布的图片的另一种解释：

设H为三个球心所在平面。设C1、C2、C3分别为球与H相交的圆。设Lij为连接Ci和Cj两个交点的直线，则三条直线L12、L23、L13相交于一个点P。设M为通过P且垂直于H的直线，则你需要求得的两个交点位于直线M上；因此，你只需要将M与任意一个球相交即可。

- David Lehavi

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Andrew Wagner · Accepted Answer

这是我刚从维基百科文章中翻译的Python答案。不需要算法，有一个闭合形式解决方案。

import numpy                                             
from numpy import sqrt, dot, cross                       
from numpy.linalg import norm                            

# Find the intersection of three spheres                 
# P1,P2,P3 are the centers, r1,r2,r3 are the radii       
# Implementaton based on Wikipedia Trilateration article.                              
def trilaterate(P1,P2,P3,r1,r2,r3):                      
    temp1 = P2-P1                                        
    e_x = temp1/norm(temp1)                              
    temp2 = P3-P1                                        
    i = dot(e_x,temp2)                                   
    temp3 = temp2 - i*e_x                                
    e_y = temp3/norm(temp3)                              
    e_z = cross(e_x,e_y)                                 
    d = norm(P2-P1)                                      
    j = dot(e_y,temp2)                                   
    x = (r1*r1 - r2*r2 + d*d) / (2*d)                    
    y = (r1*r1 - r3*r3 -2*i*x + i*i + j*j) / (2*j)       
    temp4 = r1*r1 - x*x - y*y                            
    if temp4<0:                                          
        raise Exception("The three spheres do not intersect!");
    z = sqrt(temp4)                                      
    p_12_a = P1 + x*e_x + y*e_y + z*e_z                  
    p_12_b = P1 + x*e_x + y*e_y - z*e_z                  
    return p_12_a,p_12_b