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