MATLAB等曲面函数的Python/Numpy等价

Question

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

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

在python中，我发现了一个巧妙的方法，其工作原理如下：

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等深面函数那样返回面和顶点，我需要这些面和顶点。

任何帮助都将不胜感激。

Answer 1

虽然它不是您的目标库，但是构建在VTK上的PyVista可以帮助您轻松地完成这个任务(免责声明:我是开发人员之一)。由于您在评论中似乎乐于接受基于PyVista的解决方案，下面是您如何做到的：

1. 定义网格，通常为您的数据定义一个StructuredGrid，尽管您的示例中的等距网格甚至可以用于UniformGrid，
2. 用contour滤波器计算其同构面，
3. 使用包含同构面的网格的.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)获得一个2d数组，其中每一行对应于一个三角形面(第一列为常数3)。

Answer 2

在python/numpy中，MATLAB的“等表面”函数没有等效函数。你需要移动魔方从山-图像返回脸和顶点。

立方体