python - MATLAB 等值面函数的 Python/Numpy 等价物

Question

任何人都可以在 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 等值面函数那样返回面和顶点，我需要这些面和顶点。

任何帮助将不胜感激。

score 4 · Accepted Answer

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

定义一个网格，通常作为StructuredGrid您的数据类型，尽管您示例中的等距网格甚至可以与 a 一起使用UniformGrid，
contour使用过滤器计算其等值面，
使用包含等值面的网格的方法另存为.stl文件。save

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

python - MATLAB 等值面函数的 Python/Numpy 等价物

1 回答 1

Related

Reference