Adding LineModel3D for RANSAC, unit test and example.

doc/examples/plot_ransac3D.py

@@ -0,0 +1,38 @@
+"""
+=========================================
+Robust 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 LineModel3D, 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, LineModel3D, 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, LineModel3D, CircleModel, EllipseModel, ransac
 from .block import block_reduce
 from ._label import label
 
@@ -19,6 +19,7 @@ __all__ = ['find_contours',
            'approximate_polygon',
            'subdivide_polygon',
            'LineModel',
+           'LineModel3D',
            'CircleModel',
            'EllipseModel',
            'ransac',

skimage/measure/fit.py

@@ -156,6 +156,67 @@ class LineModel(BaseModel):
         return (dist - x * math.cos(theta)) / math.sin(theta)
 
 
+class LineModel3D(BaseModel):
+    """Total least squares estimator for 3D lines.
+
+    Lines are defined by point and a unit vector (direction).
+
+    Attributes
+    ----------
+    params : tuple
+        Line model parameters in the following order `point`, `direction`.
+    """
+
+    def estimate(self, data):
+        """Estimate line model from data.
+
+        Parameters
+        ----------
+        data : (N, 3) array
+            N points with ``(x, y, z)`` coordinates, respectively.
+
+        Returns
+        -------
+        success : bool
+            True, if model estimation succeeds.
+        """
+        X0 = data.mean(axis=0)
+
+        if data.shape[0] == 2:  # well determined
+            u = data[1] - data[0]
+            u /= np.linalg.norm(u)
+        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, 3) array
+            N points with ``(x, y, z)`` coordinates, respectively.
+
+        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)
+
+
 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, LineModel3D, CircleModel, EllipseModel, ransac
 from skimage.transform import AffineTransform
 from skimage.measure.fit import _dynamic_max_trials
 from skimage._shared._warnings import expected_warnings
@@ -53,6 +53,45 @@ def test_line_model_under_determined():
     data = np.empty((1, 2))
     assert_raises(ValueError, LineModel().estimate, data)
 
+def test_line_model3D_estimate():
+    # generate original data without noise
+    model0 = LineModel3D()
+    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 = LineModel3D()
+    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_model3D_residuals():
+    model = LineModel3D()
+    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_model3D_under_determined():
+    data = np.empty((1, 3))
+    assert_raises(ValueError, LineModel().estimate, data)
+
 
 def test_circle_model_invalid_input():
     assert_raises(ValueError, CircleModel().estimate, np.empty((5, 3)))