[ENH] Update to dipy 1.13 - DO NOT MERGE -

.python-version

@@ -1,2 +1,2 @@
 3.12
->=3.10, < 3.13
+>=3.11, < 3.13

pyproject.toml

@@ -9,11 +9,11 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "scilpy"
-version = "2.2.2"
+version = "2.3.0"
 description = "Scilpy: diffusion MRI tools and utilities"
 authors = [{ name = "SCIL Team" }]
 readme = "README.md"
-requires-python = ">=3.10, <3.13"
+requires-python = ">=3.11, <3.13"
 license-files = ["LICENSE"]
 classifiers = [
     "Development Status :: 3 - Alpha",
@@ -31,7 +31,7 @@ dependencies = [
 "cvxpy==1.6.*",
 "cycler==0.12.*",
 "deepdiff>=8.1,<8.7",
-"dipy==1.11.*",
+"dipy @ git+https://github.com/Stany87/dipy.git@e109cb5",
 "dmri-amico==2.1.*",
 "dmri-commit>=2.4.2,<2.5.0",
 "docopt==0.6.*",

src/scilpy/cli/scil_bundle_diameter.py

@@ -149,7 +149,7 @@ def _prepare_bundle_view(args, sft, data_labels, centroid_smooth, radius,
                                          opacity=args.opacity,
                                          colors=coloring)
 
-    slice_actor = actor.slicer(data_labels, np.eye(4))
+    slice_actor = actor.slicer(data_labels, affine=np.eye(4))
     slice_actor.opacity(0.0)
     return tube_actor, streamlines_actor, slice_actor
 

src/scilpy/cli/scil_denoising_nlmeans.py

@@ -78,6 +78,9 @@ def _build_arg_parser():
                    help="Path to a binary mask. Only the data inside the mask "
                         "will be denoised. If not provided, only non-zero "
                         "voxels will be denoised.")
+    p.add_argument('--algorithm',
+                   choices=['blockwise','classic'], default='blockwise',
+                   help='Algorithm to use for denoising. [%(default)s]')
     p.add_argument('--gaussian', action='store_true',
                    help="If you know that your data contains gaussian noise, "
                         "use this option. Otherwise, Rician is assumed.")
@@ -187,7 +190,10 @@ def main():
     if args.sigma is not None:
         logging.info('User supplied noise standard deviation is {}'
                      .format(args.sigma))
-        sigma = np.ones(vol_data.shape[:3]) * args.sigma
+        if nb_volumes > 1:
+            sigma = np.full(nb_volumes, args.sigma, dtype=np.float32)
+        else:
+            sigma = float(args.sigma)
     elif args.basic_sigma:
         if args.mask_sigma:
             mask_sigma = get_data_as_mask(nib.load(args.mask_sigma))
@@ -206,27 +212,26 @@ def main():
         sigma = np.median(sigma)  # Managing 4D data.
         logging.info('The median noise is: {}'.format(sigma))
 
-        # Broadcast the single value to a whole 3D volume for nlmeans
-        sigma = np.ones(vol_data.shape) * sigma
+        if nb_volumes > 1:
+            sigma = np.full(nb_volumes, sigma, dtype=np.float32)
     else:  # --piesno
         logging.info("Computing sigma: one value per slice.")
         sigma, mask_noise = estimate_piesno_sigma(vol_data, args.number_coils)
-
         if args.save_piesno_mask:
             logging.info("Saving resulting Piesno noise mask in {}"
                          .format(args.save_piesno_mask))
             nib.save(nib.Nifti1Image(mask_noise, vol.affine,
                                      header=vol.header),
                      args.save_piesno_mask)
 
-        # Broadcast the values per slice to a whole 3D volume for nlmeans
+        # Keep a 3D sigma map (one value per slice) for PIESNO.
         sigma = np.ones(vol_data.shape[:3]) * sigma[None, None, :]
 
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore", category=DeprecationWarning)
-        data_denoised = nlmeans(
-            vol_data, sigma, mask=mask_denoise, rician=not args.gaussian,
-            num_threads=args.nbr_processes)
+    data_denoised = nlmeans(vol_data, sigma, 
+                            mask=mask_denoise,
+                            rician=not args.gaussian,
+                            method=args.algorithm,
+                            num_threads=args.nbr_processes)
 
     # Saving
     nib.save(nib.Nifti1Image(data_denoised, vol.affine, header=vol.header),

src/scilpy/cli/scil_tracking_local.py

@@ -66,6 +66,7 @@
 from dipy.tracking import utils as track_utils
 from dipy.tracking.local_tracking import LocalTracking
 from dipy.tracking.stopping_criterion import BinaryStoppingCriterion
+from dipy.tracking.tracker import eudx_tracking
 from scilpy.io.image import get_data_as_mask
 from scilpy.io.utils import (add_sphere_arg, add_verbose_arg,
                              assert_headers_compatible, assert_inputs_exist,
@@ -235,23 +236,44 @@ def main():
         # LocalTracking.maxlen is actually the maximum length
         # per direction, we need to filter post-tracking.
         max_steps_per_direction = int(args.max_length / args.step_size)
+        stopping_criterion = BinaryStoppingCriterion(mask_data)
 
         logging.info("Starting CPU local tracking.")
-        streamlines_generator = LocalTracking(
-            get_direction_getter(
-                args.in_odf, args.algo, args.sphere,
-                args.sub_sphere, args.theta, sh_basis,
-                voxel_size, args.sf_threshold, args.sh_to_pmf,
-                args.probe_length, args.probe_radius,
-                args.probe_quality, args.probe_count,
-                args.support_exponent, is_legacy=is_legacy),
-            BinaryStoppingCriterion(mask_data),
-            seeds, np.eye(4),
-            step_size=vox_step_size, max_cross=1,
-            maxlen=max_steps_per_direction,
-            fixedstep=True, return_all=True,
-            random_seed=args.seed,
-            save_seeds=True)
+        if args.algo == 'eudx':
+            streamlines_generator = eudx_tracking(
+                seeds,
+                stopping_criterion,
+                np.eye(4),
+                pam=get_direction_getter(
+                    args.in_odf, args.algo, args.sphere,
+                    args.sub_sphere, args.theta, sh_basis,
+                    voxel_size, args.sf_threshold, args.sh_to_pmf,
+                    args.probe_length, args.probe_radius,
+                    args.probe_quality, args.probe_count,
+                    args.support_exponent, is_legacy=is_legacy),
+                max_cross=1,
+                max_len=max_steps_per_direction,
+                step_size=vox_step_size,
+                max_angle=get_theta(args.theta, args.algo),
+                random_seed=args.seed if args.seed is not None else 0,
+                return_all=True,
+                save_seeds=True)
+        else:
+            streamlines_generator = LocalTracking(
+                get_direction_getter(
+                    args.in_odf, args.algo, args.sphere,
+                    args.sub_sphere, args.theta, sh_basis,
+                    voxel_size, args.sf_threshold, args.sh_to_pmf,
+                    args.probe_length, args.probe_radius,
+                    args.probe_quality, args.probe_count,
+                    args.support_exponent, is_legacy=is_legacy),
+                stopping_criterion,
+                seeds, np.eye(4),
+                step_size=vox_step_size, max_cross=1,
+                maxlen=max_steps_per_direction,
+                fixedstep=True, return_all=True,
+                random_seed=args.seed,
+                save_seeds=True)
 
     else:  # GPU tracking
         # we'll make our streamlines twice as long,

src/scilpy/cli/scil_tractogram_cut_streamlines.py

@@ -193,7 +193,7 @@ def main():
         new_sft = StatefulTractogram.from_sft(
             compressed_strs, sft, data_per_streamline=sft.data_per_streamline)
 
-    save_tractogram(new_sft, args.out_tractogram, args.no_empty)
+    save_tractogram(new_sft, args.out_tractogram, bbox_valid_check=args.no_empty)
 
 
 if __name__ == "__main__":

src/scilpy/cli/scil_viz_bundle_screenshot_mosaic.py

@@ -259,7 +259,7 @@ def main():
             opacity = args.opacity_background
 
             # Structural data
-            slice_actor = actor.slicer(data, affine, value_range)
+            slice_actor = actor.slicer(data, affine=affine, value_range=value_range)
             slice_actor.opacity(opacity)
             ren.add(slice_actor)
 

src/scilpy/cli/scil_viz_tractogram_seeds_3d.py

@@ -97,7 +97,7 @@ def main():
     scene.background(tuple(map(int, args.background)))
 
     seedroi_actor = actor.contour_from_label(
-        seed_map_data, seed_map_affine, color=colors)
+        seed_map_data, affine=seed_map_affine, color=colors)
     scene.add(seedroi_actor)
 
     # Load tractogram as tubes or lines, with color if specified

src/scilpy/cli/tests/test_denoising_nlmeans.py

@@ -37,6 +37,17 @@ def test_execution_basic_3d(script_runner, monkeypatch):
     assert ret.success
 
 
+def test_execution_basic_3d_classic_algo(script_runner, monkeypatch):
+    monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
+    in_img = os.path.join(SCILPY_HOME, 'others', 't1_resample.nii.gz')
+    ret = script_runner.run(['scil_denoising_nlmeans', in_img,
+                             't1_denoised1.nii.gz', '--processes', '1',
+                             '--basic_sigma', '--number_coils', 0,
+                             '--algorithm', 'classic',
+                             '--gaussian'])
+    assert ret.success
+
+
 def test_execution_basic_4d_mask(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
     in_img = os.path.join(SCILPY_HOME, 'processing', 'dwi_crop_1000.nii.gz')

src/scilpy/cli/tests/test_fibertube_compute_density.py

@@ -27,8 +27,9 @@ def init_data():
     }
     mask_img = nib.Nifti1Image(mask, affine, header, extra)
 
-    sft_fibertubes = StatefulTractogram(streamlines, mask_img, Space.VOX,
-                                        Origin.NIFTI)
+    sft_fibertubes = StatefulTractogram(streamlines, mask_img,
+                                        Space.VOX,
+                                        origin=Origin.NIFTI)
     sft_fibertubes.data_per_streamline = {
         "diameters": np.array([0.2, 0.01])
     }

src/scilpy/cli/tests/test_tractogram_convert.py

@@ -19,10 +19,10 @@ def test_help_option(script_runner):
 
 def test_execution_surface_vtk_fib(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
-    in_fib = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
-                          'gyri_fanning.fib')
+    in_trk = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
+                          'gyri_fanning.trk')
     in_fa = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
                          'fa.nii.gz')
-    ret = script_runner.run(['scil_tractogram_convert', in_fib,
-                            'gyri_fanning.trk', '--reference', in_fa])
+    ret = script_runner.run(['scil_tractogram_convert', in_trk,
+                            'gyri_fanning converted.fib', '--reference', in_fa])
     assert ret.success

src/scilpy/cli/tests/test_tractogram_flip.py

@@ -19,10 +19,10 @@ def test_help_option(script_runner):
 
 def test_execution_surface_vtk_fib(script_runner, monkeypatch):
     monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
-    in_fib = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
-                          'gyri_fanning.fib')
+    in_trk = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
+                          'gyri_fanning.trk')
     in_fa = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
                          'fa.nii.gz')
-    ret = script_runner.run(['scil_tractogram_flip', in_fib,
-                             'gyri_fanning.tck', 'x', '--reference', in_fa])
+    ret = script_runner.run(['scil_tractogram_flip', in_trk,
+                             'gyri_fanning_fliped.tck', 'x', '--reference', in_fa])
     assert ret.success

src/scilpy/denoise/tests/test_asym_filtering.py

@@ -55,7 +55,9 @@ def test_cosine_filtering():
     sharpness = 1.0
 
     sh_order, full_basis = get_sh_order_and_fullness(in_sh.shape[-1])
-    out = cosine_filtering(in_sh, sh_order, sh_basis, full_basis, legacy,
-                           sharpness, sphere_str, sigma_spatial)
+    out = cosine_filtering(in_sh, sh_order=sh_order,
+                           sh_basis=sh_basis,
+                           in_full_basis=full_basis, is_legacy=legacy,
+                           dot_sharpness=sharpness, sphere_str=sphere_str, sigma=sigma_spatial)
 
     assert np.allclose(out, fodf_3x3_order8_descoteaux07_filtered_cosine)

src/scilpy/io/utils.py

@@ -175,13 +175,17 @@ def assert_gradients_filenames_valid(parser, filename_list, input_is_fsl):
 
 
 def add_json_args(parser):
-    g1 = parser.add_argument_group(title='Json options')
-    g1.add_argument('--indent',
-                    type=int, default=2,
-                    help='Indent for json pretty print.')
-    g1.add_argument('--sort_keys',
-                    action='store_true',
-                    help='Sort keys in output json.')
+    if isinstance(parser, argparse._ArgumentGroup):
+        target = parser
+    else:
+        target = parser.add_argument_group(title='Json options')
+
+    target.add_argument('--indent',
+                        type=int, default=2,
+                        help='Indent for json pretty print.')
+    target.add_argument('--sort_keys',
+                        action='store_true',
+                        help='Sort keys in output json.')
 
 
 def add_processes_arg(parser):

src/scilpy/reconst/bingham.py

@@ -173,7 +173,7 @@ def bingham_fit_sh(sh, max_lobes=5, abs_th=0.,
     shape = sh.shape
 
     sphere = get_sphere(name='symmetric724').subdivide(n=2)
-    B_mat = sh_to_sf_matrix(sphere, order,
+    B_mat = sh_to_sf_matrix(sphere, sh_order_max=order,
                             full_basis=full_basis,
                             return_inv=False)
 

src/scilpy/reconst/fodf.py

@@ -69,7 +69,8 @@ def get_ventricles_max_fodf(data, fa, md, zoom, sh_basis,
 
     order = find_order_from_nb_coeff(data)
     sphere = get_sphere(name='repulsion100')
-    b_matrix, _ = sh_to_sf_matrix(sphere, order, sh_basis, legacy=is_legacy)
+    b_matrix, _ = sh_to_sf_matrix(sphere, sh_order_max=order,
+                                  basis_type=sh_basis, legacy=is_legacy)
     out_mask = np.zeros(data.shape[:-1])
 
     if mask is None:

src/scilpy/reconst/sh.py

@@ -107,7 +107,8 @@ def compute_sh_coefficients(dwi, gradient_table,
         sphere = Sphere(xyz=bvecs)
 
     # Fit SH
-    sh = sf_to_sh(weights, sphere, sh_order, basis_type, smooth=smooth,
+    sh = sf_to_sh(weights, sphere, sh_order_max=sh_order,
+                  basis_type=basis_type, smooth=smooth,
                   legacy=is_legacy)
 
     # Apply mask
@@ -294,7 +295,7 @@ def peaks_from_sh(shm_coeff, sphere, mask=None, relative_peak_threshold=0.5,
     """
     sh_order = order_from_ncoef(shm_coeff.shape[-1],
                                 full_basis=full_basis)
-    B, _ = sh_to_sf_matrix(sphere=sphere, sh_order=sh_order,
+    B, _ = sh_to_sf_matrix(sphere=sphere, sh_order_max=sh_order,
                            basis_type=sh_basis_type,
                            full_basis=full_basis, legacy=is_legacy)
 
@@ -447,7 +448,7 @@ def maps_from_sh(shm_coeff, peak_values, peak_indices, sphere,
         nufo_map, afd_max, afd_sum, rgb_map, gfa, qa
     """
     sh_order = order_from_ncoef(shm_coeff.shape[-1])
-    B, _ = sh_to_sf_matrix(sphere=sphere, sh_order=sh_order,
+    B, _ = sh_to_sf_matrix(sphere=sphere, sh_order_max=sh_order,
                            basis_type=sh_basis_type)
 
     data_shape = shm_coeff.shape
@@ -613,10 +614,10 @@ def convert_sh_basis(shm_coeff, sphere, mask=None,
         return shm_coeff
 
     sh_order = order_from_ncoef(shm_coeff.shape[-1])
-    B_in, _ = sh_to_sf_matrix(sphere=sphere, sh_order=sh_order,
+    B_in, _ = sh_to_sf_matrix(sphere=sphere, sh_order_max=sh_order,
                               basis_type=input_basis,
                               legacy=is_input_legacy)
-    _, invB_out = sh_to_sf_matrix(sphere=sphere, sh_order=sh_order,
+    _, invB_out = sh_to_sf_matrix(sphere=sphere, sh_order_max=sh_order,
                                   basis_type=output_basis,
                                   legacy=is_output_legacy)
 
@@ -724,7 +725,8 @@ def convert_sh_to_sf(shm_coeff, sphere, mask=None, dtype="float32",
 
     sh_order = order_from_ncoef(shm_coeff.shape[-1],
                                 full_basis=input_full_basis)
-    B_in, _ = sh_to_sf_matrix(sphere, sh_order, basis_type=input_basis,
+    B_in, _ = sh_to_sf_matrix(sphere, sh_order_max=sh_order,
+                              basis_type=input_basis,
                               full_basis=input_full_basis,
                               legacy=is_input_legacy)
     B_in = B_in.astype(dtype)

src/scilpy/reconst/tests/test_fodf.py

@@ -30,7 +30,9 @@ def test_get_ventricles_max_fodf():
     # Reconstruct SF values same as in method.
     order = find_order_from_nb_coeff(fodf_3x3_order8_descoteaux07)
     sphere = get_sphere(name='repulsion100')
-    b_matrix, _ = sh_to_sf_matrix(sphere, order, sh_basis, legacy=True)
+    b_matrix, _ = sh_to_sf_matrix(sphere, sh_order_max=order,
+                                  basis_type=sh_basis,
+                                  legacy=True)
 
     sf1 = np.dot(fodf_3x3_order8_descoteaux07[1, 0, 0], b_matrix)
     sf2 = np.dot(fodf_3x3_order8_descoteaux07[1, 1, 0], b_matrix)
@@ -56,7 +58,8 @@ def test_get_ventricles_max_fodf_median():
     # Reconstruct SF values same as in method.
     order = find_order_from_nb_coeff(fodf_3x3_order8_descoteaux07)
     sphere = get_sphere(name='repulsion100')
-    b_matrix, _ = sh_to_sf_matrix(sphere, order, sh_basis, legacy=True)
+    b_matrix, _ = sh_to_sf_matrix(sphere, sh_order_max=order,
+                                  basis_type=sh_basis, legacy=True)
 
     sf1 = np.dot(fodf_3x3_order8_descoteaux07[1, 0, 0], b_matrix)
     sf2 = np.dot(fodf_3x3_order8_descoteaux07[1, 1, 0], b_matrix)

src/scilpy/tracking/propagator.py

@@ -402,7 +402,8 @@ def __init__(self, datavolume, step_size,
             get_sh_order_and_fullness(self.datavolume.nb_coeffs)
         self.basis = basis
         self.is_legacy = is_legacy
-        self.B = sh_to_sf_matrix(self.sphere, sh_order, self.basis,
+        self.B = sh_to_sf_matrix(self.sphere, sh_order_max=sh_order,
+                                 basis_type=self.basis,
                                  smooth=0.006, return_inv=False,
                                  full_basis=full_basis, legacy=self.is_legacy)
 

src/scilpy/tracking/tracker.py

@@ -696,7 +696,8 @@ def _track(self):
         cl_manager.add_input_buffer('vertices', self.sphere.vertices)
 
         sh_order = find_order_from_nb_coeff(self.sh)
-        B_mat = sh_to_sf_matrix(self.sphere, sh_order, self.sh_basis,
+        B_mat = sh_to_sf_matrix(self.sphere, sh_order_max=sh_order,
+                                basis_type=self.sh_basis,
                                 return_inv=False, legacy=self.is_legacy)
         cl_manager.add_input_buffer('b_matrix', B_mat)