From a606a53875da213e337a814f07675f1e322be807 Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Tue, 1 Dec 2015 12:17:22 +0000
Subject: [PATCH 1/8] Adding LineModel3D for RANSAC, unit test and example.

---
 doc/examples/plot_ransac3D.py     | 38 +++++++++++++++++++
 skimage/measure/__init__.py       |  3 +-
 skimage/measure/fit.py            | 61 +++++++++++++++++++++++++++++++
 skimage/measure/tests/test_fit.py | 41 ++++++++++++++++++++-
 4 files changed, 141 insertions(+), 2 deletions(-)
 create mode 100644 doc/examples/plot_ransac3D.py

diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py
new file mode 100644
index 00000000..d3d60605
--- /dev/null
+++ b/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()
diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index 9731d6da..81805292 100755
--- a/skimage/measure/__init__.py
+++ b/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',
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index f8e65117..49d55d0f 100644
--- a/skimage/measure/fit.py
+++ b/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.
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index 7f98c971..32963cb0 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/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)))

From c1c0d18ba194097789fd3b9e7715901e4a4d1c82 Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Tue, 1 Dec 2015 14:32:25 +0000
Subject: [PATCH 2/8] reformatting code to wrap at 80 characters per line

---
 doc/examples/plot_ransac3D.py     |  6 ++++--
 skimage/measure/fit.py            | 10 ++++++----
 skimage/measure/tests/test_fit.py | 16 ++++++++++------
 3 files changed, 20 insertions(+), 12 deletions(-)

diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py
index d3d60605..af2fa373 100644
--- a/doc/examples/plot_ransac3D.py
+++ b/doc/examples/plot_ransac3D.py
@@ -32,7 +32,9 @@ 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.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()
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index 49d55d0f..9849c996 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -188,7 +188,8 @@ class LineModel3D(BaseModel):
         elif data.shape[0] > 2:  # over-determined
             data = data - X0
             # first principal component
-            # Note: without full_matrices=False Python dies with joblib parallel_for.
+            # 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
@@ -200,8 +201,8 @@ class LineModel3D(BaseModel):
 
     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.
+        For each point the shortest distance to the line is returned.
+        It is obtained by projecting the data onto the line.
 
         Parameters
         ----------
@@ -214,7 +215,8 @@ class LineModel3D(BaseModel):
             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)
+        return np.linalg.norm((data - X0) -
+                              np.dot(data - X0, u)[..., np.newaxis] * u, axis=1)
 
 
 class CircleModel(BaseModel):
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index 32963cb0..b9c77031 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/skimage/measure/tests/test_fit.py
@@ -56,10 +56,12 @@ def test_line_model_under_determined():
 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]
+    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)
@@ -70,9 +72,11 @@ def test_line_model3D_estimate():
     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
+    # 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)
+    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:

From 5bb206233d5cb91fb6059394cea44ff8af4efdab Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Tue, 1 Dec 2015 15:18:57 +0000
Subject: [PATCH 3/8] fixed typo and added blank line

---
 skimage/measure/fit.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index 9849c996..4ada55b2 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -159,7 +159,7 @@ class LineModel(BaseModel):
 class LineModel3D(BaseModel):
     """Total least squares estimator for 3D lines.
 
-    Lines are defined by point and a unit vector (direction).
+    Lines are defined by a point and a unit vector (direction).
 
     Attributes
     ----------
@@ -201,6 +201,7 @@ class LineModel3D(BaseModel):
 
     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.
 

From 040a53456d933f8587dbdcfb093b3bc03db2c5e9 Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Wed, 2 Dec 2015 08:58:16 +0000
Subject: [PATCH 4/8] Merged LineModel3D with LineModel in order to form a
 unified model.

 - model.params holds the legacy params (2D representation),
 - model.new_params holds the true ND line representation.

The proposed LineModel behaves identically as the implementation in master in the 2D case. In the 3D case, it behaves similarly to the previously proposed LineModel3D, as long as new_params is used instead of params. This implementation thus takes advantage that the 2D case is a special case of the ND general model.
---
 doc/examples/plot_ransac3D.py     |   4 +-
 skimage/measure/__init__.py       |   3 +-
 skimage/measure/fit.py            | 212 ++++++++++++++----------------
 skimage/measure/tests/test_fit.py |  40 +++---
 4 files changed, 122 insertions(+), 137 deletions(-)

diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py
index af2fa373..93cce22e 100644
--- a/doc/examples/plot_ransac3D.py
+++ b/doc/examples/plot_ransac3D.py
@@ -10,7 +10,7 @@ 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
+from skimage.measure import LineModel, ransac
 
 np.random.seed(seed=1)
 
@@ -26,7 +26,7 @@ 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,
+model_robust, inliers = ransac(xyz, LineModel, min_samples=2,
                                residual_threshold=1, max_trials=1000)
 outliers = inliers == False
 
diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index 81805292..9731d6da 100755
--- a/skimage/measure/__init__.py
+++ b/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, LineModel3D, CircleModel, EllipseModel, ransac
+from .fit import LineModel, CircleModel, EllipseModel, ransac
 from .block import block_reduce
 from ._label import label
 
@@ -19,7 +19,6 @@ __all__ = ['find_contours',
            'approximate_polygon',
            'subdivide_polygon',
            'LineModel',
-           'LineModel3D',
            'CircleModel',
            'EllipseModel',
            'ransac',
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index 4ada55b2..b4829613 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -9,10 +9,16 @@ 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):
         self.params = None
+        self.new_params = None
 
     @property
     def _params(self):
@@ -22,60 +28,61 @@ class BaseModel(object):
 
 
 class LineModel(BaseModel):
+    """Total least squares estimator for ND lines.
 
-    """Total least squares estimator for 2D lines.
-
-    Lines are parameterized using polar coordinates as functional model::
-
-        dist = x * cos(theta) + y * sin(theta)
-
-    This parameterization is able to model vertical lines in contrast to the
-    standard line model ``y = a*x + b``.
-
-    This estimator minimizes the squared distances from all points to the
-    line::
-
-        min{ sum((dist - x_i * cos(theta) + y_i * sin(theta))**2) }
-
-    A minimum number of 2 points is required to solve for the parameters.
+    Lines are defined by a point and a unit vector (direction).
 
     Attributes
     ----------
     params : tuple
-        Line model parameters in the following order `dist`, `theta`.
-
+        2D line model parameters in the following order `dist`, `theta`.
+        If dim > 2, these parameters correspond to the projection of the line
+        into the space spanned by the first two axes.
+        These parameters correspond to the functional model:
+            dist = x * cos(theta) + y * sin(theta)
+    new_params : tuple
+        ND line model parameters in the following order `X0`, `direction`.
+        These parameters correspond to the vector equation
+            X = X0 + lambda * direction
     """
 
     def estimate(self, data):
-        """Estimate line model from data using total least squares.
+        """Estimate line model from data.
 
         Parameters
         ----------
-        data : (N, 2) array
-            N points with ``(x, y)`` coordinates, respectively.
+        data : (N, dim) array
+            N points in a space of dimensionality dim >= 2.
 
         Returns
         -------
         success : bool
             True, if model estimation succeeds.
-
         """
 
-        _check_data_dim(data, dim=2)
+        _check_data_atleast_2D(data)
 
         X0 = data.mean(axis=0)
 
         if data.shape[0] == 2:  # well determined
-            theta = np.arctan2(data[1, 1] - data[0, 1],
-                               data[1, 0] - data[0, 0])
+            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
-            _, _, v = np.linalg.svd(data)
-            theta = np.arctan2(v[0, 1], v[0, 0])
+            # 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.new_params = (X0, u)
+
+        # legacy LineModel (2D case)
+        theta = np.arctan2(u[1], u[0])
         # angle perpendicular to line angle
         theta = (theta + np.pi / 2) % np.pi
         # line always passes through mean
@@ -89,37 +96,84 @@ class LineModel(BaseModel):
         """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, 2) array
-            N points with ``(x, y)`` coordinates, respectively.
+        data : (N, dim) array
+            N points in a space of dimension dim.
 
         Returns
         -------
         residuals : (N, ) array
             Residual for each data point.
-
         """
+        if self.new_params is None:
+            self.new_params = self._params_from_polar(self.params)
 
-        _check_data_dim(data, dim=2)
+        X0, u = self.new_params
+        return np.linalg.norm((data - X0) -
+                              np.dot(data - X0, u)[..., np.newaxis] * u, axis=1)
 
-        dist, theta = self.params
+    def _params_from_polar(self, params):
+        (dist, theta) = params
+        u = np.array([math.cos(theta - np.pi / 2), math.sin(theta - np.pi / 2)])
+        if math.cos(theta) == 0:
+            X0 = np.array([0, dist / math.sin(theta)])
+        else:
+            X0 = np.array([dist / math.cos(theta), 0])
+        return X0, u
 
-        x = data[:, 0]
-        y = data[:, 1]
+    def predict(self, x, axis=0, params=None, new_params=None):
+        """Predict intersection of the estimated line model with a hyperplane
+        orthogonal to a given axis.
 
-        return dist - (x * math.cos(theta) + y * math.sin(theta))
+        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 (`dist`, `theta`).
+        new_params : (2, ) array, optional
+            Optional custom parameter set in the form (`X0`, `direction`).
 
-    def predict_x(self, y, params=None):
+        Returns
+        -------
+        y : array
+            Predicted coordinates.
+
+        If the line is parallel to the given axis, a ValueError is raised.
+        """
+        if new_params is None:
+            if params is None and self.new_params is not None:
+                new_params = self.new_params
+            else:
+                new_params = self._params_from_polar(params or self.params)
+
+        X0, u = new_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 using the estimated model.
 
+        Alias for predict(y, axis=1)[:, 0].
+
         Parameters
         ----------
         y : array
             y-coordinates.
         params : (2, ) array, optional
-            Optional custom parameter set.
+            Optional custom parameter set in the form (`dist`, `theta`).
+        new_params : (2, ) array, optional
+            Optional custom parameter set in the form (`X0`, `direction`).
 
         Returns
         -------
@@ -127,21 +181,22 @@ class LineModel(BaseModel):
             Predicted x-coordinates.
 
         """
+        return self.predict(y, axis=1, params=params,
+                            new_params=new_params)[:, 0]
 
-        if params is None:
-            params = self.params
-        dist, theta = params
-        return (dist - y * math.sin(theta)) / math.cos(theta)
-
-    def predict_y(self, x, params=None):
+    def predict_y(self, x, params=None, new_params=None):
         """Predict y-coordinates using the estimated model.
 
+        Alias for predict(x, axis=1)[:, 1].
+
         Parameters
         ----------
         x : array
             x-coordinates.
         params : (2, ) array, optional
-            Optional custom parameter set.
+            Optional custom parameter set in the form (`dist`, `theta`).
+        new_params : (2, ) array, optional
+            Optional custom parameter set in the form (`X0`, `direction`).
 
         Returns
         -------
@@ -149,75 +204,8 @@ class LineModel(BaseModel):
             Predicted y-coordinates.
 
         """
-
-        if params is None:
-            params = self.params
-        dist, theta = params
-        return (dist - x * math.cos(theta)) / math.sin(theta)
-
-
-class LineModel3D(BaseModel):
-    """Total least squares estimator for 3D lines.
-
-    Lines are defined by a 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)
+        return self.predict(x, axis=0, params=params,
+                            new_params=new_params)[:, 1]
 
 
 class CircleModel(BaseModel):
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index b9c77031..ab518126 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/skimage/measure/tests/test_fit.py
@@ -1,13 +1,13 @@
 import numpy as np
 from numpy.testing import assert_equal, assert_raises, assert_almost_equal
-from skimage.measure import LineModel, LineModel3D, CircleModel, EllipseModel, ransac
+from skimage.measure import LineModel, CircleModel, EllipseModel, ransac
 from skimage.transform import AffineTransform
 from skimage.measure.fit import _dynamic_max_trials
 from skimage._shared._warnings import expected_warnings
 
 
 def test_line_model_invalid_input():
-    assert_raises(ValueError, LineModel().estimate, np.empty((5, 3)))
+    assert_raises(ValueError, LineModel().estimate, np.empty((5, 1)))
 
 
 def test_line_model_predict():
@@ -42,61 +42,59 @@ def test_line_model_residuals():
     model = LineModel()
     model.params = (0, 0)
     assert_equal(abs(model.residuals(np.array([[0, 0]]))), 0)
-    assert_equal(abs(model.residuals(np.array([[0, 10]]))), 0)
+    assert_almost_equal(abs(model.residuals(np.array([[0, 10]]))), 0)
     assert_equal(abs(model.residuals(np.array([[10, 0]]))), 10)
+    model = LineModel()
     model.params = (5, np.pi / 4)
-    assert_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
+    assert_almost_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
     assert_almost_equal(abs(model.residuals(np.array([[np.sqrt(50), 0]]))), 0)
 
 
 def test_line_model_under_determined():
     data = np.empty((1, 2))
     assert_raises(ValueError, LineModel().estimate, data)
+    data = np.empty((1, 3))
+    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))
+    model0 = LineModel()
+    model0.new_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])
+    data0 = (model0.new_params[0] +
+             10 * np.arange(-100,100)[...,np.newaxis] * model0.new_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 = LineModel()
     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)
+    assert_almost_equal(np.linalg.norm(np.cross(model0.new_params[1],
+                                                model_est.new_params[1])), 0, 1)
     # test if origins are aligned with the direction
-    a = model_est.params[0] - model0.params[0]
+    a = model_est.new_params[0] - model0.new_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)
+    assert_almost_equal(np.linalg.norm(np.cross(model0.new_params[1], a)), 0, 1)
 
 
 def test_line_model3D_residuals():
-    model = LineModel3D()
-    model.params = (np.array([0,0,0]), np.array([0,0,1]))
+    model = LineModel()
+    model.new_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)))
 

From 6a2961fdee6fb3324f2a5492f91551d55069355d Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Thu, 3 Dec 2015 22:03:10 +0000
Subject: [PATCH 5/8] we now have the old LineModel with params (dist, theta),
 and the new LineModelND with params (origin, direction).

---
 doc/examples/plot_ransac3D.py     |   4 +-
 skimage/measure/__init__.py       |   3 +-
 skimage/measure/fit.py            | 206 ++++++++++++++++++++++--------
 skimage/measure/tests/test_fit.py |  53 +++++---
 4 files changed, 192 insertions(+), 74 deletions(-)

diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py
index 93cce22e..b7816e47 100644
--- a/doc/examples/plot_ransac3D.py
+++ b/doc/examples/plot_ransac3D.py
@@ -10,7 +10,7 @@ the RANSAC algorithm.
 import numpy as np
 from matplotlib import pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
-from skimage.measure import LineModel, ransac
+from skimage.measure import LineModelND, ransac
 
 np.random.seed(seed=1)
 
@@ -26,7 +26,7 @@ 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, LineModel, min_samples=2,
+model_robust, inliers = ransac(xyz, LineModelND, min_samples=2,
                                residual_threshold=1, max_trials=1000)
 outliers = inliers == False
 
diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index 9731d6da..9e8df953 100755
--- a/skimage/measure/__init__.py
+++ b/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',
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index b4829613..f4a8df37 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -18,7 +18,6 @@ class BaseModel(object):
 
     def __init__(self):
         self.params = None
-        self.new_params = None
 
     @property
     def _params(self):
@@ -28,22 +27,151 @@ class BaseModel(object):
 
 
 class LineModel(BaseModel):
-    """Total least squares estimator for ND lines.
 
-    Lines are defined by a point and a unit vector (direction).
+    """Total least squares estimator for 2D lines.
+
+    Lines are parameterized using polar coordinates as functional model::
+
+        dist = x * cos(theta) + y * sin(theta)
+
+    This parameterization is able to model vertical lines in contrast to the
+    standard line model ``y = a*x + b``.
+
+    This estimator minimizes the squared distances from all points to the
+    line::
+
+        min{ sum((dist - x_i * cos(theta) + y_i * sin(theta))**2) }
+
+    A minimum number of 2 points is required to solve for the parameters.
 
     Attributes
     ----------
     params : tuple
-        2D line model parameters in the following order `dist`, `theta`.
-        If dim > 2, these parameters correspond to the projection of the line
-        into the space spanned by the first two axes.
-        These parameters correspond to the functional model:
-            dist = x * cos(theta) + y * sin(theta)
-    new_params : tuple
-        ND line model parameters in the following order `X0`, `direction`.
+        Line model parameters in the following order `dist`, `theta`.
+
+    """
+
+    def estimate(self, data):
+        """Estimate line model from data using total least squares.
+
+        Parameters
+        ----------
+        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)
+
+        X0 = data.mean(axis=0)
+
+        if data.shape[0] == 2:  # well determined
+            theta = np.arctan2(data[1, 1] - data[0, 1],
+                               data[1, 0] - data[0, 0])
+        elif data.shape[0] > 2:  # over-determined
+            data = data - X0
+            # first principal component
+            _, _, v = np.linalg.svd(data)
+            theta = np.arctan2(v[0, 1], v[0, 0])
+        else:  # under-determined
+            raise ValueError('At least 2 input points needed.')
+
+        # angle perpendicular to line angle
+        theta = (theta + np.pi / 2) % np.pi
+        # line always passes through mean
+        dist = X0[0] * math.cos(theta) + X0[1] * math.sin(theta)
+
+        self.params = (dist, theta)
+
+        return True
+
+    def residuals(self, data):
+        """Determine residuals of data to model.
+
+        For each point the shortest distance to the line is returned.
+
+        Parameters
+        ----------
+        data : (N, 2) array
+            N points with ``(x, y)`` coordinates, respectively.
+
+        Returns
+        -------
+        residuals : (N, ) array
+            Residual for each data point.
+
+        """
+
+        _check_data_dim(data, dim=2)
+
+        dist, theta = self.params
+
+        x = data[:, 0]
+        y = data[:, 1]
+
+        return dist - (x * math.cos(theta) + y * math.sin(theta))
+
+    def predict_x(self, y, params=None):
+        """Predict x-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        y : array
+            y-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        x : array
+            Predicted x-coordinates.
+
+        """
+
+        if params is None:
+            params = self.params
+        dist, theta = params
+        return (dist - y * math.sin(theta)) / math.cos(theta)
+
+    def predict_y(self, x, params=None):
+        """Predict y-coordinates using the estimated model.
+
+        Parameters
+        ----------
+        x : array
+            x-coordinates.
+        params : (2, ) array, optional
+            Optional custom parameter set.
+
+        Returns
+        -------
+        y : array
+            Predicted y-coordinates.
+
+        """
+
+        if params is None:
+            params = self.params
+        dist, theta = params
+        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).
+
+    Attributes
+    ----------
+    params : tuple
+        Line model parameters in the following order `origin`, `direction`.
         These parameters correspond to the vector equation
-            X = X0 + lambda * direction
+            X = origin + lambda * direction
     """
 
     def estimate(self, data):
@@ -79,16 +207,7 @@ class LineModel(BaseModel):
         else:  # under-determined
             raise ValueError('At least 2 input points needed.')
 
-        self.new_params = (X0, u)
-
-        # legacy LineModel (2D case)
-        theta = np.arctan2(u[1], u[0])
-        # angle perpendicular to line angle
-        theta = (theta + np.pi / 2) % np.pi
-        # line always passes through mean
-        dist = X0[0] * math.cos(theta) + X0[1] * math.sin(theta)
-
-        self.params = (dist, theta)
+        self.params = (X0, u)
 
         return True
 
@@ -108,23 +227,12 @@ class LineModel(BaseModel):
         residuals : (N, ) array
             Residual for each data point.
         """
-        if self.new_params is None:
-            self.new_params = self._params_from_polar(self.params)
 
-        X0, u = self.new_params
+        X0, u = self.params
         return np.linalg.norm((data - X0) -
                               np.dot(data - X0, u)[..., np.newaxis] * u, axis=1)
 
-    def _params_from_polar(self, params):
-        (dist, theta) = params
-        u = np.array([math.cos(theta - np.pi / 2), math.sin(theta - np.pi / 2)])
-        if math.cos(theta) == 0:
-            X0 = np.array([0, dist / math.sin(theta)])
-        else:
-            X0 = np.array([dist / math.cos(theta), 0])
-        return X0, u
-
-    def predict(self, x, axis=0, params=None, new_params=None):
+    def predict(self, x, axis=0, params=None):
         """Predict intersection of the estimated line model with a hyperplane
         orthogonal to a given axis.
 
@@ -146,13 +254,11 @@ class LineModel(BaseModel):
 
         If the line is parallel to the given axis, a ValueError is raised.
         """
-        if new_params is None:
-            if params is None and self.new_params is not None:
-                new_params = self.new_params
-            else:
-                new_params = self._params_from_polar(params or self.params)
 
-        X0, u = new_params
+        if params is None:
+            params = self.params
+
+        X0, u = params
 
         if u[axis] == 0:
             # line parallel to axis
@@ -162,7 +268,7 @@ class LineModel(BaseModel):
         return X0 + l[..., np.newaxis] * u
 
     def predict_x(self, y, params=None, new_params=None):
-        """Predict x-coordinates using the estimated model.
+        """Predict x-coordinates for 2D lines using the estimated model.
 
         Alias for predict(y, axis=1)[:, 0].
 
@@ -171,9 +277,7 @@ class LineModel(BaseModel):
         y : array
             y-coordinates.
         params : (2, ) array, optional
-            Optional custom parameter set in the form (`dist`, `theta`).
-        new_params : (2, ) array, optional
-            Optional custom parameter set in the form (`X0`, `direction`).
+            Optional custom parameter set in the form (`origin`, `direction`).
 
         Returns
         -------
@@ -181,11 +285,10 @@ class LineModel(BaseModel):
             Predicted x-coordinates.
 
         """
-        return self.predict(y, axis=1, params=params,
-                            new_params=new_params)[:, 0]
+        return self.predict(y, axis=1, params=params)[:, 0]
 
-    def predict_y(self, x, params=None, new_params=None):
-        """Predict y-coordinates using the estimated model.
+    def predict_y(self, x, params=None):
+        """Predict y-coordinates  for 2D lines using the estimated model.
 
         Alias for predict(x, axis=1)[:, 1].
 
@@ -194,9 +297,7 @@ class LineModel(BaseModel):
         x : array
             x-coordinates.
         params : (2, ) array, optional
-            Optional custom parameter set in the form (`dist`, `theta`).
-        new_params : (2, ) array, optional
-            Optional custom parameter set in the form (`X0`, `direction`).
+            Optional custom parameter set in the form (`origin`, `direction`).
 
         Returns
         -------
@@ -204,8 +305,7 @@ class LineModel(BaseModel):
             Predicted y-coordinates.
 
         """
-        return self.predict(x, axis=0, params=params,
-                            new_params=new_params)[:, 1]
+        return self.predict(x, axis=0, params=params)[:, 1]
 
 
 class CircleModel(BaseModel):
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index ab518126..86791bd7 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/skimage/measure/tests/test_fit.py
@@ -1,13 +1,13 @@
 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
 
 
 def test_line_model_invalid_input():
-    assert_raises(ValueError, LineModel().estimate, np.empty((5, 1)))
+    assert_raises(ValueError, LineModel().estimate, np.empty((5, 3)))
 
 
 def test_line_model_predict():
@@ -42,11 +42,10 @@ def test_line_model_residuals():
     model = LineModel()
     model.params = (0, 0)
     assert_equal(abs(model.residuals(np.array([[0, 0]]))), 0)
-    assert_almost_equal(abs(model.residuals(np.array([[0, 10]]))), 0)
+    assert_equal(abs(model.residuals(np.array([[0, 10]]))), 0)
     assert_equal(abs(model.residuals(np.array([[10, 0]]))), 10)
-    model = LineModel()
     model.params = (5, np.pi / 4)
-    assert_almost_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
+    assert_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
     assert_almost_equal(abs(model.residuals(np.array([[np.sqrt(50), 0]]))), 0)
 
 
@@ -56,45 +55,63 @@ def test_line_model_under_determined():
     data = np.empty((1, 3))
     assert_raises(ValueError, LineModel().estimate, data)
 
-def test_line_model3D_estimate():
+
+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 = LineModel()
-    model0.new_params = (np.array([0,0,0], dtype='float'),
+    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.new_params[0] +
-             10 * np.arange(-100,100)[...,np.newaxis] * model0.new_params[1])
+    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 = LineModel()
+    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.new_params[1],
-                                                model_est.new_params[1])), 0, 1)
+    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.new_params[0] - model0.new_params[0]
+    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.new_params[1], a)), 0, 1)
+    assert_almost_equal(np.linalg.norm(np.cross(model0.params[1], a)), 0, 1)
 
 
-def test_line_model3D_residuals():
-    model = LineModel()
-    model.new_params = (np.array([0,0,0]), np.array([0,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)))
 

From 61b7c20e033838f88f6cc6a4b7e8b7cc900dc50e Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Thu, 3 Dec 2015 22:16:56 +0000
Subject: [PATCH 6/8] cleaned tests and documentation, added deprecation note
 in doc for LineModel.

---
 skimage/measure/fit.py            | 6 +++---
 skimage/measure/tests/test_fit.py | 2 --
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index f4a8df37..c6097c57 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -44,6 +44,8 @@ class LineModel(BaseModel):
 
     A minimum number of 2 points is required to solve for the parameters.
 
+    **Deprecated, use ``LineModelND`` instead.**
+
     Attributes
     ----------
     params : tuple
@@ -243,9 +245,7 @@ class LineModelND(BaseModel):
         axis : int
             axis orthogonal to the hyperplane intersecting the line.
         params : (2, ) array, optional
-            Optional custom parameter set in the form (`dist`, `theta`).
-        new_params : (2, ) array, optional
-            Optional custom parameter set in the form (`X0`, `direction`).
+            Optional custom parameter set in the form (`origin`, `direction`).
 
         Returns
         -------
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index 86791bd7..af34cc29 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/skimage/measure/tests/test_fit.py
@@ -52,8 +52,6 @@ def test_line_model_residuals():
 def test_line_model_under_determined():
     data = np.empty((1, 2))
     assert_raises(ValueError, LineModel().estimate, data)
-    data = np.empty((1, 3))
-    assert_raises(ValueError, LineModel().estimate, data)
 
 
 def test_line_modelND_invalid_input():

From 6f9a55c9873f7596487bf6c1ff837acc7172d8ee Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Fri, 4 Dec 2015 08:48:08 +0000
Subject: [PATCH 7/8] added deprecation warning for LineModel and added note in
 TODO.txt

---
 TODO.txt               | 4 ++++
 skimage/measure/fit.py | 6 ++++++
 2 files changed, 10 insertions(+)

diff --git a/TODO.txt b/TODO.txt
index 2b890b89..02bd9ead 100644
--- a/TODO.txt
+++ b/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
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index c6097c57..bb31928c 100644
--- a/skimage/measure/fit.py
+++ b/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):
@@ -53,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.
 

From 0bdf4251c6f8971864f75c629ac3e3bb6df501de Mon Sep 17 00:00:00 2001
From: Kevin Keraudren <kevin@klarismo.com>
Date: Fri, 4 Dec 2015 22:25:13 +0000
Subject: [PATCH 8/8] documentation updates; added "3D" to example title

---
 doc/examples/plot_ransac3D.py |  6 +++---
 skimage/measure/fit.py        | 18 ++++++++++++------
 2 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py
index b7816e47..644cf154 100644
--- a/doc/examples/plot_ransac3D.py
+++ b/doc/examples/plot_ransac3D.py
@@ -1,7 +1,7 @@
 """
-=========================================
-Robust line model estimation using RANSAC
-=========================================
+============================================
+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.
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index bb31928c..5ab66251 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -45,7 +45,7 @@ class LineModel(BaseModel):
 
     A minimum number of 2 points is required to solve for the parameters.
 
-    **Deprecated, use ``LineModelND`` instead.**
+    **Deprecated class**. Use ``LineModelND`` instead.
 
     Attributes
     ----------
@@ -172,14 +172,16 @@ class LineModel(BaseModel):
 class LineModelND(BaseModel):
     """Total least squares estimator for N-dimensional lines.
 
-    Lines are defined by a point (origin) and a unit vector (direction).
+    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`.
-        These parameters correspond to the vector equation
-            X = origin + lambda * direction
+
     """
 
     def estimate(self, data):
@@ -276,7 +278,9 @@ class LineModelND(BaseModel):
     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].
+        Alias for::
+
+            predict(y, axis=1)[:, 0]
 
         Parameters
         ----------
@@ -296,7 +300,9 @@ class LineModelND(BaseModel):
     def predict_y(self, x, params=None):
         """Predict y-coordinates  for 2D lines using the estimated model.
 
-        Alias for predict(x, axis=1)[:, 1].
+        Alias for::
+
+            predict(x, axis=0)[:, 1]
 
         Parameters
         ----------