Apply deformation to the 3D-CT image volume from mesh displacement

isuruwi · September 20, 2022, 11:51am

I have a reference mesh and a deformed version of it (they are unstrcuted grid saved in vtk file format). The reference mesh was acquired from a 3D-CT volume.

I would want to apply the same deformation to the reference 3D-CT volume so that it align with the deformed mesh configuration.

I just wrote a code as shown below to apply the deformation to the reference 3D-CT image volume from mesh displacement field.

import nibabel
import vtk
import numpy as np
from vtk.util import numpy_support

def getSurfaceMeshPoints(mesh_path):
    meshReader = vtk.vtkUnstructuredGridReader()
    meshReader.SetFileName(mesh_path)
    meshReader.ReadAllVectorsOn()
    meshReader.ReadAllScalarsOn()
    meshReader.Update()

    meshGeometryFilter = vtk.vtkGeometryFilter()
    meshGeometryFilter.SetInputData(meshReader.GetOutput())
    meshGeometryFilter.Update()

    polyDataForMesh = meshGeometryFilter.GetOutput()

    meshSurfaceSourcePoints = polyDataForMesh.GetPoints()

    #return meshSurfaceSourcePoints
    return meshReader.GetOutput().GetPoints()


def splineTransformationToRefMeshFollwoedByImage():
    ref_image_volume_path = "./sample_data/reference_CT.nii.gz"
    ref_volumetric_mesh_path = "./sample_data/reference_mesh.vtk"
    predicted_deformed_mesh_path = "./sample_data/deformed_mesh.vtk"
    deformed_image_volume_path = "./resampled_deformed_img.nii.gz"
    
    reference_img = nibabel.load(ref_image_volume_path)
    reference_img_data = np.asanyarray(reference_img.dataobj)
    reference_img_data_shape = reference_img_data.shape
    width, height, depth = reference_img_data.shape
    
    voxel_space_matrix = reference_img.affine[:3, :3]
    space_x, space_y, space_z = [voxel_space_matrix[i][i] for i in range(len(voxel_space_matrix))]
    
    print(space_x, space_y, space_z)

    origin_x, origin_y, origin_z = reference_img.affine[:3, 3]
    print(origin_x, origin_y, origin_z)
    
    vtk_data = numpy_support.numpy_to_vtk(num_array=reference_img_data.transpose(2, 1, 0).ravel(), deep=True, array_type=vtk.VTK_FLOAT)
    
    vtkImage = vtk.vtkImageData()
    vtkImage.SetDimensions((width, height, depth))
    vtkImage.GetPointData().SetScalars(vtk_data)
    vtkImage.SetOrigin(origin_x, origin_y, origin_z)
    vtkImage.SetSpacing(space_x, space_y, space_z)
    vtkImage.SetExtent(0, width - 1, 0, height - 1, 0, depth - 1)
    
    refMeshSurfaceSourcePoints = getSurfaceMeshPoints(ref_volumetric_mesh_path)
    deformedMeshSurfaceTargetPoints = getSurfaceMeshPoints(predicted_deformed_mesh_path)
    
    splineTransform = vtk.vtkThinPlateSplineTransform()
    splineTransform.SetSourceLandmarks(refMeshSurfaceSourcePoints)
    splineTransform.SetTargetLandmarks(deformedMeshSurfaceTargetPoints)
    splineTransform.SetBasisToR2LogR()
    
    splineTransform.Inverse()
    reslice = vtk.vtkImageReslice()
    reslice.SetInputData(vtkImage)
    reslice.SetResliceTransform(splineTransform)
    reslice.SetOutputDimensionality(3)
    reslice.SetInterpolationModeToLinear()
    
    spacing = vtkImage.GetSpacing()
    print('spacing:', spacing)
    origin = vtkImage.GetOrigin()
    print('origin:', origin)
    reslice.SetOutputSpacing(spacing[0], spacing[1], spacing[2])
    reslice.SetOutputOrigin(origin[0], origin[1], origin[2])
    
    reslice.SetOutputExtent(vtkImage.GetExtent())
    
    reslice.Update()
    
    writer = vtk.vtkNIFTIImageWriter()
    writer.SetFileName(deformed_image_volume_path)
    writer.SetInputData(reslice.GetOutput())
    writer.Write()

I would want to sanity check my method since I noticed that the reference image orintation and the deformed 3D-CT iamge oriantation is not perfectly match each other. The header inforamtion outputs are below.

Output for reference 3D-CT image:

change_ref_nib_path = './sample_data/reference_CT.nii.gz'

change_ref_nib_img = nibabel.load(change_ref_nib_path)
img_data_shape = change_ref_nib_img.shape
print(img_data_shape)
# check orientation of the 3D CT
x_nib, y_nib, z_nib = nibabel.aff2axcodes(change_ref_nib_img.affine)
print(x_nib, y_nib, z_nib)
# print(change_210_nib_img.header)
print(change_ref_nib_img.affine)

(512, 512, 105)
L A S
[[ -0.9765625 0. 0. 250. ]
[ -0. 0.9765625 0. -280.9234314]
[ 0. -0. 2. -68. ]
[ 0. 0. 0. 1. ]]

Output for deformed 3D-CT image:

resampled_deformed_img_nib_path = './resampled_deformed_img.nii.gz'

resampled_deformed_img_nib_img = nibabel.load(resampled_deformed_img_nib_path)
resampled_deformed_img_img_data_shape = resampled_deformed_img_nib_img.shape
print(resampled_deformed_img_img_data_shape)
# check orientation of the 3D CT
x_nib, y_nib, z_nib = nibabel.aff2axcodes(resampled_deformed_img_nib_img.affine)
print(x_nib, y_nib, z_nib)
# print(change_210_nib_img.header)
print(resampled_deformed_img_nib_img.affine)

(512, 512, 105)
R A S
[[ 0.9765625 0. 0. 250. ]
[ 0. 0.9765625 0. -280.9234314]
[ 0. 0. 2. -68. ]
[ 0. 0. 0. 1. ]]

Can you please help me to resolve this issue ?

I have herewith shared a sample reference mesh, deformed mesh and corresponding reference CT volume in the link given below for your further consideration.

dgobbi · September 20, 2022, 10:34pm

It might be better to read the data with vtkNIFTIImageReader.

The vtkNIFTIImageReader converts the image volume to RAS when it reads it, because some VTK filters, mappers, and writers don’t work properly with LAS. The vtkNIFTIImageWriter can convert the volume back to LAS when it saves it.

The workflow would be something like the following. First, read the image:

reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(ref_image_volume_path)
reader.Update()
imageMatrix = reader.GetSFormMatrix()
if imageMatrix is None:
    imageMatrix = reader.GetQFormMatrix()

The imageMatrix contains the offset (converted from LAS to RAS) as well as the orientation.

Next, read the meshes and transform them into the same coords as the vtkImageData:

refMeshSurfaceSourcePoints = getSurfaceMeshPoints(ref_volumetric_mesh_path)
deformedMeshSurfaceTargetPoints = getSurfaceMeshPoints(predicted_deformed_mesh_path)

pointsToImage = vtk.vtkTransform()
pointsToImage.Concatenate(imageMatrix)
pointsToImage.Inverse()

newSourcePoints = vtk.vtkPoints()
pointsToImage.TransformPoints(refMeshSurfaceSourcePoints, newSourcePoints)
newTargetPoints = vtk.vtkPoints()
pointsToImage.TransformPoints(deformedMeshSurfaceTargetPoints, newTargetPoints)

Now the newSourcePoints and newTargetPoints are in the same coordinate space as the vtkImageData, so they can be used to build the reslice transform.

Finally, you can write the image. Use SetQFac(-1) to write as LAS instead of RAS:

writer = vtk.vtkNIFTIImageWriter()
writer.SetFileName(deformed_image_volume_path)
writer.SetQFac(reader.GetQFac())
writer.SetQFormMatrix(imageMatrix)
writer.SetSFormMatrix(imageMatrix)
writer.SetInputData(reslice.GetOutput())
writer.Write()

As a sanity check, the Bounds of all the data sets should be similar:

print("input bounds", reader.GetOutput().GetBounds())
print("output bounds", reslice.GetOutput().GetBounds())
print("source bounds", newSourcePoints.GetBounds())
print("target bounds", newTargetPoints.GetBounds())

All of the above code is completely untested, but it gives an overview of how I might transform a nifti file with VTK. This kind of code is very tricky, and I usually use other packages to do things like this, for example the ANTs package has WarpImageMultiTransform for warping NIFTI files.