Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

to_nibabel() produces wrong transforms #52

Closed
lmoesch opened this issue Mar 25, 2019 · 26 comments
Closed

to_nibabel() produces wrong transforms #52

lmoesch opened this issue Mar 25, 2019 · 26 comments

Comments

@lmoesch
Copy link

lmoesch commented Mar 25, 2019

While using transforms on the MNI152 template I noticed that the output of to_nibabel() was flipped along the vertical plane. While investigating the issue, I noticed, that this is caused by a wrong extraction of the q-form parameters (direction, origin) from ITK.

    ants_mni = ants.image_read('./MNI152.nii.gz')
    nii_mni = nib.load('./MNI152.nii.gz')

    ants_mni = ants_mni.to_nibabel()

    print(ants_mni.get_qform())
    print(nii_mni.get_qform())

results in:

[[  2.   0.   0. -90.]
 [  0.  -2.   0. 126.]
 [  0.   0.   2. -72.]
 [  0.   0.   0.   1.]]
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]

Using the ITK internal functions for exporting to nifti file format, everything works as expected:

    ants_mni = ants.image_read('./MNI152.nii.gz')
    ants.image_write(ants_mni, './ANTS_MNI152.nii.gz')

    ants_nii = nib.load('./ANTS_MNI152.nii.gz')
    nii_mni = nib.load('./MNI152.nii.gz')

    print(ants_nii.get_qform())
    print(nii_mni.get_qform())
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]

Although this could be used a a quick work around, like in the from_nibabel() implementation, a Nifti/Nibabel class export, should no rely on using hard-drive exports, when all the necessary information could be retrieved while loading the file initially.

@stnava
Copy link
Member

stnava commented Mar 25, 2019

very nice find. a pull request would be helpful. it's probably just changing some indexing somewhere in the function get_qform.

@lmoesch
Copy link
Author

lmoesch commented Mar 25, 2019

As for the pull request, I'm not really sure where to fix this issue (prolly somewhere in the itkNiftiImageIO.cxx). As far as I understand, the direction and origin are parsed from the q_form without paying respect to quatern_b, quatern_c and quatern_d.

For my current code I wrote a wrapper to replace ants.from_nibabel()

def nii_to_ants(nii_img : nib.Nifti1Image):
    q_form = nii_img.get_qform()
    spacing = nii_img.header['pixdim'][1:4]
    
    origin = q_form[:3,3]
    direction = q_form[:3,:3] / spacing

    return ants.from_numpy(data=nii_img.get_data().astype(np.float), origin=origin.tolist(), spacing=spacing.tolist(), direction=direction)

This however assumes 3D data and Nifti1 headers.

@stnava
Copy link
Member

stnava commented Mar 25, 2019

indeed, i am not sure where you would fix this either. when i search:

https://github.com/ANTsX/ANTsPy/search?q=get_qform&unscoped_q=get_qform

i dont find any code in ANTsPy. thoughts?

@lmoesch
Copy link
Author

lmoesch commented Mar 25, 2019

When I understand the code correctly, the direction, origin and spacing is derived from the ITK function itk::ImageIOBase::ReadImageInformation() in ANTsPy/ants/lib/LOCAL_antsImageHeaderInfo.cxx So this has to be an issue in the used ITK version.

@stnava
Copy link
Member

stnava commented Mar 25, 2019

a reproducible / self-contained example of this problem:

import ants
import nibabel as nib
fn = ants.get_data( "ch2" )
ants_mni = ants.image_read( fn )
nii_mni = nib.load( fn )
ants_mni = ants_mni.to_nibabel()
print(ants_mni.get_qform())
print(nii_mni.get_qform())

@stnava
Copy link
Member

stnava commented Mar 25, 2019

the problem may be resolvable by alterations to this bit of code:

array_data = image.numpy()
affine = np.hstack([image.direction*np.diag(image.spacing),np.array(image.origin).reshape(3,1)])
affine = np.vstack([affine, np.array([0,0,0,1.])])
new_img = nib.Nifti1Image(array_data, affine)

@lmoesch
Copy link
Author

lmoesch commented Mar 25, 2019

The problem is, that the metadata already mismatches upon loading the image via ants.image_read().

    nii_img = nib.load('MNI152.nii.gz')
    ants_img = ants.from_nibabel(nii_img) # This is the same as ants_img = ants.image_read('MNI152.nii.gz') 

    print(nii_img.header.get_qform()[:3,3].tolist())
    print(ants_img.origin) 
(90.0, -126.0, -72.0)
(-90.0, 126.0, -72.0)

@stnava
Copy link
Member

stnava commented Mar 25, 2019

this particular difference may be the difference between LPS vs whatever nibabel uses.

@lmoesch
Copy link
Author

lmoesch commented Mar 25, 2019

This issue is not nibabel related, since the qform is saved within the nifti header and should be parsed the same way in all implementations. I checked the values within SPM and MRIcroN and they are the same as the nibabel output, so there has to be an issue by extracting the values via the ReadImageInformation().

As reading and writing through ITK methods results in correct qforms, my guess is that ITK stores the quatern_b, quatern_c and quatern_d values from the header as well to correct the directions.

@stnava
Copy link
Member

stnava commented Mar 25, 2019

can you share the image you are using ( upload to issue ) so we can be on the same page.

@lmoesch
Copy link
Author

lmoesch commented Mar 25, 2019

@stnava
Copy link
Member

stnava commented Mar 25, 2019

upon further reflection, i think that this is just an image orientation difference - ie a difference in coordinate systems (itk uses LPS). @cookpa - can you briefly comment and / or suggest a relevant link ?

@cookpa
Copy link
Member

cookpa commented Mar 25, 2019

ITK uses LPS coordinates and NIFTI RAS. This is handled in the I/O, when ITK converts qform into a direction cosine matrix or writes a cosine matrix back out as a qform.

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx

See methods SetImageIOOrientationFromNIfTI and SetNIfTIOrientationFromImageIO.

According to this page

https://nipy.org/nibabel/coordinate_systems.html

NiBabel always wants RAS coordinates so I think it should be safe to apply the same conversion steps in the to_nibabel and from_nibabel functions.

Some more on this:

https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#a23ebc76de3b60bb1eeabf66cd4dabc48

This is the base ITK class that handles voxel to physical space conversion. It references this page

https://itk.org/Wiki/Proposals:Orientation#Some_notes_on_the_DICOM_convention_and_current_ITK_usage

@stnava
Copy link
Member

stnava commented Mar 26, 2019

more discussion here: https://discourse.itk.org/t/nifti-orientation-issues/431/2 it doesnt "convert images" ... rather, the reference frame is interpreted as LPS. the issue here is, i think, that the ANTsPy "to_nibabel" function is probably wrong.

@lmoesch
Copy link
Author

lmoesch commented Mar 26, 2019

To make a long story short: ITK assumes NIFTI images to be in RAS space and simply flips the first two axes upon loading (code). Therefore it is sufficient to do the same upon to_nibabel().

However, I tried to utilize ants.reorient_image2() in order to achieve this and watched again a strange behavior. When loading an NIFTI image in RAS orientation

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

this results in a direction

[[-1. 0. 0.]
 [0. -1. 0.]
 [0. 0. 1.]]

as the first axes are flipped. But when applying ants.reorient_image2(orientation='RAS') to it, it results in

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. -1.]]

where the Z-Axes is flipped, which looks like a bug within itkOrientImageFilter.hxx to me.

@stnava
Copy link
Member

stnava commented Mar 26, 2019

i do not think that is likely. behavior looks correct here:

>>> mni = ants.image_read(ants.get_data('mni'))

>>> mni
ANTsImage (RPI)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (-90.0, 126.0, -72.0)
	 Direction  : [ 1.  0.  0.  0. -1.  0.  0.  0.  1.]

>>> mniRAS
ANTsImage (RAS)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (-90.0, -91.0, 109.0)
	 Direction  : [ 1.  0.  0.  0.  1.  0.  0.  0. -1.]

>>> mniLPS
ANTsImage (LPS)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (91.0, 126.0, 109.0)
	 Direction  : [-1.  0.  0.  0. -1.  0.  0.  0. -1.]

@stnava
Copy link
Member

stnava commented Mar 26, 2019

however, your comment:

"To make a long story short: ITK assumes NIFTI images to be in RAS space and simply flips the first two axes upon loading (code). Therefore it is sufficient to do the same upon to_nibabel()."

is probably right. i guess the question is what happens in 2D and 4D.

@lmoesch
Copy link
Author

lmoesch commented Mar 26, 2019

i do not think that is likely. behavior looks correct here:

>>> mni = ants.image_read(ants.get_data('mni'))

>>> mni
ANTsImage (RPI)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (-90.0, 126.0, -72.0)
	 Direction  : [ 1.  0.  0.  0. -1.  0.  0.  0.  1.]

>>> mniRAS
ANTsImage (RAS)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (-90.0, -91.0, 109.0)
	 Direction  : [ 1.  0.  0.  0.  1.  0.  0.  0. -1.]

>>> mniLPS
ANTsImage (LPS)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (182, 218, 182)
	 Spacing    : (1.0, 1.0, 1.0)
	 Origin     : (91.0, 126.0, 109.0)
	 Direction  : [-1.  0.  0.  0. -1.  0.  0.  0. -1.]

According to nibabel documentation the z-axis value should be positive in RAS, or am I missing a point here?

@lmoesch
Copy link
Author

lmoesch commented Mar 26, 2019

[...] the question is what happens in 2D and 4D.

let me answer with the code from ITK:

const int max_defined_orientation_dims = ( dims > 3 ) ? 3 : dims;

It flips the first 2 dimensions in every case. And every dimension above the third is simply discarded.

@stnava
Copy link
Member

stnava commented Mar 26, 2019

ok - here is a fix:

>>> import ants
>>> import nibabel as nib
>>> fn = "avg152T1_brain.nii.gz"
>>> ants_img = ants.image_read( fn )
>>> nii_mni = nib.load( fn )
>>> ants_mni = ants_img.to_nibabel()
>>> print(ants_mni.get_qform())
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]
>>> print(nii_mni.get_qform())
[[  -2.    0.    0.   90.]
 [   0.    2.    0. -126.]
 [   0.    0.    2.  -72.]
 [   0.    0.    0.    1.]]
>>> ants.from_nibabel( nii_mni )
ANTsImage (RPI)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (91, 109, 91)
	 Spacing    : (2.0, 2.0, 2.0)
	 Origin     : (-90.0, 126.0, -72.0)
	 Direction  : [ 1.  0.  0.  0. -1.  0.  0.  0.  1.]

>>> ants_img 
ANTsImage (RPI)
	 Pixel Type : float (float32)
	 Components : 1
	 Dimensions : (91, 109, 91)
	 Spacing    : (2.0, 2.0, 2.0)
	 Origin     : (-90.0, 126.0, -72.0)
	 Direction  : [ 1.  0.  0.  0. -1.  0.  0.  0.  1.]

@stnava
Copy link
Member

stnava commented Mar 26, 2019

b15e286

@stnava
Copy link
Member

stnava commented Mar 26, 2019

apologies for pushing this to master but it seems to be correct .... agree?

@lmoesch
Copy link
Author

lmoesch commented Mar 26, 2019

I dont think the last-ind has to be flipped. But I'll check into this.

@lmoesch
Copy link
Author

lmoesch commented Mar 27, 2019

I'll checked into the fix and as far as I see it is not quite correct. ITK changes the signs of the first two elements of all vectors. Fortunately, this simplifies the whole fix to:

affine[:2,:] *= -1

UPDATE: The way the affine is created, you fix works nevertheless. When considering using the real transformation in the future, tho above approach is the way to go.

@stnava
Copy link
Member

stnava commented Mar 27, 2019

thx - i am leaving this open until a test is added to help prevent this from coming up in the future.

@stnava
Copy link
Member

stnava commented Mar 27, 2019

consider 1e73b82

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants