{ "cells": [ { "cell_type": "markdown", "id": "da96a3c0", "metadata": {}, "source": [ "# Aims volume orientation manipulation and nibabel\n", "\n", "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.\n", "\n", "For a complete documentation, please read:\n", "\n", "* https://brainvisa.info/cartodata/doxygen/cartovolumes.html#volume_orient\n", "* https://brainvisa.info/aimsdata/user_doc/coordinates_systems.html\n", "\n", "## Imports\n", "\n", "First of all, let's import the needed modules." ] }, { "cell_type": "code", "execution_count": 1, "id": "4b4d79ab", "metadata": {}, "outputs": [], "source": [ "from soma import aims\n", "import nibabel\n", "import numpy as np\n", "import tempfile\n", "import os" ] }, { "cell_type": "markdown", "id": "7cd6dd28", "metadata": {}, "source": [ "## load a volume using both aims and nibabel\n", "\n", "This will allow to compare" ] }, { "cell_type": "code", "execution_count": 2, "id": "71f66030", "metadata": {}, "outputs": [], "source": [ "image_path = aims.carto.Paths.findResourceFile('anatomical_templates/MNI152_T1_1mm.nii.gz')\n", "avol = aims.read(image_path)\n", "nvol = nibabel.load(image_path)" ] }, { "cell_type": "code", "execution_count": 3, "id": "698b8a58", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "aims shape: (182, 218, 182, 1)\n", "nibabel shape: (182, 218, 182)\n", "aims voxel (100, 102, 62): 5933.0\n", "nibabel: 6464.0\n" ] } ], "source": [ "print('aims shape:', avol.shape)\n", "print('nibabel shape:', nvol.shape)\n", "print('aims voxel (100, 102, 62):', avol[100, 102, 62, 0])\n", "print('nibabel: ', nvol.get_fdata()[100, 102, 62])" ] }, { "cell_type": "markdown", "id": "a6b7f0dd", "metadata": {}, "source": [ "As we see, voxels values are **not the same**: axes are actually flipped.\n", "Aims reads data in LPI orientation.\n", "Nibabel reads data in disk storage orientation (which thus depnds on how data has been saved).\n", "Here, in LAST orientation, the same voxel (100, 102, 62) in LPI should be found at (100, 217-102=115, 181-62=119):" ] }, { "cell_type": "code", "execution_count": 4, "id": "efcfa2d9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nibabel (100, 115, 119): 5933.0\n", "aims (100, 115, 119): 6464.0\n" ] } ], "source": [ "print('nibabel (100, 115, 119):', nvol.get_fdata()[100, 115, 119])\n", "print('aims (100, 115, 119):', avol[100, 115, 119, 0])" ] }, { "cell_type": "markdown", "id": "40258cc0", "metadata": {}, "source": [ "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).\n", "\n", "\n", "## Check AIMS orientation and disk storage orientation" ] }, { "cell_type": "code", "execution_count": 5, "id": "57bad4e6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "disk storage orientation: LAST\n", "AIMS volume orientation: LPIT\n" ] } ], "source": [ "storage_orient_v = avol.storageLayoutOrientation()\n", "storage_orient = avol.referential().orientationStr(storage_orient_v)\n", "print('disk storage orientation:', storage_orient)\n", "print('AIMS volume orientation:', avol.referential().orientationStr())" ] }, { "cell_type": "markdown", "id": "f92f1c9b", "metadata": {}, "source": [ "## Building an AIMS volume from nibabel data\n", "\n", "We must set the correct nibabel orientation on the volume, which is storage orientation." ] }, { "cell_type": "code", "execution_count": 6, "id": "67693634", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "aims from nibabel:\n", "orientation: LAST\n", "shape: (182, 218, 182, 1)\n" ] } ], "source": [ "avol2 = aims.Volume(nvol.get_fdata())\n", "avol2.referential().setOrientation(storage_orient)\n", "print('aims from nibabel:\\norientation:', avol2.referential().orientationStr())\n", "print('shape:', avol2.shape)" ] }, { "cell_type": "code", "execution_count": 7, "id": "751ed595", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "aims from nibabel: 6464.0 == 6464.0 (nibabel)\n", " != 5933.0 (aims)\n" ] } ], "source": [ "print('aims from nibabel:', avol2[100, 102, 62, 0], \"==\", nvol.get_fdata()[100, 102, 62], '(nibabel)')\n", "print(' !=', avol[100, 102, 62, 0], '(aims)')" ] }, { "cell_type": "markdown", "id": "2e1f029a", "metadata": {}, "source": [ "As we are in storage orientation, we should use indices as in nibabel, which actually works !\n", "\n", "### Flipping to another orientation\n", "\n", "If we want to go back to the default AIMS orientation (LPI) without copying the voxels data, it is possible:" ] }, { "cell_type": "code", "execution_count": 8, "id": "c64a19b6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "in LPI orientation: 5933.0 == 5933.0 (aims)\n" ] } ], "source": [ "avol2.flipToOrientation('LPI')\n", "print('in LPI orientation:', avol2[100, 102, 62, 0], '==', avol[100, 102, 62, 0], '(aims)')" ] }, { "cell_type": "markdown", "id": "2ca50c16", "metadata": {}, "source": [ "## Building a Nibabel volume from AIMS data\n", "\n", "First let's see what is wrong when doing things too quickly..." ] }, { "cell_type": "code", "execution_count": 9, "id": "ec3751fc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nibabel from aims (LPI):\n", "shape: (182, 218, 182, 1)\n", "5933.0 !WARNING! matches aims but is wrong in nibabel orientation\n" ] } ], "source": [ "nvol2 = nibabel.Nifti1Image(avol.np, None)\n", "print('nibabel from aims (LPI):\\nshape:', nvol2.shape)\n", "print(nvol2.get_fdata()[100, 102, 62, 0], '!WARNING! matches aims but is wrong in nibabel orientation')" ] }, { "cell_type": "markdown", "id": "77748d4e", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 10, "id": "b7d59e67", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "saving to: /tmp/aims_volume_orientidg2tnle.nii\n" ] } ], "source": [ "tfile = tempfile.mkstemp(prefix='aims_volume_orient', suffix='.nii')\n", "tfilename = tfile[1]\n", "os.close(tfile[0])\n", "print('saving to:', tfilename)" ] }, { "cell_type": "code", "execution_count": 11, "id": "0a419849", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "after nibabel save/reload: [5933.] != 6464.0 (nibabel)\n" ] } ], "source": [ "nibabel.save(nvol2, tfilename)\n", "nvol3 = nibabel.load(tfilename)\n", "print('after nibabel save/reload:', nvol3.get_fdata()[100, 102, 62], '!=', \n", " nvol.get_fdata()[100, 102, 62], '(nibabel)')" ] }, { "cell_type": "markdown", "id": "ddef43f4", "metadata": {}, "source": [ "The reloaded image is not like the original one: it is flipped since the orientation does not match the expected default one.\n", "Alternatively, providing an affine to nibabel.Nifti1Image could fix it. Let's use the aims orientation manipulation instead.\n", "\n", "Cleanup the temporary file, by the way..." ] }, { "cell_type": "code", "execution_count": 12, "id": "0f27fd0d", "metadata": {}, "outputs": [], "source": [ "os.unlink(tfilename)" ] }, { "cell_type": "markdown", "id": "d85678e6", "metadata": {}, "source": [ "If we flip the AIMS volume to storage orientation before buildig the nibabel image from it, thigs are going right." ] }, { "cell_type": "code", "execution_count": 13, "id": "b6da44ee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nibabel from aims (storage):\n", "shape: (182, 218, 182, 1)\n", "6464.0 == 6464.0 (nibabel) *CORRECT!*\n" ] } ], "source": [ "avol.flipToOrientation(storage_orient)\n", "nvol2 = nibabel.Nifti1Image(avol.np, None)\n", "print('nibabel from aims (storage):\\nshape:', nvol2.shape)\n", "print(nvol2.get_fdata()[100, 102, 62, 0], '==', nvol.get_fdata()[100, 102, 62], '(nibabel) *CORRECT!*')" ] }, { "cell_type": "markdown", "id": "8167b1e5", "metadata": {}, "source": [ "## Loading in another orientation\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 14, "id": "b1763c82", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "storage orientation: LAST\n", "indices orientation: LPIT\n", "strides: (4, -728, -158704, 28884128)\n" ] } ], "source": [ "avol2 = aims.read(image_path, options={'orientation': 'storage'})\n", "storage_orient_v = avol.storageLayoutOrientation()\n", "storage_orient = avol.referential().orientationStr(storage_orient_v)\n", "print('storage orientation:', storage_orient)\n", "print('indices orientation:', avol2.referential().orientationStr())\n", "print('strides:', avol2.np.strides)" ] }, { "cell_type": "markdown", "id": "12e059dc", "metadata": {}, "source": [ "Flipping to storage orientation will thus restore the \"default\" strides with 1st axis fastest" ] }, { "cell_type": "code", "execution_count": 15, "id": "a34610ff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "now indices orientation are: LAST\n", "and strides: (4, 728, 158704, 28884128)\n" ] } ], "source": [ "avol2.flipToOrientation(storage_orient)\n", "print('now indices orientation are:', avol2.referential().orientationStr())\n", "print('and strides:', avol2.np.strides)" ] }, { "cell_type": "code", "execution_count": 16, "id": "f8ae1ac9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nibabel strides: (2, 364, 79352)\n" ] } ], "source": [ "print('nibabel strides:', np.asanyarray(nvol.dataobj).strides)" ] }, { "cell_type": "markdown", "id": "42731146", "metadata": {}, "source": [ "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)." ] }, { "cell_type": "code", "execution_count": 17, "id": "840560fc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "nibabel: 6464.0 == 6464.0 (aims storage oriented)\n" ] } ], "source": [ "nvol4 = nibabel.Nifti1Image(avol2.np, None)\n", "print('nibabel:', nvol4.get_fdata()[100, 102, 62, 0], '==', avol2[100, 102, 62, 0], '(aims storage oriented)')" ] }, { "cell_type": "markdown", "id": "76593606", "metadata": {}, "source": [ "## Reorienting memory layout\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 18, "id": "5d02c2e8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "storage orient: LAST\n", "strides: (4, 728, 158704, 28884128)\n", "voxel (100, 102, 62): 6464.0 : same as nibabel\n" ] } ], "source": [ "avol = aims.read(image_path)\n", "# the 2nd param of flipToOrientation() specifies the memory layout\n", "avol.flipToOrientation(storage_orient, storage_orient)\n", "print('storage orient:', avol.referential().orientationStr())\n", "print('strides:', avol.np.strides)\n", "print('voxel (100, 102, 62):', avol[100, 102, 62, 0], ': same as nibabel')" ] }, { "cell_type": "markdown", "id": "cbff8378", "metadata": {}, "source": [ "## Details about headers\n", "\n", "Volume headers are updated when a volume is flipped, in order to adapt transformations to external referentials, and from storage orientation." ] }, { "cell_type": "code", "execution_count": 19, "id": "d5bcdf9e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LPI header:\n", "\n", "{ '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' }\n", "\n", "especially:\n", "transformations: [\n", "[[ -1. 0. 0. 90.]\n", " [ 0. -1. 0. 91.]\n", " [ 0. 0. -1. 109.]\n", " [ 0. 0. 0. 1.]], \n", "[[ -1. 0. 0. 90.]\n", " [ 0. -1. 0. 91.]\n", " [ 0. 0. -1. 109.]\n", " [ 0. 0. 0. 1.]]]\n", "storage_to_memory: \n", "[[ 1. 0. 0. 0.]\n", " [ 0. -1. 0. 217.]\n", " [ 0. 0. -1. 181.]\n", " [ 0. 0. 0. 1.]]\n" ] } ], "source": [ "avol = aims.read(image_path)\n", "print('LPI header:\\n')\n", "print(avol.header())\n", "print('\\nespecially:')\n", "print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])\n", "print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))" ] }, { "cell_type": "code", "execution_count": 20, "id": "57d4f521", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SAR header:\n", "{ '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' }\n", "\n", "especially:\n", "transformations: [\n", "[[ 0. 0. 1. -91.]\n", " [ 0. 1. 0. -126.]\n", " [ 1. 0. 0. -72.]\n", " [ 0. 0. 0. 1.]], \n", "[[ 0. 0. 1. -91.]\n", " [ 0. 1. 0. -126.]\n", " [ 1. 0. 0. -72.]\n", " [ 0. 0. 0. 1.]]]\n", "storage_to_memory: \n", "[[ 0. 0. 1. 0.]\n", " [ 0. 1. 0. 0.]\n", " [ -1. 0. 0. 181.]\n", " [ 0. 0. 0. 1.]]\n" ] } ], "source": [ "avol.flipToOrientation('SAR')\n", "print('SAR header:')\n", "print(avol.header())\n", "print('\\nespecially:')\n", "print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])\n", "print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))" ] }, { "cell_type": "markdown", "id": "3bc84f4b", "metadata": {}, "source": [ "Hopefully, if we get back to initial LPI orientation, the initial header is restored:" ] }, { "cell_type": "code", "execution_count": 21, "id": "32a5b30e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SAR header:\n", "{ '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' }\n", "\n", "especially:\n", "transformations: [\n", "[[ -1. 0. 0. 90.]\n", " [ 0. -1. 0. 91.]\n", " [ 0. 0. -1. 109.]\n", " [ 0. 0. 0. 1.]], \n", "[[ -1. 0. 0. 90.]\n", " [ 0. -1. 0. 91.]\n", " [ 0. 0. -1. 109.]\n", " [ 0. 0. 0. 1.]]]\n", "storage_to_memory: \n", "[[ 1. 0. 0. 0.]\n", " [ 0. -1. 0. 217.]\n", " [ 0. 0. -1. 181.]\n", " [ 0. 0. 0. 1.]]\n" ] } ], "source": [ "avol.flipToOrientation('LPI')\n", "print('SAR header:')\n", "print(avol.header())\n", "print('\\nespecially:')\n", "print('transformations:', [aims.AffineTransformationBase(t) for t in avol.header()['transformations']])\n", "print('storage_to_memory:', aims.AffineTransformationBase(avol.header()['storage_to_memory']))" ] }, { "cell_type": "markdown", "id": "4fd1d7bb", "metadata": {}, "source": [ "All this ensure that the volume data is not changed in space in regards to other coordinates systems: only the internal referential is changed." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.0" } }, "nbformat": 4, "nbformat_minor": 5 }