python - 如何绘制由 python 中的向量给出的结构的表面？

Question

我想绘制由笛卡尔坐标 x、y、z 中的 3D 向量给出的数据表面。数据不能用光滑函数表示。

因此，首先我们使用函数 eq_points(N_count, r) 生成一些虚拟数据，该函数返回一个数组 points，其中包含每个点的 x、y、z 坐标我们物体的表面。量 omega 是立体角，现在不感兴趣。

#credit to Markus Deserno from MPI
#https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
def eq_points(N_count, r):
    points = []
    a = 4*np.pi*r**2/N_count
    d = np.sqrt(a)
    M_theta = int(np.pi/d)
    d_theta = np.pi/M_theta
    d_phi = a/d_theta
    for m in range(M_theta):
        theta = np.pi*(m+0.5)/M_theta
        M_phi = int(2*np.pi*np.sin(theta)/d_phi)
        for n in range(M_phi):
            phi = 2*np.pi*n/M_phi

            points.append(np.array([r*np.sin(theta)*np.cos(phi),
                                    r*np.sin(theta)*np.sin(phi),
                                    r*np.cos(theta)]))

    omega = 4*np.pi/N_count

    return np.array(points), omega

#starting plotting sequence

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

points, omega = eq_points(400, 1.)

ax.scatter(points[:,0], points[:,1], points[:,2])
ax.scatter(0., 0., 0., c="r")
ax.set_xlabel(r'$x$ axis')
ax.set_ylabel(r'$y$ axis')
ax.set_zlabel(r'$Z$ axis')

plt.savefig("./sphere.png", format="png", dpi=300)
plt.clf()

结果是一个如下图所示的球体。 蓝色点标记来自 points 数组的数据，而红色点是原点。

我想要这样的东西取自 here .然而，mplot3d 教程中的数据始终是平滑函数的结果。除了我用于球形图的 ax.scatter() 函数。

所以最后我的目标是绘制一些仅显示其表面的数据。该数据是通过改变到每个蓝点原点的径向距离产生的。此外，还需要确保每个点都与表面接触。绘制的表面如何here例如在plot_surface()中详细构造？一些实际的实时数据如下所示:

Answer 1

我建议找到船体，然后绘制单纯形(即形成船体的三角形)。确保适本地更新 x、y、z 限制。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy.spatial import ConvexHull

N = 1000
pts = np.random.randn(N, 3)

# exclude outliers
# obviously, this is data dependent
cutoff = 3.
is_outlier = np.any(np.abs(pts) > cutoff, axis=1)
pts = pts[~is_outlier]

# plot points
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(pts[:,0], pts[:,1], pts[:,2])

ax.set_xlim(-(cutoff +1), cutoff+1)
ax.set_ylim(-(cutoff +1), cutoff+1)
ax.set_zlim(-(cutoff +1), cutoff+1)

# get and plot hull
hull = ConvexHull(pts)

fig = plt.figure()
ax = Axes3D(fig)
vertices = [pts[s] for s in hull.simplices]
triangles = Poly3DCollection(vertices, edgecolor='k')
ax.add_collection3d(triangles)

ax.set_xlim(-(cutoff +1), cutoff+1)
ax.set_ylim(-(cutoff +1), cutoff+1)
ax.set_zlim(-(cutoff +1), cutoff+1)
plt.show()