Use an external SPM Displacement Field in SimpleITK to translate a point from the source MRI to the output MRI: Annotations not matching

Beginner Questions
,
1

Hello community,
I have made point annotations in original MRIs and used SPM CAT12 to register the MRIs. Now I want to transform the points into the registered MRIs using SimpleITK. SPM CAT12 gives me a DisplacementField for this purpose. In the example below, I want the pixel marked with 10000 to be found at the same location in the registered MRI. But this is not the case.

Here is a minimal example: I have stored the data in the cloud here: OneDrive

import torchio as tio
import SimpleITK as sitk

tio_org = tio.ScalarImage('../data/tmp/test/IXI013-HH-1212-T1.nii.gz')
# Point should be at the same location in sitk_wmi-Space ...
tio_org.data[:, tio_org.data.shape[1] // 2, tio_org.data.shape[2] // 2, tio_org.data.shape[3] // 2] = 10000
sitk_org = tio_org.as_sitk()

tio_wmi = tio.ScalarImage('../data/tmp/test/wmiIXI013-HH-1212-T1.nii')
sitk_wmi = tio_wmi.as_sitk()

point_in_org_mri = (sitk_org.GetSize()[0] // 2, sitk_org.GetSize()[1] // 2, sitk_org.GetSize()[2] // 2)
_x, _y, _z = sitk_org.TransformIndexToPhysicalPoint(point_in_org_mri)

sitk_deformation_field = sitk.ReadImage('../data/tmp/test/y_IXI013-HH-1212-T1.nii')

np_deformation_field = sitk.GetArrayFromImage(sitk_deformation_field)
np_deformation_field = np_deformation_field.squeeze()
np_deformation_field = np_deformation_field.astype(np.float64)

sitk_deformation_field_2 = sitk.GetImageFromArray(np.transpose(np_deformation_field, (1, 2, 3, 0)),  isVector=True)
sitk_deformation_field_2.SetOrigin(sitk_deformation_field.GetOrigin())
sitk_deformation_field_2.SetSpacing(sitk_deformation_field.GetSpacing())

direction_5d = sitk_deformation_field.GetDirection()
direction_3d = (direction_5d[0], direction_5d[1], direction_5d[2],
                direction_5d[5], direction_5d[6], direction_5d[7],
                direction_5d[10], direction_5d[11], direction_5d[12])
sitk_deformation_field_2.SetDirection(direction_3d)
# sitk.WriteImage(sitk_deformation_field_2, '../data/tmp/test/mmm.nii')

displ_transform = sitk.DisplacementFieldTransform(sitk_deformation_field_2)

_x, _y, _z = displ_transform.TransformPoint((_x, _y, _z))
x, y, z = sitk_wmi.TransformPhysicalPointToIndex((_x, _y, _z))
# Points are at a different location.
sitk_wmi.SetPixel((x, y, z), 10000)

tio_wmi = tio.ScalarImage.from_sitk(sitk_wmi)
tio_wmi.save('../data/tmp/test/out.nii')
tio_org.save('../data/tmp/test/in.nii')

Where am I going wrong?

2

I solved it. I can’t say where my thinking went wrong.
It works now with nibabel… It’s not very elegant, though, and I’d be very grateful for suggestions for improvement. But it works now.

import nibabel as nib
import numpy as np

# Step 1: Load the deformation field
deformation_file = '../data/tmp/test/y_IXI013-HH-1212-T1.nii'
deformation_img = nib.load(deformation_file)
deformations_data = np.squeeze(deformation_img.get_fdata())

# Step 2: Load the source MRI image (MRI A)
mri_a_file = '../data/tmp/test/IXI013-HH-1212-T1.nii.gz'
mri_a_img = nib.load(mri_a_file)
mri_a_data = mri_a_img.get_fdata()
mri_a_affine = mri_a_img.affine

# Step 3: Define the point in MRI A voxel coordinates
point_a_voxel = np.array([128, 170, 75])

# Step 4: Convert the point to mm coordinates in MRI A space
point_a_mm = nib.affines.apply_affine(mri_a_affine, point_a_voxel)

# Step 5: Find the corresponding voxel in the deformation field
def find_closest_voxel(point_mm, deformations_data, tpm_affine):
    tpm_size = deformations_data.shape[:3]
    closest_indices = None
    min_distance = np.inf

    for x in tqdm.tqdm(range(tpm_size[0])):
        for y in range(tpm_size[1]):
            for z in range(tpm_size[2]):
                mapped_mm = deformations_data[x, y, z, :]

                distance = np.linalg.norm(point_mm - mapped_mm)
                if distance < min_distance:
                    min_distance = distance
                    closest_indices = np.array([x, y, z])

    return closest_indices

tpm_file = '../data/tmp/test/wmiIXI013-HH-1212-T1.nii.gz'
tpm_img = nib.load(tpm_file)
tpm_affine = tpm_img.affine

# Find the closest voxel indices in the TPM space
closest_voxel_indices = find_closest_voxel(point_a_mm, deformations_data, tpm_affine)

print("Point in MRI A (voxel coordinates):", point_a_voxel)
print("Point in MRI A (mm coordinates):", point_a_mm)
print("Closest voxel indices in TPM:", closest_voxel_indices)
3

I did not take a close look at this, but a frequent problem is that transforming points requires the inverse of the transform that is used for resampling the image. Could this be the source of your problem/confusion?