Python/Numpy 等效于 MATLAB 等值面函数

Question

任何人都可以建议 python/numpy 中 MATLAB 的“等值面”函数的等效函数。 MATLAB 等值面 returns 面和顶点。我需要面和顶点来创建 .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()

此方法的问题是它只绘制等值面，而不是 return 面和顶点，就像 MATLAB 等值面函数那样，我需要这些面和顶点。

如有任何帮助，我们将不胜感激。

Answer 1

尽管它不在您的目标库中，PyVista 基于 VTK 构建可以帮助您轻松地做到这一点。由于您在评论中似乎接受了基于 PyVista 的解决方案，因此您可以这样做：

定义一个mesh，一般为StructuredGrid for your kind of data, although the equidistant grid in your example could even work with a 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,)

现在面孔的后勤工作有点棘手。它们都被编码在一维整数数组中。在 1d 数组中，你总是有一个整数 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）。

Python/Numpy 等效于 MATLAB 等值面函数

Python/Numpy equivalent of MATLAB isosurface functions

python

matlab

numpy

isosurface