How to use a deformationField to deform a new volume?


Thank you for working on this package. I am using GitHub - SuperElastix/SimpleElastix: Multi-lingual medical image registration library which has now been integrated into SimpleITK.

I have generated a deformation field in SimpleItk (see I now have a nii volume with 3 dimensions. How can i then take this and apply it to a new image?
Here is my current unsuccesful attempt.

import nibabel as nib
from scipy.ndimage import map_coordinates
from skimage.transform import resize
import numpy as np
import SimpleITK as sitk
import nibabel as nib
import nrrd

# Load the image to be deformed
image_basepath = r"reoriented_data\\"
template_path = rf"{image_basepath}/average_template_10_reoriented.nii.gz"
# Load the nii.gz file
template_img = sitk.ReadImage(template_path)
template_arr = sitk.GetArrayFromImage(template_img)
# Load the deformation field
deformation_basepath = r"deformation_matrices\original\forward\\"
deformation_path = rf"{deformation_basepath}p56_to_p28_deformation_field.nii"
deformation_img = sitk.ReadImage(deformation_path)
deform_arr = sitk.GetArrayFromImage(deformation_img)
voxel_size = deformation_img.GetSpacing()
#since indices needs the number of dimensions first we transpose
deform_arr = np.transpose(deform_arr, (3, 0, 1, 2))
#the voxels are isotropic so we can use the same voxel size for all dimensions
indices = (deform_arr / voxel_size[0]) + np.indices(deform_arr.shape[1:])
temlpate_resized = resize(template_arr, deform_arr.shape[1:], order=0)
deformed_volume_data = map_coordinates(temlpate_resized, indices, order=0, prefilter=False)

The original run in elastix produces a good result, but when applying the deformation matrices to the same template it does not reproduce. Any help you can provide would be most appreciated.