mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-17 11:32:45 +08:00
we now have the old LineModel with params (dist, theta), and the new LineModelND with params (origin, direction).
This commit is contained in:
Kevin Keraudren
@@ -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
|
||||
|
||||
|
||||
@@ -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',
|
||||
|
||||
+153
-53
@@ -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):
|
||||
|
||||
@@ -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)))
|
||||
|
||||
|
||||
Reference in New Issue
Block a user