尽管PyVista不在您的目标库中，但基于VTK构建的PyVista可以帮助您轻松完成这项工作(免责声明：我是开发人员之一)。既然你在评论中似乎接受了基于PyVista的解决方案，下面是你会怎么做：
1.为您的数据类型定义网格，通常为StructuredGrid，尽管示例中的等距网格甚至可以与UniformGrid一起使用，
1.用contour滤镜计算其等值面，
1.使用包含等值面的网格的save方法另存为.stl文件。

import numpy as np
import pyvista as pv

# generate data grid for computing the values

X, Y, Z = np.mgrid[-5:5:40j, -5:5:40j, -5:5:40j]
values = X**2 * 0.5 + Y**2 + Z**2 * 2

# create a structured grid

# (for this simple example we could've used an unstructured grid too)

# note the fortran-order call to ravel()!

mesh = pv.StructuredGrid(X, Y, Z)
mesh.point_arrays['values'] = values.ravel(order='F')  # also the active scalars

# compute 3 isosurfaces

isos = mesh.contour(isosurfaces=3, rng=[10, 40])

# or: mesh.contour(isosurfaces=np.linspace(10, 40, 3)) etc.

# plot them interactively if you want to

isos.plot(opacity=0.7)

# save to stl

isos.save('isosurfaces.stl')

交互绘图如下所示：

颜色对应于从标量数组中选取并由标量条指示的等值线。
如果我们从文件加载回网格，我们将得到结构，但不会得到标量：

loaded = pv.read('isosurfaces.stl')
loaded.plot(opacity=0.7)

缺少标量的原因是无法将数据数组导出为.stl文件：

>>> isos  # original isosurface mesh
PolyData (0x7fa7245a2220)
  N Cells:  26664
  N Points: 13656
  X Bounds: -4.470e+00, 4.470e+00
  Y Bounds: -5.000e+00, 5.000e+00
  Z Bounds: -5.000e+00, 5.000e+00
  N Arrays: 3

>>> isos.point_arrays
pyvista DataSetAttributes
Association: POINT
Contains keys:
    values
    Normals

>>> isos.cell_arrays
pyvista DataSetAttributes
Association: CELL
Contains keys:
    Normals

>>> loaded  # read back from .stl file
PolyData (0x7fa7118e7d00)
  N Cells:  26664
  N Points: 13656
  X Bounds: -4.470e+00, 4.470e+00
  Y Bounds: -5.000e+00, 5.000e+00
  Z Bounds: -5.000e+00, 5.000e+00
  N Arrays: 0

虽然原始等值面的每个都有绑定到它们的等值(提供了第一张图中看到的颜色Map)，以及点和单元法线(由于某种原因调用.save()来计算)，但在后一种情况下没有数据。
尽管如此，因为您正在寻找顶点和面，所以这样做应该很好。如果需要，还可以在PyVista一侧访问这些内容，因为等值面网格是PolyData对象：

>>> isos.n_points, isos.n_cells
(13656, 26664)

>>> isos.points.shape  # each row is a point
(13656, 3)

>>> isos.faces
array([    3,     0,    45, ..., 13529, 13531, 13530])

>>> isos.faces.shape
(106656,)

现在，这些面孔的后勤工作有点棘手。它们都以一维整数数组的形式进行编码。在一维数组中，您总是有一个整数n告诉您给定面的大小，然后n与Points数组中的点相对应的从零开始的索引。上述等值面完全由三角形组成：

>>> isos.faces[::4]  # [3 i1 i2 i3] quadruples encode faces
array([3, 3, 3, ..., 3, 3, 3])

>>> isos.is_all_triangles()
True

这就是为什么你会看到

>>> isos.faces.size == 4 * isos.n_cells
True

你可以做isos.faces.reshape(-1, 4)得到一个二维数组，其中每一行对应一个三角面(第一列是常量3)。

在PYTHON/NumPy中没有与matlab的“isosSurface”函数等价的函数。你需要来自Scikit-Image的行进立方体来返回面和顶点。
Https://scikit-image.org/docs/stable/api/skimage.measure.html?highlight=marching%20cube#skimage.measure.marching_cubes