# Use an external SPM Displacement Field in SimpleITK to translate a point from the source MRI to the output MRI: Annotations not matching
**URL:** https://discourse.itk.org/t/use-an-external-spm-displacement-field-in-simpleitk-to-translate-a-point-from-the-source-mri-to-the-output-mri-annotations-not-matching/7518
**Category:** Beginner Questions
**Tags:** displacement-field, simpleitk
**Created:** [May 10, 2025, 6:26am UTC](https://discourse.itk.org/t/use-an-external-spm-displacement-field-in-simpleitk-to-translate-a-point-from-the-source-mri-to-the-output-mri-annotations-not-matching/7518 "2025-05-10T06:26:44Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@miaomiaomiao](https://discourse.itk.org/u/miaomiaomiao)
#### Post date: [May 10, 2025, 6:26am UTC](https://discourse.itk.org/t/use-an-external-spm-displacement-field-in-simpleitk-to-translate-a-point-from-the-source-mri-to-the-output-mri-annotations-not-matching/7518/1 "2025-05-10T06:26:45Z")
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](https://1drv.ms/f/c/ca49f3641f1c26f4/EoCR-s9PNG5NlXCPVirjmX8BtHzKKC7TSXah7NBS4a1htA?e=auuMHC)
```auto
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?
---
### Author:  [@miaomiaomiao](https://discourse.itk.org/u/miaomiaomiao)
#### Post date: [May 13, 2025, 8:37am UTC](https://discourse.itk.org/t/use-an-external-spm-displacement-field-in-simpleitk-to-translate-a-point-from-the-source-mri-to-the-output-mri-annotations-not-matching/7518/2 "2025-05-13T08:37:28Z")
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.
```auto
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)
```
---
### Author:  [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [May 13, 2025, 1:22pm UTC](https://discourse.itk.org/t/use-an-external-spm-displacement-field-in-simpleitk-to-translate-a-point-from-the-source-mri-to-the-output-mri-annotations-not-matching/7518/3 "2025-05-13T13:22:57Z")
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?