MATLAB 等值面函数的 Python/Numpy 等价物答案

【问题标题】：Python/Numpy equivalent of MATLAB isosurface functionsMATLAB 等值面函数的 Python/Numpy 等价物
【发布时间】：2021-08-18 14:20:44
【问题描述】：

任何人都可以在 python/numpy 中建议 MATLAB 的“isosurface”函数的等效函数。 MATLAB 等值面返回面和顶点。我需要面和顶点来创建 .stl 文件。 MATLAB 等值面函数如下所示：

[f,v] = isosurface(X,Y,Z,V,isovalue)

在 python 中，我在 plotly 中找到了这种方法，其工作原理如下：

import plotly.graph_objects as go
import numpy as np

X, Y, Z = np.mgrid[-5:5:40j, -5:5:40j, -5:5:40j]

# ellipsoid
values = X * X * 0.5 + Y * Y + Z * Z * 2

fig = go.Figure(data=go.Isosurface(
    x=X.flatten(),
    y=Y.flatten(),
    z=Z.flatten(),
    value=values.flatten(),
    isomin=10,
    isomax=40,
    caps=dict(x_show=False, y_show=False)
    ))
fig.show()

这种方法的问题在于它只绘制等值面，而不像 MATLAB 等值面函数那样返回面和顶点，我需要这些面和顶点。

任何帮助将不胜感激。

【问题讨论】：

我认为 numpy 或 scipy 不容易支持这一点。如果您愿意接受更深奥的解决方案，基于 VTK 的 pyvista 可以轻松做到这一点。 1. 根据您的数据创建一个网格，2. 使用 surfaces = mesh.contour(isovalues) 对其进行轮廓化，以及 3. surfaces.save('output.stl') 应该保存整个网格，包括面和顶点。
感谢您的帮助。您能否详细说明如何创建网格？

标签： python numpy matlab isosurface

【解决方案1】：

虽然它不在您的目标库中，但基于 VTK 构建的 PyVista 可以帮助您轻松完成此操作。由于您似乎接受 cmets 中基于 PyVista 的解决方案，因此您可以这样做：

为您的数据类型定义一个网格，通常为 StructuredGrid，尽管您示例中的等距网格甚至可以与 UniformGrid 一起使用，
使用contour 过滤器计算其等值面，
使用包含等值面的网格的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

虽然每个原始等值面都有绑定到它们的等值（提供第一张图中看到的颜色映射），以及点和单元法线（由于某种原因通过调用 .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 与点数组中的点相对应的从零开始的索引。上述等值面完全由三角形组成：

>>> 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）。

【讨论】：