coregister with sitk Elastix

[romainVala]

🚀 Feature
I would like to make a PR to propose, a coregistration transform
but here I want to report an issue with Euler to affine matrix problem

Motivation
Cf #428

I open a new issue, to see if someone can help with orientation issue,

Pitch

I try to apply a given affine to a image, (with given translation and Euler Angle) and then use Elastix to coregister the 2 images and retrieve the same affine (or the same euler transform)

I put here the code because I can not find the solution. @fepegar I would appreciate any insight (I am getting sick of euler to affine transform ! ) ...

it works ok for single translation or single rotation but not when I have a mix of both

Note that in order to make it run You need to install SimpleElastix, (in to of SimpleItk (because, both package have the same name ...)

import SimpleITK as sitk
import torchio as tio
import numpy as np
np.set_printoptions(2)


def spm_matrix(P, order=1, set_ZXY=False, rotation_center=None):
    """
    FORMAT [A] = spm_matrix(P )
    P(0)  - x translation
    P(1)  - y translation
    P(2)  - z translation
    P(3)  - x rotation around x in degree
    P(4)  - y rotation around y in degree
    P(5)  - z rotation around z in degree
    P(6)  - x scaling
    P(7)  - y scaling
    P(8)  - z scaling
    P(9) - x affine
    P(10) - y affine
    P(11) - z affine

    order - application order of transformations. if order (the Default): T*R*Z*S if order==0 S*Z*R*T
    """

    #[P[3], P[4], P[5]] = [P[3] * 180 / np.pi, P[4] * 180 / np.pi, P[5] * 180 / np.pi]  # degre to radian
    P[3], P[4], P[5] = P[3]*np.pi/180, P[4]*np.pi/180, P[5]*np.pi/180 #degre to radian

    T = np.array([[1,0,0,P[0]],[0,1,0,P[1]],[0,0,1,P[2]],[0,0,0,1]])
    R1 =  np.array([[1,0,0,0],
                    [0,np.cos(P[3]),np.sin(P[3]),0],
                    [0,-np.sin(P[3]),np.cos(P[3]),0],
                    [0,0,0,1]])
    R2 =  np.array([[np.cos(P[4]),0,-np.sin(P[4]),0],  #sing change compare to spm to match tio Affine
                    [0,1,0,0],
                    [np.sin(P[4]),0,np.cos(P[4]),0],
                    [0,0,0,1]])
    R3 =  np.array([[np.cos(P[5]),np.sin(P[5]),0,0],
                    [-np.sin(P[5]),np.cos(P[5]),0,0],
                    [0,0,1,0],
                    [0,0,0,1]])

    #R = R3.dot(R2.dot(R1)) #fsl convention (with a sign difference)
    if set_ZXY:
        R = (R2.dot(R1)).dot(R3)
    else:
        R = (R1.dot(R2)).dot(R3)

    Z = np.array([[P[6],0,0,0],[0,P[7],0,0],[0,0,P[8],0],[0,0,0,1]])
    S = np.array([[1,P[9],P[10],0],[0,1,P[11],0],[0,0,1,0],[0,0,0,1]])
    if order == 0:
        A = S.dot(Z.dot(R.dot(T)))
    else:
        A = T.dot(R.dot(Z.dot(S)))

    if rotation_center is not None:
        C = np.eye(4); C[:3,3] = np.array(rotation_center)
        A = np.dot(C, np.dot(A, np.linalg.inv(C)))

    return A

def ras_to_lps_vector(triplet: np.ndarray):
    return np.array((-1, -1, 1), dtype=float) * np.asarray(triplet)

def ras_to_lps_affine(affine):
    FLIPXY_44 = np.diag([-1, -1, 1, 1])
    affine = np.dot(affine, FLIPXY_44)
    affine = np.dot(FLIPXY_44, affine)
    return affine


def euler_to_affine(Euler_angle, Translation, v_center, degree_to_radian=True, make_RAS=False,
                    set_ZYX=False):

    if degree_to_radian:
        Euler_angle = np.deg2rad(Euler_angle)
    if make_RAS: #RAS to LPS   not very usefull here since, we invert it again of the affine at the end :
        Euler_angle = ras_to_lps_vector(Euler_angle)
        Translation = ras_to_lps_vector(Translation)
        #v_center = ras_to_lps(v_center)  suppose to be given already in lps ... !

    rigid_euler = sitk.Euler3DTransform(v_center,Euler_angle[0],Euler_angle[1],Euler_angle[2],Translation)
    if set_ZYX:
        rigid_euler.SetComputeZYX(True) #default is ZXY !!! arggg !!!!

    A1 = np.asarray(rigid_euler.GetMatrix()).reshape(3,3)
    A1 = np.linalg.inv(A1)   # I do not understand why I have to invert here to get same spm convention
                            # may be passiv versus active euler definition???
    c1 = np.asarray(rigid_euler.GetCenter())
    t1 = np.asarray(rigid_euler.GetTranslation())
    affine = np.eye(4)
    affine[:3,:3] = A1
    affine[:3,3] = t1+c1 - np.matmul(A1,c1)

    """LPS to RAS"""
    if make_RAS :
        affine = ras_to_lps_affine(affine)
        #affine = np.linalg.inv(affine)  #may be not to do here it is not applyed ...

    #affine = tio.io._from_itk_convention(affine) #this include the inv
    return affine

def get_affine_from_Elastixtransfo(elastixImageFilter, do_print=True):
    tp = elastixImageFilter.GetTransformParameterMap()[0]

    v_center = np.array(np.double(tp['CenterOfRotationPoint'])) #* np.array((-1, -1, 1), dtype=float)

    Euler_angle = np.array(np.double(tp['TransformParameters'][:3]))
    Translation = np.array(np.double(tp['TransformParameters'][3:]))
    affine = euler_to_affine(Euler_angle, Translation, v_center,
                             degree_to_radian=False, make_RAS=False)


    affine = ras_to_lps_affine(affine)
    affine = np.linalg.inv(affine)  #because itk store matrix for ref to src

    if do_print:
        print(f"Found Rot  {np.rad2deg(Euler_angle)} ")
        print(f"Trans      {Translation}")
        print(f"vcenter    {v_center}")

    return affine

def ElastixRegister(img_src, img_ref):
    # inputs are tio Images
    img1 = img_src.data # if img_src.data.shape[0] == 1 else no4C
    img2 = img_ref.data # if img_ref.data.shape[0] == 1 else no4C

    i1, i2 = tio.io.nib_to_sitk(img1, img_src.affine), tio.io.nib_to_sitk(img2.data, img_ref.affine)
    elastixImageFilter = sitk.ElastixImageFilter()
    elastixImageFilter.SetFixedImage(i2);
    elastixImageFilter.SetMovingImage(i1)
    elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap("rigid"))
    elastixImageFilter.LogToConsoleOff()
    elastixImageFilter.Execute()
    return get_affine_from_Elastixtransfo(elastixImageFilter), elastixImageFilter


#tes single affine
sdata = tio.datasets.Colin27()
sdata.pop('head'); sdata.pop('brain')


#sdata.t1.save('orig.nii')
P = np.zeros(6)
for k in range(0,4):
    #P = np.zeros(6)
    P[k] = 10
    Euler_angle, Translation = P[3:] , P[:3] ;

    taff = tio.Affine([1,1,1],Euler_angle, Translation, center='image');  img_center_ras = sdata.t1.get_center()
    #taff = tio.Affine([1,1,1],Euler_angle, Translation, center='origin'); img_center_ras=[0,0,0]

    img_center_lps = ras_to_lps_vector(img_center_ras)

    affitk = euler_to_affine(Euler_angle, Translation, img_center_ras, make_RAS=False)
    #affitk = euler_to_affine(Euler_angle, Translation, img_center_lps, make_RAS=True)  #identical !

    aff = spm_matrix(np.hstack([P, [1, 1, 1, 0, 0, 0]]), order=1, set_ZXY=True, rotation_center= img_center_ras)

    if 1==1:
        srot = taff(sdata)
        sname = f'tio_A_ind{k}.nii'
        srot.t1.save(sname)

        aff_elastix, elastix_trf = ElastixRegister(sdata.t1, srot.t1)
        aff_elastix[abs(aff_elastix)<0.0001]=0
        if np.allclose(affitk, aff_elastix, atol=1e-1):
            print(f'k:{k} almost equal ')
        else:
            #print(f'k:{k} transform tioAff (then elastix) \n {affitk} \n Elastix {aff_elastix}  ')
            print(f'k:{k} tioAff - elastix) \n {affitk - aff_elastix}  ')

    if 0==1: #compare euler to affine itk and spm  -> this works in all conditions
        if np.allclose(affitk, aff, atol=1e-6):
            print(f'k:{k} spm==tioAff  ')
        else:
            print(f'k:{k} transform tioAff (then spm) \n {affitk} \n {aff}  ')
            print(f'k:{k} spm - tioAff ) \n {aff - affitk}  ')

[fepegar]

here I want to report an issue with Euler to affine matrix problem

What is the problem? Is it related to TorchIO?

because, both package have the same name

Wow! That sounds like a bad idea...

---

At the moment I'm pretty busy writing my thesis, so I'm not sure when I will be able to look into this.

[romainVala]

This issue is related to #693, in the sense that I try to fully characterize the given affine one want to set with tio.Affine transform. For instance it may be precise in the doc the exact euler convention that is chosen when the affine is define from 3 rotation and 3 translation. This is the purpose of the euler_to_affine function I propose. I realize that itk is using the following rotation composition Ry * Rx * Rz (I also had to add an inversion just for the rotation that I am not sure to correctly interpret)

Ok, as for the recent PR #712 , on could argue that it does not really matter for a use of RandomAffine transform since we want random affine no matter to get a lelft/right inversion ...

It will be really important for torchio if one want to add a coregistration transform as you suggested here #428. Indeed this must be resolved, because a simple test to setup would be to applied a given Affine (with tio.Affine) and retrieve the exact affine you have applied. (or applied a given set of rotation+translation, and retrieve this set of rotation+translation, but because the euleur representation is not invertible, we have to go through the 4*4 matrix representation (Cf euler_to_affine)

For instance in the code example here, the last tested transform is a translation of [10, 10, 10] and a rotation of [10, 0, 0] applied with center='image' (because I did not find how to change elastix convention that also report the euler parameter with rotation applied at image center). The elastix report an Euler rigid transform of

Found Rot  [ 1.00e+01 -4.12e-03  1.71e-03] 
Trans      [10.   11.59 -8.11]
vcenter    [ 0. 18. 18.]

Since it is in LPS, the translation should be read as [-10 -11.59 -8.11] and since itk is storing the reference to source transfo one should add an extra inversion to find what we set. ... [10, 10, 10 ] ... well it is almost it ... !!! argg !!!

If we do only rotation and translation it is find but because we have a "mixt" affine here, we need to go with affine representation.

In[56]: affitk
Out[56]: 
array([[ 1.  ,  0.  ,  0.  , 10.  ],
       [ 0.  ,  0.98,  0.17,  6.6 ],
       [ 0.  , -0.17,  0.98,  7.15],
       [ 0.  ,  0.  ,  0.  ,  1.  ]])
In[57]: aff_elastix
Out[57]: 
array([[ 1.00e+00, -2.93e-05,  6.67e-05,  1.00e+01],
       [ 1.73e-05,  9.85e-01,  1.74e-01,  9.42e+00],
       [-7.08e-05, -1.74e-01,  9.85e-01,  3.12e+00],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  1.00e+00]])

Note that the direct itk affine we applied with tio.Affine is not [10, 10, 10] in the translation part, because the rotation center is not at [0,0,0] ...

I understand you do not have so much time, but since you made all the code to translation from nibable to itk, may be you would notice an obvious mistake ... Although this may not be obvious (I try a few days ... without success ! )

But ok, I will consider to make a working example without torchio, and try to ask to SimpleElastix folks ...

[romainVala]

Ok, I just found a workaround:
if I change
aff_elastix, elastix_trf = ElastixRegister(sdata.t1, srot.t1)
by the other way around
aff_elastix, elastix_trf = ElastixRegister( srot.t1, sdata.t1)

it gives me the correct result !

Found Rot  [-1.00e+01 -2.39e-03 -2.17e-03] 
Trans      [ -9.99 -10.    10.  ]
vcenter    [ 0. 18. 18.]

and, then :
aff_elastix == np.linalg.inv(affitk)

but I must still miss something since looking for the inverse transformation should give the inverse affine matrix ... arg ...

by the way I found the explanation why itk is storing the transformation from reference to the moving image
https://elastix.lumc.nl/download/elastix_manual_v4.7.pdf

With the transformation defined as it is, resampling is quite simple: 
loop over all voxels x in the fixed image domain ΩF , compute its mapped position y = Tμ(x), 
interpolate the moving image at y, and fill in this value at x in the output image.

[romainVala]

Just to put my numerical output (if it can helps ... ) :

        aff_elastix, elastix_trf = ElastixRegister(sdata.t1, srot.t1)
        aff_elastix2, elastix_trf2 = ElastixRegister( srot.t1, sdata.t1)
In[78]: affitk
Out[78]: 
array([[ 1.  ,  0.  ,  0.  , 10.  ],
       [ 0.  ,  0.98,  0.17,  6.6 ],
       [ 0.  , -0.17,  0.98,  7.15],
       [ 0.  ,  0.  ,  0.  ,  1.  ]])
In[75]: aff_elastix
Out[75]: 
array([[ 1.00e+00, -2.93e-05,  6.67e-05,  1.00e+01],
       [ 1.73e-05,  9.85e-01,  1.74e-01,  9.42e+00],
       [-7.08e-05, -1.74e-01,  9.85e-01,  3.12e+00],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  1.00e+00]])
In[77]: np.linalg.inv(aff_elastix2)
Out[77]: 
array([[ 1.00e+00, -3.07e-05, -4.11e-05,  9.99e+00],
       [ 3.74e-05,  9.85e-01,  1.74e-01,  6.59e+00],
       [ 3.51e-05, -1.74e-01,  9.85e-01,  7.15e+00],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  1.00e+00]])

inv(all_elastix2) is the same as the applied affine from tio.Affine, but why it is not the same as aff_elastix ... ?