diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index e0118388..7de90291 100644 --- a/skimage/measure/fit.py +++ b/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 diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 3ef7a0f0..32d9ef22 100644 --- a/skimage/measure/tests/test_fit.py +++ b/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)