Aims volume orientation manipulation and nibabel¶
In Aims 5.2, indices axes orientation, memory layout and storage (disk) layout are known and can be manipulated and changed. This allows to use the axes we want, and convert from/to nibabel objects.
For a complete documentation, please read:
https://brainvisa.info/cartodata/doxygen/cartovolumes.html#volume_orient
https://brainvisa.info/aimsdata/user_doc/coordinates_systems.html
Imports¶
First of all, let’s import the needed modules.
[1]:
from soma import aims
import nibabel
import numpy as np
import tempfile
import os
load a volume using both aims and nibabel¶
This will allow to compare
[2]:
image_path = aims.carto.Paths.findResourceFile('anatomical_templates/MNI152_T1_1mm.nii.gz')
avol = aims.read(image_path)
nvol = nibabel.load(image_path)
[3]:
print('aims shape:', avol.shape)
print('nibabel shape:', nvol.shape)
print('aims voxel (100, 102, 62):', avol[100, 102, 62, 0])
print('nibabel: ', nvol.get_fdata()[100, 102, 62])
aims shape: (182, 218, 182, 1)
nibabel shape: (182, 218, 182)
aims voxel (100, 102, 62): 5933.0
nibabel: 6464.0
As we see, voxels values are not the same: axes are actually flipped. Aims reads data in LPI orientation. Nibabel reads data in disk storage orientation (which thus depnds on how data has been saved). Here, in LAST orientation, the same voxel (100, 102, 62) in LPI should be found at (100, 217-102=115, 181-62=119):
[4]:
print('nibabel (100, 115, 119):', nvol.get_fdata()[100, 115, 119])
print('aims (100, 115, 119):', avol[100, 115, 119, 0])
nibabel (100, 115, 119): 5933.0
aims (100, 115, 119): 6464.0
which, for the nipy volume, is the same as voxel (100, 102, 62) in AIMS LPI, and symetrically, for aims, the same as the nibabel volume at (100, 102, 62).
Check AIMS orientation and disk storage orientation¶
[5]:
storage_orient_v = avol.storageLayoutOrientation()
storage_orient = avol.referential().orientationStr(storage_orient_v)
print('disk storage orientation:', storage_orient)
print('AIMS volume orientation:', avol.referential().orientationStr())
disk storage orientation: LAST
AIMS volume orientation: LPIT
Building an AIMS volume from nibabel data¶
We must set the correct nibabel orientation on the volume, which is storage orientation.
[6]:
avol2 = aims.Volume(nvol.get_fdata())
avol2.referential().setOrientation(storage_orient)
print('aims from nibabel:\norientation:', avol2.referential().orientationStr())
print('shape:', avol2.shape)
aims from nibabel:
orientation: LAST
shape: (182, 218, 182, 1)
[7]:
print('aims from nibabel:', avol2[100, 102, 62, 0], "==", nvol.get_fdata()[100, 102, 62], '(nibabel)')
print(' !=', avol[100, 102, 62, 0], '(aims)')
aims from nibabel: 6464.0 == 6464.0 (nibabel)
!= 5933.0 (aims)
As we are in storage orientation, we should use indices as in nibabel, which actually works !
Flipping to another orientation¶
If we want to go back to the default AIMS orientation (LPI) without copying the voxels data, it is possible:
[8]:
avol2.flipToOrientation('LPI')
print('in LPI orientation:', avol2[100, 102, 62, 0], '==', avol[100, 102, 62, 0], '(aims)')
in LPI orientation: 5933.0 == 5933.0 (aims)
Building a Nibabel volume from AIMS data¶
First let’s see what is wrong when doing things too quickly…
[9]:
nvol2 = nibabel.Nifti1Image(avol.np, None)
print('nibabel from aims (LPI):\nshape:', nvol2.shape)
print(nvol2.get_fdata()[100, 102, 62, 0], '!WARNING! matches aims but is wrong in nibabel orientation')
nibabel from aims (LPI):
shape: (182, 218, 182, 1)
5933.0 !WARNING! matches aims but is wrong in nibabel orientation
Actually the LPI-oriented array is used directly in nibabel. This is not wrong in itself, but nibabel does not know it is LPI, so saving the image will result in a flipped image.
[10]:
tfile = tempfile.mkstemp(prefix='aims_volume_orient', suffix='.nii')
tfilename = tfile[1]
os.close(tfile[0])
print('saving to:', tfilename)
saving to: /tmp/aims_volume_orientidg2tnle.nii
[11]:
nibabel.save(nvol2, tfilename)
nvol3 = nibabel.load(tfilename)
print('after nibabel save/reload:', nvol3.get_fdata()[100, 102, 62], '!=',
nvol.get_fdata()[100, 102, 62], '(nibabel)')
after nibabel save/reload: [5933.] != 6464.0 (nibabel)
The reloaded image is not like the original one: it is flipped since the orientation does not match the expected default one. Alternatively, providing an affine to nibabel.Nifti1Image could fix it. Let’s use the aims orientation manipulation instead.
Cleanup the temporary file, by the way…
[12]:
os.unlink(tfilename)
If we flip the AIMS volume to storage orientation before buildig the nibabel image from it, thigs are going right.
[13]:
avol.flipToOrientation(storage_orient)
nvol2 = nibabel.Nifti1Image(avol.np, None)
print('nibabel from aims (storage):\nshape:', nvol2.shape)
print(nvol2.get_fdata()[100, 102, 62, 0], '==', nvol.get_fdata()[100, 102, 62], '(nibabel) *CORRECT!*')
nibabel from aims (storage):
shape: (182, 218, 182, 1)
6464.0 == 6464.0 (nibabel) *CORRECT!*
Loading in another orientation¶
Aims allows to load a volume in another memory orientation, and especially the storage layout orientation. The memory layout will then match the asked orientation, but the indices will still be left to the default LPI:
[14]:
avol2 = aims.read(image_path, options={'orientation': 'storage'})
storage_orient_v = avol.storageLayoutOrientation()
storage_orient = avol.referential().orientationStr(storage_orient_v)
print('storage orientation:', storage_orient)
print('indices orientation:', avol2.referential().orientationStr())
print('strides:', avol2.np.strides)
storage orientation: LAST
indices orientation: LPIT
strides: (4, -728, -158704, 28884128)
Flipping to storage orientation will thus restore the “default” strides with 1st axis fastest
[15]:
avol2.flipToOrientation(storage_orient)
print('now indices orientation are:', avol2.referential().orientationStr())
print('and strides:', avol2.np.strides)
now indices orientation are: LAST
and strides: (4, 728, 158704, 28884128)
[16]:
print('nibabel strides:', np.asanyarray(nvol.dataobj).strides)
nibabel strides: (2, 364, 79352)
which demonstrates the same memory layout orientation, left aside that the data type is not the same (16 bit ints for nibabel, 32 bit floats for aims since nibabel does not apply scale factors in dataobj, and get_fdata() would return 64 bit doubles).
[17]:
nvol4 = nibabel.Nifti1Image(avol2.np, None)
print('nibabel:', nvol4.get_fdata()[100, 102, 62, 0], '==', avol2[100, 102, 62, 0], '(aims storage oriented)')
nibabel: 6464.0 == 6464.0 (aims storage oriented)
Reorienting memory layout¶
It is also possible top flip Aims volumes in memory, changing the memory layout. In this case voxels data will be actually copied and moved. The operation will (for now) require to temporarily duplicate the full data block in memory - this may be improved in the future.
[18]:
avol = aims.read(image_path)
# the 2nd param of flipToOrientation() specifies the memory layout
avol.flipToOrientation(storage_orient, storage_orient)
print('storage orient:', avol.referential().orientationStr())
print('strides:', avol.np.strides)
print('voxel (100, 102, 62):', avol[100, 102, 62, 0], ': same as nibabel')
storage orient: LAST
strides: (4, 728, 158704, 28884128)
voxel (100, 102, 62): 6464.0 : same as nibabel
Details about headers¶
Volume headers are updated when a volume is flipped, in order to adapt transformations to external referentials, and from storage orientation.
[19]:
avol = aims.read(image_path)
print('LPI header:\n')
print(avol.header())
print('\nespecially:')
print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])
print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))
LPI header:
{ 'volume_dimension' : [ 182, 218, 182, 1 ], 'sizeX' : 182, 'sizeY' : 218, 'sizeZ' : 182, 'sizeT' : 1, 'referential' : '49e6b349-b115-211a-c8b9-20d0ece9846d', 'disk_data_type' : 'S16', 'bits_allocated' : 16, 'scale_factor' : 1, 'scale_offset' : 0, 'data_type' : 'FLOAT', 'scale_factor_applied' : 0, 'possible_data_types' : [ 'FLOAT', 'S16', 'DOUBLE' ], 'cal_min' : 3000, 'cal_max' : 8000, 'freq_dim' : 0, 'phase_dim' : 0, 'slice_dim' : 0, 'slice_code' : 0, 'slice_start' : 0, 'slice_end' : 0, 'slice_duration' : 0, 'storage_to_memory' : [ 1, 0, 0, 0, 0, -1, 0, 217, 0, 0, -1, 181, 0, 0, 0, 1 ], 'voxel_size' : [ 1, 1, 1 ], 'tr' : 1, 'referentials' : [ 'Talairach-MNI template-SPM', 'Talairach-MNI template-SPM' ], 'transformations' : [ [ -1, 0, 0, 90, 0, -1, 0, 91, 0, 0, -1, 109, 0, 0, 0, 1 ], [ -1, 0, 0, 90, 0, -1, 0, 91, 0, 0, -1, 109, 0, 0, 0, 1 ] ], 'toffset' : 0, 'xyz_units' : 2, 'time_units' : 8, 'descrip' : 'FSL3.3', 'aux_file' : '', 'nifti_type' : 1, 'object_type' : 'Volume', 'uuid' : 'b4cc9054-f49d-cd20-8e95-ac17e247b769', 'file_type' : 'NIFTI-1' }
especially:
transformations: [<soma.aims.AffineTransformationBase object at 0x7167d20df7f0>
[[ -1. 0. 0. 90.]
[ 0. -1. 0. 91.]
[ 0. 0. -1. 109.]
[ 0. 0. 0. 1.]], <soma.aims.AffineTransformationBase object at 0x7167e7b67760>
[[ -1. 0. 0. 90.]
[ 0. -1. 0. 91.]
[ 0. 0. -1. 109.]
[ 0. 0. 0. 1.]]]
storage_to_memory: <soma.aims.AffineTransformationBase object at 0x7167e7b67760>
[[ 1. 0. 0. 0.]
[ 0. -1. 0. 217.]
[ 0. 0. -1. 181.]
[ 0. 0. 0. 1.]]
[20]:
avol.flipToOrientation('SAR')
print('SAR header:')
print(avol.header())
print('\nespecially:')
print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])
print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))
SAR header:
{ 'volume_dimension' : [ 182, 218, 182, 1 ], 'sizeX' : 182, 'sizeY' : 218, 'sizeZ' : 182, 'sizeT' : 1, 'referential' : '96275a96-cc6f-11f0-9838-c0a810ba0cb6', 'disk_data_type' : 'S16', 'bits_allocated' : 16, 'scale_factor' : 1, 'scale_offset' : 0, 'data_type' : 'FLOAT', 'scale_factor_applied' : 0, 'possible_data_types' : [ 'FLOAT', 'S16', 'DOUBLE' ], 'cal_min' : 3000, 'cal_max' : 8000, 'freq_dim' : 0, 'phase_dim' : 0, 'slice_dim' : 0, 'slice_code' : 0, 'slice_start' : 0, 'slice_end' : 0, 'slice_duration' : 0, 'storage_to_memory' : [ 0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, 181, 0, 0, 0, 1 ], 'voxel_size' : [ 1, 1, 1, 1 ], 'tr' : 1, 'referentials' : [ 'Talairach-MNI template-SPM', 'Talairach-MNI template-SPM' ], 'transformations' : [ [ 0, 0, 1, -91, 0, 1, 0, -126, 1, 0, 0, -72, 0, 0, 0, 1 ], [ 0, 0, 1, -91, 0, 1, 0, -126, 1, 0, 0, -72, 0, 0, 0, 1 ] ], 'toffset' : 0, 'xyz_units' : 2, 'time_units' : 8, 'descrip' : 'FSL3.3', 'aux_file' : '', 'nifti_type' : 1, 'object_type' : 'Volume', 'uuid' : 'b4cc9054-f49d-cd20-8e95-ac17e247b769', 'file_type' : 'NIFTI-1' }
especially:
transformations: [<soma.aims.AffineTransformationBase object at 0x7167d20df9a0>
[[ 0. 0. 1. -91.]
[ 0. 1. 0. -126.]
[ 1. 0. 0. -72.]
[ 0. 0. 0. 1.]], <soma.aims.AffineTransformationBase object at 0x7167d20df910>
[[ 0. 0. 1. -91.]
[ 0. 1. 0. -126.]
[ 1. 0. 0. -72.]
[ 0. 0. 0. 1.]]]
storage_to_memory: <soma.aims.AffineTransformationBase object at 0x7167e7b676d0>
[[ 0. 0. 1. 0.]
[ 0. 1. 0. 0.]
[ -1. 0. 0. 181.]
[ 0. 0. 0. 1.]]
Hopefully, if we get back to initial LPI orientation, the initial header is restored:
[21]:
avol.flipToOrientation('LPI')
print('SAR header:')
print(avol.header())
print('\nespecially:')
print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])
print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))
SAR header:
{ 'volume_dimension' : [ 182, 218, 182, 1 ], 'sizeX' : 182, 'sizeY' : 218, 'sizeZ' : 182, 'sizeT' : 1, 'referential' : '6ae542c6-cc6f-11f0-9838-c0a810ba0cb6', 'disk_data_type' : 'S16', 'bits_allocated' : 16, 'scale_factor' : 1, 'scale_offset' : 0, 'data_type' : 'FLOAT', 'scale_factor_applied' : 0, 'possible_data_types' : [ 'FLOAT', 'S16', 'DOUBLE' ], 'cal_min' : 3000, 'cal_max' : 8000, 'freq_dim' : 0, 'phase_dim' : 0, 'slice_dim' : 0, 'slice_code' : 0, 'slice_start' : 0, 'slice_end' : 0, 'slice_duration' : 0, 'storage_to_memory' : [ 1, 0, 0, 0, 0, -1, 0, 217, 0, 0, -1, 181, 0, 0, 0, 1 ], 'voxel_size' : [ 1, 1, 1, 1 ], 'tr' : 1, 'referentials' : [ 'Talairach-MNI template-SPM', 'Talairach-MNI template-SPM' ], 'transformations' : [ [ -1, 0, 0, 90, 0, -1, 0, 91, 0, 0, -1, 109, 0, 0, 0, 1 ], [ -1, 0, 0, 90, 0, -1, 0, 91, 0, 0, -1, 109, 0, 0, 0, 1 ] ], 'toffset' : 0, 'xyz_units' : 2, 'time_units' : 8, 'descrip' : 'FSL3.3', 'aux_file' : '', 'nifti_type' : 1, 'object_type' : 'Volume', 'uuid' : 'b4cc9054-f49d-cd20-8e95-ac17e247b769', 'file_type' : 'NIFTI-1' }
especially:
transformations: [<soma.aims.AffineTransformationBase object at 0x7167d20df9a0>
[[ -1. 0. 0. 90.]
[ 0. -1. 0. 91.]
[ 0. 0. -1. 109.]
[ 0. 0. 0. 1.]], <soma.aims.AffineTransformationBase object at 0x7167d20df880>
[[ -1. 0. 0. 90.]
[ 0. -1. 0. 91.]
[ 0. 0. -1. 109.]
[ 0. 0. 0. 1.]]]
storage_to_memory: <soma.aims.AffineTransformationBase object at 0x7167d20df9a0>
[[ 1. 0. 0. 0.]
[ 0. -1. 0. 217.]
[ 0. 0. -1. 181.]
[ 0. 0. 0. 1.]]
All this ensure that the volume data is not changed in space in regards to other coordinates systems: only the internal referential is changed.