Manual registration with SimpleElastix

Question

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

bogovicj · Accepted Answer · 2020-11-03T12:05:45.533

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:

transforms the center to the origin
applies the matrix T you have
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:

the axes are flipped relative to elastix expectations
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

Another question, when performing the registration of two 3D stack (two volume from MRI scanner), can we follow the same process you described for each slice of the two volumes ?(ie., I register slice0_vol_moved and slice0_vol_fixed and so on) — jeanluc, Nov 12 '20 at 09:54
Depends what you do. If you find a 3d transformation, then what you described (applying a 2d transform per slice) will not be correct. If instead, you find a series of 2d transforms from slice-to-slice, then it will be correct. Doing 3d transforms is more common since usually the slices in the two image are not in correspondence before registration. — bogovicj, Nov 12 '20 at 13:54

Manual registration with SimpleElastix

1 Answers1

Linked