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

Question

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

因此，首先我们使用函数生成一些虚拟数据，该函数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数组中的数据，而红点是原点。

我想得到这样的东西 取自这里。但是 mplot3d 教程中的数据始终是平滑函数的结果。除了ax.scatter()我用于球体图的函数。

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

Answer 1

用新规范解决了所有点都接触表面的问题。假设用户如示例中所示设置角度，则可以通过计算由单位球体上具有相同角度的点形成的外壳的单纯形来预先计算形成构成表面的单纯形的点的索引如感兴趣的数据集中一样。然后我们可以使用这些指数来获得感兴趣的表面。

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

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

def eq_points_with_random_radius(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
            rr = r * np.random.rand()
            points.append(np.array([rr*np.sin(theta)*np.cos(phi),
                                    rr*np.sin(theta)*np.sin(phi),
                                    rr*np.cos(theta)]))

    omega = 4*np.pi/N_count

    return np.array(points), omega


N = 400
pts, _ = eq_points(N, 1.)
pts_rescaled, _ = eq_points_with_random_radius(N, 1.)
extremum = 2.

# plot points
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(pts_rescaled[:,0], pts_rescaled[:,1], pts_rescaled[:,2])
ax.set_xlim(-extremum, extremum)
ax.set_ylim(-extremum, extremum)
ax.set_zlim(-extremum, extremum)

# get indices of simplices making up the surface using points on unit sphere;
# index into rescaled points  
hull = ConvexHull(pts)
vertices = [pts_rescaled[s] for s in hull.simplices]

fig = plt.figure()
ax = Axes3D(fig)
triangles = Poly3DCollection(vertices, edgecolor='k')
ax.add_collection3d(triangles)
ax.set_xlim(-extremum, extremum)
ax.set_ylim(-extremum, extremum)
ax.set_zlim(-extremum, extremum)
plt.show()