Fix random seed in tracking

dipy/tracking/local/localtracking.py

@@ -1,3 +1,5 @@
+import random
+
 import numpy as np
 
 from dipy.tracking.local.localtrack import local_tracker, pft_tracker
@@ -36,7 +38,7 @@ def _get_voxel_size(affine):
 
     def __init__(self, direction_getter, tissue_classifier, seeds, affine,
                  step_size, max_cross=None, maxlen=500, fixedstep=True,
-                 return_all=True):
+                 return_all=True, random_seed=None):
         """Creates streamlines by using local fiber-tracking.
 
         Parameters
@@ -69,6 +71,9 @@ def __init__(self, direction_getter, tissue_classifier, seeds, affine,
         return_all : bool
             If true, return all generated streamlines, otherwise only
             streamlines reaching end points or exiting the image.
+        random_seed : int
+            The seed for the random seed generator (numpy.random.seed and
+            random.seed).
         """
 
         self.direction_getter = direction_getter
@@ -88,6 +93,7 @@ def __init__(self, direction_getter, tissue_classifier, seeds, affine,
         self.max_cross = max_cross
         self.max_length = maxlen
         self.return_all = return_all
+        self.random_seed = random_seed
 
     def _tracker(self, seed, first_step, streamline):
         return local_tracker(self.direction_getter,
@@ -116,6 +122,12 @@ def _generate_streamlines(self):
         B = F.copy()
         for s in self.seeds:
             s = np.dot(lin, s) + offset
+            # Set the random seed in numpy and random
+            if self.random_seed is not None:
+                s_random_seed = hash(np.abs((np.sum(s)) + self.random_seed)) \
+                    % (2**32 - 1)
+                random.seed(s_random_seed)
+                np.random.seed(s_random_seed)
             directions = self.direction_getter.initial_direction(s)
             if directions.size == 0 and self.return_all:
                 # only the seed position
@@ -146,7 +158,8 @@ class ParticleFilteringTracking(LocalTracking):
     def __init__(self, direction_getter, tissue_classifier, seeds, affine,
                  step_size, max_cross=None, maxlen=500,
                  pft_back_tracking_dist=2, pft_front_tracking_dist=1,
-                 pft_max_trial=20, particle_count=15, return_all=True):
+                 pft_max_trial=20, particle_count=15, return_all=True,
+                 random_seed=None):
         r"""A streamline generator using the particle filtering tractography
         method [1]_.
 
@@ -192,6 +205,9 @@ def __init__(self, direction_getter, tissue_classifier, seeds, affine,
         return_all : bool
             If true, return all generated streamlines, otherwise only
             streamlines reaching end points or exiting the image.
+        random_seed : int
+            The seed for the random seed generator (numpy.random.seed and
+            random.seed).
 
         References
         ----------
@@ -239,7 +255,8 @@ def __init__(self, direction_getter, tissue_classifier, seeds, affine,
                                                         max_cross,
                                                         maxlen,
                                                         True,
-                                                        return_all)
+                                                        return_all,
+                                                        random_seed)
 
     def _tracker(self, seed, first_step, streamline):
         return pft_tracker(self.direction_getter,

dipy/tracking/local/tests/test_tracking.py

@@ -241,6 +241,15 @@ def allclose(x, y):
     # number of seeds places
     npt.assert_(np.array([len(streamlines) == len(seeds)]))
 
+    # Test reproducibility
+    tracking_1 = Streamlines(LocalTracking(dg, tc, seeds, np.eye(4),
+                                           0.5,
+                                           random_seed=0)).data
+    tracking_2 = Streamlines(LocalTracking(dg, tc, seeds, np.eye(4),
+                                           0.5,
+                                           random_seed=0)).data
+    npt.assert_equal(tracking_1, tracking_2)
+
 
 def test_particle_filtering_tractography():
     """This tests that the ParticleFilteringTracking produces
@@ -374,6 +383,15 @@ def test_particle_filtering_tractography():
         lambda: ParticleFilteringTracking(dg, tc, seeds, np.eye(4), step_size,
                                           particle_count=-1))
 
+    # Test reproducibility
+    tracking_1 = Streamlines(ParticleFilteringTracking(dg, tc, seeds, np.eye(4),
+                                                       step_size,
+                                                       random_seed=0)).data
+    tracking_2 = Streamlines(ParticleFilteringTracking(dg, tc, seeds, np.eye(4),
+                                                       step_size,
+                                                       random_seed=0)).data
+    npt.assert_equal(tracking_1, tracking_2)
+
 
 def test_maximum_deterministic_tracker():
     """This tests that the Maximum Deterministic Direction Getter plays nice

dipy/tracking/tests/test_utils.py

@@ -537,6 +537,24 @@ def test_random_seeds_from_mask():
     assert_equal(100, len(seeds))
     assert_true(np.all((seeds > 1.5) & (seeds < 2.5)))
 
+    mask = np.zeros((15, 15, 15))
+    mask[2:14, 2:14, 2:14] = 1
+    seeds_npv_2 = random_seeds_from_mask(mask, seeds_count=2,
+                                         seed_count_per_voxel=True,
+                                         random_seed=0)[:150]
+    seeds_npv_3 = random_seeds_from_mask(mask, seeds_count=3,
+                                         seed_count_per_voxel=True,
+                                         random_seed=0)[:150]
+    assert_true(np.all(seeds_npv_2 == seeds_npv_3))
+
+    seeds_nt_150 = random_seeds_from_mask(mask, seeds_count=150,
+                                          seed_count_per_voxel=False,
+                                          random_seed=0)[:150]
+    seeds_nt_500 = random_seeds_from_mask(mask, seeds_count=500,
+                                          seed_count_per_voxel=False,
+                                          random_seed=0)[:150]
+    assert_true(np.all(seeds_nt_150 == seeds_nt_500))
+
 
 def test_connectivity_matrix_shape():
     # Labels: z-planes have labels 0,1,2

dipy/tracking/utils.py

@@ -413,16 +413,15 @@ def seeds_from_mask(mask, density=[1, 1, 1], voxel_size=None, affine=None):
 
 
 def random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True,
-                           affine=None):
+                           affine=None, random_seed=None):
     """Creates randomly placed seeds for fiber tracking from a binary mask.
 
     Seeds points are placed randomly distributed in voxels of ``mask``
     which are ``True``.
     If ``seed_count_per_voxel`` is ``True``, this function is
     similar to ``seeds_from_mask()``, with the difference that instead of
     evenly distributing the seeds, it randomly places the seeds within the
-    voxels specified by the ``mask``. The initial random conditions can be set
-    using ``numpy.random.seed(...)``, prior to calling this function.
+    voxels specified by the ``mask``.
 
     Parameters
     ----------
@@ -439,6 +438,8 @@ def random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True,
         The mapping between voxel indices and the point space for seeds. A
         seed point at the center the voxel ``[i, j, k]`` will be represented as
         ``[x, y, z]`` where ``[x, y, z, 1] == np.dot(affine, [i, j, k , 1])``.
+    random_seed : int
+        The seed for the random seed generator (numpy.random.seed).
 
     See Also
     --------
@@ -453,28 +454,37 @@ def random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True,
     --------
     >>> mask = np.zeros((3,3,3), 'bool')
     >>> mask[0,0,0] = 1
-    >>> np.random.seed(1)
-    >>> random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True)
-    array([[-0.082978  ,  0.22032449, -0.49988563]])
-    >>> random_seeds_from_mask(mask, seeds_count=6, seed_count_per_voxel=True)
-    array([[-0.19766743, -0.35324411, -0.40766141],
-           [-0.31373979, -0.15443927, -0.10323253],
-           [ 0.03881673, -0.08080549,  0.1852195 ],
-           [-0.29554775,  0.37811744, -0.47261241],
-           [ 0.17046751, -0.0826952 ,  0.05868983],
-           [-0.35961306, -0.30189851,  0.30074457]])
+    >>> random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True,
+    ... random_seed=1)
+    array([[-0.0640051 , -0.47407377,  0.04966248]])
+    >>> random_seeds_from_mask(mask, seeds_count=6, seed_count_per_voxel=True,
+    ... random_seed=1)
+    array([[-0.0640051 , -0.47407377,  0.04966248],
+           [ 0.0507979 ,  0.20814782, -0.20909526],
+           [ 0.46702984,  0.04723225,  0.47268436],
+           [-0.27800683,  0.37073231, -0.29328084],
+           [ 0.39286015, -0.16802019,  0.32122912],
+           [-0.42369171,  0.27991879, -0.06159077]])
     >>> mask[0,1,2] = 1
-    >>> random_seeds_from_mask(mask, seeds_count=2, seed_count_per_voxel=True)
-    array([[ 0.46826158, -0.18657582,  0.19232262],
-           [ 0.37638915,  0.39460666, -0.41495579],
-           [-0.46094522,  0.66983042,  2.3781425 ],
-           [-0.40165317,  0.92110763,  2.45788953]])
+    >>> random_seeds_from_mask(mask, seeds_count=2, seed_count_per_voxel=True,
+    ... random_seed=1)
+    array([[-0.0640051 , -0.47407377,  0.04966248],
+           [-0.27800683,  1.37073231,  1.70671916],
+           [ 0.0507979 ,  0.20814782, -0.20909526],
+           [-0.48962585,  1.00187459,  1.99577329]])
     """
     mask = np.array(mask, dtype=bool, copy=False, ndmin=3)
     if mask.ndim != 3:
         raise ValueError('mask cannot be more than 3d')
 
-    where = np.argwhere(mask)
+    # Randomize the voxels
+    np.random.seed(random_seed)
+    shape = mask.shape
+    mask = mask.flatten()
+    indices = np.arange(len(mask))
+    np.random.shuffle(indices)
+
+    where = [np.unravel_index(i, shape) for i in indices if mask[i] == 1]
     num_voxels = len(where)
 
     if not seed_count_per_voxel:
@@ -483,16 +493,23 @@ def random_seeds_from_mask(mask, seeds_count=1, seed_count_per_voxel=True,
     else:
         seeds_per_voxel = seeds_count
 
-    # Generate as many random triplets as the number of seeds needed
-    grid = np.random.random([seeds_per_voxel * num_voxels, 3])
-    # Repeat elements of 'where' so that it can be added to grid
-    where = np.repeat(where, seeds_per_voxel, axis=0)
-    seeds = where + grid - .5
+    seeds = []
+    for i in range(1, seeds_per_voxel + 1):
+        for s in where:
+            # Set the random seed with the current seed, the current value of
+            # seeds per voxel and the global random seed.
+            if random_seed is not None:
+                s_random_seed = hash((np.sum(s) + 1) * i + random_seed) \
+                           % (2**32 - 1)
+                np.random.seed(s_random_seed)
+            # Generate random triplet
+            grid = np.random.random(3)
+            seed = s + grid - .5
+            seeds.append(seed)
     seeds = asarray(seeds)
 
     if not seed_count_per_voxel:
-        # Randomize the seeds and select the requested amount
-        np.random.shuffle(seeds)
+        # Select the requested amount
         seeds = seeds[:seeds_count]
 
     # Apply the spatial transform