diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index e0118388..d9dd11a8 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,44 @@ 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 +563,18 @@ 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 in RANSAC. This requires to generate at least + N samples (iterations): + + N >= log(1 - probability) / log(1 - e**m) + + where the probability (confidence) is typically set to high value such + as 0.99 and e is the current fraction of inliers w.r.t. the total + number of samples. Returns ------- @@ -616,13 +665,13 @@ def ransac(data, model_class, min_samples, residual_threshold, # 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 +709,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..af259777 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,38 @@ 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_deprecated_params_attribute(): model = LineModel() model.params = (10, 1)