Merge pull request #1411 from ahojnnes/rans-fix

Fix RANSAC for invalid model estimation and confidence corner case

skimage/measure/fit.py

@@ -4,9 +4,6 @@ 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)
@@ -57,6 +54,11 @@ class LineModel(BaseModel):
         data : (N, 2) array
             N points with ``(x, y)`` coordinates, respectively.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         _check_data_dim(data, dim=2)
@@ -81,6 +83,8 @@ class LineModel(BaseModel):
 
         self.params = (dist, theta)
 
+        return True
+
     def residuals(self, data):
         """Determine residuals of data to model.
 
@@ -182,6 +186,11 @@ class CircleModel(BaseModel):
         data : (N, 2) array
             N points with ``(x, y)`` coordinates, respectively.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         _check_data_dim(data, dim=2)
@@ -217,6 +226,8 @@ class CircleModel(BaseModel):
 
         self.params = params
 
+        return True
+
     def residuals(self, data):
         """Determine residuals of data to model.
 
@@ -313,6 +324,11 @@ class EllipseModel(BaseModel):
         data : (N, 2) array
             N points with ``(x, y)`` coordinates, respectively.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         _check_data_dim(data, dim=2)
@@ -373,6 +389,8 @@ class EllipseModel(BaseModel):
 
         self.params = params[:5]
 
+        return True
+
     def residuals(self, data):
         """Determine residuals of data to model.
 
@@ -471,7 +489,6 @@ class EllipseModel(BaseModel):
 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
@@ -482,21 +499,31 @@ def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
         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:
+    if n_inliers == 0:
         return float('inf')
-    return abs(float(np.ceil(np.log(nom) / np.log(denom))))
+
+    nom = 1 - probability
+    if nom == 0:
+        return float('inf')
+
+    inlier_ratio = n_inliers / float(n_samples)
+    denom = 1 - inlier_ratio ** min_samples
+    if denom == 0:
+        return 1
+    elif denom == 1:
+        return float('inf')
+
+    nom = np.log(nom)
+    denom = np.log(denom)
+    if denom == 0:
+        return 0
+
+    return int(np.ceil(nom / denom))
 
 
 def ransac(data, model_class, min_samples, residual_threshold,
@@ -542,9 +569,11 @@ def ransac(data, model_class, min_samples, residual_threshold,
     model_class : object
         Object with the following object methods:
 
-         * ``estimate(*data)``
+         * ``success = estimate(*data)``
          * ``residuals(*data)``
 
+        where `success` indicates whether the model estimation succeeded
+        (`True` or `None` for success, `False` for failure).
     min_samples : int
         The minimum number of data points to fit a model to.
     residual_threshold : float
@@ -612,6 +641,7 @@ def ransac(data, model_class, min_samples, residual_threshold,
 
     >>> model = EllipseModel()
     >>> model.estimate(data)
+    True
     >>> model.params # doctest: +SKIP
     array([ -3.30354146e+03,  -2.87791160e+03,   5.59062118e+03,
              7.84365066e+00,   7.19203152e-01])
@@ -688,7 +718,12 @@ def ransac(data, model_class, min_samples, residual_threshold,
 
         # estimate model for current random sample set
         sample_model = model_class()
-        sample_model.estimate(*samples)
+
+        success = sample_model.estimate(*samples)
+
+        if success is not None:  # backwards compatibility
+            if not success:
+                continue
 
         # check if estimated model is valid
         if is_model_valid is not None and not is_model_valid(sample_model,

skimage/measure/tests/test_fit.py

@@ -233,9 +233,9 @@ def test_ransac_dynamic_max_trials():
     # 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'))
+    # e = 0%, min_samples = 5
+    assert_equal(_dynamic_max_trials(1, 100, 5, 0), 0)
+    assert_equal(_dynamic_max_trials(1, 100, 5, 1), float('inf'))
 
 
 def test_ransac_invalid_input():

skimage/transform/_geometric.py

@@ -263,6 +263,11 @@ class ProjectiveTransform(GeometricTransform):
         dst : (N, 2) array
             Destination coordinates.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         try:
@@ -270,7 +275,7 @@ class ProjectiveTransform(GeometricTransform):
             dst_matrix, dst = _center_and_normalize_points(dst)
         except ZeroDivisionError:
             self.params = np.nan * np.empty((3, 3))
-            return
+            return False
 
         xs = src[:, 0]
         ys = src[:, 1]
@@ -309,6 +314,8 @@ class ProjectiveTransform(GeometricTransform):
 
         self.params = H
 
+        return True
+
     def __add__(self, other):
         """Combine this transformation with another.
 
@@ -459,6 +466,11 @@ class PiecewiseAffineTransform(GeometricTransform):
         dst : (N, 2) array
             Destination coordinates.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         # forward piecewise affine
@@ -481,6 +493,8 @@ class PiecewiseAffineTransform(GeometricTransform):
             affine.estimate(dst[tri, :], src[tri, :])
             self.inverse_affines.append(affine)
 
+        return True
+
     def __call__(self, coords):
         """Apply forward transformation.
 
@@ -658,6 +672,11 @@ class SimilarityTransform(ProjectiveTransform):
         dst : (N, 2) array
             Destination coordinates.
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
 
         try:
@@ -665,7 +684,7 @@ class SimilarityTransform(ProjectiveTransform):
             dst_matrix, dst = _center_and_normalize_points(dst)
         except ZeroDivisionError:
             self.params = np.nan * np.empty((3, 3))
-            return
+            return False
 
         xs = src[:, 0]
         ys = src[:, 1]
@@ -699,6 +718,7 @@ class SimilarityTransform(ProjectiveTransform):
 
         self.params = S
 
+        return True
 
     @property
     def scale(self):
@@ -798,6 +818,11 @@ class PolynomialTransform(GeometricTransform):
         order : int, optional
             Polynomial order (number of coefficients is order + 1).
 
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+
         """
         xs = src[:, 0]
         ys = src[:, 1]
@@ -828,6 +853,8 @@ class PolynomialTransform(GeometricTransform):
 
         self.params = params.reshape((2, u // 2))
 
+        return True
+
     def __call__(self, coords):
         """Apply forward transformation.