diff --git a/src/scilpy/cli/scil_dti_convert_tensors.py b/src/scilpy/cli/scil_dti_convert_tensors.py index d41b25a1f..5c2a2feae 100755 --- a/src/scilpy/cli/scil_dti_convert_tensors.py +++ b/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__": diff --git a/src/scilpy/cli/scil_sh_convert.py b/src/scilpy/cli/scil_sh_convert.py index eca4dd0b5..cbc4f8431 100755 --- a/src/scilpy/cli/scil_sh_convert.py +++ b/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,12 +50,13 @@ 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, @@ -63,8 +64,12 @@ def main(): 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__": diff --git a/src/scilpy/cli/scil_sh_fusion.py b/src/scilpy/cli/scil_sh_fusion.py index 9a8523910..97dbc8b9e 100755 --- a/src/scilpy/cli/scil_sh_fusion.py +++ b/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__': diff --git a/src/scilpy/cli/scil_sh_to_aodf.py b/src/scilpy/cli/scil_sh_to_aodf.py index 46306d198..000902e51 100755 --- a/src/scilpy/cli/scil_sh_to_aodf.py +++ b/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__': diff --git a/src/scilpy/cli/scil_sh_to_rish.py b/src/scilpy/cli/scil_sh_to_rish.py index a3733966c..3ae4a9c60 100755 --- a/src/scilpy/cli/scil_sh_to_rish.py +++ b/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__': diff --git a/src/scilpy/cli/tests/test_sh_to_aodf.py b/src/scilpy/cli/tests/test_sh_to_aodf.py index 6d0c198e7..28c21581e 100644 --- a/src/scilpy/cli/tests/test_sh_to_aodf.py +++ b/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)