python - 在 3D numpy 矩阵中绘制给定值的表面轮廓

Question

我有三个 3D 网格矩阵(X、Y、Z)对应于 xyz 坐标空间。

我还有一个 3D Numpy 矩阵 A，其中 A[i,j,k] 包含一个与点 (x,y ,z) 其中 x=X[i,j,k]、y=Y[i,j,k] 和 z= Z[i,j,k]。 A 中的浮点值是连续的(即 A 的相邻元素之间的值变化通常很小)。

有没有办法使用 Matplotlib 或任何其他基于 Python 的图形包绘制对应于 A 中给定浮点值的表面？例如，如果给定一个值 2.34，我有兴趣在 2.34 的任何位置获取矩阵 A 的绘制轮廓曲面(加上或减去一些tolerance)出现了吗？

到目前为止，我已经能够恢复 A 中所有值的 xyz 坐标，这些值在目标值的某个公差范围内，然后使用 this 制作 3D 散点图。 (下面的代码)。也许还有一种从这些点绘制曲面的方法？

def clean (A, t, dt):
    # function for making A binary for t+-dt
    # t is the target value I want in the matrix A with tolerance dt
    new_A = np.copy(A)
    new_A[np.logical_and(new_A > t-dt, new_A < t+dt)] = -1
    new_A[new_A != -1] = 0
    new_A[new_A == -1] = 1
    return (new_A)

def get_surface (X, Y, Z, new_A):
    x_vals = []
    y_vals = []
    z_vals = []

    # Retrieve (x,y,z) coordinates of surface
    for i in range(new_A.shape[0]):
        for j in range(new_A.shape[1]):
            for k in range(new_A.shape[2]):
                if new_A[i,j,k] == 1.0:
                    x_vals.append(X[i,j,k])
                    y_vals.append(Y[i,j,k])
                    z_vals.append(Z[i,j,k])

    return (np.array(x_vals), np.array(y_vals), np.array(z_vals))

cleaned_A = clean (A, t=2.5, dt=0.001)
x_f, y_f, z_f = get_surface (X, Y, Z, cleaned_A )

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect='equal')
ax.scatter(x_f, y_f, z_f, color='g', s=1)

我也尝试过 ax.plot_trisurf(x_f,y_f,z_f)，但这给了我连接不佳的情节。我猜我的数组中值的排序可能会影响这一点，在这种情况下，是否有一个包可以使用点的随机排序来做某种 3D 插值曲面图(例如通过最小化表面积或类似的东西) ？)

我感兴趣的物体大致是球形的(即每个 (x,y) 有两个 z)。我似乎找不到有人在封闭的 3D 表面上进行三角测量的任何工作示例，但也许我没有找对地方。

Answer 1

经过大量的挖掘，我认为我已经找到了一个可行的解决方案(至少对于一个球体——当我尝试一个球体的变形时会更新我的答案)。非常感谢帮助我思考正确道路的评论。我基本上使用 ConvexHull 从 scipy.spatial 进行三角剖分:

from matplotlib.tri import Triangulation
from scipy.spatial import ConvexHull

def clean (A, t, dt):
    # function for making A binary for t+-dt
    # t is the target value I want in the matrix A with tolerance dt
    new_A = np.copy(A)
    new_A[np.logical_and(new_A > t-dt, new_A < t+dt)] = -1
    new_A[new_A != -1] = 0
    new_A[new_A == -1] = 1
    return (new_A)

def get_surface (X, Y, Z, new_A):
    x_vals = []
    y_vals = []
    z_vals = []

    # Retrieve (x,y,z) coordinates of surface
    for i in range(new_A.shape[0]):
        for j in range(new_A.shape[1]):
            for k in range(new_A.shape[2]):
                if new_A[i,j,k] == 1.0:
                    x_vals.append(X[i,j,k])
                    y_vals.append(Y[i,j,k])
                    z_vals.append(Z[i,j,k])

    return (np.array(x_vals), np.array(y_vals), np.array(z_vals))

cleaned_A = clean (A, t=2.5, dt=0.001)
x_f, y_f, z_f = get_surface (X, Y, Z, cleaned_A )

Xs = np.vstack((x_f, y_f, z_f)).T
hull = ConvexHull(Xs)
x, y, z = Xs.T

tri = Triangulation(x, y, triangles=hull.simplices)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d', aspect='equal')
ax.plot_trisurf(tri, z, color='g', alpha=0.1)