Merge pull request #1176 from ahojnnes/ransac

Add stop_probability for RANSAC

skimage/measure/fit.py

@@ -4,6 +4,9 @@ import numpy as np
 from scipy import optimize
 
 
+_EPSILON = np.spacing(1)
+
+
 def _check_data_dim(data, dim):
     if data.ndim != 2 or data.shape[1] != dim:
         raise ValueError('Input data must have shape (N, %d).' % dim)
@@ -465,9 +468,41 @@ class EllipseModel(BaseModel):
         return np.concatenate((x[..., None], y[..., None]), axis=t.ndim)
 
 
+def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
+    """Determine number trials such that at least one outlier-free subset is
+    sampled for the given inlier/outlier ratio.
+
+    Parameters
+    ----------
+    n_inliers : int
+        Number of inliers in the data.
+    n_samples : int
+        Total number of samples in the data.
+    min_samples : int
+        Minimum number of samples chosen randomly from original data.
+    probability : float
+        Probability (confidence) that one outlier-free sample is generated.
+
+    Returns
+    -------
+    trials : int
+        Number of trials.
+
+    """
+    inlier_ratio = n_inliers / float(n_samples)
+    nom = max(_EPSILON, 1 - probability)
+    denom = max(_EPSILON, 1 - inlier_ratio ** min_samples)
+    if nom == 1:
+        return 0
+    if denom == 1:
+        return float('inf')
+    return abs(float(np.ceil(np.log(nom) / np.log(denom))))
+
+
 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):
+           max_trials=100, stop_sample_num=np.inf, stop_residuals_sum=0,
+           stop_probability=1):
     """Fit a model to data with the RANSAC (random sample consensus) algorithm.
 
     RANSAC is an iterative algorithm for the robust estimation of parameters
@@ -525,7 +560,19 @@ def ransac(data, model_class, min_samples, residual_threshold,
     stop_sample_num : int, optional
         Stop iteration if at least this number of inliers are found.
     stop_residuals_sum : float, optional
-        Stop iteration if sum of residuals is less equal than this threshold.
+        Stop iteration if sum of residuals is less than or equal to this
+        threshold.
+    stop_probability : float in range [0, 1], optional
+        RANSAC iteration stops if at least one outlier-free set of the
+        training data is sampled with ``probability >= stop_probability``,
+        depending on the current best model's inlier ratio and the number
+        of trials. This requires to generate at least N samples (trials):
+
+            N >= log(1 - probability) / log(1 - e**m)
+
+        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.
 
     Returns
     -------
@@ -610,19 +657,28 @@ def ransac(data, model_class, min_samples, residual_threshold,
     best_inlier_residuals_sum = np.inf
     best_inliers = None
 
+    if min_samples < 0:
+        raise ValueError("`min_samples` must be greater than zero")
+
+    if max_trials < 0:
+        raise ValueError("`max_trials` must be greater than zero")
+
+    if stop_probability < 0 or stop_probability > 1:
+        raise ValueError("`stop_probability` must be in range [0, 1]")
+
     if not isinstance(data, list) and not isinstance(data, tuple):
         data = [data]
 
     # make sure data is list and not tuple, so it can be modified below
     data = list(data)
     # number of samples
-    N = data[0].shape[0]
+    num_samples = data[0].shape[0]
 
-    for _ in range(max_trials):
+    for num_trials in range(max_trials):
 
         # choose random sample set
         samples = []
-        random_idxs = np.random.randint(0, N, min_samples)
+        random_idxs = np.random.randint(0, num_samples, min_samples)
         for d in data:
             samples.append(d[random_idxs])
 
@@ -660,6 +716,9 @@ def ransac(data, model_class, min_samples, residual_threshold,
             if (
                 best_inlier_num >= stop_sample_num
                 or best_inlier_residuals_sum <= stop_residuals_sum
+                or num_trials
+                    >= _dynamic_max_trials(best_inlier_num, num_samples,
+                                           min_samples, stop_probability)
             ):
                 break
 

skimage/measure/tests/test_fit.py

@@ -2,6 +2,7 @@ import numpy as np
 from numpy.testing import assert_equal, assert_raises, assert_almost_equal
 from skimage.measure import LineModel, CircleModel, EllipseModel, ransac
 from skimage.transform import AffineTransform
+from skimage.measure.fit import _dynamic_max_trials
 
 
 def test_line_model_invalid_input():
@@ -204,6 +205,51 @@ def test_ransac_is_model_valid():
     assert_equal(inliers, None)
 
 
+def test_ransac_dynamic_max_trials():
+    # Numbers hand-calculated and confirmed on page 119 (Table 4.3) in
+    #   Hartley, R.~I. and Zisserman, A., 2004,
+    #   Multiple View Geometry in Computer Vision, Second Edition,
+    #   Cambridge University Press, ISBN: 0521540518
+
+    # e = 0%, min_samples = X
+    assert_equal(_dynamic_max_trials(100, 100, 2, 0.99), 1)
+
+    # e = 5%, min_samples = 2
+    assert_equal(_dynamic_max_trials(95, 100, 2, 0.99), 2)
+    # e = 10%, min_samples = 2
+    assert_equal(_dynamic_max_trials(90, 100, 2, 0.99), 3)
+    # e = 30%, min_samples = 2
+    assert_equal(_dynamic_max_trials(70, 100, 2, 0.99), 7)
+    # e = 50%, min_samples = 2
+    assert_equal(_dynamic_max_trials(50, 100, 2, 0.99), 17)
+
+    # e = 5%, min_samples = 8
+    assert_equal(_dynamic_max_trials(95, 100, 8, 0.99), 5)
+    # e = 10%, min_samples = 8
+    assert_equal(_dynamic_max_trials(90, 100, 8, 0.99), 9)
+    # e = 30%, min_samples = 8
+    assert_equal(_dynamic_max_trials(70, 100, 8, 0.99), 78)
+    # e = 50%, min_samples = 8
+    assert_equal(_dynamic_max_trials(50, 100, 8, 0.99), 1177)
+
+    # e = 0%, min_samples = 10
+    assert_equal(_dynamic_max_trials(1, 100, 10, 0), 0)
+    assert_equal(_dynamic_max_trials(1, 100, 10, 1), float('inf'))
+
+
+def test_ransac_invalid_input():
+    assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=-1,
+                                     residual_threshold=0)
+    assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2,
+                                     residual_threshold=0, max_trials=-1)
+    assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2,
+                                     residual_threshold=0,
+                                     stop_probability=-1)
+    assert_raises(ValueError, ransac, np.zeros((10, 2)), None, min_samples=2,
+                                     residual_threshold=0,
+                                     stop_probability=1.01)
+
+
 def test_deprecated_params_attribute():
     model = LineModel()
     model.params = (10, 1)