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