Merge pull request #1795 from kevin-keraudren/ransac-linemodel3D

Adding LineModel3D for RANSAC, unit test and example.

TODO.txt

@@ -8,6 +8,10 @@ Version 0.14
 * Remove deprecated ``skimage.restoration.nl_means_denoising``.
 * Remove deprecated ``skimage.filters.gaussian_filter``.
 * Remove deprecated ``skimage.filters.gabor_filter``.
+* Remove deprecated ``skimage.measure.LineModel`` and
+  add an alias LineModel = LineModelND. While the deprecated LineModel has for
+  parameters `(dist, theta)`, LineModelND has the more general parameters
+  `(origin, direction)`.
 
 
 Version 0.13

doc/examples/plot_ransac3D.py

@@ -0,0 +1,40 @@
+"""
+============================================
+Robust 3D line model estimation using RANSAC
+============================================
+
+In this example we see how to robustly fit a 3D line model to faulty data using
+the RANSAC algorithm.
+
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from skimage.measure import LineModelND, ransac
+
+np.random.seed(seed=1)
+
+# generate coordinates of line
+point = np.array([0, 0, 0], dtype='float')
+direction = np.array([1, 1, 1], dtype='float') / np.sqrt(3)
+xyz = point + 10 * np.arange(-100, 100)[..., np.newaxis] * direction
+
+# add gaussian noise to coordinates
+noise = np.random.normal(size=xyz.shape)
+xyz += 0.5 * noise
+xyz[::2] += 20 * noise[::2]
+xyz[::4] += 100 * noise[::4]
+
+# robustly fit line only using inlier data with RANSAC algorithm
+model_robust, inliers = ransac(xyz, LineModelND, min_samples=2,
+                               residual_threshold=1, max_trials=1000)
+outliers = inliers == False
+
+fig = plt.figure()
+ax = fig.add_subplot(111, projection='3d')
+ax.scatter(xyz[inliers][:, 0], xyz[inliers][:, 1], xyz[inliers][:, 2], c='b',
+           marker='o', label='Inlier data')
+ax.scatter(xyz[outliers][:, 0], xyz[outliers][:, 1], xyz[outliers][:, 2], c='r',
+           marker='o', label='Outlier data')
+ax.legend(loc='lower left')
+plt.show()

skimage/measure/__init__.py

@@ -7,7 +7,7 @@ from ._polygon import approximate_polygon, subdivide_polygon
 from ._pnpoly import points_in_poly, grid_points_in_poly
 from ._moments import moments, moments_central, moments_normalized, moments_hu
 from .profile import profile_line
-from .fit import LineModel, CircleModel, EllipseModel, ransac
+from .fit import LineModel, LineModelND, CircleModel, EllipseModel, ransac
 from .block import block_reduce
 from ._label import label
 
@@ -19,6 +19,7 @@ __all__ = ['find_contours',
            'approximate_polygon',
            'subdivide_polygon',
            'LineModel',
+           'LineModelND',
            'CircleModel',
            'EllipseModel',
            'ransac',

skimage/measure/fit.py

@@ -2,6 +2,7 @@ import math
 import warnings
 import numpy as np
 from scipy import optimize
+from .._shared.utils import skimage_deprecation
 
 
 def _check_data_dim(data, dim):
@@ -9,6 +10,11 @@ def _check_data_dim(data, dim):
         raise ValueError('Input data must have shape (N, %d).' % dim)
 
 
+def _check_data_atleast_2D(data):
+    if data.ndim < 2 or data.shape[1] < 2:
+        raise ValueError('Input data must be at least 2D.')
+
+
 class BaseModel(object):
 
     def __init__(self):
@@ -39,6 +45,8 @@ class LineModel(BaseModel):
 
     A minimum number of 2 points is required to solve for the parameters.
 
+    **Deprecated class**. Use ``LineModelND`` instead.
+
     Attributes
     ----------
     params : tuple
@@ -46,6 +54,11 @@ class LineModel(BaseModel):
 
     """
 
+    def __init__(self):
+        self.params = None
+        warnings.warn(skimage_deprecation('`LineModel` is deprecated, '
+                      'use `LineModelND` instead.'))
+
     def estimate(self, data):
         """Estimate line model from data using total least squares.
 
@@ -156,6 +169,157 @@ class LineModel(BaseModel):
         return (dist - x * math.cos(theta)) / math.sin(theta)
 
 
+class LineModelND(BaseModel):
+    """Total least squares estimator for N-dimensional lines.
+
+    Lines are defined by a point (origin) and a unit vector (direction)
+    according to the following vector equation::
+
+        X = origin + lambda * direction
+
+    Attributes
+    ----------
+    params : tuple
+        Line model parameters in the following order `origin`, `direction`.
+
+    """
+
+    def estimate(self, data):
+        """Estimate line model from data.
+
+        Parameters
+        ----------
+        data : (N, dim) array
+            N points in a space of dimensionality dim >= 2.
+
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+        """
+
+        _check_data_atleast_2D(data)
+
+        X0 = data.mean(axis=0)
+
+        if data.shape[0] == 2:  # well determined
+            u = data[1] - data[0]
+            norm = np.linalg.norm(u)
+            if norm > 0:
+                u /= norm
+        elif data.shape[0] > 2:  # over-determined
+            data = data - X0
+            # first principal component
+            # Note: without full_matrices=False Python dies with joblib
+            # parallel_for.
+            _, _, u = np.linalg.svd(data, full_matrices=False)
+            u = u[0]
+        else:  # under-determined
+            raise ValueError('At least 2 input points needed.')
+
+        self.params = (X0, u)
+
+        return True
+
+    def residuals(self, data):
+        """Determine residuals of data to model.
+
+        For each point the shortest distance to the line is returned.
+        It is obtained by projecting the data onto the line.
+
+        Parameters
+        ----------
+        data : (N, dim) array
+            N points in a space of dimension dim.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for each data point.
+        """
+
+        X0, u = self.params
+        return np.linalg.norm((data - X0) -
+                              np.dot(data - X0, u)[..., np.newaxis] * u, axis=1)
+
+    def predict(self, x, axis=0, params=None):
+        """Predict intersection of the estimated line model with a hyperplane
+        orthogonal to a given axis.
+
+        Parameters
+        ----------
+        x : array
+            coordinates along an axis.
+        axis : int
+            axis orthogonal to the hyperplane intersecting the line.
+        params : (2, ) array, optional
+            Optional custom parameter set in the form (`origin`, `direction`).
+
+        Returns
+        -------
+        y : array
+            Predicted coordinates.
+
+        If the line is parallel to the given axis, a ValueError is raised.
+        """
+
+        if params is None:
+            params = self.params
+
+        X0, u = params
+
+        if u[axis] == 0:
+            # line parallel to axis
+            raise ValueError('Line parallel to axis %s' % axis)
+
+        l = (x - X0[axis]) / u[axis]
+        return X0 + l[..., np.newaxis] * u
+
+    def predict_x(self, y, params=None, new_params=None):
+        """Predict x-coordinates for 2D lines using the estimated model.
+
+        Alias for::
+
+            predict(y, axis=1)[:, 0]
+
+        Parameters
+        ----------
+        y : array
+            y-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set in the form (`origin`, `direction`).
+
+        Returns
+        -------
+        x : array
+            Predicted x-coordinates.
+
+        """
+        return self.predict(y, axis=1, params=params)[:, 0]
+
+    def predict_y(self, x, params=None):
+        """Predict y-coordinates  for 2D lines using the estimated model.
+
+        Alias for::
+
+            predict(x, axis=0)[:, 1]
+
+        Parameters
+        ----------
+        x : array
+            x-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set in the form (`origin`, `direction`).
+
+        Returns
+        -------
+        y : array
+            Predicted y-coordinates.
+
+        """
+        return self.predict(x, axis=0, params=params)[:, 1]
+
+
 class CircleModel(BaseModel):
 
     """Total least squares estimator for 2D circles.

skimage/measure/tests/test_fit.py

@@ -1,6 +1,6 @@
 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.measure import LineModel, LineModelND, CircleModel, EllipseModel, ransac
 from skimage.transform import AffineTransform
 from skimage.measure.fit import _dynamic_max_trials
 from skimage._shared._warnings import expected_warnings
@@ -54,6 +54,62 @@ def test_line_model_under_determined():
     assert_raises(ValueError, LineModel().estimate, data)
 
 
+def test_line_modelND_invalid_input():
+    assert_raises(ValueError, LineModelND().estimate, np.empty((5, 1)))
+
+
+def test_line_modelND_predict():
+    model = LineModelND()
+    model.params = (np.array([0,0]), np.array([0.2,0.98]))
+    x = np.arange(-10, 10)
+    y = model.predict_y(x)
+    assert_almost_equal(x, model.predict_x(y))
+
+
+def test_line_modelND_estimate():
+    # generate original data without noise
+    model0 = LineModelND()
+    model0.params = (np.array([0,0,0], dtype='float'),
+                         np.array([1,1,1], dtype='float')/np.sqrt(3))
+    # we scale the unit vector with a factor 10 when generating points on the
+    # line in order to compensate for the scale of the random noise
+    data0 = (model0.params[0] +
+             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)
+
+    # estimate parameters of noisy data
+    model_est = LineModelND()
+    model_est.estimate(data)
+
+    # test whether estimated parameters are correct
+    # we use the following geometric property: two aligned vectors have
+    # a cross-product equal to zero
+    # test if direction vectors are aligned
+    assert_almost_equal(np.linalg.norm(np.cross(model0.params[1],
+                                                model_est.params[1])), 0, 1)
+    # test if origins are aligned with the direction
+    a = model_est.params[0] - model0.params[0]
+    if np.linalg.norm(a) > 0:
+        a /= np.linalg.norm(a)
+    assert_almost_equal(np.linalg.norm(np.cross(model0.params[1], a)), 0, 1)
+
+
+def test_line_modelND_residuals():
+    model = LineModelND()
+    model.params = (np.array([0,0,0]), np.array([0,0,1]))
+    assert_equal(abs(model.residuals(np.array([[0, 0,0]]))), 0)
+    assert_equal(abs(model.residuals(np.array([[0,0,1]]))), 0)
+    assert_equal(abs(model.residuals(np.array([[10, 0,0]]))), 10)
+
+
+def test_line_modelND_under_determined():
+    data = np.empty((1, 3))
+    assert_raises(ValueError, LineModelND().estimate, data)
+
+
 def test_circle_model_invalid_input():
     assert_raises(ValueError, CircleModel().estimate, np.empty((5, 3)))