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)))