VTK调用vtkCellDataToPointData后Polydata点数据为空问题

我需要将同时包含点数据与单元数据的vtkPolyData对象全部转换为点数据格式。该vtkPolyData初始通过读取STL文件生成，后续我在其上计算了属于点数据的高斯曲率，以及属于单元数据的标量梯度。
我使用vtkCellDataToPointData滤波器完成所有数据到点数据的转换，当前实现代码如下：

import vtk

# 计算梯度
ocurv = vtk.vtkCurvatures()
ocurv.SetInputData(odata)
ocurv.SetCurvatureTypeToGaussian()
ocurv.Update()
ocurv = ocurv.GetOutput()
ocurv.GetPointData().SetActiveScalars('Gauss_Curvature')
ocurvdiff = vtk.vtkCellDerivatives()
ocurvdiff.SetInputData(ocurv)
ocurvdiff.SetVectorModeToComputeGradient()
ocurvdiff.Update()
ocurvdiff = ocurvdiff.GetOutput()

# 单元数据转点数据
grad = vtk.vtkCellDataToPointData()
grad.ProcessAllArraysOn()
grad.SetInputData(ocurvdiff)
grad.Update()
grad = grad.GetOutput()

# 提取点上的梯度数据
grad_data = vtk_to_numpy(grad.GetPointData().GetArray(2))

上述代码中odata为读取STL文件实例化得到的vtkPolyData对象。

异常现象 #

当前获取到的grad_data为空，调试输出如下：

>>> grad.GetPointData().GetArrayName(2) 
'ScalarGradient'
>>> vtk_to_numpy(grad.GetPointData().GetArray(2))
array([], shape=(0, 3), dtype=float64)

将ocurvdiff写入VTP文件时，可正常识别名为ScalarGradient的向量数据，其Mag/x/y/z分量均以单元数据形式存在；但经过单元转点数据滤波后，将grad写入VTP文件时该数据完全丢失。
已确认代码中调用了Update()触发管道更新，待转换数据为向量类型而非标量类型，仍出现该问题。

该问题是VTK中vtkCellDataToPointData滤波器的常见兼容问题，核心原因有两点：

• 部分VTK版本中ProcessAllArraysOn()开关对多分量向量数组的识别存在缺陷，不会主动转换未被标记为活动标量/向量的单元数组，仅会创建空的同名数组装入点数据
• 通过数组索引取值的方式不稳定，滤波过程中数组顺序可能发生偏移，导致取到未被正确填充的空数组

可按优先级选择以下方案修复：

方案1：修正vtkCellDataToPointData配置 #

显式指定待转换的数组，调整滤波参数，不要依赖索引取值：

grad = vtk.vtkCellDataToPointData()
grad.SetInputData(ocurvdiff)
# 显式指定要转换的单元向量数组
grad.SetInputArrayToProcess(0, 0, 0, 
    vtk.VTK_SCALAR_MODE_USE_CELL_FIELD_DATA, 
    'ScalarGradient')
grad.SetContributingCellOptionToAll() # 所有邻接单元参与点值计算
grad.ProcessAllArraysOn()
grad.Update()
grad = grad.GetOutput()
# 直接通过数组名获取数据，不要用索引
grad_data = vtk_to_numpy(grad.GetPointData().GetArray('ScalarGradient'))

方案2：绕开单元转点流程，直接输出点数据梯度 #

vtkCellDerivatives输出的是单元中心的梯度值，本身就需要插值到点上，可直接替换为vtkGradientFilter，该滤波器可直接输出点数据格式的梯度，完全规避单元转点的异常问题：

# 替换原有的vtkCellDerivatives计算逻辑
ocurv = vtk.vtkCurvatures()
ocurv.SetInputData(odata)
ocurv.SetCurvatureTypeToGaussian()
ocurv.Update()
ocurv = ocurv.GetOutput()

# 直接计算点数据上的梯度，输出结果为点数据
grad_filter = vtk.vtkGradientFilter()
grad_filter.SetInputData(ocurv)
grad_filter.SetInputArrayToProcess(0, 0, 0,
    vtk.VTK_SCALAR_MODE_USE_POINT_FIELD_DATA,
    'Gauss_Curvature')
grad_filter.SetResultArrayName('ScalarGradient')
grad_filter.Update()
grad = grad_filter.GetOutput()

# 直接提取点数据梯度
grad_data = vtk_to_numpy(grad.GetPointData().GetArray('ScalarGradient'))

该方案计算效率更高，也不会出现数组丢失的问题，是更推荐的实现方式。

内容的提问来源于stack exchange，提问作者Derek Eden