From 040a53456d933f8587dbdcfb093b3bc03db2c5e9 Mon Sep 17 00:00:00 2001 From: Kevin Keraudren Date: Wed, 2 Dec 2015 08:58:16 +0000 Subject: [PATCH] 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)))