英文:
Visualizing the coalescence of 2 spheres
问题
我正在尝试通过Matplotlib可视化两个球体的合并。
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
r = 10
Init = [0,0,0] ; Final = [5,5,5]
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi:20j]
x = Init[0] + r*np.cos(u) * np.sin(v)
y = Init[1] + r*np.sin(u) * np.sin(v)
z = Init[2] + r*np.cos(v)
ax.plot_surface(x, y, z, color='r', alpha=0.5, linewidth=0)
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi:20j]
x = Final[0] + r*np.cos(u) * np.sin(v)
y = Final[1] + r*np.sin(u) * np.sin(v)
z = Final[2] + r*np.cos(v)
ax.plot_surface(x, y, z, color='r', alpha=0.5, linewidth=0)
plt.show()
我得到了两个球体,它们像预期的那样紧密排列。但是,我想删除两个球体之间的重叠部分,以便我得到一个仅由两个球体不重叠部分组成的表面。
Matplotlib中是否有任何功能可以绘制这个?
英文:
I am trying to visualize the coalescence of 2 spheres through Matplotlib.
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
r = 10
Init = [0,0,0] ; Final = [5,5,5]
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi:20j]
x = Init[0] + r*np.cos(u) * np.sin(v)
y = Init[1] + r*np.sin(u) * np.sin(v)
z = Init[2] + r*np.cos(v)
ax.plot_surface(x, y, z, color = 'r', alpha = 0.5, linewidth=0)
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi:20j]
x = Final[0] + r*np.cos(u) * np.sin(v)
y = Final[1] + r*np.sin(u) * np.sin(v)
z = Final[2] + r*np.cos(v)
ax.plot_surface(x, y, z, color = 'r', alpha = 0.5, linewidth=0)
plt.show()
I am getting 2 spheres closely situated as expected. However, I would like to delete the overlapping (shared) portion between the two spheres, so I get a surface composed of just the non-overlapping parts of the spheres.
Is there any feature in Matplotlib that can plot this?
答案1
得分: 0
这是该过程的步骤:
- 球体位于空间中。设
d
为连接它们中心的段C1C2
的距离。 - 让我们考虑一个参考坐标系位于第一个球的中心,正 x 轴沿着连接它们中心的段的方向。现在,第一个球的中心位于原点,第二个球的中心位于
[d, 0, 0]
。 - 计算两个球之间的交点圆。
- 选择参数化方式,使得球体与它们的“极点”相交于 x 轴。这样可以轻松在新的参考坐标系中绘制直到交点圆的球冠。
- 计算段
C1C2
在原始坐标系中与新参考坐标系的 x 轴的方向。实际上,你需要两次旋转来将C1C2
对准新坐标系的 x 轴。 - 最后,将球体的点旋转回它们的原始方向。这是通过旋转矩阵来实现的。
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
# 球体的中心
c1 = np.array([-1, -2, -3])
c2 = np.array([1, 2, 3])
# 绘制中心和连接它们的段
centers = np.stack([c1, c2])
ax.scatter(centers[:, 0], centers[:, 1], centers[:, 2])
ax.plot(centers[:, 0], centers[:, 1], centers[:, 2])
r1 = 5
r2 = 3
# 两个球的交点是一个圆
# 计算这个圆的半径。
# https://mathworld.wolfram.com/Sphere-SphereIntersection.html
# 中心之间的距离
d = np.linalg.norm(c1 - c2)
# 圆的半径
a = np.sqrt((-d + r2 - r1) * (-d - r2 + r1) * (-d + r2 + r1) * (d + r2 + r1)) / (2 * d)
# 我们想要绘制被截断在交点圆上的球冠。
# 这些是参数化的极限角度
alpha = np.arcsin(a / r2)
beta = np.arcsin(a / r1)
# 参数化,使得“极点”相交于 x 轴
u, v = np.mgrid[0:2 * np.pi:30j, beta:np.pi:20j]
x1 = r1 * np.cos(v)
y1 = r1 * np.cos(u) * np.sin(v)
z1 = r1 * np.sin(u) * np.sin(v)
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi-alpha:20j]
x2 = d + r2 * np.cos(v)
y2 = r2 * np.cos(u) * np.sin(v)
z2 = r2 * np.sin(u) * np.sin(v)
# 让我们计算两个旋转角度,将连接中心的段对准 x 轴
c2new = c2 - c1
theta = np.arctan2(c2new[1], c2new[0])
phi = np.arctan2(c2new[2], np.sqrt(c2new[0]**2 + c2new[1]**2))
# 旋转矩阵。4x4 因为它们更容易应用于坐标矩阵
Ry = lambda t: np.array([
[np.cos(t), 0, np.sin(t), 0],
[0, 1, 0, 0],
[-np.sin(t), 0, np.cos(t), 0],
[0, 0, 0, 1]
])
Rz = lambda t: np.array([
[np.cos(t), np.sin(t), 0, 0],
[-np.sin(t), np.cos(t), 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]
])
# 首先,围绕 z 轴旋转 theta,然后围绕 y 轴旋转 phi
rot = Ry(phi) @ Rz(theta)
def rotate_points(x, y, z):
shape = x.shape
# 添加一列1,以便我们可以将结果为3xN的数组与旋转矩阵相乘
coord = np.stack([t.flatten() for t in [x, y, z, np.ones_like(x)]])
coord = rot.T @ coord
x, y, z, _ = [t.reshape(shape) for t in coord]
return x, y, z
x1, y1, z1 = rotate_points(x1, y1, z1)
x2, y2, z2 = rotate_points(x2, y2, z2)
# 对球体应用必要的偏移并绘制它们
ax.plot_surface(x1 + c1[0], y1 + c1[1], z1 + c1[2], color = 'r', alpha = 0.5, linewidth=0)
ax.plot_surface(x2 + c1[0], y2 + c1[1], z2 + c1[2], color = 'g', alpha = 0.5, linewidth=0)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.view_init(30, -20)
plt.show()
英文:
This is the procedure:
- the spheres are located in space. Let
d
be the distance of the segmentC1C2
connecting their centers. - let's consider a reference frame located at the center of the first sphere, with the positive x-axis directed along the segment connecting their centers. Now, the center of the first sphere is located at the origin, the center of the second sphere is located at
[d, 0, 0]
. - Compute the intersection circle between the two spheres.
- chose a parameterization such that the sphere is plotted with their "poles" intersecting the x-axis. This make it easy to draw a sphere cap up to the intersection circle (in the new reference frame).
- Compute the orientation of the segment
C1C2
in the original frame with the x-axis of the new reference frame. Essentially, you need two rotations to alignC1C2
to the x-axis of the new frame. - Finally, rotate the points of the spheres back to their original orientation. This is done with rotation matrices.
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
# centers of the spheres
c1 = np.array([-1, -2, -3])
c2 = np.array([1, 2, 3])
# draw centers and a segment connecting them
centers = np.stack([c1, c2])
ax.scatter(centers[:, 0], centers[:, 1], centers[:, 2])
ax.plot(centers[:, 0], centers[:, 1], centers[:, 2])
r1 = 5
r2 = 3
# intersection of two sphere is a circle
# compute the radius of this circle.
# https://mathworld.wolfram.com/Sphere-SphereIntersection.html
# distance between centers
d = np.linalg.norm(c1 - c2)
# radius of the circle
a = np.sqrt((-d + r2 - r1) * (-d - r2 + r1) * (-d + r2 + r1) * (d + r2 + r1)) / (2 * d)
# we want to draw sphere caps that are cut at the intersecting circle.
# these are the limiting angles for the parameterization
alpha = np.arcsin(a / r2)
beta = np.arcsin(a / r1)
# parameterization such that the "poles" intersect the x-axis
u, v = np.mgrid[0:2 * np.pi:30j, beta:np.pi:20j]
x1 = r1 * np.cos(v)
y1 = r1 * np.cos(u) * np.sin(v)
z1 = r1 * np.sin(u) * np.sin(v)
u, v = np.mgrid[0:2 * np.pi:30j, 0:np.pi-alpha:20j]
x2 = d + r2 * np.cos(v)
y2 = r2 * np.cos(u) * np.sin(v)
z2 = r2 * np.sin(u) * np.sin(v)
# let's compute the two rotation angles that aligns the
# segment connecting the centers to the x-axis
c2new = c2 - c1
theta = np.arctan2(c2new[1], c2new[0])
phi = np.arctan2(c2new[2], np.sqrt(c2new[0]**2 + c2new[1]**2))
# rotations matrices. 4x4 because they are easier to apply
# to a matrix of coordinates
Ry = lambda t: np.array([
[np.cos(t), 0, np.sin(t), 0],
[0, 1, 0, 0],
[-np.sin(t), 0, np.cos(t), 0],
[0, 0, 0, 1]
])
Rz = lambda t: np.array([
[np.cos(t), np.sin(t), 0, 0],
[-np.sin(t), np.cos(t), 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]
])
# first, rotate by theta around the z-axis,
# then rotate by phi around y-axis
rot = Ry(phi) @ Rz(theta)
def rotate_points(x, y, z):
shape = x.shape
# add a column of 1s so that we can multiply the
# resulting 3xN array with the rotation matrix
coord = np.stack([t.flatten() for t in [x, y, z, np.ones_like(x)]])
coord = rot.T @ coord
x, y, z, _ = [t.reshape(shape) for t in coord]
return x, y, z
x1, y1, z1 = rotate_points(x1, y1, z1)
x2, y2, z2 = rotate_points(x2, y2, z2)
# apply the necessary offset to the spheres and plot them
ax.plot_surface(x1 + c1[0], y1 + c1[1], z1 + c1[2], color = 'r', alpha = 0.5, linewidth=0)
ax.plot_surface(x2 + c1[0], y2 + c1[1], z2 + c1[2], color = 'g', alpha = 0.5, linewidth=0)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.view_init(30, -20)
plt.show()
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论