From 70ac214f89be19ec16153099aac6e3718de4f6ea Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Wed, 4 Oct 2023 10:35:23 -0400
Subject: [PATCH 01/10] Add option to use exact same seed every time

---
 scilpy/tracking/seed.py                    | 28 +++++++++++++---------
 scilpy/tracking/tracker.py                 | 19 +++++++++------
 scripts/scil_compute_local_tracking_dev.py |  9 ++++++-
 3 files changed, 37 insertions(+), 19 deletions(-)

diff --git a/scilpy/tracking/seed.py b/scilpy/tracking/seed.py
index b55a2f29a..d71434707 100644
--- a/scilpy/tracking/seed.py
+++ b/scilpy/tracking/seed.py
@@ -6,7 +6,7 @@
 from dipy.io.stateful_tractogram import Space, Origin
 
 
-class SeedGenerator(object):
+class SeedGenerator:
     """
     Class to get seeding positions.
 
@@ -19,7 +19,8 @@ class SeedGenerator(object):
     in the range x = [0, 3], y = [3, 6], z = [6, 9].
     """
     def __init__(self, data, voxres,
-                 space=Space('vox'), origin=Origin('center')):
+                 space=Space('vox'), origin=Origin('center'),
+                 randomize_positions=True):
         """
         Parameters
         ----------
@@ -28,8 +29,12 @@ def __init__(self, data, voxres,
             to find all voxels with values > 0, but will not be kept in memory.
         voxres: np.array(3,)
             The pixel resolution, ex, using img.header.get_zooms()[:3].
+        randomize_positions: bool
+            By default, seed position is moved randomly inside the voxel. Set to
+            False to have all seeds centered at the middle of the voxel.
         """
         self.voxres = voxres
+        self.randomize_positions = randomize_positions
 
         self.origin = origin
         self.space = space
@@ -71,16 +76,17 @@ def get_next_pos(self, random_generator, indices, which_seed):
         ind = which_seed % len_seeds
         x, y, z = self.seeds_vox[indices[ind]]
 
-        # Subvoxel initial positioning. Right now x, y, z are in vox space,
-        # origin=corner, so between 0 and 1.
-        r_x = random_generator.uniform(0, 1)
-        r_y = random_generator.uniform(0, 1)
-        r_z = random_generator.uniform(0, 1)
+        if self.randomize_positions:
+            # Subvoxel initial positioning. Right now x, y, z are in vox space,
+            # origin=corner, so between 0 and 1.
+            r_x = random_generator.uniform(0, 1)
+            r_y = random_generator.uniform(0, 1)
+            r_z = random_generator.uniform(0, 1)
 
-        # Moving inside the voxel
-        x += r_x
-        y += r_y
-        z += r_z
+            # Moving inside the voxel
+            x += r_x
+            y += r_y
+            z += r_z
 
         if self.origin == Origin('center'):
             # Bound [0, 0, 0] is now [-0.5, -0.5, -0.5]
diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py
index 32dbae453..87e3f6226 100644
--- a/scilpy/tracking/tracker.py
+++ b/scilpy/tracking/tracker.py
@@ -329,6 +329,18 @@ def _get_streamlines(self, chunk_id, lock=None):
             seed = self.seed_generator.get_next_pos(
                 random_generator, indices, first_seed_of_chunk + s)
 
+            # Setting the random value.
+            # Previous usage (and usage in Dipy) is to set the random seed
+            # based on the (real) seed position. However, in the case where we
+            # like to have exactly the same seed more than once, this will lead
+            # to exactly the same line, even in probabilistic tracking.
+            # Changing to seed position + seed number.
+            # toDo See numpy's doc: np.random.seed:
+            #  This is a convenience, legacy function.
+            #  The best practice is to not reseed a BitGenerator, rather to
+            #  recreate a new one. This method is here for legacy reasons.
+            np.random.seed(np.uint32(hash((seed + (s, s, s), self.rng_seed))))
+
             # Forward and backward tracking
             line = self._get_line_both_directions(seed)
 
@@ -375,13 +387,6 @@ def _get_line_both_directions(self, seeding_pos):
         line: list of 3D positions
             The generated streamline for seeding_pos.
         """
-
-        # toDo See numpy's doc: np.random.seed:
-        #  This is a convenience, legacy function.
-        #  The best practice is to not reseed a BitGenerator, rather to
-        #  recreate a new one. This method is here for legacy reasons.
-        np.random.seed(np.uint32(hash((seeding_pos, self.rng_seed))))
-
         # Forward
         line = [np.asarray(seeding_pos)]
         tracking_info = self.propagator.prepare_forward(seeding_pos)
diff --git a/scripts/scil_compute_local_tracking_dev.py b/scripts/scil_compute_local_tracking_dev.py
index 4157631a6..20cb048c3 100755
--- a/scripts/scil_compute_local_tracking_dev.py
+++ b/scripts/scil_compute_local_tracking_dev.py
@@ -120,6 +120,11 @@ def _build_arg_parser():
              "direction \n(based on the propagator's definition of invalid; "
              "ex when \nangle is too sharp of sh_threshold not reached) are "
              "never added.")
+    track_g.add_argument(
+        "--do_not_randomize_seed_positions", action="store_true",
+        help="By default, seed position is moved randomly inside the voxel. "
+             "Use this option to have all seeds centered at the middle of the "
+             "voxel.")
 
     add_seeding_options(p)
 
@@ -194,8 +199,10 @@ def main():
                       'seeding mask.'.format(args.in_seed))
 
     seed_res = seed_img.header.get_zooms()[:3]
+    randomize_positions = not args.do_not_randomize_seed_positions
     seed_generator = SeedGenerator(seed_data, seed_res,
-                                   space=our_space, origin=our_origin)
+                                   space=our_space, origin=our_origin,
+                                   randomize_positions=randomize_positions)
     if args.npv:
         # toDo. This will not really produce n seeds per voxel, only true
         #  in average.

From 85f2f60c4b55f88df5b784e103500f626a673e35 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Wed, 4 Oct 2023 10:39:21 -0400
Subject: [PATCH 02/10] Modify option discard_last_point to discard by default

---
 scripts/scil_compute_local_tracking_dev.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/scripts/scil_compute_local_tracking_dev.py b/scripts/scil_compute_local_tracking_dev.py
index 20cb048c3..8b828330b 100755
--- a/scripts/scil_compute_local_tracking_dev.py
+++ b/scripts/scil_compute_local_tracking_dev.py
@@ -113,13 +113,15 @@ def _build_arg_parser():
                          help="Mask interpolation: nearest-neighbor or "
                               "trilinear. [%(default)s]")
     track_g.add_argument(
-        '--discard_last_out_point', action='store_true',
-        help="If set, discard the last point (once out of the tracking mask) "
-             "of \nthe streamline. Default: append them. This is the default "
+        '--keep_last_out_point', action='store_true',
+        help="If set, keep the last point (once out of the tracking mask) "
+             "of \nthe streamline. Default: discard them. This is the default "
              " in \nDipy too. Note that points obtained after an invalid "
              "direction \n(based on the propagator's definition of invalid; "
              "ex when \nangle is too sharp of sh_threshold not reached) are "
-             "never added.")
+             "never added.\n"
+             "REMARK: Our results (our endpoints) seem to differ from dipy's. "
+             "This should be investigated.")
     track_g.add_argument(
         "--do_not_randomize_seed_positions", action="store_true",
         help="By default, seed position is moved randomly inside the voxel. "
@@ -243,7 +245,6 @@ def main():
         space=our_space, origin=our_origin)
 
     logging.debug("Instantiating tracker.")
-    append_last_point = not args.discard_last_out_point
     tracker = Tracker(propagator, mask, seed_generator, nbr_seeds, min_nbr_pts,
                       max_nbr_pts, args.max_invalid_nb_points,
                       compression_th=args.compress,
@@ -251,7 +252,8 @@ def main():
                       save_seeds=args.save_seeds,
                       mmap_mode='r+', rng_seed=args.rng_seed,
                       track_forward_only=args.forward_only,
-                      skip=args.skip, append_last_point=append_last_point,
+                      skip=args.skip,
+                      append_last_point=args.keep_last_out_point,
                       verbose=args.verbose)
 
     start = time.time()

From 9e4b1143f6543bfd0f6459339dda84b9b21785b1 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Wed, 4 Oct 2023 11:16:39 -0400
Subject: [PATCH 03/10] Changing 3 calls to uniform --> 1 call

---
 scilpy/tracking/seed.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/scilpy/tracking/seed.py b/scilpy/tracking/seed.py
index d71434707..99fdcd139 100644
--- a/scilpy/tracking/seed.py
+++ b/scilpy/tracking/seed.py
@@ -24,10 +24,10 @@ def __init__(self, data, voxres,
         """
         Parameters
         ----------
-        data: np.array
+        data: np.ndarray
             The data, ex, loaded from nibabel img.get_fdata(). It will be used
             to find all voxels with values > 0, but will not be kept in memory.
-        voxres: np.array(3,)
+        voxres: np.ndarray(3,)
             The pixel resolution, ex, using img.header.get_zooms()[:3].
         randomize_positions: bool
             By default, seed position is moved randomly inside the voxel. Set to
@@ -79,9 +79,7 @@ def get_next_pos(self, random_generator, indices, which_seed):
         if self.randomize_positions:
             # Subvoxel initial positioning. Right now x, y, z are in vox space,
             # origin=corner, so between 0 and 1.
-            r_x = random_generator.uniform(0, 1)
-            r_y = random_generator.uniform(0, 1)
-            r_z = random_generator.uniform(0, 1)
+            r_x, r_y, r_z = random_generator.uniform(0, 1, size=3)
 
             # Moving inside the voxel
             x += r_x

From 1196a5e100f81a3673e49156c430604af18836f3 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Wed, 4 Oct 2023 11:19:58 -0400
Subject: [PATCH 04/10] Include process number into the rng to avoid same value
 at each process

---
 scilpy/tracking/tracker.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py
index 87e3f6226..4119bd3cc 100644
--- a/scilpy/tracking/tracker.py
+++ b/scilpy/tracking/tracker.py
@@ -335,11 +335,15 @@ def _get_streamlines(self, chunk_id, lock=None):
             # like to have exactly the same seed more than once, this will lead
             # to exactly the same line, even in probabilistic tracking.
             # Changing to seed position + seed number.
+            # Then in the case of multiprocessing, adding also a fraction based
+            # on current process ID.
             # toDo See numpy's doc: np.random.seed:
             #  This is a convenience, legacy function.
             #  The best practice is to not reseed a BitGenerator, rather to
             #  recreate a new one. This method is here for legacy reasons.
-            np.random.seed(np.uint32(hash((seed + (s, s, s), self.rng_seed))))
+            eps = s + chunk_id / (self.nbr_processes + 1)
+            np.random.seed(np.uint32(hash((seed + (eps, eps, eps),
+                                           self.rng_seed))))
 
             # Forward and backward tracking
             line = self._get_line_both_directions(seed)

From 12abe8dc66b04c28b6f8175dbc8243f2847c9b6e Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Wed, 4 Oct 2023 11:41:17 -0400
Subject: [PATCH 05/10] Manage toDo: remove deprecated usage of seed

---
 scilpy/tracking/propagator.py | 16 ++++++++++++----
 scilpy/tracking/tools.py      |  6 ++++--
 scilpy/tracking/tracker.py    | 15 ++++++---------
 3 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/scilpy/tracking/propagator.py b/scilpy/tracking/propagator.py
index bd4ddddd4..a26c91608 100644
--- a/scilpy/tracking/propagator.py
+++ b/scilpy/tracking/propagator.py
@@ -67,6 +67,8 @@ def __init__(self, datavolume, step_size, rk_order, space, origin):
         # By default, normalizing directions. Adding option for child classes.
         self.normalize_directions = True
 
+        self.line_rng_generator = None   # Will be reset at each new streamline.
+
     def reset_data(self, new_data=None):
         """
         Reset data before starting a new process. In current implementation,
@@ -81,7 +83,7 @@ def reset_data(self, new_data=None):
         """
         self.datavolume.data = new_data
 
-    def prepare_forward(self, seeding_pos):
+    def prepare_forward(self, seeding_pos, random_generator):
         """
         Prepare information necessary at the first point of the
         streamline for forward propagation: v_in and any other information
@@ -92,6 +94,7 @@ def prepare_forward(self, seeding_pos):
         seeding_pos: tuple(x,y,z)
             The seeding position. Important, position must be in the same space
             and origin as self.space, self.origin!
+        random_generator: numpy Generator.
 
         Returns
         -------
@@ -100,6 +103,8 @@ def prepare_forward(self, seeding_pos):
             Return PropagationStatus.ERROR if no good tracking direction can be
             set at current seeding position.
         """
+        # To be defined by child classes.
+        # Should set self.line_rng_generator = random_generator
         raise NotImplementedError
 
     def prepare_backward(self, line, forward_dir):
@@ -421,7 +426,7 @@ def _get_sf(self, pos):
             sf /= sf_max
         return sf
 
-    def prepare_forward(self, seeding_pos):
+    def prepare_forward(self, seeding_pos, random_generator):
         """
         Prepare information necessary at the first point of the
         streamline for forward propagation: v_in and any other information
@@ -437,6 +442,7 @@ def prepare_forward(self, seeding_pos):
         seeding_pos: tuple(x,y,z)
             The seeding position. Important, position must be in the same space
             and origin as self.space, self.origin!
+        random_generator: numpy Generator
 
         Returns
         -------
@@ -452,9 +458,10 @@ def prepare_forward(self, seeding_pos):
         # "more probable" peak.
         sf = self._get_sf(seeding_pos)
         sf[sf < self.sf_threshold_init] = 0
+        self.line_rng_generator = random_generator
 
         if np.sum(sf) > 0:
-            ind = sample_distribution(sf)
+            ind = sample_distribution(sf, self.line_rng_generator)
             return TrackingDirection(self.dirs[ind], ind)
 
         # Else: sf at current position is smaller than acceptable threshold in
@@ -485,7 +492,8 @@ def _sample_next_direction(self, pos, v_in):
 
             # Sampling one.
             if np.sum(sf) > 0:
-                v_out = directions[sample_distribution(sf)]
+                v_out = directions[sample_distribution(sf,
+                                                       self.line_rng_generator)]
             else:
                 return None
         elif self.algo == 'det':
diff --git a/scilpy/tracking/tools.py b/scilpy/tracking/tools.py
index b6eae1d4b..90cddea6d 100644
--- a/scilpy/tracking/tools.py
+++ b/scilpy/tracking/tools.py
@@ -15,12 +15,13 @@ def get_theta(requested_theta, tracking_type):
     return theta
 
 
-def sample_distribution(dist):
+def sample_distribution(dist, random_generator: np.random.Generator):
     """
     Parameters
     ----------
     dist: numpy.array
         The empirical distribution to sample from.
+    random_generator: numpy Generator
 
     Return
     ------
@@ -30,4 +31,5 @@ def sample_distribution(dist):
     cdf = dist.cumsum()
     if cdf[-1] == 0:
         return None
-    return cdf.searchsorted(np.random.random() * cdf[-1])
+
+    return cdf.searchsorted(random_generator.random() * cdf[-1])
diff --git a/scilpy/tracking/tracker.py b/scilpy/tracking/tracker.py
index 4119bd3cc..d81a7c902 100644
--- a/scilpy/tracking/tracker.py
+++ b/scilpy/tracking/tracker.py
@@ -337,16 +337,12 @@ def _get_streamlines(self, chunk_id, lock=None):
             # Changing to seed position + seed number.
             # Then in the case of multiprocessing, adding also a fraction based
             # on current process ID.
-            # toDo See numpy's doc: np.random.seed:
-            #  This is a convenience, legacy function.
-            #  The best practice is to not reseed a BitGenerator, rather to
-            #  recreate a new one. This method is here for legacy reasons.
             eps = s + chunk_id / (self.nbr_processes + 1)
-            np.random.seed(np.uint32(hash((seed + (eps, eps, eps),
-                                           self.rng_seed))))
+            line_generator = np.random.default_rng(
+                np.uint32(hash((seed + (eps, eps, eps), self.rng_seed))))
 
             # Forward and backward tracking
-            line = self._get_line_both_directions(seed)
+            line = self._get_line_both_directions(seed, line_generator)
 
             if line is not None:
                 streamline = np.array(line, dtype='float32')
@@ -376,7 +372,7 @@ def _get_streamlines(self, chunk_id, lock=None):
                 p.close()
         return streamlines, seeds
 
-    def _get_line_both_directions(self, seeding_pos):
+    def _get_line_both_directions(self, seeding_pos, line_generator):
         """
         Generate a streamline from an initial position following the tracking
         parameters.
@@ -393,7 +389,8 @@ def _get_line_both_directions(self, seeding_pos):
         """
         # Forward
         line = [np.asarray(seeding_pos)]
-        tracking_info = self.propagator.prepare_forward(seeding_pos)
+        tracking_info = self.propagator.prepare_forward(seeding_pos,
+                                                        line_generator)
         if tracking_info == PropagationStatus.ERROR:
             # No good tracking direction can be found at seeding position.
             return None

From e27a5cdd917729b150fc3d379bd5f59f25d5da69 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Fri, 27 Oct 2023 09:14:50 -0400
Subject: [PATCH 06/10] Bring changes to get_n_pos too

---
 scilpy/tracking/seed.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/scilpy/tracking/seed.py b/scilpy/tracking/seed.py
index 99fdcd139..dc38a23c1 100644
--- a/scilpy/tracking/seed.py
+++ b/scilpy/tracking/seed.py
@@ -130,9 +130,14 @@ def get_next_n_pos(self, random_generator, indices, which_seeds):
         # Sub-voxel initial positioning
         # Prepare sub-voxel random movement now (faster out of loop)
         n = len(which_seeds)
-        r_x = random_generator.uniform(0, 1, size=n)
-        r_y = random_generator.uniform(0, 1, size=n)
-        r_z = random_generator.uniform(0, 1, size=n)
+        if self.randomize_positions:
+            r_x = random_generator.uniform(0, 1, size=n)
+            r_y = random_generator.uniform(0, 1, size=n)
+            r_z = random_generator.uniform(0, 1, size=n)
+        else:
+            r_x = np.zeros(n)
+            r_y = np.zeros(n)
+            r_z = np.zeros(n)
 
         seeds = []
         # Looping. toDo, see if can be done faster.

From 79141a84fe995c83a0c8727c31fa5650c083c9df Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Fri, 27 Oct 2023 09:16:49 -0400
Subject: [PATCH 07/10] Move toDo out of arparser help

---
 scripts/scil_compute_local_tracking_dev.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/scripts/scil_compute_local_tracking_dev.py b/scripts/scil_compute_local_tracking_dev.py
index 8b828330b..ab800d2be 100755
--- a/scripts/scil_compute_local_tracking_dev.py
+++ b/scripts/scil_compute_local_tracking_dev.py
@@ -119,9 +119,9 @@ def _build_arg_parser():
              " in \nDipy too. Note that points obtained after an invalid "
              "direction \n(based on the propagator's definition of invalid; "
              "ex when \nangle is too sharp of sh_threshold not reached) are "
-             "never added.\n"
-             "REMARK: Our results (our endpoints) seem to differ from dipy's. "
-             "This should be investigated.")
+             "never added.\n")
+    # ToDo Our results (our endpoints) seem to differ from dipy's, with or
+    #  witout option. This should be investigated.
     track_g.add_argument(
         "--do_not_randomize_seed_positions", action="store_true",
         help="By default, seed position is moved randomly inside the voxel. "

From e064d9c9e6f6a8e74223a672280b0b4b07adedf2 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Fri, 27 Oct 2023 09:18:30 -0400
Subject: [PATCH 08/10] Rename long param name

---
 scripts/scil_compute_local_tracking_dev.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/scripts/scil_compute_local_tracking_dev.py b/scripts/scil_compute_local_tracking_dev.py
index ab800d2be..9f11cf32f 100755
--- a/scripts/scil_compute_local_tracking_dev.py
+++ b/scripts/scil_compute_local_tracking_dev.py
@@ -123,8 +123,8 @@ def _build_arg_parser():
     # ToDo Our results (our endpoints) seem to differ from dipy's, with or
     #  witout option. This should be investigated.
     track_g.add_argument(
-        "--do_not_randomize_seed_positions", action="store_true",
-        help="By default, seed position is moved randomly inside the voxel. "
+        "--no_seed_displacement", action="store_true",
+        help="By default, seed position is moved randomly inside the voxel.\n"
              "Use this option to have all seeds centered at the middle of the "
              "voxel.")
 
@@ -201,7 +201,7 @@ def main():
                       'seeding mask.'.format(args.in_seed))
 
     seed_res = seed_img.header.get_zooms()[:3]
-    randomize_positions = not args.do_not_randomize_seed_positions
+    randomize_positions = not args.no_seed_displacement
     seed_generator = SeedGenerator(seed_data, seed_res,
                                    space=our_space, origin=our_origin,
                                    randomize_positions=randomize_positions)

From 79817ec37327ea8caddfe7851a1f82aee2bcea79 Mon Sep 17 00:00:00 2001
From: Charles Poirier <charles.poirier@usherbrooke.ca>
Date: Mon, 30 Oct 2023 15:34:19 -0400
Subject: [PATCH 09/10] Combined Emma + Charles: changing for option n_repeats.
 Adding doc. Adding unit test.

---
 scilpy/tracking/seed.py                    | 178 ++++++++++++++-------
 scilpy/tracking/tests/notes.txt            |   5 +
 scilpy/tracking/tests/test_seed.py         |  45 ++++++
 scripts/scil_compute_local_tracking_dev.py |  17 +-
 4 files changed, 175 insertions(+), 70 deletions(-)
 create mode 100644 scilpy/tracking/tests/notes.txt
 create mode 100644 scilpy/tracking/tests/test_seed.py

diff --git a/scilpy/tracking/seed.py b/scilpy/tracking/seed.py
index dc38a23c1..c7c0159ef 100644
--- a/scilpy/tracking/seed.py
+++ b/scilpy/tracking/seed.py
@@ -1,6 +1,4 @@
 # -*- coding: utf-8 -*-
-import logging
-
 import numpy as np
 
 from dipy.io.stateful_tractogram import Space, Origin
@@ -19,8 +17,7 @@ class SeedGenerator:
     in the range x = [0, 3], y = [3, 6], z = [6, 9].
     """
     def __init__(self, data, voxres,
-                 space=Space('vox'), origin=Origin('center'),
-                 randomize_positions=True):
+                 space=Space('vox'), origin=Origin('center'), n_repeats=1):
         """
         Parameters
         ----------
@@ -29,62 +26,79 @@ def __init__(self, data, voxres,
             to find all voxels with values > 0, but will not be kept in memory.
         voxres: np.ndarray(3,)
             The pixel resolution, ex, using img.header.get_zooms()[:3].
-        randomize_positions: bool
-            By default, seed position is moved randomly inside the voxel. Set to
-            False to have all seeds centered at the middle of the voxel.
+        n_repeats: int
+            Number of times a same seed position is returned.
+            If used, we supposed that calls to either get_next_pos or
+            get_next_n_pos are used sequentially. Not verified.
         """
         self.voxres = voxres
-        self.randomize_positions = randomize_positions
-
+        self.n_repeats = n_repeats
         self.origin = origin
         self.space = space
         if space == Space.RASMM:
             raise NotImplementedError("We do not support rasmm space.")
+        elif space not in [Space.VOX, Space.VOXMM]:
+            raise ValueError("Space should be a choice of Dipy Space.")
+        if origin not in [Origin.NIFTI, Origin.TRACKVIS]:
+            raise ValueError("Origin should be a choice of Dipy Origin.")
 
         # self.seed are all the voxels where a seed could be placed
         # (voxel space, origin=corner, int numbers).
-        self.seeds_vox = np.array(np.where(np.squeeze(data) > 0),
-                                  dtype=float).transpose()
+        self.seeds_vox_corner = np.array(np.where(np.squeeze(data) > 0),
+                                         dtype=float).transpose()
 
-        if len(self.seeds_vox) == 0:
-            logging.warning("There are no positive voxels in the seeding "
-                            "mask!")
+        if len(self.seeds_vox_corner) == 0:
+            raise ValueError("There are no positive voxels in the seeding "
+                             "mask!")
 
-    def get_next_pos(self, random_generator, indices, which_seed):
+        # We use this to remember last offset if n_repeats > 1:
+        self.previous_offset = None
+
+    def get_next_pos(self, random_generator, shuffled_indices, which_seed):
         """
         Generate the next seed position (Space=voxmm, origin=corner).
+        See self.init()_generator to get the generator and shuffled_indices.
+
+        To be used with self.n_repeats, we suppose that sequential
+        get_next_pos calls are used with sequentials values of which_seed.
 
         Parameters
         ----------
-        random_generator : numpy random generator
+        random_generator: numpy random generator
             Initialized numpy number generator.
-        indices : List
+        shuffled_indices: np.array
             Indices of current seeding map.
-        which_seed : int
+        which_seed: int
             Seed number to be processed.
+            (which_seed // self.n_repeats corresponds to the index of the
+            chosen seed in the flattened seeding mask).
 
         Return
         ------
         seed_pos: tuple
             Position of next seed expressed in mm.
         """
-        len_seeds = len(self.seeds_vox)
-        if len_seeds == 0:
-            return []
+        nb_seed_voxels = len(self.seeds_vox_corner)
 
         # Voxel selection from the seeding mask
-        ind = which_seed % len_seeds
-        x, y, z = self.seeds_vox[indices[ind]]
+        ind = (which_seed // self.n_repeats) % nb_seed_voxels
+        x, y, z = self.seeds_vox_corner[shuffled_indices[ind]]
 
-        if self.randomize_positions:
+        if which_seed % self.n_repeats == 0:
             # Subvoxel initial positioning. Right now x, y, z are in vox space,
             # origin=corner, so between 0 and 1.
             r_x, r_y, r_z = random_generator.uniform(0, 1, size=3)
+            self.previous_offset = (r_x, r_y, r_z)
+        else:
+            # Supposing that calls to get_next_pos are used correctly:
+            # previous_offset should already exist and correspond to the
+            # correct offset.
+            r_x, r_y, r_z = self.previous_offset
 
-            # Moving inside the voxel
-            x += r_x
-            y += r_y
-            z += r_z
+        # Moving inside the voxel
+        x += r_x
+        y += r_y
+        z += r_z
 
         if self.origin == Origin('center'):
             # Bound [0, 0, 0] is now [-0.5, -0.5, -0.5]
@@ -99,50 +113,81 @@ def get_next_pos(self, random_generator, indices, which_seed):
         else:
             raise NotImplementedError("We do not support rasmm space.")
 
-    def get_next_n_pos(self, random_generator, indices, which_seeds):
+    def get_next_n_pos(self, random_generator, shuffled_indices,
+                       which_seed_start, n):
         """
         Generate the next n seed positions. Intended for GPU usage.
+        Equivalent to:
+            for s in range(which_seed_start, which_seed_start + nb_seeds):
+                self.get_next_pos(..., s)
+        See description of get_next_pos for more information.
+
+        To be used with self.n_repeats, we suppose that sequential
+        get_next_n_pos calls are used with sequential values of
+        which_seed_start (with steps of nb_seeds).
 
         Parameters
         ----------
-        random_generator : numpy random generator
+        random_generator: numpy random generator
             Initialized numpy number generator.
-        indices : numpy array
+        shuffled_indices: np.array
             Indices of current seeding map.
-        which_seeds : numpy array
-            Seed numbers (i.e. IDs) to be processed.
+        which_seed_start: int
+            First seed numbers to be processed.
+            (which_seed_start // self.n_repeats corresponds to the index of the
+            chosen seed in the flattened seeding mask).
+        n: int
+            Number of seeds to get.
 
         Return
         ------
-        seed_pos: List[List]
+        seeds: List[List]
             Positions of next seeds expressed seed_generator's space and
             origin.
         """
-
-        len_seeds = len(self.seeds_vox)
-
-        if len_seeds == 0:
-            return []
+        nb_seed_voxels = len(self.seeds_vox_corner)
+        which_seeds = np.arange(which_seed_start, which_seed_start + n)
 
         # Voxel selection from the seeding mask
-        inds = which_seeds % len_seeds
+        # Same seed is re-used n_repeats times
+        inds = (which_seeds // self.n_repeats) % nb_seed_voxels
 
-        # Sub-voxel initial positioning
         # Prepare sub-voxel random movement now (faster out of loop)
-        n = len(which_seeds)
-        if self.randomize_positions:
-            r_x = random_generator.uniform(0, 1, size=n)
-            r_y = random_generator.uniform(0, 1, size=n)
-            r_z = random_generator.uniform(0, 1, size=n)
-        else:
-            r_x = np.zeros(n)
-            r_y = np.zeros(n)
-            r_z = np.zeros(n)
+        r_x = np.zeros((n,))
+        r_y = np.zeros((n,))
+        r_z = np.zeros((n,))
+
+        # Find where which_seeds % self.n_repeats == 0
+        # Note. If where_new_offsets[0] is False, supposing that calls to
+        # get_next_n_pos are used correctly: previous_offset should already
+        # exist and correspond to the correct offset.
+        where_new_offsets = ~(which_seeds % self.n_repeats).astype(bool)
+        ind_first = np.where(where_new_offsets)[0][0]
+        if not where_new_offsets[0]:
+            assert self.previous_offset is not None
+
+            # First continuing previous_offset.
+            r_x[0:ind_first] = self.previous_offset[0]
+            r_y[0:ind_first] = self.previous_offset[1]
+            r_z[0:ind_first] = self.previous_offset[2]
+
+        # Then starting and repeating new offsets.
+        nb_new_offsets = np.sum(where_new_offsets)
+        new_offsets_x = random_generator.uniform(0, 1, size=nb_new_offsets)
+        new_offsets_y = random_generator.uniform(0, 1, size=nb_new_offsets)
+        new_offsets_z = random_generator.uniform(0, 1, size=nb_new_offsets)
+        nb_r = n - ind_first
+        r_x[ind_first:] = np.repeat(new_offsets_x, self.n_repeats)[:nb_r]
+        r_y[ind_first:] = np.repeat(new_offsets_y, self.n_repeats)[:nb_r]
+        r_z[ind_first:] = np.repeat(new_offsets_z, self.n_repeats)[:nb_r]
+
+        # Save previous offset for next batch
+        self.previous_offset = (r_x[-1], r_y[-1], r_z[-1])
 
         seeds = []
         # Looping. toDo, see if can be done faster.
         for i in range(len(which_seeds)):
-            x, y, z = self.seeds_vox[indices[inds[i]]]
+            x, y, z = self.seeds_vox_corner[shuffled_indices[inds[i]]]
 
             # Moving inside the voxel
             x += r_x[i]
@@ -167,29 +212,40 @@ def get_next_n_pos(self, random_generator, indices, which_seeds):
 
         return seeds
 
-    def init_generator(self, random_initial_value, first_seed_of_chunk):
+    def init_generator(self, rng_seed, numbers_to_skip):
         """
-        Initialize numpy number generator according to user's parameter
-        and indexes from the seeding map.
+        Initialize a numpy number generator according to user's parameters.
+        Returns also the suffled index of all voxels.
+
+        The values are not stored in this classed, but are returned to the
+        user, who should add them as arguments in the methods
+        self.get_next_pos()
+        self.get_next_n_pos()
+        The use of this is that with multiprocessing, each process may have its
+        own generator, with less risk of using the wrong one when they are
+        managed by the user.
 
         Parameters
         ----------
-        random_initial_value : int
+        rng_seed : int
             The "seed" for the random generator.
-        first_seed_of_chunk : int
-            Number of seeds to skip (skip parameter + multi-processor skip).
+        numbers_to_skip : int
+            Number of seeds (i.e. voxels) to skip. Useful if you want to
+            continue sampling from the same generator as in a first experiment
+            (with a fixed rng_seed).
 
         Return
         ------
         random_generator : numpy random generator
             Initialized numpy number generator.
         indices : ndarray
-            Indices of current seeding map.
+            Shuffled indices of current seeding map, shuffled with current
+            generator.
         """
-        random_generator = np.random.RandomState(random_initial_value)
+        random_generator = np.random.RandomState(rng_seed)
 
         # 1. Initializing seeding maps indices (shuffling in-place)
-        indices = np.arange(len(self.seeds_vox))
+        indices = np.arange(len(self.seeds_vox_corner))
         random_generator.shuffle(indices)
 
         # 2. Initializing the random generator
@@ -198,7 +254,7 @@ def init_generator(self, random_initial_value, first_seed_of_chunk):
         # process (i.e this chunk)'s set of random numbers. Producing only
         # 100000 at the time to prevent RAM overuse.
         # (Multiplying by 3 for x,y,z)
-        random_numbers_to_skip = first_seed_of_chunk * 3
+        random_numbers_to_skip = numbers_to_skip * 3
         # toDo: see if 100000 is ok, and if we can create something not
         #  hard-coded
         while random_numbers_to_skip > 100000:
diff --git a/scilpy/tracking/tests/notes.txt b/scilpy/tracking/tests/notes.txt
new file mode 100644
index 000000000..0ecedff87
--- /dev/null
+++ b/scilpy/tracking/tests/notes.txt
@@ -0,0 +1,5 @@
+
+
+We will not test the tracker / propagator : too big to be tested, and only used
+in scil_compute_local_tracking_dev, which is intented for developping and
+testing new parameters.
\ No newline at end of file
diff --git a/scilpy/tracking/tests/test_seed.py b/scilpy/tracking/tests/test_seed.py
new file mode 100644
index 000000000..97436f1c8
--- /dev/null
+++ b/scilpy/tracking/tests/test_seed.py
@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+from dipy.io.stateful_tractogram import Space, Origin
+import numpy as np
+
+from scilpy.tracking.seed import SeedGenerator
+
+
+def test_seed_generation():
+
+    mask = np.zeros((5, 5, 5))
+
+    # Two seeds:
+    mask[1, 1, 1] = 1
+    mask[4, 3, 2] = 1
+
+    # With n_repeats: it is our role to get seeds sequentially.
+    generator = SeedGenerator(mask, voxres=[1, 1, 1], space=Space('vox'),
+                              origin=Origin('corner'), n_repeats=3)
+
+    rng_generator, shuffled_indices = generator.init_generator(
+        rng_seed=1, numbers_to_skip=0)
+
+    # Generating 3 seeds, should be the same
+    seed0 = generator.get_next_pos(rng_generator, shuffled_indices, 0)
+    seed1 = generator.get_next_pos(rng_generator, shuffled_indices, 1)
+    seed2 = generator.get_next_pos(rng_generator, shuffled_indices, 2)
+
+    assert np.array_equal(seed0, seed1)
+    assert np.array_equal(seed0, seed2)
+    assert np.array_equal(np.floor(seed0), [1, 1, 1])
+
+    # Generating more seeds, should be the same, in the other voxel
+    seed3 = generator.get_next_pos(rng_generator, shuffled_indices, 3)
+    seed4 = generator.get_next_pos(rng_generator, shuffled_indices, 4)
+    _ = generator.get_next_pos(rng_generator, shuffled_indices, 5)
+    assert np.array_equal(seed3, seed4)
+    assert np.array_equal(np.floor(seed3), [4, 3, 2])
+
+    # Generating n=4 seeds at once, back to voxel 1
+    seeds = generator.get_next_n_pos(rng_generator, shuffled_indices, 6, 4)
+    assert np.array_equal(seeds[0], seeds[1])
+    assert np.array_equal(seeds[0], seeds[2])
+    assert np.array_equal(np.floor(seeds[0]), [1, 1, 1])
+    assert np.array_equal(np.floor(seeds[3]), [4, 3, 2])
+
diff --git a/scripts/scil_compute_local_tracking_dev.py b/scripts/scil_compute_local_tracking_dev.py
index 9f11cf32f..ad3457c38 100755
--- a/scripts/scil_compute_local_tracking_dev.py
+++ b/scripts/scil_compute_local_tracking_dev.py
@@ -123,10 +123,10 @@ def _build_arg_parser():
     # ToDo Our results (our endpoints) seem to differ from dipy's, with or
     #  witout option. This should be investigated.
     track_g.add_argument(
-        "--no_seed_displacement", action="store_true",
-        help="By default, seed position is moved randomly inside the voxel.\n"
-             "Use this option to have all seeds centered at the middle of the "
-             "voxel.")
+        "--n_repeats_per_seed", type=int, default=1,
+        help="By default, each seed position is used only once. This option\n"
+             "allows for tracking from the exact same seed n_repeats_per_seed"
+             "\ntimes. [%(default)s]")
 
     add_seeding_options(p)
 
@@ -201,20 +201,19 @@ def main():
                       'seeding mask.'.format(args.in_seed))
 
     seed_res = seed_img.header.get_zooms()[:3]
-    randomize_positions = not args.no_seed_displacement
     seed_generator = SeedGenerator(seed_data, seed_res,
                                    space=our_space, origin=our_origin,
-                                   randomize_positions=randomize_positions)
+                                   n_repeats=args.n_repeats_per_seed)
     if args.npv:
         # toDo. This will not really produce n seeds per voxel, only true
         #  in average.
-        nbr_seeds = len(seed_generator.seeds_vox) * args.npv
+        nbr_seeds = len(seed_generator.seeds_vox_corner) * args.npv
     elif args.nt:
         nbr_seeds = args.nt
     else:
         # Setting npv = 1.
-        nbr_seeds = len(seed_generator.seeds_vox)
-    if len(seed_generator.seeds_vox) == 0:
+        nbr_seeds = len(seed_generator.seeds_vox_corner)
+    if len(seed_generator.seeds_vox_corner) == 0:
         parser.error('Seed mask "{}" does not have any voxel with value > 0.'
                      .format(args.in_seed))
 

From ae7d874ae6212e700e6614d88c0b196bf5bcaab2 Mon Sep 17 00:00:00 2001
From: EmmaRenauld <emmanuelle.renauld@usherbrooke.ca>
Date: Thu, 9 Nov 2023 08:53:31 -0500
Subject: [PATCH 10/10] Change notes.txt to dummpy test_propagator,
 test_tracker

---
 scilpy/tracking/tests/notes.txt          |  5 -----
 scilpy/tracking/tests/test_propagator.py | 10 ++++++++++
 scilpy/tracking/tests/test_tracker.py    | 10 ++++++++++
 3 files changed, 20 insertions(+), 5 deletions(-)
 delete mode 100644 scilpy/tracking/tests/notes.txt
 create mode 100644 scilpy/tracking/tests/test_propagator.py
 create mode 100644 scilpy/tracking/tests/test_tracker.py

diff --git a/scilpy/tracking/tests/notes.txt b/scilpy/tracking/tests/notes.txt
deleted file mode 100644
index 0ecedff87..000000000
--- a/scilpy/tracking/tests/notes.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-
-
-We will not test the tracker / propagator : too big to be tested, and only used
-in scil_compute_local_tracking_dev, which is intented for developping and
-testing new parameters.
\ No newline at end of file
diff --git a/scilpy/tracking/tests/test_propagator.py b/scilpy/tracking/tests/test_propagator.py
new file mode 100644
index 000000000..bba328ea0
--- /dev/null
+++ b/scilpy/tracking/tests/test_propagator.py
@@ -0,0 +1,10 @@
+# -*- coding: utf-8 -*-
+
+
+def test_class_propagator():
+    """
+    We will not test the tracker / propagator : too big to be tested, and only
+    used in scil_compute_local_tracking_dev, which is intented for developping
+    and testing new parameters.
+    """
+    pass
diff --git a/scilpy/tracking/tests/test_tracker.py b/scilpy/tracking/tests/test_tracker.py
new file mode 100644
index 000000000..b0a524842
--- /dev/null
+++ b/scilpy/tracking/tests/test_tracker.py
@@ -0,0 +1,10 @@
+# -*- coding: utf-8 -*-
+
+
+def test_class_tracker():
+    """
+    We will not test the tracker / propagator : too big to be tested, and only
+    used in scil_compute_local_tracking_dev, which is intented for developping
+    and testing new parameters.
+    """
+    pass