2

Given 3D MRI scans A, B, and C I want to perform an affine (co)registration of B onto A, take the transformation affine matrix of the registration and apply it on C.

My problem is that the affine matrix of the registration transform has the wrong signs. Maybe due to wrong orientation?

The TransformParameters contain 12 values of which the first 9 are the rotation matrix in row-major order and the last 3 are the translation values.

TransformParameters = [R1, R2, R3, R4, R5, R6, R7, R8, R9, Tx, Ty, Tz]

registration_affine = [[R1, R2, R3, Tx],
                       [R4, R5, R6, Ty],
                       [R7, R8, R9, Tz],
                       [0,  0,  0,  1 ]]

I know that ITK holds images in LPS orientation and nibabel in RAS. So I tried to apply a change with respect to the orientation difference to the transform_affine but this did not work out.

I can not get the same registration output as ITK, below I will show some number examples and my minimal code example.

To test this I applied an affine transformation to an existing image. The inverse of that transformation matrix is the true affine the registration could find.

array([[ 1.02800583,  0.11462834, -0.11426342, -0.43383606],
       [ 0.11462834,  1.02800583, -0.11426342,  0.47954143],
       [-0.11426342, -0.11426342,  1.02285268, -0.20457054],
       [ 0.        ,  0.        ,  0.        ,  1.        ]])

But the affine constructed as explained above, yields:

array([[ 1.02757335,  0.11459412,  0.11448339,  0.23000557],
       [ 0.11410441,  1.02746452,  0.11413955, -0.20848751],
       [ 0.11398788,  0.11411115,  1.02255042, -0.04884404],
       [ 0.        ,  0.        ,  0.        ,  1.        ]])

You can see that the values are quite close but that only the signs are wrong. In fact, if I manually set the same signs as in the "true" matrix, the transformation matrix is good.

In the ITK loader of MONAI I found code that suggested to do the following to convert an ITK affine to a nibabel affine:

np.diag([-1, -1, 1, 1]) @ registration_affine

If I use nibabels ornt_transform methods to get the ornt transform from LPS to RAS, this returns [-1, -1, 1] and matches what is done in the ITK loader of MONAI.

But applying this to the affine from above does not actually yield the correct signs (only in the translation bit):

array([[-1.02757335, -0.11459412, -0.11448339, -0.23000557],
       [-0.11410441, -1.02746452, -0.11413955,  0.20848751],
       [ 0.11398788,  0.11411115,  1.02255042, -0.04884404],
       [ 0.        ,  0.        ,  0.        ,  1.        ]])

So I am a bit stuck here.

Here a complete minimal code example to run what I'm doing / trying to do. See below also the example data and package versions.

import nibabel
import numpy as np
from monai.transforms import Affine
from nibabel import Nifti1Image
import itk

# Import Images
moving_image = itk.imread('moving_2mm.nii.gz', itk.F)
fixed_image = itk.imread('fixed_2mm.nii.gz', itk.F)

# Import Default Parameter Map
parameter_object = itk.ParameterObject.New()
affine_parameter_map = parameter_object.GetDefaultParameterMap('affine', 4)
affine_parameter_map['FinalBSplineInterpolationOrder'] = ['1']
parameter_object.AddParameterMap(affine_parameter_map)

# Call registration function
result_image, result_transform_parameters = itk.elastix_registration_method(
    fixed_image, moving_image, parameter_object=parameter_object)
parameter_map = result_transform_parameters.GetParameterMap(0)
transform_parameters = np.array(parameter_map['TransformParameters'], dtype=float)

itk.imwrite(result_image, 'reg_itk.nii.gz', compression=True)

# Convert ITK params to affine matrix
rotation = transform_parameters[:9].reshape(3, 3)
translation = transform_parameters[-3:][..., np.newaxis]
reg_affine: np.ndarray = np.append(rotation, translation, axis=1)  # type: ignore
reg_affine = np.append(reg_affine, [[0, 0, 0, 1]], axis=0)  # type: ignore

# Apply affine transform matrix via MONAI
moving_image_ni: Nifti1Image = nibabel.load('moving_2mm.nii.gz')
fixed_image_ni: Nifti1Image = nibabel.load('fixed_2mm.nii.gz')
moving_image_np: np.ndarray = moving_image_ni.get_fdata()  # type: ignore

LPS = nibabel.orientations.axcodes2ornt(('L', 'P', 'S'))
RAS = nibabel.orientations.axcodes2ornt(('R', 'A', 'S'))
ornt_transform = nibabel.orientations.ornt_transform(LPS, RAS)[:, -1]  # type: ignore

affine_transform = Affine(affine=np.diag([*ornt_transform, 1]) @ reg_affine, image_only=False)
out_img, out_affine = affine_transform(moving_image_np[np.newaxis, ...])
reg_monai = np.squeeze(out_img)

out = Nifti1Image(reg_monai, fixed_image_ni.affine, header=fixed_image_ni.header)

nibabel.save(out, 'reg_monai.nii.gz')

Input data:

Output data:

Package versions:

itk-elastix==0.12.0
monai==0.8.0
nibabel==3.1.1
numpy==1.19.2

I did ask this question before on the ITKElastix project on GitHub #145 but could not resolve my issues. Thanks to dzenanz and mstaring who tried to help over there.

Spenhouet
  • 6,556
  • 12
  • 51
  • 76

2 Answers2

2

After a lot of trying and discussing with my team we came to a realization on what is going on.

We already established how to read the ITK TransformParameters with the first 9 numbers being part of the rotation matrix read in row-major order and the last three part of the translation matrix.

rot00, rot01, rot02, rot10, rot11, rot12, rot20, rot21, rot22, tx, ty, tz = parameter_map['TransformParameters']
affine = np.array([
    [rot00, rot01, rot02, tx],
    [rot10, rot11, rot12, ty],
    [rot20, rot21, rot22, tz],
    [    0,     0,     0,  1],
], dtype=np.float32)  # yapf: disable

We already knew that nibabel has images in RAS orientation and ITK in LPS orientation.

We also already knew that if we want to change the orientation of an image, we need to flip the respective axis. LPS to RAS means flipping L->R and P->A. So a flip of the first two axis. Flip hereby is represented via -1 and no flip by 1. So a flip of the first two axis can be described by [-1, -1, 1]. We can construct an affine transformation matrix for this flip with np.diag([-1, -1, 1, 1]) (the last 1 is only for computation purposes). The affine transformation matrix to flip between LPS and RAS is therefore:

flip_LPS_RAS = np.array([[-1,  0,  0,  0],
                         [ 0, -1,  0,  0],
                         [ 0,  0,  1,  0],
                         [ 0,  0,  0,  1]])

Note that this flip works both ways. LPS -> RAS and RAS -> LPS.

If you have a 3D image and the respective affine matrix, you can flip the axis in this image by applying the 'flip_LPS_RAS'. If you want to calculate the new affine of this image, you can do:

flipped = flip_LPS_RAS @ image_affine

Since we have the foundation layed out, let's now look at what we failed to figure out.

We were aware that the affine matrix of the registration transform is based on an image in the LPS orientation and we were aware that the nibabel image is in RAS. The thought process was that we need to convert the transformation affine from LPS orientation into RAS orientation similar to the image reorientation mentioned above. So we applied the flip_LPS_RAS affine to the registration affine. What we got wrong is that this does not make the affine transformation a RAS oriented transformation.

The thing is that the registration affine expects to be applied to an image in LPS orientation and outputs an image in LPS orientation. Let's recap, that the image affine has RAS orientation and that the registration affine expects to be applied to an image in LPS orientation. Now it becomes easier to see that to apply the registration transform on an image in RAS orientation we first need to change the orientation of the image to LPS and after the registration back to RAS.

image -> flip_LPS_RAS -> registration_lps -> flip_LPS_RAS

We are only interested in the affine matrix for the registration transform, so let's ignore the image in the above chain of transformations. Writing this affine transformation chain in code:

registration_ras = flip_LPS_RAS @ registration_lps @ flip_LPS_RAS

This will yield an affine matrix which takes in a RAS oriented image, changes it to LPS orientation, performs the registration in LPS orientation and then changes the orientation back to RAS - giving us a single affine transformation which performs the ITK registration on an RAS oriented image.

From the minimal code example above, the following should now work:

affine_transform = Affine(affine=registration_ras, image_only=True)
out_img = affine_transform(moving_image_np[np.newaxis, ...])
Spenhouet
  • 6,556
  • 12
  • 51
  • 76
-1

Taking a look at this diff, you might be more interested in the old way of doing it. It directly constructs an ITK transform from 4x4 matrix.

But beware, I think there is a bug somewhere in this code. I added this recently and it decreased test accuracy, which makes me believe there is a bug somewhere in there.

Dženan
  • 3,329
  • 3
  • 31
  • 44