affine transform from voxelmorph not compatible with ITK affine transform

[jmanjon]

Hi

I used your spatial_transformer within a CNN to estimate and affine transform bettwen 2 volumes.
spatial_transformer = SpatialTransformer_with_disp(add_identity=False,shift_center=True,name='transformer')

everything works well and the tranformed image looks great.

However, when i save the obtained affine transform in itk and later apply it the transformation does not work as expected. Specifically i used this code:

#save transform
affine_transform = sitk.AffineTransform(3)
output_file_path = "Affine.txt"
matrix = np.array(t[0:3, 0:3]).flatten().tolist()
translation = np.array([t[0, 3], t[1, 3], t[2, 3]]).tolist()
affine_transform.SetMatrix(matrix)
affine_transform.SetTranslation(translation)
sitk.WriteTransform(affine_transform, output_file_path)

#apply transform
input_image = sitk.ReadImage(path+os.sep+T1[i])
affine_transform = sitk.ReadTransform(output_file_path)
output_image = sitk.Resample(input_image, affine_transform)
sitk.WriteImage(output_image, respath+os.sep+"itk_"+T1[i])

can you tell me what i am missing?

thanks!!

[mu40]

Hi @jmanjon, I do not recognize SpatialTransformer_with_disp. In general, VoxelMorph with shift_center=True works in a centered voxel coordinate system, assumed identical for the fixed and moving image. In contrast, ITK operates in (LPS) world space. Before using ITK tools, you will likely need to convert your affine transform to that space.

Calling the matrix that transforms fixed to moving centered voxel coordinates (or, identically, the moving to the fixed image) trans_cen, first convert it to a transform trans_ind operating on voxel indices. For 3D images of shape shape:

# Transform converting zero-based indices to centered voxel coordinates.
ind_to_cen = np.eye(4)
ind_to_cen[:-1, -1] = -0.5 * (shape - 1)

# Transform converting centered voxel coordinates to zero-based indices.
cen_to_ind = np.eye(4)
cen_to_ind[:-1, -1] = +0.5 * (shape - 1)

trans_ind = cen_to_ind @ trans_cen @ ind_to_cen

Now you can convert the index transform trans_ind to a physical transform trans_lps using the voxel-to-world matrices fix_to_lps and mov_to_lps stored with your original images.

lps_to_fix = np.linalg.inv(fix_to_lps)
trans_lps = mov_to_lps @ trans_ind @ lps_to_fix

For reference, VoxelMorph and NiBabel return voxel-to-RAS instead of voxel-to-LPS matrices when loading images.