如何用mplot3D等显示3D数组的等值面3D图？

Question

我有一个三维的numpy数组。我想在matplotlib中显示这个数组的等值面图（更准确地说，是通过在样本点之间插值来定义的三维标量场的等值面）。

matplotlib的mplot3D部分提供了很好的三维绘图支持，但（就我所见）它的接口没有简单的方法可以直接将一个三维标量值数组显示为等值面。不过，它确实支持显示一组多边形，所以我想我可以实现“行进立方体”算法来生成这些多边形。

看起来很可能已经有一个适合scipy的“行进立方体”算法被实现过，只是我还没有找到，或者我错过了一些简单的方法来做到这一点。或者，我也欢迎任何关于其他工具的建议，这些工具可以方便地用于可视化3D数组数据，并且可以在Python/numpy/scipy的环境中使用。

Answer 1

如果你想在matplotlib中保存你的图像（我觉得这比mayavi简单得多，更容易制作出适合发表的高质量图片），你可以使用skimage库中的marching_cubes函数，然后用下面的代码在matplotlib中绘制结果：

mpl_toolkits.mplot3d.art3d.Poly3DCollection

就像上面链接中展示的那样。matplotlib在渲染等值面方面表现得相当不错。这里有一个我用真实的断层扫描数据制作的例子：

Answer 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。下面是“新”的代码片段。

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()

Answer 3

我想进一步解释一下我之前的评论，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()