Improve orientation robustness in SH/DTI scripts

src/scilpy/cli/scil_dti_convert_tensors.py

@@ -9,9 +9,10 @@
 import argparse
 import logging
 
-import nibabel as nib
 import numpy as np
 
+from scilpy.io.stateful_image import StatefulImage
+
 from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg,
                              assert_inputs_exist, assert_outputs_exist)
 from scilpy.io.tensor import (supported_tensor_formats,
@@ -52,14 +53,14 @@ def main():
     assert_inputs_exist(parser, args.in_file)
     assert_outputs_exist(parser, args, args.out_file)
 
-    in_tensors_img = nib.load(args.in_file)
-    in_tensors = in_tensors_img.get_fdata(dtype=np.float32)
+    in_tensors_simg = StatefulImage.load(args.in_file)
+    in_tensors = in_tensors_simg.get_fdata(dtype=np.float32)
 
     out_tensors = convert_tensor_format(in_tensors, args.in_format,
                                         args.out_format)
-    out_tensors_img = nib.Nifti1Image(
-        out_tensors.astype(np.float32), in_tensors_img.affine)
-    nib.save(out_tensors_img, args.out_file)
+    StatefulImage.create_from(
+        out_tensors.astype(np.float32), in_tensors_simg,
+        is_orientation=False).save(args.out_file)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_sh_convert.py

@@ -12,9 +12,9 @@
 import logging
 
 from dipy.data import get_sphere
-import nibabel as nib
 import numpy as np
 
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.reconst.sh import convert_sh_basis
 from scilpy.io.utils import (add_overwrite_arg, add_sh_basis_args,
                              add_processes_arg, add_verbose_arg,
@@ -50,21 +50,26 @@ def main():
     assert_outputs_exist(parser, args, args.out_sh)
 
     sphere = get_sphere(name='repulsion724').subdivide(n=1)
-    img = nib.load(args.in_sh)
-    data = img.get_fdata(dtype=np.float32)
-
     in_sh_basis, is_in_legacy, out_sh_basis, is_out_legacy \
         = parse_sh_basis_arg(args)
 
+    sh_simg = StatefulImage.load(args.in_sh, is_orientation=True,
+                                 sh_basis=in_sh_basis, is_legacy=is_in_legacy)
+    data = sh_simg.get_fdata(dtype=np.float32)
+
     new_data = convert_sh_basis(data, sphere,
                                 input_basis=in_sh_basis,
                                 output_basis=out_sh_basis,
                                 is_input_legacy=is_in_legacy,
                                 is_output_legacy=is_out_legacy,
                                 nbr_processes=args.nbr_processes)
 
-    nib.save(nib.Nifti1Image(new_data, img.affine, header=img.header),
-             args.out_sh)
+    out_simg = StatefulImage.create_from(
+        new_data.astype(np.float32), sh_simg,
+        is_orientation=True,
+        sh_basis=out_sh_basis,
+        is_legacy=is_out_legacy)
+    out_simg.save(args.out_sh)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_sh_fusion.py

@@ -28,12 +28,12 @@
 import argparse
 import logging
 
-import nibabel as nib
 import numpy as np
 
-from scilpy.io.image import assert_same_resolution
-from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
-                             assert_outputs_exist, add_verbose_arg)
+from scilpy.io.stateful_image import StatefulImage
+from scilpy.io.utils import (add_overwrite_arg, assert_headers_compatible,
+                             assert_inputs_exist, assert_outputs_exist,
+                             add_verbose_arg)
 from scilpy.version import version_string
 
 
@@ -60,20 +60,22 @@ def main():
 
     assert_inputs_exist(parser, args.in_shs)
     assert_outputs_exist(parser, args, args.out_sh)
-    assert_same_resolution(args.in_shs)
+    assert_headers_compatible(parser, args.in_shs)
 
-    first_im = nib.load(args.in_shs[0])
-    out_coeffs = first_im.get_fdata(dtype=np.float32)
+    first_simg = StatefulImage.load(args.in_shs[0], is_orientation=True)
+    out_coeffs = first_simg.get_fdata(dtype=np.float32)
 
     for sh_file in args.in_shs[1:]:
-        im_dat = nib.load(sh_file).get_fdata(dtype=np.float32)
+        curr_simg = StatefulImage.load(sh_file, is_orientation=True)
+        curr_simg.reorient(first_simg.axcodes)
+        im_dat = curr_simg.get_fdata(dtype=np.float32)
 
         out_coeffs = np.where(np.abs(im_dat) > np.abs(out_coeffs),
                               im_dat, out_coeffs)
 
-    nib.save(nib.Nifti1Image(out_coeffs, first_im.affine,
-                             header=first_im.header),
-             args.out_sh)
+    StatefulImage.create_from(
+        out_coeffs.astype(np.float32), first_simg,
+        is_orientation=True).save(args.out_sh)
 
 
 if __name__ == '__main__':

src/scilpy/cli/scil_sh_to_aodf.py

@@ -28,9 +28,11 @@
 import argparse
 import logging
 import time
-import nibabel as nib
+
 import numpy as np
 
+from scilpy.io.stateful_image import StatefulImage
+
 from dipy.data import SPHERE_FILES
 from dipy.reconst.shm import sph_harm_ind_list
 from scilpy.reconst.utils import get_sh_order_and_fullness
@@ -137,11 +139,11 @@ def main():
     assert_inputs_exist(parser, args.in_sh)
 
     # Prepare data
-    sh_img = nib.load(args.in_sh)
-    data = sh_img.get_fdata(dtype=np.float32)
-
-    sh_order, full_basis = get_sh_order_and_fullness(data.shape[-1])
     sh_basis, is_legacy = parse_sh_basis_arg(args)
+    sh_simg = StatefulImage.load(args.in_sh, is_orientation=True,
+                                 sh_basis=sh_basis, is_legacy=is_legacy)
+    data = sh_simg.get_fdata(dtype=np.float32)
+    sh_order, full_basis = get_sh_order_and_fullness(data.shape[-1])
 
     t0 = time.perf_counter()
     logging.info('Filtering SH image.')
@@ -176,13 +178,17 @@ def main():
     logging.info('Elapsed time (s): {0}'.format(t1 - t0))
 
     logging.info('Saving filtered SH to file {0}.'.format(args.out_sh))
-    nib.save(nib.Nifti1Image(asym_sh, sh_img.affine), args.out_sh)
+    StatefulImage.create_from(
+        asym_sh.astype(np.float32), sh_simg,
+        is_orientation=True).save(args.out_sh)
 
     if args.out_sym:
         _, orders = sph_harm_ind_list(sh_order, full_basis=True)
         logging.info('Saving symmetric SH to file {0}.'.format(args.out_sym))
-        nib.save(nib.Nifti1Image(asym_sh[..., orders % 2 == 0], sh_img.affine),
-                 args.out_sym)
+        sym_sh = asym_sh[..., orders % 2 == 0].astype(np.float32)
+        StatefulImage.create_from(
+            sym_sh, sh_simg,
+            is_orientation=True).save(args.out_sym)
 
 
 if __name__ == '__main__':

src/scilpy/cli/scil_sh_to_rish.py

@@ -24,10 +24,10 @@
 import logging
 
 from dipy.reconst.shm import order_from_ncoef, sph_harm_ind_list
-import nibabel as nib
 import numpy as np
 
 from scilpy.io.image import get_data_as_mask
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist,
                              assert_outputs_exist, add_verbose_arg,
                              assert_headers_compatible)
@@ -68,10 +68,13 @@ def main():
     assert_headers_compatible(parser, args.in_sh, optional=args.mask)
 
     # Load data
-    sh_img = nib.load(args.in_sh)
-    sh = sh_img.get_fdata(dtype=np.float32)
-    mask = get_data_as_mask(nib.load(args.mask),
-                            dtype=bool) if args.mask else None
+    sh_simg = StatefulImage.load(args.in_sh, is_orientation=True)
+    sh = sh_simg.get_fdata(dtype=np.float32)
+    mask = None
+    if args.mask:
+        mask_simg = StatefulImage.load(args.mask)
+        mask_simg.reorient(sh_simg.axcodes)
+        mask = get_data_as_mask(mask_simg, dtype=bool)
 
     # Precompute output filenames to check if they exist
     sh_order = order_from_ncoef(sh.shape[-1], full_basis=args.full_basis)
@@ -89,7 +92,9 @@ def main():
     # Save each RISH feature as a separate file
     for i, fname in enumerate(output_fnames):
         logging.info("Saving {}".format(fname))
-        nib.save(nib.Nifti1Image(rish[..., i], sh_img.affine), fname)
+        StatefulImage.create_from(
+            rish[..., i].astype(np.float32), sh_simg,
+            is_orientation=False).save(fname)
 
 
 if __name__ == '__main__':

src/scilpy/cli/tests/test_sh_to_aodf.py

@@ -8,8 +8,8 @@
 import numpy as np
 import pytest
 
-from scilpy.io.dvc import pull_test_case_package
 from scilpy.gpuparallel.opencl_utils import have_opencl
+from scilpy.io.dvc import pull_test_case_package
 
 # If they already exist, this only takes 5 seconds (check md5sum)
 test_data_root = pull_test_case_package("aodf")
@@ -127,3 +127,49 @@ def test_cosine_method(script_runner, in_fodf, out_fodf, monkeypatch):
     test_fodf = nib.load(out_fodf)
 
     assert np.allclose(ret_fodf.get_fdata(), test_fodf.get_fdata())
+
+
+def test_non_ras_sh_to_aodf(script_runner, monkeypatch):
+    monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
+
+    in_fodf = os.path.join(test_data_root, 'fodf_descoteaux07_sub.nii.gz')
+    img = nib.load(in_fodf)
+    # Use a small crop in the middle of the volume for fast execution
+    size = 7
+    start = [s // 2 - size // 2 for s in img.shape[:3]]
+    img_crop = img.slicer[start[0]:start[0] + size,
+                          start[1]:start[1] + size,
+                          start[2]:start[2] + size]
+    data_crop = img_crop.get_fdata(dtype=np.float32)
+    aff_ras = img_crop.affine.copy()
+    nib.save(nib.Nifti1Image(data_crop, aff_ras), 'sh_ras.nii.gz')
+
+    aff_las = aff_ras.copy()
+    aff_las[0, 0] = -aff_las[0, 0]
+    aff_las[0, 3] = aff_las[0, 3] + (data_crop.shape[0] - 1) * aff_ras[0, 0]
+    data_las = data_crop[::-1].copy()
+    nib.save(nib.Nifti1Image(data_las, aff_las), 'sh_las.nii.gz')
+
+    # Run on RAS
+    ret_ras = script_runner.run([
+        'scil_sh_to_aodf', 'sh_ras.nii.gz', 'out_ras.nii.gz',
+        '--sphere', 'repulsion100', '--method', 'cosine', '-f',
+        '--sh_basis', 'descoteaux07_legacy'])
+    assert ret_ras.success
+
+    # Run on LAS
+    ret_las = script_runner.run([
+        'scil_sh_to_aodf', 'sh_las.nii.gz', 'out_las.nii.gz',
+        '--sphere', 'repulsion100', '--method', 'cosine', '-f',
+        '--sh_basis', 'descoteaux07_legacy'])
+    assert ret_las.success
+
+    # Verify orientation is preserved
+    out_l = nib.load('out_las.nii.gz')
+    assert nib.orientations.aff2axcodes(out_l.affine) == ('L', 'A', 'S')
+
+    out_r = nib.load('out_ras.nii.gz')
+    assert nib.orientations.aff2axcodes(out_r.affine) == ('R', 'A', 'S')
+
+    # Compare output values
+    assert np.allclose(out_r.get_fdata(), out_l.get_fdata()[::-1], atol=1e-5)