将标量映射到网格,将结果导出为 vtk 文件
Mapping scalars to a mesh, exporting the result as a vtk file
我有:
- 具有标量值和定义的原点、间距和方向的 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
一定意味着我做错了什么...
事实证明我所做的一切都很好,除了三角形质心计算应该是 centroid = np_pts.mean(axis=0)
而不是 centroid = np_pts.mean(axis=1)
我有:
- 具有标量值和定义的原点、间距和方向的 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
一定意味着我做错了什么...
事实证明我所做的一切都很好,除了三角形质心计算应该是 centroid = np_pts.mean(axis=0)
而不是 centroid = np_pts.mean(axis=1)