wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-10 23:57:34 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #440 from ahojnnes/fitting

Browse Source 
ENH: Line, Circle, Ellipse total least squares fitting and RANSAC algorithm.
This commit is contained in:
Stefan van der Walt
2013-05-07 04:34:54 -07:00
parent 4d7905951b 219d621725
commit 2eae93124a
6 changed files with 1085 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_matching.py
+140
View File
@@ -0,0 +1,140 @@
"""
============================
Robust matching using RANSAC
============================
In this simplified example we first generate two synthetic images as if they
were taken from different view points.
In the next step we find interest points in both images and find
correspondencies based on a weighted sum of squared differences of a small
neighbourhood around them. Note, that this measure is only robust towards
linear radiometric and not geometric distortions and is thus only usable with
slight view point changes.
After finding the correspondencies we end up having a set of source and
destination coordinates which can be used to estimate the geometric
transformation between both images. However, many of the correspondencies are
faulty and simply estimating the parameter set with all coordinates is not
sufficient. Therefore, the RANSAC algorithm is used on top of the normal model
to robustly estimate the parameter set by detecting outliers.
"""
import numpy as np
from matplotlib import pyplot as plt
from skimage import data
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.transform import warp, AffineTransform
from skimage.exposure import rescale_intensity
from skimage.color import rgb2gray
from skimage.measure import ransac
# generate synthetic checkerboard image and add gradient for the later matching
checkerboard = data.checkerboard()
img_orig = np.zeros(list(checkerboard.shape) + [3])
img_orig[..., 0] = checkerboard
gradient_r, gradient_c = np.mgrid[0:img_orig.shape[0], 0:img_orig.shape[1]]
img_orig[..., 1] = gradient_r
img_orig[..., 2] = gradient_c
img_orig = rescale_intensity(img_orig)
img_orig_gray = rgb2gray(img_orig)
# warp synthetic image
tform = AffineTransform(scale=(0.9, 0.9), rotation=0.2, translation=(20, -10))
img_warped = warp(img_orig, tform.inverse, output_shape=(200, 200))
img_warped_gray = rgb2gray(img_warped)
# extract corners using Harris' corner measure
coords_orig = corner_peaks(corner_harris(img_orig_gray), threshold_rel=0.001,
min_distance=5)
coords_warped = corner_peaks(corner_harris(img_warped_gray),
threshold_rel=0.001, min_distance=5)
# determine subpixel corner position
coords_orig_subpix = corner_subpix(img_orig_gray, coords_orig, window_size=10)
coords_warped_subpix = corner_subpix(img_warped_gray, coords_warped,
window_size=10)
def gaussian_weights(window_ext, sigma=1):
y, x = np.mgrid[-window_ext:window_ext+1, -window_ext:window_ext+1]
g = np.zeros(y.shape, dtype=np.double)
g[:] = np.exp(-0.5 * (x**2 / sigma**2 + y**2 / sigma**2))
g /= 2 * np.pi * sigma * sigma
return g
def match_corner(coord, window_ext=5):
r, c = np.round(coord)
window_orig = img_orig[r-window_ext:r+window_ext+1,
c-window_ext:c+window_ext+1, :]
# weight pixels depending on distance to center pixel
weights = gaussian_weights(window_ext, 3)
weights = np.dstack((weights, weights, weights))
# compute sum of squared differences to all corners in warped image
SSDs = []
for cr, cc in coords_warped:
window_warped = img_warped[cr-window_ext:cr+window_ext+1,
cc-window_ext:cc+window_ext+1, :]
SSD = np.sum(weights * (window_orig - window_warped)**2)
SSDs.append(SSD)
# use corner with minimum SSD as correspondency
min_idx = np.argmin(SSDs)
return coords_warped_subpix[min_idx]
# find correspondencies using simple weighted sum of squared differences
src = []
dst = []
for coord in coords_orig_subpix:
src.append(coord)
dst.append(match_corner(coord))
src = np.array(src)
dst = np.array(dst)
# estimate affine transform model using all coordinates
model = AffineTransform()
model.estimate(src, dst)
# robustly estimate affine transform model with RANSAC
model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3,
residual_threshold=2, max_trials=100)
outliers = inliers == False
# compare "true" and estimated transform parameters
print tform.scale, tform.translation, tform.rotation
print model.scale, model.translation, model.rotation
print model_robust.scale, model_robust.translation, model_robust.rotation
# visualize correspondencies
img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1)
fig, ax = plt.subplots(nrows=2, ncols=1)
plt.gray()
ax[0].imshow(img_combined, interpolation='nearest')
ax[0].axis('off')
ax[0].axis((0, 400, 200, 0))
ax[0].set_title('Correct correspondencies')
ax[1].imshow(img_combined, interpolation='nearest')
ax[1].axis('off')
ax[1].axis((0, 400, 200, 0))
ax[1].set_title('Faulty correspondencies')
for ax_idx, (m, color) in enumerate(((inliers, 'g'), (outliers, 'r'))):
ax[ax_idx].plot((src[m, 1], dst[m, 1] + 200), (src[m, 0], dst[m, 0]), '-',
color=color)
ax[ax_idx].plot(src[m, 1], src[m, 0], '.', markersize=10, color=color)
ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10,
color=color)
plt.show()
doc/examples/plot_ransac.py
+55
View File
@@ -0,0 +1,55 @@
"""
=========================================
Robust line model estimation using RANSAC
=========================================
In this example we see how to robustly fit a line model to faulty data using
the RANSAC algorithm.
"""
import numpy as np
from matplotlib import pyplot as plt
from skimage.measure import LineModel, ransac
np.random.seed(seed=1)
# generate coordinates of line
x = np.arange(-200, 200)
y = 0.2 * x + 20
data = np.column_stack([x, y])
# add faulty data
faulty = np.array(30 * [(180, -100)])
faulty += 5 * np.random.normal(size=faulty.shape)
data[:faulty.shape[0]] = faulty
# add gaussian noise to coordinates
noise = np.random.normal(size=data.shape)
data += 0.5 * noise
data[::2] += 5 * noise[::2]
data[::4] += 20 * noise[::4]
# fit line using all data
model = LineModel()
model.estimate(data)
# robustly fit line only using inlier data with RANSAC algorithm
model_robust, inliers = ransac(data, LineModel, min_samples=2,
residual_threshold=1, max_trials=1000)
outliers = inliers == False
# generate coordinates of estimated models
line_x = np.arange(-250, 250)
line_y = model.predict_y(line_x)
line_y_robust = model_robust.predict_y(line_x)
plt.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
label='Inlier data')
plt.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
label='Outlier data')
plt.plot(line_x, line_y, '-k', label='Line model from all data')
plt.plot(line_x, line_y_robust, '-b', label='Robust line model')
plt.legend(loc='lower left')
plt.show()
skimage/measure/__init__.py
+7 -1
View File
@@ -2,10 +2,16 @@ from .find_contours import find_contours
from ._regionprops import regionprops, perimeter
from ._structural_similarity import structural_similarity
from ._polygon import approximate_polygon, subdivide_polygon
from .fit import LineModel, CircleModel, EllipseModel, ransac
__all__ = ['find_contours',
'regionprops',
'perimeter',
'structural_similarity',
'approximate_polygon',
'subdivide_polygon']
'subdivide_polygon',
'LineModel',
'CircleModel',
'EllipseModel',
'ransac']
skimage/measure/fit.py
+653
View File
@@ -0,0 +1,653 @@
import math
import numpy as np
from scipy import optimize
def _check_data_dim(data, dim):
if data.ndim != 2 or data.shape[1] != dim:
raise ValueError('Input data must have shape (N, %d).' % dim)
class BaseModel(object):
def __init__(self):
self._params = None
class LineModel(BaseModel):
"""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) }
The ``_params`` attribute contains the parameters in the following order::
dist, theta
A minimum number of 2 points is required to solve for the parameters.
"""
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.
"""
_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)
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 CircleModel(BaseModel):
"""Total least squares estimator for 2D circles.
The functional model of the circle is::
r**2 = (x - xc)**2 + (y - yc)**2
This estimator minimizes the squared distances from all points to the
circle::
min{ sum((r - sqrt((x_i - xc)**2 + (y_i - yc)**2))**2) }
The ``_params`` attribute contains the parameters in the following order::
xc, yc, r
A minimum number of 3 points is required to solve for the parameters.
"""
def estimate(self, data):
"""Estimate circle model from data using total least squares.
Parameters
----------
data : (N, 2) array
N points with ``(x, y)`` coordinates, respectively.
"""
_check_data_dim(data, dim=2)
x = data[:, 0]
y = data[:, 1]
# pre-allocate jacobian for all iterations
A = np.zeros((3, data.shape[0]), dtype=np.double)
# same for all iterations: r
A[2, :] = -1
def dist(xc, yc):
return np.sqrt((x - xc)**2 + (y - yc)**2)
def fun(params):
xc, yc, r = params
return dist(xc, yc) - r
def Dfun(params):
xc, yc, r = params
d = dist(xc, yc)
A[0, :] = -(x - xc) / d
A[1, :] = -(y - yc) / d
# same for all iterations, so not changed in each iteration
#A[2, :] = -1
return A
xc0 = x.mean()
yc0 = y.mean()
r0 = dist(xc0, yc0).mean()
params0 = (xc0, yc0, r0)
params, _ = optimize.leastsq(fun, params0, Dfun=Dfun, col_deriv=True)
self._params = params
def residuals(self, data):
"""Determine residuals of data to model.
For each point the shortest distance to the circle 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)
xc, yc, r = self._params
x = data[:, 0]
y = data[:, 1]
return r - np.sqrt((x - xc)**2 + (y - yc)**2)
def predict_xy(self, t, params=None):
"""Predict x- and y-coordinates using the estimated model.
Parameters
----------
t : array
Angles in circle in radians. Angles start to count from positive
x-axis to positive y-axis in a right-handed system.
params : (3, ) array, optional
Optional custom parameter set.
Returns
-------
xy : (..., 2) array
Predicted x- and y-coordinates.
"""
if params is None:
params = self._params
xc, yc, r = params
x = xc + r * np.cos(t)
y = yc + r * np.sin(t)
return np.concatenate((x[..., None], y[..., None]), axis=t.ndim)
class EllipseModel(BaseModel):
"""Total least squares estimator for 2D ellipses.
The functional model of the ellipse is::
xt = xc + a*cos(theta)*cos(t) - b*sin(theta)*sin(t)
yt = yc + a*sin(theta)*cos(t) + b*cos(theta)*sin(t)
d = sqrt((x - xt)**2 + (y - yt)**2)
where ``(xt, yt)`` is the closest point on the ellipse to ``(x, y)``. Thus
d is the shortest distance from the point to the ellipse.
This estimator minimizes the squared distances from all points to the
ellipse::
min{ sum(d_i**2) } = min{ sum((x_i - xt)**2 + (y_i - yt)**2) }
Thus you have ``2 * N`` equations (x_i, y_i) for ``N + 5`` unknowns (t_i,
xc, yc, a, b, theta), which gives you an effective redundancy of ``N - 5``.
The ``_params`` attribute contains the parameters in the following order::
xc, yc, a, b, theta
A minimum number of 5 points is required to solve for the parameters.
"""
def estimate(self, data):
"""Estimate circle model from data using total least squares.
Parameters
----------
data : (N, 2) array
N points with ``(x, y)`` coordinates, respectively.
"""
_check_data_dim(data, dim=2)
x = data[:, 0]
y = data[:, 1]
N = data.shape[0]
# pre-allocate jacobian for all iterations
A = np.zeros((N + 5, 2 * N), dtype=np.double)
# same for all iterations: xc, yc
A[0, :N] = -1
A[1, N:] = -1
diag_idxs = np.diag_indices(N)
def fun(params):
xyt = self.predict_xy(params[5:], params[:5])
fx = x - xyt[:, 0]
fy = y - xyt[:, 1]
return np.append(fx, fy)
def Dfun(params):
xc, yc, a, b, theta = params[:5]
t = params[5:]
ct = np.cos(t)
st = np.sin(t)
ctheta = math.cos(theta)
stheta = math.sin(theta)
# derivatives for fx, fy in the following order:
# xc, yc, a, b, theta, t_i
# fx
A[2, :N] = - ctheta * ct
A[3, :N] = stheta * st
A[4, :N] = a * stheta * ct + b * ctheta * st
A[5:, :N][diag_idxs] = a * ctheta * st + b * stheta * ct
# fy
A[2, N:] = - stheta * ct
A[3, N:] = - ctheta * st
A[4, N:] = - a * ctheta * ct + b * stheta * st
A[5:, N:][diag_idxs] = a * stheta * st - b * ctheta * ct
return A
# initial guess of parameters using a circle model
params0 = np.empty((N + 5, ), dtype=np.double)
xc0 = x.mean()
yc0 = y.mean()
r0 = np.sqrt((x - xc0)**2 + (y - yc0)**2).mean()
params0[:5] = (xc0, yc0, r0, 0, 0)
params0[5:] = np.arctan2(y - yc0, x - xc0)
params, _ = optimize.leastsq(fun, params0, Dfun=Dfun, col_deriv=True)
self._params = params[:5]
def residuals(self, data):
"""Determine residuals of data to model.
For each point the shortest distance to the ellipse 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)
xc, yc, a, b, theta = self._params
ctheta = math.cos(theta)
stheta = math.sin(theta)
x = data[:, 0]
y = data[:, 1]
N = data.shape[0]
def fun(t, xi, yi):
ct = math.cos(t)
st = math.sin(t)
xt = xc + a * ctheta * ct - b * stheta * st
yt = yc + a * stheta * ct + b * ctheta * st
return (xi - xt)**2 + (yi - yt)**2
# def Dfun(t, xi, yi):
# ct = math.cos(t)
# st = math.sin(t)
# xt = xc + a * ctheta * ct - b * stheta * st
# yt = yc + a * stheta * ct + b * ctheta * st
# dfx_t = - 2 * (xi - xt) * (- a * ctheta * st
# - b * stheta * ct)
# dfy_t = - 2 * (yi - yt) * (- a * stheta * st
# + b * ctheta * ct)
# return [dfx_t + dfy_t]
residuals = np.empty((N, ), dtype=np.double)
# initial guess for parameter t of closest point on ellipse
t0 = np.arctan2(y - yc, x - xc) - theta
# determine shortest distance to ellipse for each point
for i in range(N):
xi = x[i]
yi = y[i]
# faster without Dfun, because of the python overhead
t, _ = optimize.leastsq(fun, t0[i], args=(xi, yi))
residuals[i] = np.sqrt(fun(t, xi, yi))
return residuals
def predict_xy(self, t, params=None):
"""Predict x- and y-coordinates using the estimated model.
Parameters
----------
t : array
Angles in circle in radians. Angles start to count from positive
x-axis to positive y-axis in a right-handed system.
params : (5, ) array, optional
Optional custom parameter set.
Returns
-------
xy : (..., 2) array
Predicted x- and y-coordinates.
"""
if params is None:
params = self._params
xc, yc, a, b, theta = params
ct = np.cos(t)
st = np.sin(t)
ctheta = math.cos(theta)
stheta = math.sin(theta)
x = xc + a * ctheta * ct - b * stheta * st
y = yc + a * stheta * ct + b * ctheta * st
return np.concatenate((x[..., None], y[..., None]), axis=t.ndim)
def ransac(data, model_class, min_samples, residual_threshold,
is_data_valid=None, is_model_valid=None,
max_trials=100, stop_sample_num=np.inf, stop_residuals_sum=0):
"""Fit a model to data with the RANSAC (random sample consensus) algorithm.
RANSAC is an iterative algorithm for the robust estimation of parameters
from a subset of inliers from the complete data set. Each iteration
performs the following tasks:
1. Select `min_samples` random samples from the original data and check
whether the set of data is valid (see `is_data_valid`).
2. Estimate a model to the random subset
(`model_cls.estimate(*data[random_subset]`) and check whether the
estimated model is valid (see `is_model_valid`).
3. Classify all data as inliers or outliers by calculating the residuals
to the estimated model (`model_cls.residuals(*data)`) - all data samples
with residuals smaller than the `residual_threshold` are considered as
inliers.
4. Save estimated model as best model if number of inlier samples is
maximal. In case the current estimated model has the same number of
inliers, it is only considered as the best model if it has less sum of
residuals.
These steps are performed either a maximum number of times or until one of
the special stop criteria are met. The final model is estimated using all
inlier samples of the previously determined best model.
Parameters
----------
data : [list, tuple of] (N, D) array
Data set to which the model is fitted, where N is the number of data
points and D the dimensionality of the data.
If the model class requires multiple input data arrays (e.g. source and
destination coordinates of ``skimage.transform.AffineTransform``),
they can be optionally passed as tuple or list. Note, that in this case
the functions ``estimate(*data)``, ``residuals(*data)``,
``is_model_valid(model, *random_data)`` and
``is_data_valid(*random_data)`` must all take each data array as
separate arguments.
model_class : object
Object with the following object methods:
* ``estimate(*data)``
* ``residuals(*data)``
min_samples : int
The minimum number of data points to fit a model to.
residual_threshold : float
Maximum distance for a data point to be classified as an inlier.
is_data_valid : function, optional
This function is called with the randomly selected data before the
model is fitted to it: `is_data_valid(*random_data)`.
is_model_valid : function, optional
This function is called with the estimated model and the randomly
selected data: `is_model_valid(model, *random_data)`, .
max_trials : int, optional
Maximum number of iterations for random sample selection.
stop_sample_num : int, optional
Stop iteration if at least this number of inliers are found.
stop_residuals_sum : float, optional
Stop iteration if sum of residuals is less equal than this threshold.
Returns
-------
model : object
Best model with largest consensus set.
inliers : (N, ) array
Boolean mask of inliers classified as ``True``.
References
----------
.. [1] "RANSAC", Wikipedia, http://en.wikipedia.org/wiki/RANSAC
Examples
--------
Generate ellipse data without tilt and add noise:
>>> t = np.linspace(0, 2 * np.pi, 50)
>>> a = 5
>>> b = 10
>>> xc = 20
>>> yc = 30
>>> x = xc + a * np.cos(t)
>>> y = yc + b * np.sin(t)
>>> data = np.column_stack([x, y])
>>> np.random.seed(seed=1234)
>>> data += np.random.normal(size=data.shape)
Add some faulty data:
>>> data[0] = (100, 100)
>>> data[1] = (110, 120)
>>> data[2] = (120, 130)
>>> data[3] = (140, 130)
Estimate ellipse model using all available data:
>>> model = EllipseModel()
>>> model.estimate(data)
>>> model._params
array([ 4.85808595e+02, 4.51492793e+02, 1.15018491e+03,
5.52428289e+00, 7.32420126e-01])
Estimate ellipse model using RANSAC:
>>> ransac_model, inliers = ransac(data, EllipseModel, 5, 3, max_trials=50)
>>> # ransac_model._params, inliers
Should give the correct result estimated without the faulty data:
[ 20.12762373, 29.73563061, 4.81499637, 10.4743584, 0.05217117]
[ 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
Robustly estimate geometric transformation:
>>> from skimage.transform import SimilarityTransform
>>> src = 100 * np.random.random((50, 2))
>>> model0 = SimilarityTransform(scale=0.5, rotation=1,
... translation=(10, 20))
>>> dst = model0(src)
>>> dst[0] = (10000, 10000)
>>> dst[1] = (-100, 100)
>>> dst[2] = (50, 50)
>>> model, inliers = ransac((src, dst), SimilarityTransform, 2, 10)
>>> inliers
array([ 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49])
"""
best_model = None
best_inlier_num = 0
best_inlier_residuals_sum = np.inf
best_inliers = None
if not isinstance(data, list) and not isinstance(data, tuple):
data = [data]
# make sure data is list and not tuple, so it can be modified below
data = list(data)
# number of samples
N = data[0].shape[0]
for _ in range(max_trials):
# choose random sample set
samples = []
random_idxs = np.random.randint(0, N, min_samples)
for d in data:
samples.append(d[random_idxs])
# check if random sample set is valid
if is_data_valid is not None and not is_data_valid(*samples):
continue
# estimate model for current random sample set
sample_model = model_class()
sample_model.estimate(*samples)
# check if estimated model is valid
if is_model_valid is not None and not is_model_valid(sample_model,
*samples):
continue
sample_model_residuals = np.abs(sample_model.residuals(*data))
# consensus set / inliers
sample_model_inliers = sample_model_residuals < residual_threshold
sample_model_residuals_sum = np.sum(sample_model_residuals**2)
# choose as new best model if number of inliers is maximal
sample_inlier_num = np.sum(sample_model_inliers)
if (
# more inliers
sample_inlier_num > best_inlier_num
# same number of inliers but less "error" in terms of residuals
or (sample_inlier_num == best_inlier_num
and sample_model_residuals_sum < best_inlier_residuals_sum)
):
best_model = sample_model
best_inlier_num = sample_inlier_num
best_inlier_residuals_sum = sample_model_residuals_sum
best_inliers = sample_model_inliers
if (
best_inlier_num >= stop_sample_num
or best_inlier_residuals_sum <= stop_residuals_sum
):
break
# estimate final model using all inliers
if best_inliers is not None:
# select inliers for each data array
for i in range(len(data)):
data[i] = data[i][best_inliers]
best_model.estimate(*data)
return best_model, best_inliers
skimage/measure/tests/test_fit.py
+208
View File
@@ -0,0 +1,208 @@
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.transform import AffineTransform
def test_line_model_invalid_input():
assert_raises(ValueError, LineModel().estimate, np.empty((5, 3)))
def test_line_model_predict():
model = LineModel()
model._params = (10, 1)
x = np.arange(-10, 10)
y = model.predict_y(x)
assert_almost_equal(x, model.predict_x(y))
def test_line_model_estimate():
# generate original data without noise
model0 = LineModel()
model0._params = (10, 1)
x0 = np.arange(-100, 100)
y0 = model0.predict_y(x0)
data0 = np.column_stack([x0, y0])
# 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.estimate(data)
# test whether estimated parameters almost equal original parameters
assert_almost_equal(model0._params, model_est._params, 1)
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_equal(abs(model.residuals(np.array([[10, 0]]))), 10)
model._params = (5, np.pi / 4)
assert_equal(abs(model.residuals(np.array([[0, 0]]))), 5)
assert_equal(abs(model.residuals(np.array([[np.sqrt(50), 0]]))), 5)
def test_line_model_under_determined():
data = np.empty((1, 2))
assert_raises(ValueError, LineModel().estimate, data)
def test_circle_model_invalid_input():
assert_raises(ValueError, CircleModel().estimate, np.empty((5, 3)))
def test_circle_model_predict():
model = CircleModel()
r = 5
model._params = (0, 0, r)
t = np.arange(0, 2 * np.pi, np.pi / 2)
xy = np.array(((5, 0), (0, 5), (-5, 0), (0, -5)))
assert_almost_equal(xy, model.predict_xy(t))
def test_circle_model_estimate():
# generate original data without noise
model0 = CircleModel()
model0._params = (10, 12, 3)
t = np.linspace(0, 2 * np.pi, 1000)
data0 = model0.predict_xy(t)
# add gaussian noise to data
np.random.seed(1234)
data = data0 + np.random.normal(size=data0.shape)
# estimate parameters of noisy data
model_est = CircleModel()
model_est.estimate(data)
# test whether estimated parameters almost equal original parameters
assert_almost_equal(model0._params, model_est._params, 1)
def test_circle_model_residuals():
model = CircleModel()
model._params = (0, 0, 5)
assert_almost_equal(abs(model.residuals(np.array([[5, 0]]))), 0)
assert_almost_equal(abs(model.residuals(np.array([[6, 6]]))),
np.sqrt(2 * 6**2) - 5)
assert_almost_equal(abs(model.residuals(np.array([[10, 0]]))), 5)
def test_ellipse_model_invalid_input():
assert_raises(ValueError, EllipseModel().estimate, np.empty((5, 3)))
def test_ellipse_model_predict():
model = EllipseModel()
r = 5
model._params = (0, 0, 5, 10, 0)
t = np.arange(0, 2 * np.pi, np.pi / 2)
xy = np.array(((5, 0), (0, 10), (-5, 0), (0, -10)))
assert_almost_equal(xy, model.predict_xy(t))
def test_ellipse_model_estimate():
# generate original data without noise
model0 = EllipseModel()
model0._params = (10, 20, 15, 25, 0)
t = np.linspace(0, 2 * np.pi, 100)
data0 = model0.predict_xy(t)
# add gaussian noise to data
np.random.seed(1234)
data = data0 + np.random.normal(size=data0.shape)
# estimate parameters of noisy data
model_est = EllipseModel()
model_est.estimate(data)
# test whether estimated parameters almost equal original parameters
assert_almost_equal(model0._params, model_est._params, 0)
def test_line_model_residuals():
model = EllipseModel()
# vertical line through origin
model._params = (0, 0, 10, 5, 0)
assert_almost_equal(abs(model.residuals(np.array([[10, 0]]))), 0)
assert_almost_equal(abs(model.residuals(np.array([[0, 5]]))), 0)
assert_almost_equal(abs(model.residuals(np.array([[0, 10]]))), 5)
def test_ransac_shape():
np.random.seed(1)
# generate original data without noise
model0 = CircleModel()
model0._params = (10, 12, 3)
t = np.linspace(0, 2 * np.pi, 1000)
data0 = model0.predict_xy(t)
# add some faulty data
outliers = (10, 30, 200)
data0[outliers[0], :] = (1000, 1000)
data0[outliers[1], :] = (-50, 50)
data0[outliers[2], :] = (-100, -10)
# estimate parameters of corrupted data
model_est, inliers = ransac(data0, CircleModel, 3, 5)
# test whether estimated parameters equal original parameters
assert_equal(model0._params, model_est._params)
for outlier in outliers:
assert outlier not in inliers
def test_ransac_geometric():
np.random.seed(1)
# generate original data without noise
src = 100 * np.random.random((50, 2))
model0 = AffineTransform(scale=(0.5, 0.3), rotation=1,
translation=(10, 20))
dst = model0(src)
# add some faulty data
outliers = (0, 5, 20)
dst[outliers[0]] = (10000, 10000)
dst[outliers[1]] = (-100, 100)
dst[outliers[2]] = (50, 50)
# estimate parameters of corrupted data
model_est, inliers = ransac((src, dst), AffineTransform, 2, 20)
# test whether estimated parameters equal original parameters
assert_almost_equal(model0._matrix, model_est._matrix)
assert np.all(np.nonzero(inliers == False)[0] == outliers)
def test_ransac_is_data_valid():
np.random.seed(1)
is_data_valid = lambda data: data.shape[0] > 2
model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf,
is_data_valid=is_data_valid)
assert_equal(model, None)
assert_equal(inliers, None)
def test_ransac_is_model_valid():
np.random.seed(1)
def is_model_valid(model, data):
return False
model, inliers = ransac(np.empty((10, 2)), LineModel, 2, np.inf,
is_model_valid=is_model_valid)
assert_equal(model, None)
assert_equal(inliers, None)
if __name__ == "__main__":
np.testing.run_module_suite()
skimage/transform/_geometric.py
+22
View File
@@ -41,6 +41,28 @@ class GeometricTransform(object):
"""
raise NotImplementedError()
def residuals(self, src, dst):
"""Determine residuals of transformed destination coordinates.
For each transformed source coordinate the euclidean distance to the
respective destination coordinate is determined.
Parameters
----------
src : (N, 2) array
Source coordinates.
dst : (N, 2) array
Destination coordinates.
Returns
-------
residuals : (N, ) array
Residual for coordinate.
"""
return np.sqrt(np.sum((self(src) - dst)**2, axis=1))
def __add__(self, other):
"""Combine this transformation with another.