我有一个三维的numpy数组。我想在matplotlib中显示这个数组的等值面（或更严格地说，显示由插值样本点定义的三维标量场的等值面）。
matplotlib的mplot3D部分提供了很好的三维图形支持，但是（就我所知）它的API没有任何东西可以简单地获取一个三维标量值数组并显示等值面。然而，它支持显示一组多边形，因此我可以实现Marching Cubes算法来生成这样的多边形。
似乎已经在某个地方实现了一个友好的scipy Marching Cubes算法，而我还没有找到它，或者我错过了一些简单的方法。另外，我欢迎任何指向从Python/numpy/scipy世界轻松使用的其他3D数组数据可视化工具的指针。
我之前的评论是想进一步解释，matplotlib的3D绘图并不适用于像等值面这样复杂的绘图。它旨在为简单的三维绘图生成漂亮的、出版质量的矢量输出。它无法处理复杂的三维多边形，因此即使你自己使用Marching Cubes创建等值面，它也无法正确地呈现。
但是，您可以使用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()
from mayavi import mlab。 - Jean-Pat补充@DanHickstein的回答，您还可以使用trisurf来可视化在Marching Cubes阶段获得的多边形。
补充@DanHickstein的回答，您还可以使用trisurf来可视化在Marching Cubes阶段获得的多边形。
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()
如@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()
补充我的先前回答，我应该提到自那时以来已经发布了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()
以下是结果
如@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()
verts, faces, sumit, sumitelse = measure.marching_cubes(vol, 0, spacing=(0.1, 0.1, 0.1))，上面的代码就可以工作了。 - DrBwtsmarching_cubes已被重命名为marching_cubes_lewiner。 - Henri Menke如果您想保留matplotlib中的图形（在我看来比mayavi更容易生成出版品质量的图像），那么您可以使用在skimage中实现的marching_cubes函数，然后使用matplotlib绘制结果。
mpl_toolkits.mplot3d.art3d.Poly3DCollection
如上链接所示。Matplotlib在渲染等值面方面表现得非常好。以下是我根据一些真实层析成像数据制作的示例:
import numpy as np
from skimage.measure import marching_cubes_lewiner
import matplotlib.pyplot as plt
import s3dlib.surface as s3d
import s3dlib.cmap_utilities as cmu
bcmap = cmu.binary_cmap('orange', 'yellowgreen')
# ...................................................................
fig = plt.figure(figsize=plt.figaspect(.5))
minmax = (-2*np.pi, 2*np.pi)
for i,[dmn,title] in enumerate ( [ \
[2*np.pi, r'-2$\pi$ < surface < 2$\pi$'],
[ np.pi, r'-$\pi$ < central section < $\pi$'] ] ) :
x, y, z = dmn*np.mgrid[-1:1:61j, -1:1:61j, -1:1:61j]
vol = np.cos(x) + np.cos(y) + np.cos(z) # Schwarz P
verts, faces, _, _ = marching_cubes_lewiner(vol, 0.001, (0.1, 0.1, 0.1))
verts = dmn*( verts-(3,3,3) )/3
surface = s3d.Surface3DCollection( verts, faces, cmap=bcmap)
ax =fig.add_subplot(121+i, projection='3d')
ax.set(xlabel='X',ylabel='Y',zlabel='Z',title=title,
xlim=minmax,ylim=minmax,zlim=minmax )
surface.map_cmap_from_normals(direction=ax)
ax.add_collection3d(surface.shade(.3).hilite(0.7,focus=2))
fig.tight_layout()
plt.show()
