# Apply deformation to the 3D-CT image volume from mesh displacement
**URL:** https://discourse.vtk.org/t/apply-deformation-to-the-3d-ct-image-volume-from-mesh-displacement/9407
**Category:** Support
**Tags:** python, code
**Created:** [September 20, 2022, 11:51am UTC](https://discourse.vtk.org/t/apply-deformation-to-the-3d-ct-image-volume-from-mesh-displacement/9407 "2022-09-20T11:51:18Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@isuruwi](https://discourse.vtk.org/u/isuruwi)
#### Post date: [September 20, 2022, 11:51am UTC](https://discourse.vtk.org/t/apply-deformation-to-the-3d-ct-image-volume-from-mesh-displacement/9407/1 "2022-09-20T11:51:18Z")
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.
```auto
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:
```auto
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:
```auto
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.
> **[sample\_data.zip](https://drive.google.com/file/d/1MQuB0xBVVO10rK3Xb8mo1KpGWqyM8wcF/view?usp=sharing)**
>
> Google Drive file.
---
### Author:  [@dgobbi](https://discourse.vtk.org/u/dgobbi)
#### Post date: [September 20, 2022, 10:34pm UTC](https://discourse.vtk.org/t/apply-deformation-to-the-3d-ct-image-volume-from-mesh-displacement/9407/2 "2022-09-20T22:34:47Z")
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:
```python
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:
```python
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:
```python
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:
```nohighlight
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.