Get Affine Transformation Matrix from TransformParameters

[Spenhouet]

I'm using ITK to perform a affine coregistration on 3D MRI scans.
Now similar to what is suggested in https://github.com/InsightSoftwareConsortium/ITKElastix/blob/master/examples/ITK_Example08_SimpleTransformix.ipynb I would like to perform the affine transformation on another image but we are using another library for that.

So I would need a full affine matrix for the transformation that describes the affine coregistration.

I noticed that the registration returns the transform parameters:

result_image, result_transform_parameters = itk.elastix_registration_method( ...
parameter_map = transform_parameters.GetParameterMap(0)
transform_parameters = np.array(parameter_map['TransformParameters'], dtype=float)

From other GitHub issues I gathered that the last 3 values of the 12 values in transform_parameters is for the translation and the first 9 are for the rotation. Not sure if that is correct.

Based on this I started playing around and got close to the correct affine transformation matrix (the resulting images are similar but not quite right).

rotation = transform_parameters[:9].reshape(3, 3)
translation = transform_parameters[-3:][..., np.newaxis]
translation = -np.flip(translation)
reg_affine: np.ndarray = np.append(rotation, translation, axis=1) 
reg_affine = np.append(reg_affine, [[0,0 ,0,1]], axis=0) 

The translation might be right but the rotation is slightly off. What am I missing / what is wrong? Also why do I need to flip and negate the translation? Do I need to do something similar with the rotation matrix?

Would really appreciate some insight. The ITK methods are a blackbox to me.

---

Not sure what example data I can share. Here some example values for the above variables:

transform_parameters = [ 9.98573286e-01, -3.92533565e-03, -7.45617495e-03,  6.11828178e-04,
                        9.98081930e-01,  7.21948277e-03,  2.74329272e-03, -1.22292851e-02,
                        1.00763651e+00, -1.52612188e-01,  8.89476641e-01,  2.55258072e+00]

with the resulting affine matrix:

reg_affine = [[ 1.00763651e+00, -1.22292851e-02,  2.74329272e-03, -2.55258072e+00],
             [ 7.21948277e-03,  9.98081930e-01,  6.11828178e-04, -8.89476641e-01],
             [-7.45617495e-03, -3.92533565e-03,  9.98573286e-01, 1.52612188e-01],
             [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, 1.00000000e+00]

[dzenanz]

I think you are getting the directionality of the transform wrong. The transform from X (fixed) to Y (moving) is to be applied to points/pixels of Y, and it will yield corresponding points/pixels in X. This has confused me more than once.

So what you might want to do is construct the matrix without negating the translation, and then invert it. If you don't invert it, you will have a X<-Y transform. If you invert it, you will have Y<-X transform. Alternatively, you might invert just the rotation part (and that might boil down to a transpose), and negate the translation.

[Spenhouet]

Thanks for the quick reply and thanks for helping. I know this is not an issue with ITK but just my problem of trying to get this solved. So I appreciate any help.

I am applying the affine transformation matrix to Y (the moving image). Not seeing a mistake there.

But I'm really unsure about what itk returns as transform parameters. From my point of view it is just a list of values which somehow translates into an affine matrix. Is my general assumption about rotation being the first part and translation the last correct?

I did as you suggested:

rotation = transform_parameters[:9].reshape(3, 3)
translation = transform_parameters[-3:][..., np.newaxis]
reg_affine: np.ndarray = np.append(rotation, translation, axis=1)
reg_affine = np.append(reg_affine, [[0,0,0,1]], axis=0)
reg_affine = np.linalg.inv(reg_affine)

But this still does not result in the correct output.

There is also the CenterOfRotationPoint, do I need that? At least ITKs AffineTransform seems to want it as parameter.

I'm working on a full minimal example. Will share shortly.

[Spenhouet]

Here my minimal example code with input and output data:

import nibabel
from nibabel import Nifti1Image
import numpy as np
from monai.transforms import Affine
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 transform params to affine transform 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
reg_affine = np.linalg.inv(reg_affine)

# 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

affine_transform = Affine(affine=reg_affine, image_only=True)
reg_monai = np.squeeze(affine_transform(moving_image_np[np.newaxis, ...]))

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

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

Input data:
fixed_2mm.nii.gz
moving_2mm.nii.gz

Output data:
reg_itk.nii.gz
reg_monai.nii.gz

Package versions:

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

[dzenanz]

This is what ITK returns as parameters:
https://github.com/InsightSoftwareConsortium/ITK/blob/v5.3rc03/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx#L488-L511
Here is what elastix docs say:
https://github.com/SuperElastix/elastix/blob/last-itk5.2-support/Components/Transforms/AdvancedAffineTransform/elxAdvancedAffineTransform.h#L35-L36

[dzenanz]

Yes, fixed parameters (usually center of rotation) matter. They are usually initialized to all zeroes, but that might not be true in your case. Also, elastix transforms are not the same as ITK transforms and differences are possible.

[dzenanz]

nibabel library might be using RAS+ coordinate system by default. ITK uses LPS+. That is another possible source of discrepancy.

[mstaring]

Yes, you need to consider the center of rotation. For affine it reads:
T(x) = A ( x - c ) + (t + c)
with the first 9 parameters filling A in row-major order, the last 3 filling t, and the fixed parameters filling c.

[Spenhouet]

Thanks for the code links. The ITK and Elastic docs are rather sparse. The only thing I got from them is that my assumption about the first 9 and last 3 numbers was right.

I don't think it has something to do with orientation. The output is pretty close and not mirrored. I also tried changing the orientation of the image to LPS but this actually resulted in a mirrored output.

@mstaring Thanks for the math on the center of rotation. I only want to define the full affine for the transformation without applying it. Can I still apply the center of rotation?

Also in what space are the center of rotation coordinates? In the voxel space of fixed or moving or in the world space?

For the above example images and code, ITK returns this center of rotation:

[  69.46565247,   66.69592285, -109.80116272]

Using this as correction factor on the translation leads to very wrong translations. So applying it does not seem right.

The images have a voxel shape of 128 x 128 x 128.
From intuition the center refers to the moving image.
If I project the center into the world space using the affine of the moving image

center = np.array(parameter_map['CenterOfRotationPoint'], dtype=float)
points_ = np.insert(center, center.shape[-1], 1, axis=-1)
center_projected = (points_ @ moving_image_ni.affine.T)[:-1]

I get these coordinates:

array([ -63.5       ,  -63.5       , -283.10232544])

But again, the value range makes me doubt, that I need to calculate these in.

Maybe some interesting tidbit. The moving image was actually produced using this custom affine transformation:

array([[ 0.99500417, -0.09983342,  0.1       ,  0.5       ],
       [-0.09983342,  0.99500417,  0.1       , -0.5       ],
       [ 0.1       ,  0.1       ,  1.        ,  0.2       ],
       [ 0.        ,  0.        ,  0.        ,  1.        ]])

So the true transformation affine the registration should find is the inverse of that, which is:

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.        ]])

The code from my minimal example (including the inverse) yields this transformation affine:

array([[ 0.99538262, -0.09987735, -0.10029317, -0.25466545],
       [-0.09944811,  0.99546527, -0.09998216,  0.22553216],
       [-0.0998615 , -0.09995486,  1.00028445,  0.0509873 ],
       [ 0.        ,  0.        ,  0.        ,  1.        ]])

without the inverse:

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 as you can see, these are already close to the expected affine but both are slightly off and have partially wrong signs.
You can also see that from the value range, applying the center of rotation will lead to a totally different affine.

[Spenhouet]

I'm out of ideas.

If anyone knows how to make my minimal example code work with the provided data that would be great.

[Spenhouet]

Maybe it is an orientation issue after all.

If I take my transformation affine without the inverse, and manually switch all signs according to the "true" transform affine, then the results match the results of the ITK registration output.

Currently looking into how I can switch these signs based on the LPS vs. RAS difference directly on the transformation affine matrix.
Will write if I find something.

[Spenhouet]

Short update here. Still no luck.

Nibabel has that concept of "ornt" / orientation, from which I can get the flip values for the translation values:

LPS = nibabel.orientations.axcodes2ornt(('L', 'P', 'S'))
RAS = nibabel.orientations.axcodes2ornt(('R', 'A', 'S'))
ornt_trans = nibabel.orientations.ornt_transform(LPS, RAS)

Where this outputs:

array([[ 0., -1.],
       [ 1., -1.],
       [ 2.,  1.]])

If I would multiply the [-1, -1, 1] with the translation array, this would yield the correct signs but this does not address the rotation matrix. I don't know how to address that.

In the MONAI ITK loader I found this code snipped:

flip_diag = [[-1, 1], [-1, -1, 1], [-1, -1, 1, 1]][sr - 1]  # itk to nibabel affine
affine = np.diag(flip_diag) @ affine

where [-1, -1, 1, 1] matches again the [-1, -1, 1] with an additional 1 for the affine computation.

I already thought I found the solution but when I perform the computation

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

the signs of the rotation matrix do not match at all:

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.        ]])

[dzenanz]

ITK's center of rotation should be in physical (world) coordinates. This index to physical discussion might be useful to you.

Here is TorchIO code for transforming between 4x4 matrix and SimpleITK image. It handles LPS/RAS conversion.

Finally, you might need to pay attention to row-major vs column-major, AKA row-vectors vs column vectors layout of the matrices.

[Spenhouet]

@dzenanz I still don't get why the center of rotation should be important for me. The affine transformation does not really care where the origin or center is if I directly apply it to the data matrix. So I'm still not sure it has any relevance for what I'm trying to do. Not sure if it is clear what I'm trying to do?

With respect to row-major vs. column-major. I thought the rotation matrix bit is row-major. That is what the C++ code looked like, that is what the docs say and what was said here. So my reshaping should be correct (row-major).

Thanks for the thread on index to physical but I missed how that helps in my case.

Thanks for the link to TorchIOs code. What you linked there is basically what I'm already doing. See my code above. But it is not working.

Maybe to clarify that I'm working on the transformation affine and not the image affine.

I'm really stuck here. I feel I'm sooo close to the solution but don't quite get it. Missing something.
I already put in a lot of time trying to solve this.
Did open a stackoverflow question on this now since I really need a solution: https://stackoverflow.com/questions/71441883/how-to-get-transformation-affine-from-itk-registration

Maybe someone is willing to execute my code with my files and see what the issue is.

[dzenanz]

4x4 matrix implicitly rotates around origin (0,0,0). If the transform you are trying to convert into 4x4 matrix has non-zero rotation center, you need to take it into account. I think Marius gave the formula.

Your numbers and type of problem looks very similar to mine. No way around debugging it. I wish there was somebody to debug it for me when I hit a block. Best way: let it rest for a few days, then come back to it with fresh mind.

[Spenhouet]

I think Marius gave the formula.

The formula was to apply the transformation, but I don't want to do that. I just want to adjust the transformation matrix itself. And again, adding the center results in a completely different (unusable) affine.

No way around debugging it.

Not sure what there is to debug? I'm expecting I'm missing code I don't know about. Can't debug non existing code. That's why I'm asking in the first place because I don't know what is missing.
I'm already on this for 2 days. I'm not getting anywhere except if someone actually looks at the code I shared.
I'm not asking for someone to "debug" my code. Usually it is highly welcomed if one shares a minimal working example which can just be executed.