PyVista获取网格三角形索引：顶点与索引不匹配问题求解

问题描述 #

原本使用PyMCubes可便捷输出匹配的网格点与索引，但因存在奇怪变换转而使用PyVista。但提取的顶点与网格索引无法对应，独立绘制时结果异常。测试用的球体创建代码如下：

def sphere(x, y, z):
    scalar = 8 * np.ones(len(x*y*z))
    dist = (x**2 + y**2 + z**2)**(1/2)
    scalar[(dist >= 1**2) | (dist <= 0.75**2) | (z <= 0)] = 0
    return scalar

# create a uniform grid to sample the function with
n = 100
x_min, y_min, z_min = -2, -2, -2
grid = pv.ImageData(
    dimensions=(n, n, n),
    spacing=(abs(x_min) / n * 2, abs(y_min) / n * 2, abs(z_min) / n * 2),
    origin=(x_min, y_min, z_min),
)

x, y, z = grid.points.T

values = sphere(x, y, z)

fig = go.Figure([go.Scatter3d(x=x, y=y, z=z, name='inner', mode='markers', marker=dict(size=values)),])
# fig.show(renderer='browser')

mesh = grid.contour([1], values, method='marching_cubes').smooth()
dist = np.linalg.norm(mesh.points, axis=1)
mesh.plot(scalars=dist, smooth_shading=True, specular=1, opacity=0.3, cmap="plasma", show_scalar_bar=False)

triangulated = mesh.extract_surface()
triangles = triangulated.surface_indices().reshape(-1, 3)
vertices = mesh.points

提取的vertices与triangles无法对应，独立绘制显示异常。

解决方案 #

问题核心是混淆了原始网格的顶点和提取表面后的顶点，两者的索引空间完全不同。正确做法是从提取出的表面网格中同时获取顶点和对应的三角形索引：

1. 不要混用mesh.points（原始平滑后网格的顶点）和triangulated.surface_indices()（表面网格的索引），extract_surface()会生成新的精简顶点集，索引仅对应该集合。
2. 修改代码如下：

# 原有代码保持到mesh生成
mesh = grid.contour([1], values, method='marching_cubes').smooth()

# 提取表面后，从表面网格本身获取顶点和索引
triangulated = mesh.extract_surface()
# 获取表面网格的顶点
vertices = triangulated.points
# 处理faces格式：PyVista中每个面以顶点数开头（三角形为3），后续是顶点索引，需跳过第一个元素
triangles = triangulated.faces.reshape(-1, 4)[:, 1:]

关键说明 #

• PyVista的faces属性存储格式为[3, idx0, idx1, idx2, 3, idx3, idx4, idx5,...]，所以用reshape(-1,4)拆分后，取第1到3列即可得到纯三角形索引。
• surface_indices()返回的是原始网格中属于表面的顶点索引，而非表面网格自身的顶点索引，这就是导致不匹配的根本原因。
• 验证代码（以matplotlib为例）：

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(vertices[:,0], vertices[:,1], vertices[:,2], triangles=triangles, cmap='plasma')
plt.show()

内容的提问来源于stack exchange，提问作者Sterling Butters