Reorientation to Canonical RAS with gradient support

src/scilpy/cli/scil_dti_metrics.py

@@ -30,7 +30,6 @@
 
 from dipy.core.gradients import gradient_table
 import dipy.denoise.noise_estimate as ne
-from dipy.io.gradients import read_bvals_bvecs
 from dipy.reconst.dti import (TensorModel, color_fa, fractional_anisotropy,
                               geodesic_anisotropy, mean_diffusivity,
                               axial_diffusivity, norm,
@@ -41,6 +40,7 @@
 from scilpy.dwi.operations import compute_residuals, \
     compute_residuals_statistics
 from scilpy.io.image import get_data_as_mask
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.io.utils import (add_b0_thresh_arg, add_overwrite_arg,
                              add_skip_b0_check_arg, add_verbose_arg,
                              assert_inputs_exist, assert_outputs_exist,
@@ -184,19 +184,29 @@ def main():
     assert_headers_compatible(parser, args.in_dwi, args.mask)
 
     # Loading
-    img = nib.load(args.in_dwi)
-    data = img.get_fdata(dtype=np.float32)
-    affine = img.affine
-    mask = get_data_as_mask(nib.load(args.mask),
-                            dtype=bool) if args.mask else None
+    simg = StatefulImage.load(args.in_dwi)
+    simg.load_gradients(args.in_bval, args.in_bvec)
 
-    logging.info('Tensor estimation with the {} method...'.format(args.method))
-    bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec)
+    # Reorient to RAS for DIPY
+    simg.to_ras()
+
+    data = simg.get_fdata(dtype=np.float32)
+    affine = simg.affine
+    bvals = simg.bvals
+    bvecs = simg.bvecs
 
     if not is_normalized_bvecs(bvecs):
-        logging.warning('Your b-vectors do not seem normalized...')
+        logger.warning('Your b-vectors do not seem normalized...')
         bvecs = normalize_bvecs(bvecs)
 
+    mask = None
+    if args.mask:
+        mask_simg = StatefulImage.load(args.mask)
+        mask_simg.reorient(simg.axcodes)
+        mask = get_data_as_mask(mask_simg, dtype=bool)
+
+    logging.info('Tensor estimation with the {} method...'.format(args.method))
+
     # How the b0_threshold is used: gtab.b0s_mask is used
     # 1) In TensorModel in Dipy:
     #       - The S0 images used as any other image in the design matrix and in
@@ -231,7 +241,8 @@ def main():
 
         fiber_tensors = nib.Nifti1Image(
             tensor_vals_reordered.astype(np.float32), affine)
-        nib.save(fiber_tensors, args.tensor)
+        # Use StatefulImage.create_from to ensure original orientation
+        StatefulImage.create_from(fiber_tensors, simg).save(args.tensor)
 
         del tensor_vals, fiber_tensors, tensor_vals_reordered
 
@@ -240,29 +251,34 @@ def main():
         FA[np.isnan(FA)] = 0
         FA = np.clip(FA, 0, 1)
         if args.fa:
-            nib.save(nib.Nifti1Image(FA.astype(np.float32), affine), args.fa)
+            fa_img = nib.Nifti1Image(FA.astype(np.float32), affine)
+            StatefulImage.create_from(fa_img, simg).save(args.fa)
 
         if args.rgb:
             RGB = color_fa(FA, tenfit.evecs)
-            nib.save(nib.Nifti1Image(np.array(255 * RGB, 'uint8'), affine),
-                     args.rgb)
+            rgb_img = nib.Nifti1Image(np.array(255 * RGB, 'uint8'), affine)
+            StatefulImage.create_from(rgb_img, simg).save(args.rgb)
 
     if args.ga:
         GA = geodesic_anisotropy(tenfit.evals)
         GA[np.isnan(GA)] = 0
-        nib.save(nib.Nifti1Image(GA.astype(np.float32), affine), args.ga)
+        ga_img = nib.Nifti1Image(GA.astype(np.float32), affine)
+        StatefulImage.create_from(ga_img, simg).save(args.ga)
 
     if args.md:
         MD = mean_diffusivity(tenfit.evals)
-        nib.save(nib.Nifti1Image(MD.astype(np.float32), affine), args.md)
+        md_img = nib.Nifti1Image(MD.astype(np.float32), affine)
+        StatefulImage.create_from(md_img, simg).save(args.md)
 
     if args.ad:
         AD = axial_diffusivity(tenfit.evals)
-        nib.save(nib.Nifti1Image(AD.astype(np.float32), affine), args.ad)
+        ad_img = nib.Nifti1Image(AD.astype(np.float32), affine)
+        StatefulImage.create_from(ad_img, simg).save(args.ad)
 
     if args.rd:
         RD = radial_diffusivity(tenfit.evals)
-        nib.save(nib.Nifti1Image(RD.astype(np.float32), affine), args.rd)
+        rd_img = nib.Nifti1Image(RD.astype(np.float32), affine)
+        StatefulImage.create_from(rd_img, simg).save(args.rd)
 
     if args.mode:
         # Compute tensor mode
@@ -271,31 +287,37 @@ def main():
         # Since the mode computation can generate NANs when not masked,
         # we need to remove them.
         non_nan_indices = np.isfinite(inter_mode)
-        mode = np.zeros(inter_mode.shape)
-        mode[non_nan_indices] = inter_mode[non_nan_indices]
-        nib.save(nib.Nifti1Image(mode.astype(np.float32), affine), args.mode)
+        mode_data = np.zeros(inter_mode.shape)
+        mode_data[non_nan_indices] = inter_mode[non_nan_indices]
+        mode_img = nib.Nifti1Image(mode_data.astype(np.float32), affine)
+        StatefulImage.create_from(mode_img, simg).save(args.mode)
 
     if args.norm:
         NORM = norm(tenfit.quadratic_form)
-        nib.save(nib.Nifti1Image(NORM.astype(np.float32), affine), args.norm)
+        norm_img = nib.Nifti1Image(NORM.astype(np.float32), affine)
+        StatefulImage.create_from(norm_img, simg).save(args.norm)
 
     if args.evecs:
-        evecs = tenfit.evecs.astype(np.float32)
-        nib.save(nib.Nifti1Image(evecs, affine), args.evecs)
+        evecs_data = tenfit.evecs.astype(np.float32)
+        evecs_img = nib.Nifti1Image(evecs_data, affine)
+        StatefulImage.create_from(evecs_img, simg).save(args.evecs)
 
         # save individual e-vectors also
         for i in range(3):
-            nib.save(nib.Nifti1Image(evecs[..., i], affine),
-                     add_filename_suffix(args.evecs, '_v'+str(i+1)))
+            ev_img = nib.Nifti1Image(evecs_data[..., i], affine)
+            StatefulImage.create_from(ev_img, simg).save(
+                add_filename_suffix(args.evecs, '_v'+str(i+1)))
 
     if args.evals:
-        evals = tenfit.evals.astype(np.float32)
-        nib.save(nib.Nifti1Image(evals, affine), args.evals)
+        evals_data = tenfit.evals.astype(np.float32)
+        evals_img = nib.Nifti1Image(evals_data, affine)
+        StatefulImage.create_from(evals_img, simg).save(args.evals)
 
         # save individual e-values also
         for i in range(3):
-            nib.save(nib.Nifti1Image(evals[..., i], affine),
-                     add_filename_suffix(args.evals, '_e' + str(i+1)))
+            eval_img = nib.Nifti1Image(evals_data[..., i], affine)
+            StatefulImage.create_from(eval_img, simg).save(
+                add_filename_suffix(args.evals, '_e' + str(i+1)))
 
     if args.p_i_signal:
         S0 = np.mean(data[..., gtab.b0s_mask], axis=-1, keepdims=True)
@@ -305,17 +327,18 @@ def main():
         if args.mask is not None:
             pis_mask *= mask
 
-        nib.save(nib.Nifti1Image(pis_mask.astype(np.int16), affine),
-                 args.p_i_signal)
+        pis_img = nib.Nifti1Image(pis_mask.astype(np.int16), affine)
+        StatefulImage.create_from(pis_img, simg).save(args.p_i_signal)
 
     if args.pulsation:
         STD = np.std(data[..., ~gtab.b0s_mask], axis=-1)
 
         if args.mask is not None:
             STD *= mask
 
-        nib.save(nib.Nifti1Image(STD.astype(np.float32), affine),
-                 add_filename_suffix(args.pulsation, '_std_dwi'))
+        std_img = nib.Nifti1Image(STD.astype(np.float32), affine)
+        StatefulImage.create_from(std_img, simg).save(
+            add_filename_suffix(args.pulsation, '_std_dwi'))
 
         if np.sum(gtab.b0s_mask) <= 1:
             logger.info('Not enough b=0 images to output standard '
@@ -330,8 +353,9 @@ def main():
             if args.mask is not None:
                 STD *= mask
 
-            nib.save(nib.Nifti1Image(STD.astype(np.float32), affine),
-                     add_filename_suffix(args.pulsation, '_std_b0'))
+            std_b0_img = nib.Nifti1Image(STD.astype(np.float32), affine)
+            StatefulImage.create_from(std_b0_img, simg).save(
+                add_filename_suffix(args.pulsation, '_std_b0'))
 
     if args.residual:
         if mask is None:
@@ -354,7 +378,8 @@ def main():
         R, data_diff = compute_residuals(
             predicted_data=tenfit2_predict.astype(np.float32),
             real_data=data, b0s_mask=gtab.b0s_mask, mask=mask)
-        nib.save(nib.Nifti1Image(R.astype(np.float32), affine), args.residual)
+        res_img = nib.Nifti1Image(R.astype(np.float32), affine)
+        StatefulImage.create_from(res_img, simg).save(args.residual)
 
         # Each volume's residual statistics
         R_k, q1, q3, iqr, std = compute_residuals_statistics(data_diff)

src/scilpy/cli/scil_fodf_metrics.py

@@ -40,6 +40,7 @@
 from dipy.direction.peaks import reshape_peaks_for_visualization
 
 from scilpy.io.image import get_data_as_mask
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.io.utils import (add_overwrite_arg, add_sh_basis_args,
                              add_processes_arg, add_verbose_arg,
                              assert_inputs_exist, assert_outputs_exist,
@@ -138,11 +139,14 @@ def main():
     assert_headers_compatible(parser, args.in_fODF, args.mask)
 
     # Loading
-    vol = nib.load(args.in_fODF)
-    data = vol.get_fdata(dtype=np.float32)
-    affine = vol.affine
-    mask = get_data_as_mask(nib.load(args.mask),
-                            dtype=bool) if args.mask else None
+    simg = StatefulImage.load(args.in_fODF)
+    data = simg.get_fdata(dtype=np.float32)
+    affine = simg.affine
+    mask = None
+    if args.mask:
+        mask_simg = StatefulImage.load(args.mask)
+        mask_simg.reorient(simg.axcodes)
+        mask = get_data_as_mask(mask_simg, dtype=bool)
 
     sphere = get_sphere(name=args.sphere)
     sh_basis, is_legacy = parse_sh_basis_arg(args)
@@ -168,26 +172,26 @@ def main():
 
         # Save result
         if args.nufo:
-            nib.save(nib.Nifti1Image(nufo_map.astype(np.float32), affine),
-                     args.nufo)
+            nufo_img = nib.Nifti1Image(nufo_map.astype(np.float32), affine)
+            StatefulImage.create_from(nufo_img, simg).save(args.nufo)
 
         if args.afd_max:
-            nib.save(nib.Nifti1Image(afd_max.astype(np.float32), affine),
-                     args.afd_max)
+            afd_max_img = nib.Nifti1Image(afd_max.astype(np.float32), affine)
+            StatefulImage.create_from(afd_max_img, simg).save(args.afd_max)
 
         if args.afd_total:
             # this is the analytical afd total
             afd_tot = data[:, :, :, 0]
-            nib.save(nib.Nifti1Image(afd_tot.astype(np.float32), affine),
-                     args.afd_total)
+            afd_tot_img = nib.Nifti1Image(afd_tot.astype(np.float32), affine)
+            StatefulImage.create_from(afd_tot_img, simg).save(args.afd_total)
 
         if args.afd_sum:
-            nib.save(nib.Nifti1Image(afd_sum.astype(np.float32), affine),
-                     args.afd_sum)
+            afd_sum_img = nib.Nifti1Image(afd_sum.astype(np.float32), affine)
+            StatefulImage.create_from(afd_sum_img, simg).save(args.afd_sum)
 
         if args.rgb:
-            nib.save(nib.Nifti1Image(rgb_map.astype('uint8'), affine),
-                     args.rgb)
+            rgb_img = nib.Nifti1Image(rgb_map.astype('uint8'), affine)
+            StatefulImage.create_from(rgb_img, simg).save(args.rgb)
 
     if args.peaks or args.peak_values:
         if not args.abs_peaks_and_values:
@@ -196,15 +200,19 @@ def main():
                                     where=peak_values[..., 0, None] != 0)
             peak_dirs[...] *= peak_values[..., :, None]
         if args.peaks:
-            nib.save(nib.Nifti1Image(
+            peaks_img = nib.Nifti1Image(
                 reshape_peaks_for_visualization(peak_dirs),
-                affine), args.peaks)
+                affine)
+            StatefulImage.create_from(peaks_img, simg).save(args.peaks)
         if args.peak_values:
-            nib.save(nib.Nifti1Image(peak_values, vol.affine),
-                     args.peak_values)
+            peak_vals_img = nib.Nifti1Image(peak_values, affine)
+            StatefulImage.create_from(peak_vals_img, simg).save(
+                args.peak_values)
 
     if args.peak_indices:
-        nib.save(nib.Nifti1Image(peak_indices, vol.affine), args.peak_indices)
+        peak_indices_img = nib.Nifti1Image(peak_indices, affine)
+        StatefulImage.create_from(peak_indices_img, simg).save(
+            args.peak_indices)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_fodf_msmt.py

@@ -22,15 +22,14 @@
 
 from dipy.core.gradients import gradient_table, unique_bvals_tolerance
 from dipy.data import get_sphere
-from dipy.io.gradients import read_bvals_bvecs
 from dipy.reconst.mcsd import MultiShellDeconvModel, multi_shell_fiber_response
-import nibabel as nib
 import numpy as np
 
 from scilpy.gradients.bvec_bval_tools import (check_b0_threshold,
                                               normalize_bvecs,
                                               is_normalized_bvecs)
 from scilpy.io.image import get_data_as_mask
+from scilpy.io.stateful_image import StatefulImage
 from scilpy.io.utils import (add_overwrite_arg, add_processes_arg,
                              assert_inputs_exist, assert_outputs_exist,
                              add_sh_basis_args, add_skip_b0_check_arg,
@@ -132,18 +131,32 @@ def main():
     wm_frf = np.loadtxt(args.in_wm_frf)
     gm_frf = np.loadtxt(args.in_gm_frf)
     csf_frf = np.loadtxt(args.in_csf_frf)
-    vol = nib.load(args.in_dwi)
-    data = vol.get_fdata(dtype=np.float32)
-    bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec)
+
+    simg = StatefulImage.load(args.in_dwi)
+    simg.load_gradients(args.in_bval, args.in_bvec)
+
+    # Orientation standardization?
+    # Reconstruction logic (dipy/scilpy) often prefers a specific orientation or consistency.
+    # We reorient secondary inputs to match the primary one.
+    # If we want to be fully robust, we could force RAS here, but let's see.
+    # scil_frf_msmt used to_ras(), so let's be consistent.
+    simg.to_ras()
+
+    data = simg.get_fdata(dtype=np.float32)
+    bvals = simg.bvals
+    bvecs = simg.bvecs
 
     # Checking data and sh_order
     wm_frf, gm_frf, csf_frf = verify_frf_files(wm_frf, gm_frf, csf_frf)
     verify_data_vs_sh_order(data, args.sh_order)
     sh_basis, is_legacy = parse_sh_basis_arg(args)
 
     # Checking mask
-    mask = get_data_as_mask(nib.load(args.mask),
-                            dtype=bool) if args.mask else None
+    mask = None
+    if args.mask:
+        mask_simg = StatefulImage.load(args.mask)
+        mask_simg.reorient(simg.axcodes)
+        mask = get_data_as_mask(mask_simg, dtype=bool)
 
     # Checking bvals, bvecs values and loading gtab
     if not is_normalized_bvecs(bvecs):
@@ -206,8 +219,8 @@ def main():
                                     is_input_legacy=True,
                                     is_output_legacy=is_legacy,
                                     nbr_processes=args.nbr_processes)
-        nib.save(nib.Nifti1Image(wm_coeff.astype(np.float32),
-                                 vol.affine), args.wm_out_fODF)
+        res_simg = StatefulImage.from_data(wm_coeff.astype(np.float32), simg)
+        res_simg.save(args.wm_out_fODF)
 
     if args.gm_out_fODF:
         gm_coeff = shm_coeff[..., 1]
@@ -218,8 +231,8 @@ def main():
                                     is_input_legacy=True,
                                     is_output_legacy=is_legacy,
                                     nbr_processes=args.nbr_processes)
-        nib.save(nib.Nifti1Image(gm_coeff.astype(np.float32),
-                                 vol.affine), args.gm_out_fODF)
+        res_simg = StatefulImage.from_data(gm_coeff.astype(np.float32), simg)
+        res_simg.save(args.gm_out_fODF)
 
     if args.csf_out_fODF:
         csf_coeff = shm_coeff[..., 0]
@@ -230,18 +243,18 @@ def main():
                                      is_input_legacy=True,
                                      is_output_legacy=is_legacy,
                                      nbr_processes=args.nbr_processes)
-        nib.save(nib.Nifti1Image(csf_coeff.astype(np.float32),
-                                 vol.affine), args.csf_out_fODF)
+        res_simg = StatefulImage.from_data(csf_coeff.astype(np.float32), simg)
+        res_simg.save(args.csf_out_fODF)
 
     if args.vf:
-        nib.save(nib.Nifti1Image(vf.astype(np.float32),
-                                 vol.affine), args.vf)
+        res_simg = StatefulImage.from_data(vf.astype(np.float32), simg)
+        res_simg.save(args.vf)
 
     if args.vf_rgb:
         vf_rgb = vf / np.max(vf) * 255
         vf_rgb = np.clip(vf_rgb, 0, 255)
-        nib.save(nib.Nifti1Image(vf_rgb.astype(np.uint8),
-                                 vol.affine), args.vf_rgb)
+        res_simg = StatefulImage.from_data(vf_rgb.astype(np.uint8), simg)
+        res_simg.save(args.vf_rgb)
 
 
 if __name__ == "__main__":