将标量映射到网格，将结果导出为 vtk 文件

Question

我有：

• 具有标量值和定义的原点、间距和方向的 MetaImage (.mha) 文件
• 具有三角形的 vtk 网格文件，其中原点、间距和方向匹配

我正在尝试将标量值映射到网格，我想在 paraview 中可视化结果。

现在，我的输出只有 2 个不同的标量值，而它肯定应该有更多。我知道我的输入数据没有问题，因为我有一个可以进行映射的软件，结果很好。我想把它写成一个脚本来自动化它，并开始理解 VTK 的概念。

In [1]: import vtk
from vtk.util.numpy_support import numpy_to_vtk, vtk_to_numpy

output_file_name = '/tmp/test.vtk'

input_mesh = 'vtk_file.vtk'
input_mha = 'mha_file.mha'

# Read the MHA file (scalars)
reader = vtk.vtkMetaImageReader()
reader.SetFileName(input_mha)
reader.Update()
vtk_image = reader.GetOutput()
scalar_data = vtk_image.GetPointData().GetScalars()

# Read the VTK file (scalars)
reader = vtk.vtkPolyDataReader()
reader.SetFileName(input_mesh)
reader.Update()
polydata = reader.GetOutput()

In [2]: # Getting the number of triangles
n_cells = polydata.GetNumberOfCells()
# Initializing an array with the right shape
cell_data = np.zeros(n_cells)

# Iterating through the triangles and getting the corresponding scalar value
for i in range(n_cells):
    cell = polydata.GetCell(i)
    points = cell.GetPoints()    
    np_pts = np.array([points.GetPoint(i) for i in range(points.GetNumberOfPoints())])
    centroid = np_pts.mean(axis=1)
    centroid_id = vtk_image.FindPoint(centroid)
    value = scalar_data.GetValue(centroid_id)
    cell_data[i] = value

polydata.GetCellData().SetScalars(numpy_to_vtk(cell_data))

Out[2]: 0

In [3]: writer = vtk.vtkPolyDataWriter()
writer.SetFileName(output_file_name)
writer.SetInputData(polydata)
writer.SetScalarsName('a_name')
writer.Write()

Out[3]: 1

这个1作为最终输出一定意味着我做错了什么......

Answer 1

事实证明我所做的一切都很好，除了三角形质心计算应该是centroid = np_pts.mean(axis=0)而不是centroid = np_pts.mean(axis=1)