ANTsX / ANTsPy

A fast medical imaging analysis library in Python with algorithms for registration, segmentation, and more.
https://antspyx.readthedocs.io
Apache License 2.0
638 stars 162 forks source link

Results are flipped when using ants.apply_transforms() #422

Closed DimitriosGk closed 1 year ago

DimitriosGk commented 1 year ago

Hello,

I am using the ants.apply_transforms() to apply a transformation to some additional images after registration of 2 MRI volumes. Specifically in the code:

"""
Registration-------------------------------------------------------------------
"""
print('STEP: Registration (Initial affine)..................')    
tx_filename = os.path.join(output_dir,'init_affine_Transform.mat')
## Make initial Affine transform, to bring images on the same space------------
txfile = ants.affine_initializer(GRE_Magn_biasCorr_denoised, 
                                 MP2RAGE_INV2, 
                                 search_factor=20, 
                                 radian_fraction=0.6, 
                                 use_principal_axis=False, 
                                 local_search_iterations=20, 
                                 mask=None, 
                                 txfn=tx_filename)
tx = ants.read_transform(txfile, dimension=2)

## Apply transform to moving image---------------------------------------------
init_trasf_INV2 = tx.apply_to_image(MP2RAGE_INV2, 
                              reference=GRE_Magn_biasCorr_denoised, 
                                  interpolation='linear')

## Check Result visually-------------------------------------------------------
print('Post-initial affine comparison----------------------------------------')
ants.plot( GRE_Magn_biasCorr_denoised, init_trasf_INV2, overlay_alpha = 0.4 )

print('STEP: Registration (Non-linear)..................')    
## Non-linear registration of images-------------------------------------------
reg = ants.registration(fixed=GRE_Magn_biasCorr_denoised, 
                        moving=init_trasf_INV2, 
                        type_of_transform = 'SyN',
                        reg_iterations = [100,100,50] )

MP2RAGE_INV2_reg = reg['warpedmovout']
## Check Result visually-------------------------------------------------------
print('Final registration comparison-----------------------------------------')
ants.plot( GRE_Magn_biasCorr_denoised, MP2RAGE_INV2_reg, overlay_alpha = 0.4, nslices = 20 )

print('STEP: Applying transforms to the  wanted volumes..................')    
## Apply transform to other images (masks from segmentation)-------------------
fullTrasform = reg[ 'invtransforms'] 
fullTrasform.append(txfile)

MP2RAGE_UNI_biasCorr_denoised_reg = ants.apply_transforms( fixed = GRE_Magn_biasCorr_denoised,
                                                          moving = MP2RAGE_UNI_biasCorr_denoised, transformlist = fullTrasform, 
                                                          interpolator  = 'linear', whichtoinvert = [False,False,False])

The direct results of the registration (reg['warpedmovout']), both in ants.plot() and in the saved .nii files appear to be fine. The problem occurs on the volumes that I try to apply the transforms to. The final result appears to be flipped both in ants.plot() and in the saved .nii files (please see images attached, image 1 shows the registration result, Image 2 the problematic one).

1 Image 1

2 Image 2

In advance I need to say that the image used for registration and the ones I apply the transform to post-registration are coming from the same processing source and have no differences in headers that I can see (can be found here: https://drive.google.com/drive/folders/1q0zwpbnQjWR7l77UMq3rXsyxOdWGMTw_?usp=sharing). Any idea why this occurs?

Thanks a lot! Dimitri

ps. in the google drive files, the ref is Magn_echo_1_masked.nii.gz, the moving image is GTSPc004_MP2RAGE_INV2_masked.nii.gz and the rest are in the same space as the moving image and resulting transform is applied on them.

ntustison commented 1 year ago

To make it as simple as possible, please add everything one would need to simply cut-and-paste to run your example including reading of images and any other files whose provenance is ambiguous.

DimitriosGk commented 1 year ago

Hello, Please find the script (a bit simplified, but tested), in the same google drive folder. All images needed are there. When you download the files from google drive, you can just change the directory in the script to read them from the right place. Best, Dimitri

ntustison commented 1 year ago

I downloaded the directory "flipping_issues" and tried to run preProcessing_3T_Registration_PET2QSM.py after changing your_dir to the appropriate location.

However, I'm getting a "file does not exist" error. After a quick look, I see on line 17, you have a subdirectory QSM which is not in the downloaded directory:

GRE_Magn_QSM_dir = os.path.join(input_dir,'QSM/Magn_echo_1_masked.nii.gz')

Please check again and make the appropriate changes so I can just run it without having to go through the code to make the necessary changes.

DimitriosGk commented 1 year ago

Hello,

That's right, I forgot to remove the directory. I changed it and tested it now again, it shall run without interruption.

Dimitri

ntustison commented 1 year ago

I'm guessing that this well-known issue is the problem. I'll invite @cookpa to weigh in as he's more familiar with the details but it would appear that your QSM image is not amenable to ITK processing:

Screenshot 2023-01-16 at 4 10 00 AM
$ PrintHeader Magn_echo_1_masked.nii.gz
 Spacing [1, 1, 1]
 Origin [0, 0, 0]
 Direction 
1 0 0
0 1 0
0 0 1

 Size : 256 256  192

  Image Dimensions   : [256, 256, 192]
  Bounding Box       : {[0 0 0], [256 256 192]}
  Voxel Spacing      : [1, 1, 1]
  Intensity Range    : [0, 9.13813]
  Mean Intensity     : 0.67338
  Direction Cos Mtx. : 
1 0 0
0 1 0
0 0 1

  Voxel->RAS x-form  : 
  Image Metadata: 
    ITK_FileNotes = 6.0.3:b862cdd5
    ITK_original_direction of unsupported type N3itk6MatrixIdLj3ELj3EEE
    ITK_original_spacing of unsupported type NSt3__16vectorIdNS_9allocatorIdEEEE
    bitpix = 64
    cal_max = 0
    cal_min = 0
    datatype = 64
    descrip = 6.0.3:b862cdd5
    dim[0] = 3
    dim[1] = 256
    dim[2] = 256
    dim[3] = 192
    dim[4] = 1
    dim[5] = 1
    dim[6] = 1
    dim[7] = 1
    dim_info = 0
    intent_code = 0
    intent_p1 = 0
    intent_p2 = 0
    intent_p3 = 0
    nifti_type = 1
    pixdim[0] = 0
    pixdim[1] = 1
    pixdim[2] = 1
    pixdim[3] = 1
    pixdim[4] = 0
    pixdim[5] = 0
    pixdim[6] = 0
    pixdim[7] = 0
    qfac = 0
    qform_code = 0
    qform_code_name = NIFTI_XFORM_UNKNOWN
    qoffset_x = 0
    qoffset_y = 0
    qoffset_z = 0
    qto_xyz of unsupported type N3itk6MatrixIfLj4ELj4EEE
    quatern_b = 0
    quatern_c = 0
    quatern_d = 0
    scl_inter = 0
    scl_slope = 1
    sform_code = 0
    sform_code_name = NIFTI_XFORM_UNKNOWN
    slice_code = 0
    slice_duration = 0
    slice_end = 0
    slice_start = 0
    srow_x = 0 0 0 0
    srow_y = 0 0 0 0
    srow_z = 0 0 0 0
    toffset = 0
    vox_offset = 352
    xyzt_units = 10

Your first clue should be the incorrect anatomical orientation.

In contrast, your other images appear to be fine, e.g.,

PrintHeader GTSPc004_MP2RAGE_UNI_masked.nii.gz 
 Spacing [1, 1, 1]
 Origin [-111.924, -139.296, -130.044]
 Direction 
0.990292 0.138455 0.0123139
-0.137415 0.988486 -0.0633421
-0.0209421 0.0610351 0.997916

 Size : 176 256  256

  Image Dimensions   : [176, 256, 256]
  Bounding Box       : {[-111.924 -139.296 -130.044], [64.0762 116.704 125.956]}
  Voxel Spacing      : [1, 1, 1]
  Intensity Range    : [0, 4095]
  Mean Intensity     : 260.301
  Direction Cos Mtx. : 
0.990292 0.138455 0.0123139
-0.137415 0.988486 -0.0633421
-0.0209421 0.0610351 0.997916

  Voxel->RAS x-form  : 
  Image Metadata: 
    ITK_FileNotes = 6.0.3:b862cdd5
    ITK_original_direction of unsupported type N3itk6MatrixIdLj3ELj3EEE
    ITK_original_spacing of unsupported type NSt3__16vectorIdNS_9allocatorIdEEEE
    bitpix = 32
    cal_max = 0
    cal_min = 0
    datatype = 16
    descrip = 6.0.3:b862cdd5
    dim[0] = 3
    dim[1] = 176
    dim[2] = 256
    dim[3] = 256
    dim[4] = 1
    dim[5] = 1
    dim[6] = 1
    dim[7] = 1
    dim_info = 0
    intent_code = 0
    intent_p1 = 0
    intent_p2 = 0
    intent_p3 = 0
    nifti_type = 1
    pixdim[0] = 0
    pixdim[1] = 1
    pixdim[2] = 1
    pixdim[3] = 1
    pixdim[4] = 0.296
    pixdim[5] = 0
    pixdim[6] = 0
    pixdim[7] = 0
    qfac = 1
    qform_code = 1
    qform_code_name = NIFTI_XFORM_SCANNER_ANAT
    qoffset_x = 111.924
    qoffset_y = 139.296
    qoffset_z = -130.044
    qto_xyz of unsupported type N3itk6MatrixIfLj4ELj4EEE
    quatern_b = -0.00833833
    quatern_c = 0.0311853
    quatern_d = 0.997083
    scl_inter = 0
    scl_slope = 1
    sform_code = 2
    sform_code_name = NIFTI_XFORM_ALIGNED_ANAT
    slice_code = 0
    slice_duration = 0
    slice_end = 0
    slice_start = 0
    srow_x = -0.990292 -0.138454 -0.0123139 111.924
    srow_y = 0.137414 -0.988486 0.0633421 139.296
    srow_z = -0.0209421 0.0610351 0.997916 -130.044
    toffset = 0
    vox_offset = 352
    xyzt_units = 10
cookpa commented 1 year ago

Yes, if qform_code and sform_code are both 0, ITK falls back on Analyze 7.5 logic and attempts to guess the orientation. However, ITK will always write images with qform and sform set to the same rigid transform.

The reference image header should be fixed so that the image has a valid sform (preferred) or qform, and the anatomical labels look correct when the image is loaded with ITK-SNAP.

DimitriosGk commented 1 year ago

Hi, thanks for the effort. I understand the part of the qform/sform and I suppose it is easily fixable in the reference image. I still try to understand the involvement of the actual anatomical orientation of the image though. How does it affect the result? The registration seems correct, even with the reference image being non-anatomically correct labeled.

**ps. Do you think passing a valid sform in the image manually will not create further issues?

ntustison commented 1 year ago

I still try to understand the involvement of the actual anatomical orientation of the image though. How does it affect the result? The registration seems correct, even with the reference image being non-anatomically correct labeled.

From the standpoint of an ANTs registration, the behavior is not a bug given the input. If you want further clarification as to the behavior of these ambiguous header scenarios, I recommend you ask over at the ITK discussion forum.

cookpa commented 1 year ago

I think there's an additional complication in the pipeline given the different input spaces but can't elaborate until I get on a computer later

cookpa commented 1 year ago

I ran the script but there's an error here

tx = ants.read_transform(txfile, dimension=2)

read_transform doesn't take a dimension argument, and in any case we are dealing with three dimensional images.

DimitriosGk commented 1 year ago

Hello, According to the docs, https://buildmedia.readthedocs.org/media/pdf/antspy/latest/antspy.pdf, p.12, read_transform does take a dimension argument: 'ants.read_transform(filename, dimension=2, precision=’float’)'.

Although, checking the result, adding dimension=2, runs fine on my version of ANTs, but doesn't really result in any difference.

My main question is, how is this possible that the 'warpedmovout' result produces the 'right' orientation, while applying the same extracted transform to the other volumes is not working in the same way. At the same time, I did the following 'experiment': Read all images involved as numpy matrices, saved as .nii again using ants.from_numpy() and ants.image_write(). All of them end up with the same sform and still, when I load them again using ants.image_read() and run the script, it seems that the flip is still happening.

Would you have a robust solution or direction to propose for solving this issue?

cookpa commented 1 year ago

I think that document is outdated. I'm running 0.3.3 and the script stopped with an error there. What version of antspy do you have?

There might be a nicer way to get this but if you type

>>> ants.version.__version__

it should say.

I agree that it's strange how a flip is introduced. The workflow is slightly unusual in that you are resampling the moving image and doing registration on init_trasf_INV2.

Normally we recommend using the original input image for registration, and supplying the initial transform to the registration call with the initial_transform parameter. This way, there's one less transform to keep track of (it's combined into the output affine) and the registration moving space has the same header for all the inputs - currently you are registering init_trasf_INV2, which has the header of the preprocessed reference image.

Still, while that may be suboptimal, I'm not sure how a flip is introduced. I'm writing out the transforms to try to trace it.

cookpa commented 1 year ago

OK from looking at the input images, I think they are all flipped with respect to the reference. I loaded INV2 into ITK-SNAP

itksnap -g Magn_echo_1_masked.nii.gz -o GTSPc004_MP2RAGE_INV2_masked.nii.gz

and used the manual registration tool (Tools -> Registration) to align centers of mass, and manually rotate.

It looks flipped to me - I think this is obscured by the intial affine, which scales and shears to make the images more similar. Then the deformable SyN can make the overlap look decent.

image

The other images look worse because the inverse transforms are being applied. I think

fullTrasform = reg[ 'invtransforms'] 

should be

fullTrasform = reg[ 'fwdtransforms'] 

This won't solve the flip, but makes the output look consistent between the input contrasts

image

So my recommendation would be:

  1. Fix the headers to remove the flip (I'm assuming the reference image is flipped but can't be sure)
  2. Do SyN registration with moving = MP2RAGE_INV2 and inital_transform = txfile
  3. Apply transforms with transformlist = reg['fwdtransforms'], whichtoinvert = [False,False]

If these are intra-subject images, I would also consider what the actual transformation is likely to be. Is the true transform between the images rigid (motion only), rigid + affine, rigid + nonlinear, rigid + affine + nonlinear? That will affect your initialization strategy. The affine_initializer might be unnecessary once the headers are made consistent.

cookpa commented 1 year ago

FYI, the readthedocs link above might be from an outdated link.

https://antspyx.readthedocs.io/_/downloads/en/stable/pdf/

Note antspyx and not antspy. Sorry about this confusion, I'll see if we can resolve it.

DimitriosGk commented 1 year ago

Hello,

OK, thanks a lot for the suggestions. I will make sure to keep the updated docs and try to fix the headers issue effectively before moving on to the further registrations optimization. Dimitri