只是详细说明我上面的评论，matplotlib的3D绘图真的不是为了像等值面这样复杂的东西。它的目的是为非常简单的3D绘图生成漂亮的、出版质量的矢量输出。它不能处理复杂的3D多边形，所以即使你自己实现移动立方体来创建等值面，它也不能正确地渲染它。
但是，您可以使用mayavi（它的mlab API比直接使用mayavi更方便），它使用VTK来处理和可视化多维数据。
作为一个快速的例子（修改自mayavi gallery示例之一）：

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

为了补充@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()

结果如下

更新：二〇二〇年二月二十四日

正如@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()

如果您想将绘图保留在matplotlib中（在我看来，生成出版质量的图像比mayavi要容易得多），那么您可以使用skimage中实现的marking_cubes函数，然后在matplotlib中使用

mpl_toolkits.mplot3d.art3d.Poly3DCollection

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