Hi,
I’m currently working on generating a binary mask image volume from a volumetric mesh in VTK unstructured grid format. To accomplish this, I began by extracting the poly data through a geometry filter and then performed voxelization using the vtkPolyDataToImageStencil method, as illustrated in the code snippet below.
However, upon visual inspection of the resulting binary mask image volume, I observed that all voxel values appear as black. I would greatly appreciate any assistance or guidance on how to address and resolve this particular issue or am I missed any steps hrer in the code ?
ref_image_path = 'ref_image.nii.gz'
deformed_mesh_paths = glob.glob("./deformed_meshes/*.vtk")
save_image_volume_folder = "./deformed_masks"
for def_mesh_path in deformed_mesh_paths:
pred_def_mesh_name = os.path.basename(def_mesh_path).rsplit('.', 1)[0]
save_image_volume_path = os.path.join(save_image_volume_folder, pred_def_mesh_name + ".nii.gz")
imageReader = vtk.vtkNIFTIImageReader()
imageReader.SetFileName(ref_image_path)
imageReader.Update()
imageMatrix = imageReader.GetSFormMatrix()
print('imageMatrix:', imageMatrix)
if imageMatrix is None:
imageMatrix = imageReader.GetQFormMatrix()
print('imageMatrix:', imageMatrix)
meshReader = vtk.vtkUnstructuredGridReader()
meshReader.SetFileName(def_mesh_path)
meshReader.ReadAllVectorsOn()
meshReader.ReadAllScalarsOn()
meshReader.Update()
meshGeometryFilter = vtk.vtkGeometryFilter()
meshGeometryFilter.SetInputData(meshReader.GetOutput())
meshGeometryFilter.SetPassThroughCellIds(1)
meshGeometryFilter.SetPassThroughPointIds(1)
meshGeometryFilter.Update()
deformedMeshPoints = meshGeometryFilter.GetOutput().GetPoints()
pointsToImage = vtk.vtkTransform()
pointsToImage.Concatenate(imageMatrix)
pointsToImage.Inverse()
# the newdeformedMeshPoints are in the same coordinate space as the vtkImageData.
newdeformedMeshPoints = vtk.vtkPoints()
pointsToImage.TransformPoints(deformedMeshPoints, newdeformedMeshPoints)
pd = vtk.vtkPolyData()
pd.SetPoints(newdeformedMeshPoints)
origin = [0.0,0.0,0.0]
origin[0] = imageMatrix.GetElement(0, 3)
origin[1] = imageMatrix.GetElement(1, 3)
origin[2] = imageMatrix.GetElement(2, 3)
spacing = imageReader.GetOutput().GetSpacing()
dim = imageReader.GetOutput().GetDimension()
extent = imageReader.GetOutput().GetExtent()
whiteImage = vtk.vtkImageData()
whiteImage.SetSpacing(spacing)
whiteImage.SetDimensions(dim)
whiteImage.SetExtent(extent)
whiteImage.SetOrigin(origin)
whiteImage.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
# fill the image with foreground voxels:
inval = 255
outval = 0
count = whiteImage.GetNumberOfPoints()
for i in range(count):
whiteImage.GetPointData().GetScalars().SetTuple1(i, inval)
pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(pd)
pol2stenc.SetOutputOrigin(origin)
pol2stenc.SetOutputSpacing(spacing)
pol2stenc.SetOutputWholeExtent(extent)
pol2stenc.Update()
imgstenc = vtk.vtkImageStencil()
imgstenc.SetInputData(whiteImage)
imgstenc.SetStencilConnection(pol2stenc.GetOutputPort())
imgstenc.ReverseStencilOff()
imgstenc.SetBackgroundValue(outval)
imgstenc.Update()
writer = vtk.vtkNIFTIImageWriter()
writer.SetFileName(save_image_volume_path)
writer.SetQFac(-1)
writer.SetQFormMatrix(imageMatrix)
writer.SetSFormMatrix(imageMatrix)
writer.SetInputData(imgstenc.GetOutput())
writer.Write()