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/doc/examples/plot_ransac3D.py b/doc/examples/plot_ransac3D.py new file mode 100644 index 00000000..644cf154 --- /dev/null +++ b/doc/examples/plot_ransac3D.py @@ -0,0 +1,40 @@ +""" +============================================ +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. + +""" +import numpy as np +from matplotlib import pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from skimage.measure import LineModelND, 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, LineModelND, 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..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 f8e65117..5ab66251 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): @@ -9,6 +10,11 @@ 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): @@ -39,6 +45,8 @@ class LineModel(BaseModel): A minimum number of 2 points is required to solve for the parameters. + **Deprecated class**. Use ``LineModelND`` instead. + Attributes ---------- params : tuple @@ -46,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. @@ -156,6 +169,157 @@ class LineModel(BaseModel): 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) + according to the following vector equation:: + + X = origin + lambda * direction + + Attributes + ---------- + params : tuple + Line model parameters in the following order `origin`, `direction`. + + """ + + def estimate(self, data): + """Estimate line model from data. + + Parameters + ---------- + data : (N, dim) array + N points in a space of dimensionality dim >= 2. + + Returns + ------- + success : bool + True, if model estimation succeeds. + """ + + _check_data_atleast_2D(data) + + X0 = data.mean(axis=0) + + if data.shape[0] == 2: # well determined + 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 + # 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, dim) array + N points in a space of dimension dim. + + 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) + + def predict(self, x, axis=0, params=None): + """Predict intersection of the estimated line model with a hyperplane + orthogonal to a given axis. + + 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 (`origin`, `direction`). + + Returns + ------- + y : array + Predicted coordinates. + + If the line is parallel to the given axis, a ValueError is raised. + """ + + if params is None: + params = self.params + + X0, u = 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 for 2D lines using the estimated model. + + Alias for:: + + predict(y, axis=1)[:, 0] + + Parameters + ---------- + y : array + y-coordinates. + params : (2, ) array, optional + Optional custom parameter set in the form (`origin`, `direction`). + + Returns + ------- + x : array + Predicted x-coordinates. + + """ + return self.predict(y, axis=1, params=params)[:, 0] + + def predict_y(self, x, params=None): + """Predict y-coordinates for 2D lines using the estimated model. + + Alias for:: + + predict(x, axis=0)[:, 1] + + Parameters + ---------- + x : array + x-coordinates. + params : (2, ) array, optional + Optional custom parameter set in the form (`origin`, `direction`). + + Returns + ------- + y : array + Predicted y-coordinates. + + """ + return self.predict(x, axis=0, params=params)[:, 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..af34cc29 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, LineModelND, CircleModel, EllipseModel, ransac from skimage.transform import AffineTransform from skimage.measure.fit import _dynamic_max_trials from skimage._shared._warnings import expected_warnings @@ -54,6 +54,62 @@ def test_line_model_under_determined(): assert_raises(ValueError, LineModel().estimate, data) +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 = 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.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 = 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.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_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)))