Manual registration with SimpleElastix

I’m using SimpleElastix (https://simpleelastix.github.io/) for the registration (Affine) of the two 2D images (see attached) . For this I’m using this code :

import SimpleITK as sitk 
elastixImageFilter = sitk.ElastixImageFilter()
elastixImageFilter.SetFixedImage(sitk.ReadImage("fixed_image.nii"))
elastixImageFilter.SetMovingImage(sitk.ReadImage("float_image.nii"))
elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap("affine"))
resultImage=elastixImageFilter.Execute()
sitk.WriteImage(resultImage,"registred_affine.nii")

After the execution of the latter, I obtain the following TransformParameters0.txt that contains the transformation matrix :

(Transform "AffineTransform")
(NumberOfParameters 6)
(TransformParameters 0.820320 0.144798 -0.144657 0.820386 -13.106613 -11.900934)
(InitialTransformParametersFileName "NoInitialTransform")
(UseBinaryFormatForTransformationParameters "false")
(HowToCombineTransforms "Compose")

// Image specific
(FixedImageDimension 2)
(MovingImageDimension 2)
(FixedInternalImagePixelType "float")
(MovingInternalImagePixelType "float")
(Size 221 257)
(Index 0 0)
(Spacing 1.0000000000 1.0000000000)
(Origin 0.0000000000 0.0000000000)
(Direction 1.0000000000 0.0000000000 0.0000000000 1.0000000000)
(UseDirectionCosines "true")

// AdvancedAffineTransform specific
(CenterOfRotationPoint 110.0000000000 128.0000000000)

// ResampleInterpolator specific
(ResampleInterpolator "FinalBSplineInterpolator")
(FinalBSplineInterpolationOrder 3)

// Resampler specific
(Resampler "DefaultResampler")
(DefaultPixelValue 0.000000)
(ResultImageFormat "nii")
(ResultImagePixelType "float")
(CompressResultImage "false")

My aim is to use this matrix-tranformation to register the floating image and get a registrered image similar to the one obtained by SimpleElastix. For this I’m using this small script :

import SimpleITK as sitk
import numpy as np

T= np.array([[0.82, 0.144, -13.1], [-0.144, 0.82, -11.9], [0, 0, 1]] ) #matrix transformation
 
img_moved_orig = plt.imread('moved.png')
img_fixed_orig = plt.imread('fixed.png')

img_transformed = np.zeros((img_moved_orig.shape[0],img_moved_orig.shape[1])) 
for i  in range(img_moved_orig.shape[0]): 
    for j in range(img_moved_orig.shape[1]): 
        pixel_data = img_moved_orig[i, j] 
        input_coords = np.array([i, j,1]) 
        i_out, j_out, _ = T @ input_coords 
        img_transformed[int(i_out), int(j_out)] = pixel_data 

I obtain this registered image which I compare with the result of SimpleElastix (see image attached). We can observe that the scaling hasn’t been operated and there is a problem with the translation. I wonder if I missed something in the transformation matrix, since SimpleElastix provide a good registration result.

Any ideas ?

Thank you

Answer

The best and safest way to apply the transformation is with the sitk.TransformixImageFilter(), but I assume you have reasons to do it a different way. With that out of the way…

First issue: you have to take into account the center of rotation. The total matrix does the following:

1. transforms the center to the origin
2. applies the matrix T you have
3. translates the result back, like this

T = np.array([[0.82, 0.144, -13.1], [-0.144, 0.82, -11.9], [0, 0, 1]] )

center = np.array([[1, 0, 110], [0, 1, 128], [0, 0, 1]] )
center_inverse = np.array([[1, 0, -110], [0, 1, -128], [0, 0, 1]] )

total_matrix = center @ T @ center_inverse

I strongly recommend using scikit-image to do the transforming for you.

from skimage.transform import AffineTransform
from skimage.transform import warp

total_affine = AffineTransform( matrix=total_matrix )
img_moving_transformed = warp( img_moved_orig, total_affine )

If you really must do the transforming yourself, there are two things that need changing in your code:

1. the axes are flipped relative to elastix expectations
2. the transformation is from fixed coordinates to moving coordinates

img_transformed = np.zeros((img_moved_orig.shape[0],img_moved_orig.shape[1])) 
for i in range(img_moved_orig.shape[0]): 
    for j in range(img_moved_orig.shape[1]): 
        
        # j is the first dimension for the elastix transform
        j_xfm, i_xfm, _ = total_matrix @ np.array([j, i, 1]) 

        pixel_data = 0
        # notice this annoying check we have to do that skimage handles for us
        if( j_xfm >= 0 and j_xfm < img_moved_orig.shape[1] and i_xfm >=0 and i_xfm < img_moved_orig.shape[0] ):
            # transformed coordinates index the moving image
            pixel_data = img_moved_orig[int(i_xfm), int(j_xfm), 0] # "nearest-neighbor" interpolation

        # "loop" indices index the output space
        img_transformed[i, j] = pixel_data