Python：球体的交集

我对编程非常陌生，但我决定承担一个有趣的项目，因为我最近学会了如何以参数形式表示球体。当三个球体相交时，有两个不同的交点，除非它们仅在一个奇点处重叠。

球体的参数表示：

我的代码是根据答案修改的Python/matplotlib：绘制 3d 立方体、球体和向量？ https://stackoverflow.com/questions/11140163/python-matplotlib-plotting-a-3d-cube-a-sphere-and-a-vector，添加了指定 x、y 和 z 原点以及球体半径的功能。很多类似的问题都是用 C++、Java 和 C# 写的，我根本无法理解（我几乎不知道自己在做什么，所以放轻松）。

My Code:

import numpy as np

def make_sphere_x(x, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  x += radius * np.cos(u) * np.sin(v)
  return x

def make_sphere_y(y, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  y += radius * np.sin(u) * np.sin(v)
  return y

def make_sphere_z(z, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  z += radius * np.cos(v)
  return z

#x values
sphere_1_x = make_sphere_x(0, 2)
sphere_2_x = make_sphere_x(1, 3)
sphere_3_x = make_sphere_x(-1, 4)
#y values
sphere_1_y = make_sphere_y(0, 2)
sphere_2_y = make_sphere_y(1, 3)
sphere_3_y = make_sphere_y(0, 4)
#z values
sphere_1_z = make_sphere_z(0, 2)
sphere_2_z = make_sphere_z(1, 3)
sphere_3_z = make_sphere_z(-2, 4)

#intercept of x-values
intercept_x = list(filter(lambda x: x in sphere_1_x, sphere_2_x))
intercept_x = list(filter(lambda x: x in intercept_x, sphere_3_x))
print(intercept_x)

问题：

1. 显然，必须有更好的方法来找到截距。现在，代码以相等的间隔生成点，其中的间隔数是我在虚数下指定的np.mgrid。如果增加的话，交叉路口的机会应该会增加（我认为），但是当我尝试将其增加到10000j或以上，它只是吐出内存错误。

2. 数组中存在明显的间隙，即使我可以访问超级计算机并且可以将值提高到令人厌恶的值，这种方法也很可能是错误的。现在，代码结果为空集。

3. 代码效率极低，并不是说这是一个优先事项，而是人们喜欢三连的东西，对吗？

如果我在编码中犯了菜鸟错误或者在 Stack Overflow 上提问，请随意批评我。我们非常重视您的帮助。

Using scipy.optimize.fsolve如果初始猜测位于解的范围内，您可以找到给定函数的根。我用这种方法来解决你的问题，它似乎对我有用。唯一的缺点是它只为您提供一个交叉路口。要找到第二个，您必须修改初始条件，直到fsolve找到第二个根。

首先，我们通过定义每个球体的（任意）半径和中心来定义球体：

a1 = np.array([0,0,0])
r1 = .4
a2 = np.array([.3,0,0])
r2 = .5
a3 = np.array([0,.3,0])
r3 = .5

然后我们定义如何在给定角度的情况下转换回笛卡尔坐标u,v

def position(a,r,u,v):
    return a + r*np.array([np.cos(u)*np.sin(v),np.sin(u)*np.sin(v),np.cos(v)])

现在我们考虑一下我们需要求什么方程的根。对于任何交点，它都满足完美u1,v1,u2,v2,u3,v3职位position(a1,r1,u1,v1) = position(a2,r2,u2,v2) = position(a3,r3,u3,v3)是平等的。因此我们找到三个必须为零的方程，即两个位置向量的差。事实上，由于每个向量都有 3 个分量，因此我们有 9 个方程，足以确定 6 个变量。

我们找到最小化的函数为：

def f(args):
    u1,v1,u2,v2,u3,v3,_,_,_ = args
    pos1 = position(a1,r1,u1,v1)
    pos2 = position(a2,r2,u2,v2)
    pos3 = position(a3,r3,u3,v3)
    return np.array([pos1 - pos2, pos1 - pos3, pos2 - pos3]).flatten()

fsolve需要相同数量的输入和输出参数。由于我们有 9 个方程，但只有 6 个变量，我只使用了 3 个虚拟变量，以便尺寸匹配。需要展平最后一行的数组，因为fsolve仅接受一维数组。

现在可以使用 fsolve 和（相当随机的）猜测找到交集：

guess = np.array([np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,0,0,0])

x0 = fsolve(f,guess)
u1,v1,u2,v2,u3,v3,_,_,_ = x0

您可以通过将收到的角度插入到结果中来检查结果是否正确position功能。