Merge pull request #2072 from kevin-keraudren/ransac_random_state

ENH: Pass np.random.RandomState to RANSAC

skimage/_shared/utils.py

@@ -3,6 +3,7 @@ import functools
 import sys
 import numpy as np
 import types
+import numbers
 
 import six
 
@@ -191,3 +192,30 @@ def copy_func(f, name=None):
     return types.FunctionType(six.get_function_code(f),
                               six.get_function_globals(f), name or f.__name__,
                               six.get_function_defaults(f), six.get_function_closure(f))
+
+
+def check_random_state(seed):
+    """Turn seed into a `np.random.RandomState` instance.
+
+    Parameters
+    ----------
+    seed : None, int or np.random.RandomState
+           If `seed` is None, return the RandomState singleton used by `np.random`.
+           If `seed` is an int, return a new RandomState instance seeded with `seed`.
+           If `seed` is already a RandomState instance, return it.
+
+    Raises
+    ------
+    ValueError
+        If `seed` is of the wrong type.
+
+    """
+    # Function originally from scikit-learn's module sklearn.utils.validation
+    if seed is None or seed is np.random:
+        return np.random.mtrand._rand
+    if isinstance(seed, (numbers.Integral, np.integer)):
+        return np.random.RandomState(seed)
+    if isinstance(seed, np.random.RandomState):
+        return seed
+    raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
+                     ' instance' % seed)

skimage/measure/fit.py

@@ -1,7 +1,7 @@
 import math
 import numpy as np
 from scipy import optimize
-from .._shared.utils import skimage_deprecation, warn
+from .._shared.utils import check_random_state, skimage_deprecation, warn
 
 
 def _check_data_dim(data, dim):
@@ -692,7 +692,7 @@ def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
 def ransac(data, model_class, min_samples, residual_threshold,
            is_data_valid=None, is_model_valid=None,
            max_trials=100, stop_sample_num=np.inf, stop_residuals_sum=0,
-           stop_probability=1):
+           stop_probability=1, random_state=None):
     """Fit a model to data with the RANSAC (random sample consensus) algorithm.
 
     RANSAC is an iterative algorithm for the robust estimation of parameters
@@ -765,6 +765,12 @@ def ransac(data, model_class, min_samples, residual_threshold,
         where the probability (confidence) is typically set to a high value
         such as 0.99, and e is the current fraction of inliers w.r.t. the
         total number of samples.
+    random_state : int, RandomState instance or None, optional
+        If int, random_state is the seed used by the random number generator;
+        If RandomState instance, random_state is the random number generator;
+        If None, the random number generator is the RandomState instance used
+        by `np.random`.
+    
 
     Returns
     -------
@@ -849,6 +855,8 @@ def ransac(data, model_class, min_samples, residual_threshold,
     best_inlier_num = 0
     best_inlier_residuals_sum = np.inf
     best_inliers = None
+    
+    random_state = check_random_state(random_state)
 
     if min_samples < 0:
         raise ValueError("`min_samples` must be greater than zero")
@@ -871,7 +879,7 @@ def ransac(data, model_class, min_samples, residual_threshold,
 
         # choose random sample set
         samples = []
-        random_idxs = np.random.randint(0, num_samples, min_samples)
+        random_idxs = random_state.randint(0, num_samples, min_samples)
         for d in data:
             samples.append(d[random_idxs])
 

skimage/measure/tests/test_fit.py

@@ -32,7 +32,8 @@ def test_line_model_estimate():
     model_est.estimate(data)
 
     # test whether estimated parameters almost equal original parameters
-    x = np.random.rand(100, 2)
+    random_state = np.random.RandomState(1234)
+    x = random_state.rand(100, 2)
     assert_almost_equal(model0.predict(x), model_est.predict(x), 1)
 
 
@@ -75,8 +76,8 @@ def test_line_modelND_estimate():
              10 * np.arange(-100,100)[...,np.newaxis] * model0.params[1])
 
     # add gaussian noise to data
-    np.random.seed(1234)
-    data = data0 + np.random.normal(size=data0.shape)
+    random_state = np.random.RandomState(1234)
+    data = data0 + random_state.normal(size=data0.shape)
 
     # estimate parameters of noisy data
     model_est = LineModelND()
@@ -130,8 +131,8 @@ def test_circle_model_estimate():
     data0 = model0.predict_xy(t)
 
     # add gaussian noise to data
-    np.random.seed(1234)
-    data = data0 + np.random.normal(size=data0.shape)
+    random_state = np.random.RandomState(1234)
+    data = data0 + random_state.normal(size=data0.shape)
 
     # estimate parameters of noisy data
     model_est = CircleModel()
@@ -172,8 +173,8 @@ def test_ellipse_model_estimate():
     data0 = model0.predict_xy(t)
 
     # add gaussian noise to data
-    np.random.seed(1234)
-    data = data0 + np.random.normal(size=data0.shape)
+    random_state = np.random.RandomState(1234)
+    data = data0 + random_state.normal(size=data0.shape)
 
     # estimate parameters of noisy data
     model_est = EllipseModel()
@@ -193,8 +194,6 @@ def test_ellipse_model_residuals():
 
 
 def test_ransac_shape():
-    np.random.seed(1)
-
     # generate original data without noise
     model0 = CircleModel()
     model0.params = (10, 12, 3)
@@ -208,7 +207,8 @@ def test_ransac_shape():
     data0[outliers[2], :] = (-100, -10)
 
     # estimate parameters of corrupted data
-    model_est, inliers = ransac(data0, CircleModel, 3, 5)
+    model_est, inliers = ransac(data0, CircleModel, 3, 5,
+                                random_state=1)
 
     # test whether estimated parameters equal original parameters
     assert_equal(model0.params, model_est.params)
@@ -217,10 +217,10 @@ def test_ransac_shape():
 
 
 def test_ransac_geometric():
-    np.random.seed(1)
+    random_state = np.random.RandomState(1)
 
     # generate original data without noise
-    src = 100 * np.random.random((50, 2))
+    src = 100 * random_state.random_sample((50, 2))
     model0 = AffineTransform(scale=(0.5, 0.3), rotation=1,
                              translation=(10, 20))
     dst = model0(src)
@@ -232,7 +232,8 @@ def test_ransac_geometric():
     dst[outliers[2]] = (50, 50)
 
     # estimate parameters of corrupted data
-    model_est, inliers = ransac((src, dst), AffineTransform, 2, 20)
+    model_est, inliers = ransac((src, dst), AffineTransform, 2, 20,
+                                random_state=random_state)
 
     # test whether estimated parameters equal original parameters
     assert_almost_equal(model0.params, model_est.params)
@@ -240,22 +241,20 @@ def test_ransac_geometric():
 
 
 def test_ransac_is_data_valid():
-    np.random.seed(1)
 
     is_data_valid = lambda data: data.shape[0] > 2
     model, inliers = ransac(np.empty((10, 2)), LineModelND, 2, np.inf,
-                            is_data_valid=is_data_valid)
+                            is_data_valid=is_data_valid, random_state=1)
     assert_equal(model, None)
     assert_equal(inliers, None)
 
 
 def test_ransac_is_model_valid():
-    np.random.seed(1)
 
     def is_model_valid(model, data):
         return False
     model, inliers = ransac(np.empty((10, 2)), LineModelND, 2, np.inf,
-                            is_model_valid=is_model_valid)
+                            is_model_valid=is_model_valid, random_state=1)
     assert_equal(model, None)
     assert_equal(inliers, None)