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?