在 3D numpy 矩阵中为给定值绘制表面轮廓答案

【问题标题】：Plotting a surface contour for a given value in a 3D numpy matrix在 3D numpy 矩阵中为给定值绘制表面轮廓
【发布时间】：2019-08-13 00:11:37
【问题描述】：

我有对应xyz坐标空间的三个3D网格矩阵（X、Y、Z）。

我还有一个 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 的绘制轮廓表面？

到目前为止，我已经能够恢复 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)，但这给了我一个连接不良的情节。我猜我的数组中值的顺序可能会影响这一点，在这种情况下，是否有一个包可以对点进行随机排序（例如，通过最小化表面积或类似的东西） ?)

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

【问题讨论】：

你或多或少知道情节的形状吗？例如，有多少 z 值对应相同的 x,y 对？
我认为plot_trisurf 是正确的方法。但是，正如您所指出的，如果您选择的点没有形成一个漂亮、光滑的表面，它就会失败。作为第一遍，我将因此使用scipy.spatial.Delaunay 计算您的点云周围的外壳，然后绘制该三角剖分的点。
@xg.plt.py 该图将大致是一个可能有畸形的球体，但每对 (x,y) 总是两个 z 值。有没有办法在绘图中应用这些信息？
@PaulBrodersen 感谢您的提示。我对如何/为什么在plot_trisurf 之前实施 Delaunay 感到有些困惑（Delaunay 不是plot_trisurf 中的默认训练还是我遗漏了什么？）
Matplotlib 的plot_trisurf 默认确实使用了 Delaunay 三角剖分，而 Delaunay 三角剖分只考虑了 x 和 y 值进行三角剖分。因此，如果相同的 x,y 对有多个可能的值，它最终会绘制垂直或接近垂直的三角形。因此，结果不是图形的表面，而是与图形交叉的一堆三角形。

标签： python numpy matplotlib

【解决方案1】：

经过大量挖掘，我想我已经找到了一个可行的解决方案（至少对于一个球体——当我尝试球体的变形时会更新我的答案）。非常感谢帮助我思考正确道路的 cmets。我基本上使用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)

【讨论】：