如何在 matplotlib mplot3D 或类似文件中显示 3D 数组等值面的 3D 图？

Posted 2023-02-21

【中文标题】如何在 matplotlib mplot3D 或类似文件中显示 3D 数组等值面的 3D 图？

【英文标题】：How to display a 3D plot of a 3D array isosurface in matplotlib mplot3D or similar?

【发布时间】：2011-08-27 04:05:08

【问题描述】：

我有一个 3 维 numpy 数组。我想（在 matplotlib 中）显示该数组的等值面的漂亮 3D 图（或者更严格地说，显示通过在采样点之间插值定义的 3D 标量场的等值面）。

matplotlib 的 mplot3D 部分提供了很好的 3D 绘图支持，但是（据我所知）它的 API 没有任何东西可以简单地采用 3D 标量值数组并显示等值面。但是，它确实支持显示一组多边形，所以我想我可以实现行进立方体算法来生成这样的多边形。

似乎很可能已经在某个地方实现了一个对 scipy 友好的行进立方体，但我还没有找到它，或者我错过了一些简单的方法。或者，我欢迎任何指向其他工具的指针，这些工具可以从 Python/numpy/scipy 世界轻松使用可视化 3D 数组数据。

【问题讨论】：

Matplotlib 的 3D 绘图真的不适合这样的事情。 （它旨在为简单的 3D 绘图生成矢量输出，而不是完整的 3D 绘图引擎。）如果需要等值面，请使用 mayavi/mlab。

【参考方案1】：

只是为了详细说明我上面的评论，matplotlib 的 3D 绘图确实不适用于像等值面这样复杂的东西。它旨在为非常简单的 3D 绘图生成漂亮的、出版质量的矢量输出。它无法处理复杂的 3D 多边形，因此即使自己实施行进立方体来创建等值面，也无法正确渲染它。

但是，您可以改为使用mayavi（mlab API 比直接使用 mayavi 更方便），它使用VTK 来处理和可视化多维数据。

作为一个简单的示例（从 mayavi 库示例之一修改）：

import numpy as np
from enthought.mayavi import mlab

x, y, z = np.ogrid[-10:10:20j, -10:10:20j, -10:10:20j]
s = np.sin(x*y*z)/(x*y*z)

src = mlab.pipeline.scalar_field(s)
mlab.pipeline.iso_surface(src, contours=[s.min()+0.1*s.ptp(), ], opacity=0.3)
mlab.pipeline.iso_surface(src, contours=[s.max()-0.1*s.ptp(), ],)

mlab.show()

【讨论】：

完美！ apt-get install mayavi2，运行你的代码... Just Works。正是我正在寻找的。多年来我一直在想我是否不应该以某种方式使用 VTK。这看起来像是从 scipy 世界进入它的好方法。天哪，这就像发现了一个全新的星球......

还有一个 mlab contour3d 函数可以使上述内容变得更加简单：github.enthought.com/mayavi/mayavi/auto/…

取消这一点，传入特定值列表似乎在最新版本中可以完美运行，无论它值多少钱。

我刚刚看到它与一些 512^3 数组一起工作得很好。有趣的是，contour3d 的峰值内存消耗似乎大大低于上面的“管道”版本（大约 2.5GB 对 8GB；幸运的是我在一个大的 64 位系统上）。还没有尝试过使用 np.array(...,dtype=np.int16) 之类的东西做任何事情（我认为 np 数组默认为双精度）。

在 ubuntu 15.04 上，我添加了一些修改如下：from mayavi import mlab

【参考方案2】：

补充@DanHickstein 的答案，您还可以使用trisurf 来可视化在行进立方体阶段获得的多边形。

import numpy as np
from numpy import sin, cos, pi
from skimage import measure
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
   
   
def fun(x, y, z):
    return cos(x) + cos(y) + cos(z)
    
x, y, z = pi*np.mgrid[-1:1:31j, -1:1:31j, -1:1:31j]
vol = fun(x, y, z)
iso_val=0.0
verts, faces = measure.marching_cubes(vol, iso_val, spacing=(0.1, 0.1, 0.1))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2],
                cmap='Spectral', lw=1)
plt.show()

更新：2018 年 5 月 11 日

正如@DrBwts 所述，现在 marching_cubes 返回4 个值。以下代码有效。

import numpy as np
from numpy import sin, cos, pi
from skimage import measure
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def fun(x, y, z):
    return cos(x) + cos(y) + cos(z)

x, y, z = pi*np.mgrid[-1:1:31j, -1:1:31j, -1:1:31j]
vol = fun(x, y, z)
iso_val=0.0
verts, faces, _, _ = measure.marching_cubes(vol, iso_val, spacing=(0.1, 0.1, 0.1))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2],
                cmap='Spectral', lw=1)
plt.show()

更新：2020 年 2 月 2 日

除了我之前的回答，我应该提一下，从那时起PyVista 已经发布，它使这个 有点轻松的任务。

按照与之前相同的示例。

from numpy import cos, pi, mgrid
import pyvista as pv

#%% Data
x, y, z = pi*mgrid[-1:1:31j, -1:1:31j, -1:1:31j]
vol = cos(x) + cos(y) + cos(z)
grid = pv.StructuredGrid(x, y, z)
grid["vol"] = vol.flatten()
contours = grid.contour([0])

#%% Visualization
pv.set_plot_theme('document')
p = pv.Plotter()
p.add_mesh(contours, scalars=contours.points[:, 2], show_scalar_bar=False)
p.show()

结果如下

更新：2020 年 2 月 24 日

正如@HenriMenke 所说，marching_cubes 已重命名为marching_cubes_lewiner。 “新”的 sn-p 如下。

import numpy as np
from numpy import cos, pi
from skimage.measure import marching_cubes_lewiner
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x, y, z = pi*np.mgrid[-1:1:31j, -1:1:31j, -1:1:31j]
vol = cos(x) + cos(y) + cos(z)
iso_val=0.0
verts, faces, _, _ = marching_cubes_lewiner(vol, iso_val, spacing=(0.1, 0.1, 0.1))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2], cmap='Spectral',
                lw=1)
plt.show()

【讨论】：

marching_cubes 现在返回 4 个值，如果您更改为 verts, faces, sumit, sumitelse = measure.marching_cubes(vol, 0, spacing=(0.1, 0.1, 0.1))，上述代码可以工作

@AndrasDeak，我会说文档非常好并且有很多示例。话虽如此，使用起来有点不同，但并不太难。

好的，我已经仔细检查了 PyVista 文档。你是对的，它看起来非常灵活和强大，即使它的 API 看起来与我习惯的有点不同，并且文档确实充满了示例和有用的交叉链接。有时间我一定会试一试的。

好像marching_cubes 已重命名为marching_cubes_lewiner。

trisurf 的唯一问题是颜色图：它会跟随 Z 轴 (verts[:,2])，它没有 facecolors，对于定量评估可能有问题即使据说 mayavi复杂一点，它可以绘制 iso_surface 但您也可以将数据导出为 XML 并处理 Paraview 或其他 VTK 查看器

【参考方案3】：

如果你想把你的图保存在 matplotlib 中（在我看来，产生出版质量的图像比 mayavi 容易得多），那么你可以使用marching_cubes function implemented in skimage，然后在 matplotlib 中使用

mpl_toolkits.mplot3d.art3d.Poly3DCollection

如上面的链接所示。 Matplotlib 在渲染等值面方面做得很好。这是我用一些真实的断层扫描数据制作的一个例子：