diff --git a/.python-version b/.python-version
index 885016e20..4ff528a2f 100644
--- a/.python-version
+++ b/.python-version
@@ -1,2 +1,2 @@
 3.12
->=3.10, < 3.13
+>=3.11, < 3.13
diff --git a/pyproject.toml b/pyproject.toml
index 9e1e810f7..0cb5645ef 100644
--- a/pyproject.toml
+++ b/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.*",
diff --git a/src/scilpy/cli/scil_bundle_diameter.py b/src/scilpy/cli/scil_bundle_diameter.py
index b913ab55a..14d268891 100755
--- a/src/scilpy/cli/scil_bundle_diameter.py
+++ b/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
 
diff --git a/src/scilpy/cli/scil_denoising_nlmeans.py b/src/scilpy/cli/scil_denoising_nlmeans.py
index 49bc8417c..6a0bc0431 100755
--- a/src/scilpy/cli/scil_denoising_nlmeans.py
+++ b/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,12 +212,11 @@ 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))
@@ -219,14 +224,14 @@ def main():
                                      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),
diff --git a/src/scilpy/cli/scil_tracking_local.py b/src/scilpy/cli/scil_tracking_local.py
index fed1ab1d9..1d04da368 100755
--- a/src/scilpy/cli/scil_tracking_local.py
+++ b/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,
diff --git a/src/scilpy/cli/scil_tractogram_cut_streamlines.py b/src/scilpy/cli/scil_tractogram_cut_streamlines.py
index b127bde4a..df1c28b69 100755
--- a/src/scilpy/cli/scil_tractogram_cut_streamlines.py
+++ b/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__":
diff --git a/src/scilpy/cli/scil_viz_bundle_screenshot_mosaic.py b/src/scilpy/cli/scil_viz_bundle_screenshot_mosaic.py
index 70c77eecb..406017332 100755
--- a/src/scilpy/cli/scil_viz_bundle_screenshot_mosaic.py
+++ b/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)
 
diff --git a/src/scilpy/cli/scil_viz_tractogram_seeds_3d.py b/src/scilpy/cli/scil_viz_tractogram_seeds_3d.py
index c35ebfe11..781c3121f 100755
--- a/src/scilpy/cli/scil_viz_tractogram_seeds_3d.py
+++ b/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
diff --git a/src/scilpy/cli/tests/test_denoising_nlmeans.py b/src/scilpy/cli/tests/test_denoising_nlmeans.py
index 41cc2bc9f..e2c095a97 100644
--- a/src/scilpy/cli/tests/test_denoising_nlmeans.py
+++ b/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')
diff --git a/src/scilpy/cli/tests/test_fibertube_compute_density.py b/src/scilpy/cli/tests/test_fibertube_compute_density.py
index c2953fc2a..926c018e0 100644
--- a/src/scilpy/cli/tests/test_fibertube_compute_density.py
+++ b/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])
     }
diff --git a/src/scilpy/cli/tests/test_tractogram_convert.py b/src/scilpy/cli/tests/test_tractogram_convert.py
index 6805bc535..8796630e0 100644
--- a/src/scilpy/cli/tests/test_tractogram_convert.py
+++ b/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
diff --git a/src/scilpy/cli/tests/test_tractogram_flip.py b/src/scilpy/cli/tests/test_tractogram_flip.py
index 7fcfe7297..96e149fab 100644
--- a/src/scilpy/cli/tests/test_tractogram_flip.py
+++ b/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
diff --git a/src/scilpy/denoise/tests/test_asym_filtering.py b/src/scilpy/denoise/tests/test_asym_filtering.py
index b1405548f..9b1b628e0 100644
--- a/src/scilpy/denoise/tests/test_asym_filtering.py
+++ b/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)
diff --git a/src/scilpy/io/utils.py b/src/scilpy/io/utils.py
index 68056bcfd..71e6e8a19 100644
--- a/src/scilpy/io/utils.py
+++ b/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):
diff --git a/src/scilpy/reconst/bingham.py b/src/scilpy/reconst/bingham.py
index e27f50038..36c518b2b 100644
--- a/src/scilpy/reconst/bingham.py
+++ b/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)
 
diff --git a/src/scilpy/reconst/fodf.py b/src/scilpy/reconst/fodf.py
index 55d2edc4c..ed8ceb7c5 100644
--- a/src/scilpy/reconst/fodf.py
+++ b/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:
diff --git a/src/scilpy/reconst/sh.py b/src/scilpy/reconst/sh.py
index c75f0f61f..1190fdf7c 100644
--- a/src/scilpy/reconst/sh.py
+++ b/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)
diff --git a/src/scilpy/reconst/tests/test_fodf.py b/src/scilpy/reconst/tests/test_fodf.py
index ab957c8fd..89911702d 100644
--- a/src/scilpy/reconst/tests/test_fodf.py
+++ b/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)
diff --git a/src/scilpy/tracking/propagator.py b/src/scilpy/tracking/propagator.py
index 777f46c53..32f7e6b5d 100644
--- a/src/scilpy/tracking/propagator.py
+++ b/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)
 
diff --git a/src/scilpy/tracking/tracker.py b/src/scilpy/tracking/tracker.py
index aee5f2a6a..eadfc9a20 100644
--- a/src/scilpy/tracking/tracker.py
+++ b/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)
 
diff --git a/src/scilpy/tracking/utils.py b/src/scilpy/tracking/utils.py
index 551d80959..99821d872 100644
--- a/src/scilpy/tracking/utils.py
+++ b/src/scilpy/tracking/utils.py
@@ -395,8 +395,8 @@ def get_direction_getter(in_img, algo, sphere, sub_sphere, theta, sh_basis,
             peak_indices = np.full((img_shape_3d + (npeaks,)), -1,
                                    dtype='int')
             b_matrix, _ = sh_to_sf_matrix(sphere,
-                                          find_order_from_nb_coeff(img_data),
-                                          sh_basis, legacy=is_legacy)
+                                          sh_order_max=find_order_from_nb_coeff(img_data),
+                                          basis_type=sh_basis, legacy=is_legacy)
 
             for idx in np.argwhere(np.sum(img_data, axis=-1)):
                 idx = tuple(idx)
diff --git a/src/scilpy/tractanalysis/afd_along_streamlines.py b/src/scilpy/tractanalysis/afd_along_streamlines.py
index d3638ab92..555509230 100644
--- a/src/scilpy/tractanalysis/afd_along_streamlines.py
+++ b/src/scilpy/tractanalysis/afd_along_streamlines.py
@@ -91,7 +91,8 @@ def _afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
     fodf_data = fodf.get_fdata(dtype=np.float32)
     order = find_order_from_nb_coeff(fodf_data)
     sphere = get_sphere(name='repulsion724')
-    b_matrix, _ = sh_to_sf_matrix(sphere, order, fodf_basis, legacy=is_legacy)
+    b_matrix, _ = sh_to_sf_matrix(sphere, sh_order_max=order,
+                                  basis_type=fodf_basis, legacy=is_legacy)
     _, n = sph_harm_ind_list(order)
     legendre0_at_n = legendre_p_all(order, 0)[0][n]
     sphere_norm = np.linalg.norm(sphere.vertices)
diff --git a/src/scilpy/tractanalysis/reproducibility_measures.py b/src/scilpy/tractanalysis/reproducibility_measures.py
index 6a1790390..c299db683 100755
--- a/src/scilpy/tractanalysis/reproducibility_measures.py
+++ b/src/scilpy/tractanalysis/reproducibility_measures.py
@@ -684,7 +684,8 @@ def tractogram_pairwise_comparison(sft_one, sft_two, mask, nbr_cpu=1,
     sft_2.to_center()
 
     global B
-    B, _ = sh_to_sf_matrix(get_sphere(name='repulsion724'), 8, 'descoteaux07')
+    B, _ = sh_to_sf_matrix(get_sphere(name='repulsion724'), sh_order_max=8,
+                           basis_type='descoteaux07')
 
     diff_data, acc_data = _compare_tractogram_wrapper(
         mask, nbr_cpu, skip_streamlines_distance)
diff --git a/src/scilpy/tractanalysis/streamlines_metrics.pyx b/src/scilpy/tractanalysis/streamlines_metrics.pyx
index f24c15422..3d331d4d0 100644
--- a/src/scilpy/tractanalysis/streamlines_metrics.pyx
+++ b/src/scilpy/tractanalysis/streamlines_metrics.pyx
@@ -46,12 +46,12 @@ def compute_tract_counts_map(streamlines, vol_dims):
     # This array counts the number of different tracks going through each voxel.
     # Need to keep both the array and the memview on it to be able to
     # reshape and return in the end.
-    traversal_tags = np.zeros((n_voxels,), dtype=int)
-    cdef np.int_t[:] traversal_tags_v = traversal_tags
+    traversal_tags = np.zeros((n_voxels,), dtype=np.intp)
+    cdef np.intp_t[:] traversal_tags_v = traversal_tags
 
     # This array keeps track of whether the current track has already been
     # flagged in a specific voxel.
-    cdef np.int_t[:] touched_tags_v = np.zeros((n_voxels,), dtype=int)
+    cdef np.intp_t[:] touched_tags_v = np.zeros((n_voxels,), dtype=np.intp)
 
     cdef int streamlines_len = len(streamlines)
 
@@ -71,7 +71,7 @@ def compute_tract_counts_map(streamlines, vol_dims):
     cdef np.double_t[:] cur_edge = np.zeros(3, dtype=np.double)
 
     # Memview for the coordinates of the current voxel
-    cdef np.int_t[:] cur_voxel_coords = np.zeros(3, dtype=int)
+    cdef np.intp_t[:] cur_voxel_coords = np.zeros(3, dtype=np.intp)
 
     # various temporary loop and working variables
     #cdef:
diff --git a/src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py b/src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py
index f4d7d32ce..3896f219c 100644
--- a/src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py
+++ b/src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py
@@ -38,3 +38,18 @@ def test_subdivide():
 
     # streamline 3 should be split into 11 (12 points)
     assert split[2].shape[0] == 12
+
+
+def test_subdivide_float64_streamlines():
+    streamline = np.asarray([[0.0, 0.0, 0.0],
+                             [1.2, 0.5, 0.0]], dtype=np.float64)
+
+    fake_ref = nib.Nifti1Image(np.zeros((3, 3, 3)), affine=np.eye(4))
+    sft = StatefulTractogram([streamline], reference=fake_ref,
+                             origin=Origin('corner'), space=Space('vox'))
+
+    split = subdivide_streamlines_at_voxel_faces(sft.streamlines)
+
+    assert split[0].dtype == np.float32
+    assert np.allclose(split[0][0], streamline[0])
+    assert np.allclose(split[0][-1], streamline[-1])
diff --git a/src/scilpy/tractanalysis/voxel_boundary_intersection.pyx b/src/scilpy/tractanalysis/voxel_boundary_intersection.pyx
index fd8f72134..70cf12510 100644
--- a/src/scilpy/tractanalysis/voxel_boundary_intersection.pyx
+++ b/src/scilpy/tractanalysis/voxel_boundary_intersection.pyx
@@ -45,13 +45,16 @@ def subdivide_streamlines_at_voxel_faces(streamlines):
         Updated streamline coordinates with added coordinate points
         at voxel boundaries.
     """
+    streamlines_data = np.ascontiguousarray(streamlines._data,
+                                            dtype=np.float32)
+
     cdef:
         cnp.npy_intp nb_streamlines = len(streamlines._lengths)
         cnp.npy_intp at_point = 0
 
         # Multiplying by 6 is simply a heuristic to avoiding resizing too many
         # times. In my bundles tests, I had either 0 or 1 resize.
-        cnp.npy_intp max_points = (streamlines.get_data().size // 6) * 12
+        cnp.npy_intp max_points = (streamlines_data.size // 6) * 12
 
     new_array_sequence = nib.streamlines.array_sequence.ArraySequence()
     new_array_sequence._lengths.resize(nb_streamlines)
@@ -61,7 +64,7 @@ def subdivide_streamlines_at_voxel_faces(streamlines):
     cdef:
         cnp.npy_intp[:] lengths_view_in = streamlines._lengths
         cnp.npy_intp[:] offsets_view_in = streamlines._offsets
-        float[:, :] data_view_in = streamlines._data
+        float[:, :] data_view_in = streamlines_data
         cnp.npy_intp[:] lengths_view_out = new_array_sequence._lengths
         cnp.npy_intp[:] offsets_view_out = new_array_sequence._offsets
         cnp.float32_t[:] data_view_out = new_array_sequence._data
diff --git a/src/scilpy/viz/slice.py b/src/scilpy/viz/slice.py
index a53bb9d3e..c8977d2c0 100644
--- a/src/scilpy/viz/slice.py
+++ b/src/scilpy/viz/slice.py
@@ -250,18 +250,23 @@ def create_odf_slicer(sh_fodf, orientation, slice_index, sphere, sh_order,
     if nb_subdivide is not None:
         sphere = sphere.subdivide(n=nb_subdivide)
 
-    fodf = sh_to_sf(sh_fodf, sphere, sh_order, sh_basis,
+    fodf = sh_to_sf(sh_fodf, sphere, 
+                    sh_order_max=sh_order,
+                    basis_type=sh_basis,
                     full_basis=full_basis, legacy=is_legacy)
 
     fodf_var = None
     B_mat = None
     if sh_variance is not None:
-        fodf_var = sh_to_sf(sh_variance, sphere, sh_order, sh_basis,
+        fodf_var = sh_to_sf(sh_variance, sphere,
+                            sh_order_max=sh_order,
+                            basis_type=sh_basis,
                             full_basis=full_basis, legacy=is_legacy)
     else:
         fodf = sh_fodf
-        B_mat = sh_to_sf_matrix(sphere, sh_order, sh_basis,
-                                full_basis, return_inv=False)
+        B_mat = sh_to_sf_matrix(sphere, sh_order_max=sh_order,
+                                basis_type=sh_basis,
+                                full_basis=full_basis, return_inv=False)
 
     odf_actor, var_actor = create_odf_actors(fodf, sphere, scale, fodf_var,
                                              mask, radial_scale,