与 SymPy 的三球相交(三边测量)

Three spheres intersection (trilateration) with SymPy

有没有办法用 SymPy 解决三球相交(三边测量)问题? sympy.geometry 没有球体对象,因此直接方法不可行。 SymPy 能否求解非线性方程组,如 Trilateration and the Intersection of Three Spheres?

所示

是的。有不同的方法,例如:

In [21]: x, y, z = symbols('x, y, z', real=True)

In [22]: eq1 = (x-1)**2 + (y-2)**2 + (z-3)**2 - 1

In [23]: eq2 = (x-1)**2 + (y-S(5)/2)**2 + (z-3)**2 - 1

In [24]: eq3 = (x-S.Half)**2 + (y-S(5)/2)**2 + (z-3)**2 - 1

In [25]: solve([eq1, eq2, eq3], [x, y, z])
Out[25]: 
⎡⎛              √14⎞  ⎛          √14    ⎞⎤
⎢⎜3/4, 9/4, 3 - ───⎟, ⎜3/4, 9/4, ─── + 3⎟⎥
⎣⎝               4 ⎠  ⎝           4     ⎠⎦

这 post 是为了防止有人从我下面 post 编辑的内容中受益。

实际上可以更几何地解决这个问题,而不是求解二次方程组。下面的方法可以处理以中心和半径给定的三个球体或由三个二次方程表示的球体。即使将三个球体作为三个三元二次方程给出,也可以提取系数并获得三个球体的半径和中心坐标。然后,通过多次应用余弦定律,可以从几何上找到两个交点。

def norm_2(v):
    return v.dot(v)
def norm(v):
    return np.sqrt(norm_2(v))

def extract_center_and_radius(sphere_coefficients):
    '''
    sphere_equation should be 
    x^2 + y^2 + z^2 + a1*x + a2*y + a3*z + a4 = 0
    so sphere equation should be parametrized
    as a list of 4 coefficients
    [a1, a2, a3, a4]
    '''
    a = np.array(sphere_coefficients)
    center = - a[0:3] / 2
    radius = np.sqrt(center.dot(center) - a[3])
    return center, radius

def intersecton_3_spheres(sphere1, sphere2, sphere3):
    '''
    spherei = (center, radius) is a tuple or a list, 
    with first entry an np.array representing
    the coordinates of the center, 
    the second being the radius
    '''
    A, rA = sphere1
    B, rB = sphere2
    C, rC = sphere3
    AB = B-A
    AC = C-A
    l_AB = norm(AB)
    l_AC = norm(AC)
    l_BC = norm(C-B)
    x2 = (l_AB**2 + l_AC**2 - l_BC**2) / (2 * l_AB)
    y2 = np.sqrt(l_AC**2 - x2**2)
    x3 = (l_AB**2 + rA**2 - rB**2) / (2*l_AB)
    y3 = (l_AC**2 + rA**2 - rC**2) / 2
    y3 = (y3 - x2*x3) / y2
    z3 = np.sqrt(rA**2 - x3**2 - y3**2)
    n = np.cross(AB, AC)
    n = n / norm(n)
    AB = AB / l_AB
    b = np.cross(n, AB)
    return A + x3*AB + y3*b + z3*n, A + x3*AB + y3*b - z3*n

def intersecton_of_3_spheres(sphere1, sphere2, sphere3):
    '''
    the input is three lists of 4 coefficients each
    e.g. [a1, a2, a3, a4] which describe the coefficients 
    of three sphere equations 
    x^2 + y^2 + z^2 + a1*x + a2*y + a3*z + a4 = 0    
    '''
    center_radius_1 = extract_center_and_radius(sphere1)
    center_radius_2 = extract_center_and_radius(sphere2)
    center_radius_3 = extract_center_and_radius(sphere3)
    return intersecton_3_spheres(center_radius_1, center_radius_2, center_radius_3)