From e35d2399628750059ad5608003fda6718ba68855 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 27 Feb 2013 19:45:34 +0100 Subject: [PATCH 01/54] Add line estimation model --- skimage/measure/fit.py | 125 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 skimage/measure/fit.py diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py new file mode 100644 index 00000000..007a485d --- /dev/null +++ b/skimage/measure/fit.py @@ -0,0 +1,125 @@ +import numpy as np + + +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: + + 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`. + + ''' + + 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. + + ''' + + 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] * np.cos(theta) + X0[1] * np.sin(theta) + + self._params = np.array([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. + + ''' + + dist, theta = self._params + data_dists = (data[:, 0] * np.cos(theta) + data[:, 1] * np.sin(theta)) + return np.abs(dist - data_dists) + + @classmethod + def is_degenerate(cls, data): + '''Check whether set of points is degenerate. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + Returns + ------- + flag : bool + Flag indicating if data is degenerate. + + ''' + + return data.shape[0] < 2 + + def predict_x(self, y): + '''Predict x-coordinates using the estimated model. + + Parameters + ---------- + y : array + y-coordinates. + + Returns + ------- + x : array + Predicted x-coordinates. + + ''' + + dist, theta = self._params + return (dist - y * np.cos(theta)) / np.cos(theta) + + def predict_y(self, x): + '''Predict y-coordinates using the estimated model. + + Parameters + ---------- + x : array + x-coordinates. + + Returns + ------- + y : array + Predicted y-coordinates. + + ''' + + dist, theta = self._params + return (dist - x * np.cos(theta)) / np.sin(theta) From 3f650f07245a8d6d05527268ac55c4ca81874fcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 27 Feb 2013 19:45:51 +0100 Subject: [PATCH 02/54] Add RANSAC algorithm --- skimage/measure/fit.py | 68 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 007a485d..dac66e72 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -123,3 +123,71 @@ class LineModel(BaseModel): dist, theta = self._params return (dist - x * np.cos(theta)) / np.sin(theta) + + +def ransac(data, model_class, min_samples, residual_threshold, + max_trials=1000): + ''' + Fits a model to data with the RANSAC (random sample consensus) algorithm. + + Parameters + ---------- + data : (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. + model_class : object + Object with the following methods implemented: + + * `estimate(data)` + * `residuals(data)` + + min_samples : int + The minimum number of data points to fit a model. + residual_threshold : float + Maximum distance for a data point to be classified as an inlier. + max_trials : int, optional + Maximum number of iterations for random sample selection. + + Returns + ------- + model : object + Best model with largest consensus set. + inliers : (N,) array + Indices of inliers. + + ''' + + best_model = None + best_inlier_num = 0 + best_inliers = None + data_idxs = np.arange(data.shape[0]) + + for _ in range(max_trials): + + # choose random sample + sample = data[np.random.randint(0, data.shape[0], 2)] + + # check if random sample is degenerate + if model_class.is_degenerate(sample): + continue + + # create new instance of model class for current sample + sample_model = model_class() + sample_model.estimate(sample) + sample_model_residuals = sample_model.residuals(data) + # consensus set / inliers + sample_model_inliers = data_idxs[sample_model_residuals + < residual_threshold] + + # choose as new best model if number of inliers is maximal + sample_inlier_num = sample_model_inliers.shape[0] + if sample_inlier_num > best_inlier_num: + best_model = sample_model + best_inlier_num = sample_inlier_num + best_inliers = sample_model_inliers + + # estimate final model using all inliers + if best_inliers is not None: + best_model.estimate(data[best_inliers]) + + return best_model, best_inliers From 21d1e093302933b2b23661fb30655cd09a2459f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 09:34:58 +0100 Subject: [PATCH 03/54] Add circle estimator model and fix some other bugs --- skimage/measure/fit.py | 153 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 146 insertions(+), 7 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index dac66e72..f1f48009 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -1,4 +1,5 @@ import numpy as np +from scipy import optimize class BaseModel(object): @@ -11,13 +12,21 @@ class LineModel(BaseModel): '''Total least squares estimator for 2D lines. - Lines are parameterized using polar coordinates: + 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 + ''' def estimate(self, data): @@ -33,7 +42,7 @@ class LineModel(BaseModel): 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]) + 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 @@ -47,7 +56,7 @@ class LineModel(BaseModel): # line always passes through mean dist = X0[0] * np.cos(theta) + X0[1] * np.sin(theta) - self._params = np.array([dist, theta]) + self._params = (dist, theta) def residuals(self, data): '''Determine residuals of data to model. @@ -67,8 +76,11 @@ class LineModel(BaseModel): ''' dist, theta = self._params - data_dists = (data[:, 0] * np.cos(theta) + data[:, 1] * np.sin(theta)) - return np.abs(dist - data_dists) + + x = data[:, 0] + y = data[:, 1] + + return dist - (x * np.cos(theta) + y * np.sin(theta)) @classmethod def is_degenerate(cls, data): @@ -125,6 +137,133 @@ class LineModel(BaseModel): return (dist - x * np.cos(theta)) / np.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 + + ''' + + 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. + + ''' + + x = data[:, 0] + y = data[:, 1] + # pre-allocate for all iterations + A = np.empty((3, data.shape[0]), dtype=np.double) + # same for all iterations + 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 + #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 line is returned. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + Returns + ------- + residuals : (N, ) array + Residual for each data point. + + ''' + + xc, yc, r = self._params + + x = data[:, 0] + y = data[:, 1] + + return r - np.sqrt((x - xc)**2 + (y - yc)**2) + + @classmethod + def is_degenerate(cls, data): + '''Check whether set of points is degenerate. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + Returns + ------- + flag : bool + Flag indicating if data is degenerate. + + ''' + + return data.shape[0] < 2 + + def predict_xy(self, theta): + '''Predict x- and y-coordinates using the estimated model. + + Parameters + ---------- + theta : array + Angles in circle in radians. Angles start to count from positive + x-axis to positive y-axis in a right-handed system. + + Returns + ------- + x : array + Predicted x-coordinates. + y : array + Predicted y-coordinates. + + ''' + + xc, yc, r = self._params + + x = xc + r * np.cos(theta) + y = yc + r * np.sin(theta) + + return x, y + + def ransac(data, model_class, min_samples, residual_threshold, max_trials=1000): ''' @@ -165,7 +304,7 @@ def ransac(data, model_class, min_samples, residual_threshold, for _ in range(max_trials): # choose random sample - sample = data[np.random.randint(0, data.shape[0], 2)] + sample = data[np.random.randint(0, data.shape[0], min_samples)] # check if random sample is degenerate if model_class.is_degenerate(sample): @@ -176,7 +315,7 @@ def ransac(data, model_class, min_samples, residual_threshold, sample_model.estimate(sample) sample_model_residuals = sample_model.residuals(data) # consensus set / inliers - sample_model_inliers = data_idxs[sample_model_residuals + sample_model_inliers = data_idxs[np.abs(sample_model_residuals) < residual_threshold] # choose as new best model if number of inliers is maximal From 4a93f3839589f385652d9bba023e7a40e7b32902 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 17:38:38 +0100 Subject: [PATCH 04/54] Fix doc string of ransac function --- skimage/measure/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index f1f48009..fb0570e6 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -266,8 +266,7 @@ class CircleModel(BaseModel): def ransac(data, model_class, min_samples, residual_threshold, max_trials=1000): - ''' - Fits a model to data with the RANSAC (random sample consensus) algorithm. + '''Fits a model to data with the RANSAC (random sample consensus) algorithm. Parameters ---------- @@ -279,6 +278,7 @@ def ransac(data, model_class, min_samples, residual_threshold, * `estimate(data)` * `residuals(data)` + * `is_degenerate(data)` min_samples : int The minimum number of data points to fit a model. From 47c3ebac1c178b15998a30a99b184d8413c10c4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:06:39 +0100 Subject: [PATCH 05/54] Add ellipse estimator model --- skimage/measure/fit.py | 199 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 188 insertions(+), 11 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index fb0570e6..b0a8c4e7 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -27,6 +27,8 @@ class LineModel(BaseModel): dist, theta + A minimum number of 2 points is required to solve for the parameters. + ''' def estimate(self, data): @@ -100,13 +102,15 @@ class LineModel(BaseModel): return data.shape[0] < 2 - def predict_x(self, y): + 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 ------- @@ -115,16 +119,20 @@ class LineModel(BaseModel): ''' - dist, theta = self._params + if params is None: + params = self._params + dist, theta = params return (dist - y * np.cos(theta)) / np.cos(theta) - def predict_y(self, x): + 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 ------- @@ -133,7 +141,9 @@ class LineModel(BaseModel): ''' - dist, theta = self._params + if params is None: + params = self._params + dist, theta = params return (dist - x * np.cos(theta)) / np.sin(theta) @@ -154,6 +164,8 @@ class CircleModel(BaseModel): xc, yc, r + A minimum number of 3 points is required to solve for the parameters. + ''' def estimate(self, data): @@ -199,7 +211,7 @@ class CircleModel(BaseModel): def residuals(self, data): '''Determine residuals of data to model. - For each point the shortest distance to the line is returned. + For each point the shortest distance to the circle is returned. Parameters ---------- @@ -236,16 +248,176 @@ class CircleModel(BaseModel): ''' - return data.shape[0] < 2 + return data.shape[0] < 3 - def predict_xy(self, theta): + def predict_xy(self, t, params=None): '''Predict x- and y-coordinates using the estimated model. Parameters ---------- - theta : array + 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 + ------- + x : array + Predicted x-coordinates. + y : array + Predicted 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 x, y + + +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 line model from data using total least squares. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + ''' + + x = data[:, 0] + y = data[:, 1] + + N = data.shape[0] + + A = np.empty((5, 2 * N), dtype=np.double) + + def fun(params): + xt, yt = self.predict_xy(params[5:], params[:5]) + fx = x - xt + fy = y - yt + return np.append(fx, fy) + + # 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. + + ''' + + xc, yc, a, b, theta = self._params + + x = data[:, 0] + y = data[:, 1] + + N = data.shape[0] + + def fun(t, xi, yi): + xt, yt = self.predict_xy(t) + return (xi - xt)**2 + (yi - yt)**2 + + def Dfun(t, xi, yi): + xt, yt = self.predict_xy(t) + dfx = - 2 * (xi - xt) * (- a * np.cos(theta) * np.sin(t) + - b * np.sin(theta) * np.cos(t)) + dfy = - 2 * (yi - yt) * (- a * np.sin(theta) * np.sin(t) + + b * np.cos(theta) * np.cos(t)) + return dfx + dfy + + residuals = np.empty((N, ), dtype=np.double) + + for i in range(N): + xi = x[i] + yi = y[i] + t, _ = optimize.leastsq(fun, 0, args=(xi, yi), Dfun=Dfun, + col_deriv=True) + residuals[i] = np.sqrt(fun(t, xi, yi)) + + return residuals + + @classmethod + def is_degenerate(cls, data): + '''Check whether set of points is degenerate. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + Returns + ------- + flag : bool + Flag indicating if data is degenerate. + + ''' + + return data.shape[0] < 5 + + 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 ------- @@ -256,10 +428,15 @@ class CircleModel(BaseModel): ''' - xc, yc, r = self._params + if params is None: + params = self._params + xc, yc, a, b, theta = params - x = xc + r * np.cos(theta) - y = yc + r * np.sin(theta) + ct = np.cos(t) + st = np.sin(t) + + x = xc + a * np.cos(theta) * ct - b * np.sin(theta) * st + y = yc + a * np.sin(theta) * ct + b * np.cos(theta) * st return x, y From a384ca2befcdbf37744fdb4f9cfdbec58aedf89e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:14:14 +0100 Subject: [PATCH 06/54] Rename variables for better readability --- skimage/measure/fit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index b0a8c4e7..d767681b 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -373,10 +373,10 @@ class EllipseModel(BaseModel): def Dfun(t, xi, yi): xt, yt = self.predict_xy(t) - dfx = - 2 * (xi - xt) * (- a * np.cos(theta) * np.sin(t) - - b * np.sin(theta) * np.cos(t)) - dfy = - 2 * (yi - yt) * (- a * np.sin(theta) * np.sin(t) - + b * np.cos(theta) * np.cos(t)) + dfx_t = - 2 * (xi - xt) * (- a * np.cos(theta) * np.sin(t) + - b * np.sin(theta) * np.cos(t)) + dfy_t = - 2 * (yi - yt) * (- a * np.sin(theta) * np.sin(t) + + b * np.cos(theta) * np.cos(t)) return dfx + dfy residuals = np.empty((N, ), dtype=np.double) From 4dfe927059a816fa4680b0529d74cc1bbe234b59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:15:06 +0100 Subject: [PATCH 07/54] Add todo comment --- skimage/measure/fit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index d767681b..58f7f452 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -339,7 +339,9 @@ class EllipseModel(BaseModel): 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) + # TODO: add function to analytically build jacobian using Dfun for + # faster computation + params, _ = optimize.leastsq(fun, params0) self._params = params[:5] From e337c7e1e5f52359165c6ce1664a1cd07ac84314 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:22:07 +0100 Subject: [PATCH 08/54] Add example to ransac function --- skimage/measure/fit.py | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 58f7f452..af0f5aa9 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -473,6 +473,40 @@ def ransac(data, model_class, min_samples, residual_threshold, inliers : (N,) array Indices of inliers. + 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.cos(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) + >>> print model._params + + Estimate ellipse model using RANSAC: + + >>> ransac_model, _ = ransac(data, EllipseModel, 10, 3, max_trials=50) + >>> print ransac_model._params + ''' best_model = None From d45fb029dc4bf0119062a07b962dbc7fff1f300a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:39:04 +0100 Subject: [PATCH 09/54] Add imports of fit to subpackage --- skimage/measure/__init__.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 226c4235..89f707c1 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -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'] From 65f26e7e9d2722e6556f6838ba1af961ecfad631 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:39:14 +0100 Subject: [PATCH 10/54] Fix example of ransac --- skimage/measure/fit.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index af0f5aa9..b2fb3c13 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -379,7 +379,7 @@ class EllipseModel(BaseModel): - b * np.sin(theta) * np.cos(t)) dfy_t = - 2 * (yi - yt) * (- a * np.sin(theta) * np.sin(t) + b * np.cos(theta) * np.cos(t)) - return dfx + dfy + return dfx_t + dfy_t residuals = np.empty((N, ), dtype=np.double) @@ -484,7 +484,7 @@ def ransac(data, model_class, min_samples, residual_threshold, >>> xc = 20 >>> yc = 30 >>> x = xc + a * np.cos(t) - >>> y = yc + b * 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) @@ -500,12 +500,21 @@ def ransac(data, model_class, min_samples, residual_threshold, >>> model = EllipseModel() >>> model.estimate(data) - >>> print model._params + >>> 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, _ = ransac(data, EllipseModel, 10, 3, max_trials=50) - >>> print ransac_model._params + >>> ransac_model, inliers = ransac(data, EllipseModel, 5, 3, max_trials=50) + >>> # ransac_model._params, inliers + + Should give the correct result estimated without the fauly 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]) ''' From df801a887666a489bf6148490d10598af3f6367b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 19:47:54 +0100 Subject: [PATCH 11/54] Reduce default number of max trials --- skimage/measure/fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index b2fb3c13..aaccc060 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -444,7 +444,7 @@ class EllipseModel(BaseModel): def ransac(data, model_class, min_samples, residual_threshold, - max_trials=1000): + max_trials=100): '''Fits a model to data with the RANSAC (random sample consensus) algorithm. Parameters From 2c329d68ec7967fd26ea183185041a6f0a0b960a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 22:59:07 +0100 Subject: [PATCH 12/54] Add analytical jacobian for ellipse model for 100x speedup --- skimage/measure/fit.py | 43 +++++++++++++++++++++++++++++++++++------- 1 file changed, 36 insertions(+), 7 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index aaccc060..fd24f78e 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -180,9 +180,9 @@ class CircleModel(BaseModel): x = data[:, 0] y = data[:, 1] - # pre-allocate for all iterations - A = np.empty((3, data.shape[0]), dtype=np.double) - # same for all iterations + # 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): @@ -323,7 +323,13 @@ class EllipseModel(BaseModel): N = data.shape[0] - A = np.empty((5, 2 * N), dtype=np.double) + # 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): xt, yt = self.predict_xy(params[5:], params[:5]) @@ -331,6 +337,31 @@ class EllipseModel(BaseModel): fy = y - yt 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) + + # derivatives for fx, fy in the following order: + # xc, yc, a, b, theta, t_i + + # fx + A[2, :N] = - np.cos(theta) * ct + A[3, :N] = np.sin(theta) * st + A[4, :N] = a * np.sin(theta) * ct + b * np.cos(theta) * st + A[5:, :N][diag_idxs] = a * np.cos(theta) * st \ + + b * np.sin(theta) * ct + # fy + A[2, N:] = - np.sin(theta) * ct + A[3, N:] = - np.cos(theta) * st + A[4, N:] = - a * np.cos(theta) * ct + b * np.sin(theta) * st + A[5:, N:][diag_idxs] = a * np.sin(theta) * st \ + - b * np.cos(theta) * ct + + return A + # initial guess of parameters using a circle model params0 = np.empty((N + 5, ), dtype=np.double) xc0 = x.mean() @@ -339,9 +370,7 @@ class EllipseModel(BaseModel): params0[:5] = (xc0, yc0, r0, 0, 0) params0[5:] = np.arctan2(y - yc0, x - xc0) - # TODO: add function to analytically build jacobian using Dfun for - # faster computation - params, _ = optimize.leastsq(fun, params0) + params, _ = optimize.leastsq(fun, params0, Dfun=Dfun, col_deriv=True) self._params = params[:5] From c9353ce8f5c6734727e5da33ff2c1ed796dc67dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 23:11:59 +0100 Subject: [PATCH 13/54] Reduce number of function calls for speedup --- skimage/measure/fit.py | 38 ++++++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 16 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index fd24f78e..232cbde3 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -343,22 +343,22 @@ class EllipseModel(BaseModel): ct = np.cos(t) st = np.sin(t) + ctheta = np.cos(theta) + stheta = np.sin(theta) # derivatives for fx, fy in the following order: # xc, yc, a, b, theta, t_i # fx - A[2, :N] = - np.cos(theta) * ct - A[3, :N] = np.sin(theta) * st - A[4, :N] = a * np.sin(theta) * ct + b * np.cos(theta) * st - A[5:, :N][diag_idxs] = a * np.cos(theta) * st \ - + b * np.sin(theta) * ct + 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:] = - np.sin(theta) * ct - A[3, N:] = - np.cos(theta) * st - A[4, N:] = - a * np.cos(theta) * ct + b * np.sin(theta) * st - A[5:, N:][diag_idxs] = a * np.sin(theta) * st \ - - b * np.cos(theta) * ct + 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 @@ -404,10 +404,14 @@ class EllipseModel(BaseModel): def Dfun(t, xi, yi): xt, yt = self.predict_xy(t) - dfx_t = - 2 * (xi - xt) * (- a * np.cos(theta) * np.sin(t) - - b * np.sin(theta) * np.cos(t)) - dfy_t = - 2 * (yi - yt) * (- a * np.sin(theta) * np.sin(t) - + b * np.cos(theta) * np.cos(t)) + ct = np.cos(t) + st = np.sin(t) + ctheta = np.cos(theta) + stheta = np.sin(theta) + 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) @@ -465,9 +469,11 @@ class EllipseModel(BaseModel): ct = np.cos(t) st = np.sin(t) + ctheta = np.cos(theta) + stheta = np.sin(theta) - x = xc + a * np.cos(theta) * ct - b * np.sin(theta) * st - y = yc + a * np.sin(theta) * ct + b * np.cos(theta) * st + x = xc + a * ctheta * ct - b * stheta * st + y = yc + a * stheta * ct + b * ctheta * st return x, y From 444d51ceb7857e0be8918f5090298a5431db00c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 28 Feb 2013 23:18:25 +0100 Subject: [PATCH 14/54] Calculate initial guess for closest point on ellipse for speedup --- skimage/measure/fit.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 232cbde3..ed89c8b8 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -416,10 +416,14 @@ class EllipseModel(BaseModel): 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] - t, _ = optimize.leastsq(fun, 0, args=(xi, yi), Dfun=Dfun, + t, _ = optimize.leastsq(fun, t0[i], args=(xi, yi), Dfun=Dfun, col_deriv=True) residuals[i] = np.sqrt(fun(t, xi, yi)) From dd714b910c5d6487466463f5cb8f5a2260535be5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 1 Mar 2013 09:50:26 +0100 Subject: [PATCH 15/54] Replace numpy with math for scalar functions and remove Dfun from ellipse residuals for speedup --- skimage/measure/fit.py | 50 +++++++++++++++++++++++------------------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index ed89c8b8..7aee1290 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -1,3 +1,4 @@ +import math import numpy as np from scipy import optimize @@ -56,7 +57,7 @@ class LineModel(BaseModel): # angle perpendicular to line angle theta = (theta + np.pi / 2) % np.pi # line always passes through mean - dist = X0[0] * np.cos(theta) + X0[1] * np.sin(theta) + dist = X0[0] * math.cos(theta) + X0[1] * math.sin(theta) self._params = (dist, theta) @@ -82,7 +83,7 @@ class LineModel(BaseModel): x = data[:, 0] y = data[:, 1] - return dist - (x * np.cos(theta) + y * np.sin(theta)) + return dist - (x * math.cos(theta) + y * math.sin(theta)) @classmethod def is_degenerate(cls, data): @@ -122,7 +123,7 @@ class LineModel(BaseModel): if params is None: params = self._params dist, theta = params - return (dist - y * np.cos(theta)) / np.cos(theta) + return (dist - y * math.cos(theta)) / math.cos(theta) def predict_y(self, x, params=None): '''Predict y-coordinates using the estimated model. @@ -144,7 +145,7 @@ class LineModel(BaseModel): if params is None: params = self._params dist, theta = params - return (dist - x * np.cos(theta)) / np.sin(theta) + return (dist - x * math.cos(theta)) / math.sin(theta) class CircleModel(BaseModel): @@ -343,8 +344,8 @@ class EllipseModel(BaseModel): ct = np.cos(t) st = np.sin(t) - ctheta = np.cos(theta) - stheta = np.sin(theta) + ctheta = math.cos(theta) + stheta = math.sin(theta) # derivatives for fx, fy in the following order: # xc, yc, a, b, theta, t_i @@ -393,26 +394,31 @@ class EllipseModel(BaseModel): 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): - xt, yt = self.predict_xy(t) + 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): - xt, yt = self.predict_xy(t) - ct = np.cos(t) - st = np.sin(t) - ctheta = np.cos(theta) - stheta = np.sin(theta) - 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 + # 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) @@ -423,8 +429,8 @@ class EllipseModel(BaseModel): for i in range(N): xi = x[i] yi = y[i] - t, _ = optimize.leastsq(fun, t0[i], args=(xi, yi), Dfun=Dfun, - col_deriv=True) + # 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 @@ -473,8 +479,8 @@ class EllipseModel(BaseModel): ct = np.cos(t) st = np.sin(t) - ctheta = np.cos(theta) - stheta = np.sin(theta) + ctheta = math.cos(theta) + stheta = math.sin(theta) x = xc + a * ctheta * ct - b * stheta * st y = yc + a * stheta * ct + b * ctheta * st From b2fa37d337138fd48719080d44e865b49e23c279 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 1 Mar 2013 17:20:23 +0100 Subject: [PATCH 16/54] Fix typos --- skimage/measure/fit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 7aee1290..8335c069 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -170,7 +170,7 @@ class CircleModel(BaseModel): ''' def estimate(self, data): - '''Estimate line model from data using total least squares. + '''Estimate circle model from data using total least squares. Parameters ---------- @@ -310,7 +310,7 @@ class EllipseModel(BaseModel): ''' def estimate(self, data): - '''Estimate line model from data using total least squares. + '''Estimate circle model from data using total least squares. Parameters ---------- @@ -554,7 +554,7 @@ def ransac(data, model_class, min_samples, residual_threshold, >>> ransac_model, inliers = ransac(data, EllipseModel, 5, 3, max_trials=50) >>> # ransac_model._params, inliers - Should give the correct result estimated without the fauly data: + 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, From 62b25455dc8b7a1312e0401e8e5991bcacfc26f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 10:30:32 +0200 Subject: [PATCH 17/54] Add input data shape test to avoid out of memory Jacobians --- skimage/measure/fit.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 8335c069..cd9a7795 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -3,6 +3,11 @@ 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): @@ -42,6 +47,8 @@ class LineModel(BaseModel): ''' + _check_data_dim(data, dim=2) + X0 = data.mean(axis=0) if data.shape[0] == 2: # well determined @@ -78,6 +85,8 @@ class LineModel(BaseModel): ''' + _check_data_dim(data, dim=2) + dist, theta = self._params x = data[:, 0] @@ -179,6 +188,8 @@ class CircleModel(BaseModel): ''' + _check_data_dim(data, dim=2) + x = data[:, 0] y = data[:, 1] # pre-allocate jacobian for all iterations @@ -226,6 +237,8 @@ class CircleModel(BaseModel): ''' + _check_data_dim(data, dim=2) + xc, yc, r = self._params x = data[:, 0] @@ -319,6 +332,8 @@ class EllipseModel(BaseModel): ''' + _check_data_dim(data, dim=2) + x = data[:, 0] y = data[:, 1] @@ -392,6 +407,8 @@ class EllipseModel(BaseModel): ''' + _check_data_dim(data, dim=2) + xc, yc, a, b, theta = self._params ctheta = math.cos(theta) From b58f52e40e5e8dbf3557124f31362bbfaeddca00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:15:09 +0200 Subject: [PATCH 18/54] Fix bug in predict_x of line model --- skimage/measure/fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index cd9a7795..3b2d7d32 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -132,7 +132,7 @@ class LineModel(BaseModel): if params is None: params = self._params dist, theta = params - return (dist - y * math.cos(theta)) / math.cos(theta) + return (dist - y * math.sin(theta)) / math.cos(theta) def predict_y(self, x, params=None): '''Predict y-coordinates using the estimated model. From c1bb30f625abb9ecba5a92b4ded37e9e945307d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:15:47 +0200 Subject: [PATCH 19/54] Return combined xy-coordinate array with arbitrary dimensions in predict_xy --- skimage/measure/fit.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 3b2d7d32..4f3445bf 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -277,10 +277,8 @@ class CircleModel(BaseModel): Returns ------- - x : array - Predicted x-coordinates. - y : array - Predicted y-coordinates. + xy : (..., 2) array + Predicted x- and y-coordinates. ''' if params is None: @@ -290,7 +288,7 @@ class CircleModel(BaseModel): x = xc + r * np.cos(t) y = yc + r * np.sin(t) - return x, y + return np.concatenate((x[..., None], y[..., None]), axis=t.ndim) class EllipseModel(BaseModel): @@ -483,10 +481,8 @@ class EllipseModel(BaseModel): Returns ------- - x : array - Predicted x-coordinates. - y : array - Predicted y-coordinates. + xy : (..., 2) array + Predicted x- and y-coordinates. ''' @@ -502,7 +498,7 @@ class EllipseModel(BaseModel): x = xc + a * ctheta * ct - b * stheta * st y = yc + a * stheta * ct + b * ctheta * st - return x, y + return np.concatenate((x[..., None], y[..., None]), axis=t.ndim) def ransac(data, model_class, min_samples, residual_threshold, From 44ff37c8ee5cf4c3d1c95544d8316920cca0dcd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:26:03 +0200 Subject: [PATCH 20/54] Fix bug in ellipse model estimation due to changed predict_xy output --- skimage/measure/fit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 4f3445bf..4c63df92 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -346,9 +346,9 @@ class EllipseModel(BaseModel): diag_idxs = np.diag_indices(N) def fun(params): - xt, yt = self.predict_xy(params[5:], params[:5]) - fx = x - xt - fy = y - yt + xyt = self.predict_xy(params[5:], params[:5]) + fx = x - xyt[:, 0] + fy = y - xyt[:, 1] return np.append(fx, fy) def Dfun(params): From 53b7e6e41867298b1546093aaf6bdf1e3163d9ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:27:00 +0200 Subject: [PATCH 21/54] Add test cases for line model --- skimage/measure/tests/test_fit.py | 43 +++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 skimage/measure/tests/test_fit.py diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py new file mode 100644 index 00000000..137a6e44 --- /dev/null +++ b/skimage/measure/tests/test_fit.py @@ -0,0 +1,43 @@ +import numpy as np +from numpy.testing import assert_equal, assert_raises, assert_almost_equal +from skimage.measure import LineModel, CircleModel, EllipseModel + + +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_is_degenerate(): + assert_equal(LineModel().is_degenerate(np.empty((1, 2))), True) + + +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 equals original parameters + assert_almost_equal(model0._params, model_est._params, 1) + + +if __name__ == "__main__": + np.testing.run_module_suite() From e880f155a91d61a2416945d9875f93d1b97d101a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:27:25 +0200 Subject: [PATCH 22/54] Add test cases for circle model --- skimage/measure/tests/test_fit.py | 37 +++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 137a6e44..c3750789 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -39,5 +39,42 @@ def test_line_model_estimate(): assert_almost_equal(model0._params, model_est._params, 1) +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_is_degenerate(): + assert_equal(CircleModel().is_degenerate(np.empty((1, 2))), True) + + +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 equals original parameters + assert_almost_equal(model0._params, model_est._params, 1) + + if __name__ == "__main__": np.testing.run_module_suite() From 80827d3959023ba5ed9d8e617c7122ef9fc0734d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 11:34:54 +0200 Subject: [PATCH 23/54] Add test case for ransac algorithm --- skimage/measure/tests/test_fit.py | 63 +++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index c3750789..83b73840 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 +from skimage.measure import LineModel, CircleModel, EllipseModel, ransac def test_line_model_invalid_input(): @@ -35,7 +35,7 @@ def test_line_model_estimate(): model_est = LineModel() model_est.estimate(data) - # test whether estimated parameters almost equals original parameters + # test whether estimated parameters almost equal original parameters assert_almost_equal(model0._params, model_est._params, 1) @@ -72,9 +72,66 @@ def test_circle_model_estimate(): model_est = CircleModel() model_est.estimate(data) - # test whether estimated parameters almost equals original parameters + # test whether estimated parameters almost equal original parameters assert_almost_equal(model0._params, model_est._params, 1) +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_is_degenerate(): + assert_equal(EllipseModel().is_degenerate(np.empty((1, 2))), True) + + +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_ransac(): + # 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) + 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 + + if __name__ == "__main__": np.testing.run_module_suite() From aa309fd926f1ff50906b43b1b3371d0ed34cda81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 12:15:05 +0200 Subject: [PATCH 24/54] Replace triple-single-quotes with triple-double-quotes for doc strings --- skimage/measure/fit.py | 68 +++++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 4c63df92..7aa94d4c 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -16,7 +16,7 @@ class BaseModel(object): class LineModel(BaseModel): - '''Total least squares estimator for 2D lines. + """Total least squares estimator for 2D lines. Lines are parameterized using polar coordinates as functional model: @@ -35,17 +35,17 @@ class LineModel(BaseModel): 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. + """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) @@ -69,7 +69,7 @@ class LineModel(BaseModel): self._params = (dist, theta) def residuals(self, data): - '''Determine residuals of data to model. + """Determine residuals of data to model. For each point the shortest distance to the line is returned. @@ -83,7 +83,7 @@ class LineModel(BaseModel): residuals : (N, ) array Residual for each data point. - ''' + """ _check_data_dim(data, dim=2) @@ -96,7 +96,7 @@ class LineModel(BaseModel): @classmethod def is_degenerate(cls, data): - '''Check whether set of points is degenerate. + """Check whether set of points is degenerate. Parameters ---------- @@ -108,12 +108,12 @@ class LineModel(BaseModel): flag : bool Flag indicating if data is degenerate. - ''' + """ return data.shape[0] < 2 def predict_x(self, y, params=None): - '''Predict x-coordinates using the estimated model. + """Predict x-coordinates using the estimated model. Parameters ---------- @@ -127,7 +127,7 @@ class LineModel(BaseModel): x : array Predicted x-coordinates. - ''' + """ if params is None: params = self._params @@ -135,7 +135,7 @@ class LineModel(BaseModel): return (dist - y * math.sin(theta)) / math.cos(theta) def predict_y(self, x, params=None): - '''Predict y-coordinates using the estimated model. + """Predict y-coordinates using the estimated model. Parameters ---------- @@ -149,7 +149,7 @@ class LineModel(BaseModel): y : array Predicted y-coordinates. - ''' + """ if params is None: params = self._params @@ -159,7 +159,7 @@ class LineModel(BaseModel): class CircleModel(BaseModel): - '''Total least squares estimator for 2D circles. + """Total least squares estimator for 2D circles. The functional model of the circle is: @@ -176,17 +176,17 @@ class CircleModel(BaseModel): 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. + """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) @@ -221,7 +221,7 @@ class CircleModel(BaseModel): self._params = params def residuals(self, data): - '''Determine residuals of data to model. + """Determine residuals of data to model. For each point the shortest distance to the circle is returned. @@ -235,7 +235,7 @@ class CircleModel(BaseModel): residuals : (N, ) array Residual for each data point. - ''' + """ _check_data_dim(data, dim=2) @@ -248,7 +248,7 @@ class CircleModel(BaseModel): @classmethod def is_degenerate(cls, data): - '''Check whether set of points is degenerate. + """Check whether set of points is degenerate. Parameters ---------- @@ -260,12 +260,12 @@ class CircleModel(BaseModel): flag : bool Flag indicating if data is degenerate. - ''' + """ return data.shape[0] < 3 def predict_xy(self, t, params=None): - '''Predict x- and y-coordinates using the estimated model. + """Predict x- and y-coordinates using the estimated model. Parameters ---------- @@ -280,7 +280,7 @@ class CircleModel(BaseModel): xy : (..., 2) array Predicted x- and y-coordinates. - ''' + """ if params is None: params = self._params xc, yc, r = params @@ -293,7 +293,7 @@ class CircleModel(BaseModel): class EllipseModel(BaseModel): - '''Total least squares estimator for 2D ellipses. + """Total least squares estimator for 2D ellipses. The functional model of the ellipse is: @@ -318,17 +318,17 @@ class EllipseModel(BaseModel): 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. + """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) @@ -389,7 +389,7 @@ class EllipseModel(BaseModel): self._params = params[:5] def residuals(self, data): - '''Determine residuals of data to model. + """Determine residuals of data to model. For each point the shortest distance to the ellipse is returned. @@ -403,7 +403,7 @@ class EllipseModel(BaseModel): residuals : (N, ) array Residual for each data point. - ''' + """ _check_data_dim(data, dim=2) @@ -452,7 +452,7 @@ class EllipseModel(BaseModel): @classmethod def is_degenerate(cls, data): - '''Check whether set of points is degenerate. + """Check whether set of points is degenerate. Parameters ---------- @@ -464,12 +464,12 @@ class EllipseModel(BaseModel): flag : bool Flag indicating if data is degenerate. - ''' + """ return data.shape[0] < 5 def predict_xy(self, t, params=None): - '''Predict x- and y-coordinates using the estimated model. + """Predict x- and y-coordinates using the estimated model. Parameters ---------- @@ -484,7 +484,7 @@ class EllipseModel(BaseModel): xy : (..., 2) array Predicted x- and y-coordinates. - ''' + """ if params is None: params = self._params @@ -503,7 +503,7 @@ class EllipseModel(BaseModel): def ransac(data, model_class, min_samples, residual_threshold, max_trials=100): - '''Fits a model to data with the RANSAC (random sample consensus) algorithm. + """Fits a model to data with the RANSAC (random sample consensus) algorithm. Parameters ---------- @@ -574,7 +574,7 @@ def ransac(data, model_class, min_samples, residual_threshold, 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 From 5b78b2a69e1616893a5c48472cc471638ba72d0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 12:26:04 +0200 Subject: [PATCH 25/54] Add is_degenerate to geometric transforms --- skimage/transform/_geometric.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 448829de..4b2a38d5 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -41,6 +41,26 @@ class GeometricTransform(object): """ raise NotImplementedError() + @classmethod + def is_degenerate(cls, data): + """Check whether set of points is degenerate. + + Parameters + ---------- + data : (N, 2) array + N points with `(x, y)` coordinates, respectively. + + Returns + ------- + flag : bool + Flag indicating if data is degenerate. + + """ + + # by default never degenerate since system equations can be under-, + # well- and over-determined. + return False + def __add__(self, other): """Combine this transformation with another. From bf2a5a429b5853a08ab3a18b261c0deb848e59f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 12:57:30 +0200 Subject: [PATCH 26/54] Add residuals method to geometric transforms --- skimage/transform/_geometric.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 4b2a38d5..1fcb55bb 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -61,6 +61,28 @@ class GeometricTransform(object): # well- and over-determined. return False + 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. From 5da3418fff112eee05b3055c7b1f6724ccfddc60 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 13:28:10 +0200 Subject: [PATCH 27/54] Change ransac function to accept multiple input data arrays --- skimage/measure/fit.py | 44 +++++++++++++++++++++++++++++++++--------- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 7aa94d4c..ecd0f6d7 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -507,9 +507,12 @@ def ransac(data, model_class, min_samples, residual_threshold, Parameters ---------- - data : (N, D) array + 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. model_class : object Object with the following methods implemented: @@ -578,35 +581,58 @@ def ransac(data, model_class, min_samples, residual_threshold, best_model = None best_inlier_num = 0 + best_inlier_residuals_sum = np.inf best_inliers = None - data_idxs = np.arange(data.shape[0]) + + 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] + + data_idxs = np.arange(N) for _ in range(max_trials): # choose random sample - sample = data[np.random.randint(0, data.shape[0], min_samples)] + samples = [] + random_idxs = np.random.randint(0, N, min_samples) + for d in data: + samples.append(d[random_idxs]) # check if random sample is degenerate - if model_class.is_degenerate(sample): + if model_class.is_degenerate(*samples): continue # create new instance of model class for current sample sample_model = model_class() - sample_model.estimate(sample) - sample_model_residuals = sample_model.residuals(data) + sample_model.estimate(*samples) + sample_model_residuals = np.abs(sample_model.residuals(*data)) # consensus set / inliers - sample_model_inliers = data_idxs[np.abs(sample_model_residuals) + sample_model_inliers = data_idxs[sample_model_residuals < residual_threshold] + sample_model_residuals_sum = np.sum(sample_model_residuals) # choose as new best model if number of inliers is maximal sample_inlier_num = sample_model_inliers.shape[0] - if sample_inlier_num > best_inlier_num: + 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 # estimate final model using all inliers if best_inliers is not None: - best_model.estimate(data[best_inliers]) + for i in range(len(data)): + data[i] = data[i][best_inliers] + best_model.estimate(*data) return best_model, best_inliers From 087fe8891e425aec0d384b00dc856243ded338f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 13:29:10 +0200 Subject: [PATCH 28/54] Fix is_degenerate function to accept src and dst arrays separately --- skimage/transform/_geometric.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 1fcb55bb..12983e6d 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -42,13 +42,15 @@ class GeometricTransform(object): raise NotImplementedError() @classmethod - def is_degenerate(cls, data): + def is_degenerate(cls, src, dst): """Check whether set of points is degenerate. Parameters ---------- - data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + src : (N, 2) array + Source coordinates. + dst : (N, 2) array + Destination coordinates. Returns ------- From fb74a9d8ae75dc75e979c6e4b36f17c3cbafd489 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 13:59:01 +0200 Subject: [PATCH 29/54] Add new stop criteria to ransac and improve doc string --- skimage/measure/fit.py | 44 ++++++++++++++++++++++++++++----- skimage/transform/_geometric.py | 2 +- 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index ecd0f6d7..ac55534e 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -502,9 +502,28 @@ class EllipseModel(BaseModel): def ransac(data, model_class, min_samples, residual_threshold, - max_trials=100): + max_trials=100, stop_sample_num=np.inf, stop_residuals_sum=0): """Fits 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. + 2. Estimate a model to the random subset (`estimate(*data[random_subset]`). + 3. Classify all data as inliers or outliers by calculating the residuals to + the estimated model (`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 @@ -512,13 +531,16 @@ def ransac(data, model_class, min_samples, residual_threshold, 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. + can be optionally passed as tuple or list. Note, that in this case the + functions `estimate(*data)`, `residuals(*data)` and + `is_degenerate(*data)` must all take each data array as separate + arguments. model_class : object Object with the following methods implemented: - * `estimate(data)` - * `residuals(data)` - * `is_degenerate(data)` + * `estimate(*data)` + * `residuals(*data)` + * `is_degenerate(*data)` min_samples : int The minimum number of data points to fit a model. @@ -526,6 +548,10 @@ def ransac(data, model_class, min_samples, residual_threshold, Maximum distance for a data point to be classified as an inlier. 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 ------- @@ -613,7 +639,7 @@ def ransac(data, model_class, min_samples, residual_threshold, # consensus set / inliers sample_model_inliers = data_idxs[sample_model_residuals < residual_threshold] - sample_model_residuals_sum = np.sum(sample_model_residuals) + sample_model_residuals_sum = np.sum(sample_model_residuals**2) # choose as new best model if number of inliers is maximal sample_inlier_num = sample_model_inliers.shape[0] @@ -628,9 +654,15 @@ def ransac(data, model_class, min_samples, residual_threshold, 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) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 12983e6d..54b0cb61 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -83,7 +83,7 @@ class GeometricTransform(object): """ - return np.sqrt(np.sum((self(src) - dst) ** 2, axis=1)) + return np.sqrt(np.sum((self(src) - dst)**2, axis=1)) def __add__(self, other): """Combine this transformation with another. From 22224727d3be1d0d53756bfdc480eef4ae9a158b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 14:14:09 +0200 Subject: [PATCH 30/54] Add example for ransac and geometric transforms --- skimage/measure/fit.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index ac55534e..bf491688 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -603,6 +603,23 @@ def ransac(data, model_class, min_samples, residual_threshold, 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 From 82566157f65be3068ede857713b0961971e78665 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 14:18:00 +0200 Subject: [PATCH 31/54] Add test case for ransac applied to geometric transform --- skimage/measure/tests/test_fit.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 83b73840..e565df7b 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -1,6 +1,7 @@ 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(): @@ -113,12 +114,14 @@ def test_ellipse_model_estimate(): assert_almost_equal(model0._params, model_est._params, 0) -def test_ransac(): +def test_ransac_shape(): # 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) @@ -133,5 +136,27 @@ def test_ransac(): assert outlier not in inliers +def test_ransac_geometric(): + # 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[0] = (10000, 10000) + dst[1] = (-100, 100) + dst[2] = (50, 50) + + # estimate parameters of corrupted data + model_est, inliers = ransac((src, dst), AffineTransform, 2, 10) + + # test whether estimated parameters equal original parameters + assert_almost_equal(model0._matrix, model_est._matrix) + for outlier in outliers: + assert outlier not in inliers + + if __name__ == "__main__": np.testing.run_module_suite() From db496d0af9dbcdff22d74372beba0919161e9ce4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 14:22:07 +0200 Subject: [PATCH 32/54] Remove trailing brackets --- skimage/measure/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index bf491688..b0d58e28 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -598,10 +598,10 @@ def ransac(data, model_class, min_samples, residual_threshold, Should give the correct result estimated without the faulty data: - [ 20.12762373, 29.73563061, 4.81499637, 10.4743584, 0.05217117]) + [ 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]) + 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49] Robustly estimate geometric transformation: From 28bda25da88e02fd43727e07d99620e8b45f37a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 1 May 2013 18:39:16 +0200 Subject: [PATCH 33/54] Add is_degenerate to description of ransac iteration --- skimage/measure/fit.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index b0d58e28..ef81f8f9 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -509,7 +509,8 @@ def ransac(data, model_class, min_samples, residual_threshold, 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. + 1. Select `min_samples` random samples from the original data and check + whether the set of points is valid (`is_degenerate(*data)`). 2. Estimate a model to the random subset (`estimate(*data[random_subset]`). 3. Classify all data as inliers or outliers by calculating the residuals to the estimated model (`residuals(*data)`) - all data samples with From a73076157c1e13372026f514b538681f7e5eadb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 2 May 2013 18:12:23 +0200 Subject: [PATCH 34/54] Replace is_degenerate with is_model_valid and is_data_valid in ransac --- skimage/measure/fit.py | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index ef81f8f9..8fde1ca2 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -502,6 +502,7 @@ class EllipseModel(BaseModel): 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): """Fits a model to data with the RANSAC (random sample consensus) algorithm. @@ -510,11 +511,13 @@ def ransac(data, model_class, min_samples, residual_threshold, the following tasks: 1. Select `min_samples` random samples from the original data and check - whether the set of points is valid (`is_degenerate(*data)`). - 2. Estimate a model to the random subset (`estimate(*data[random_subset]`). + 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 (`residuals(*data)`) - all data samples with - residuals smaller than the `residual_threshold` are considered as + 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 @@ -537,16 +540,21 @@ def ransac(data, model_class, min_samples, residual_threshold, `is_degenerate(*data)` must all take each data array as separate arguments. model_class : object - Object with the following methods implemented: + Object with the following object methods: * `estimate(*data)` * `residuals(*data)` - * `is_degenerate(*data)` min_samples : int The minimum number of data points to fit a model. 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 @@ -640,19 +648,25 @@ def ransac(data, model_class, min_samples, residual_threshold, for _ in range(max_trials): - # choose random sample + # 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 is degenerate - if model_class.is_degenerate(*samples): + # check if random sample set is valid + if is_data_valid is not None and not is_data_valid(*samples): continue - # create new instance of model class for current sample + # 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 = data_idxs[sample_model_residuals From 8d92f5c02b3a4c6a2cde311b79c21b8f68584f80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 2 May 2013 18:13:33 +0200 Subject: [PATCH 35/54] Remove is_degenerate methods from estimation models --- skimage/measure/fit.py | 54 --------------------------------- skimage/transform/_geometric.py | 22 -------------- 2 files changed, 76 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 8fde1ca2..59750940 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -94,24 +94,6 @@ class LineModel(BaseModel): return dist - (x * math.cos(theta) + y * math.sin(theta)) - @classmethod - def is_degenerate(cls, data): - """Check whether set of points is degenerate. - - Parameters - ---------- - data : (N, 2) array - N points with `(x, y)` coordinates, respectively. - - Returns - ------- - flag : bool - Flag indicating if data is degenerate. - - """ - - return data.shape[0] < 2 - def predict_x(self, y, params=None): """Predict x-coordinates using the estimated model. @@ -246,24 +228,6 @@ class CircleModel(BaseModel): return r - np.sqrt((x - xc)**2 + (y - yc)**2) - @classmethod - def is_degenerate(cls, data): - """Check whether set of points is degenerate. - - Parameters - ---------- - data : (N, 2) array - N points with `(x, y)` coordinates, respectively. - - Returns - ------- - flag : bool - Flag indicating if data is degenerate. - - """ - - return data.shape[0] < 3 - def predict_xy(self, t, params=None): """Predict x- and y-coordinates using the estimated model. @@ -450,24 +414,6 @@ class EllipseModel(BaseModel): return residuals - @classmethod - def is_degenerate(cls, data): - """Check whether set of points is degenerate. - - Parameters - ---------- - data : (N, 2) array - N points with `(x, y)` coordinates, respectively. - - Returns - ------- - flag : bool - Flag indicating if data is degenerate. - - """ - - return data.shape[0] < 5 - def predict_xy(self, t, params=None): """Predict x- and y-coordinates using the estimated model. diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 54b0cb61..1ff8e9d5 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -41,28 +41,6 @@ class GeometricTransform(object): """ raise NotImplementedError() - @classmethod - def is_degenerate(cls, src, dst): - """Check whether set of points is degenerate. - - Parameters - ---------- - src : (N, 2) array - Source coordinates. - dst : (N, 2) array - Destination coordinates. - - Returns - ------- - flag : bool - Flag indicating if data is degenerate. - - """ - - # by default never degenerate since system equations can be under-, - # well- and over-determined. - return False - def residuals(self, src, dst): """Determine residuals of transformed destination coordinates. From c8fa3b48e557ddf3d171b19ded81fc707958e70d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 2 May 2013 18:19:39 +0200 Subject: [PATCH 36/54] Replace is_degenerate tests with is_data_valid and is_model_valid tests --- skimage/measure/tests/test_fit.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index e565df7b..7faa3c0b 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -16,10 +16,6 @@ def test_line_model_predict(): assert_almost_equal(x, model.predict_x(y)) -def test_line_model_is_degenerate(): - assert_equal(LineModel().is_degenerate(np.empty((1, 2))), True) - - def test_line_model_estimate(): # generate original data without noise model0 = LineModel() @@ -54,10 +50,6 @@ def test_circle_model_predict(): assert_almost_equal(xy, model.predict_xy(t)) -def test_circle_model_is_degenerate(): - assert_equal(CircleModel().is_degenerate(np.empty((1, 2))), True) - - def test_circle_model_estimate(): # generate original data without noise model0 = CircleModel() @@ -91,10 +83,6 @@ def test_ellipse_model_predict(): assert_almost_equal(xy, model.predict_xy(t)) -def test_ellipse_model_is_degenerate(): - assert_equal(EllipseModel().is_degenerate(np.empty((1, 2))), True) - - def test_ellipse_model_estimate(): # generate original data without noise model0 = EllipseModel() @@ -158,5 +146,22 @@ def test_ransac_geometric(): assert outlier not in inliers +def test_ransac_is_data_valid(): + 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(): + 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() From 660a6892559ab81d0d7384606dac802aaeaf8838 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 4 May 2013 11:41:25 +0200 Subject: [PATCH 37/54] Add reference for RANSAC --- skimage/measure/fit.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 59750940..3f0636c8 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -515,6 +515,10 @@ def ransac(data, model_class, min_samples, residual_threshold, inliers : (N,) array Indices of inliers. + References + ---------- + .. [1] http://en.wikipedia.org/wiki/RANSAC + Examples -------- From d857c452d31855f04d4274d3189fccf0f3f2da2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 4 May 2013 11:42:43 +0200 Subject: [PATCH 38/54] Add example matching script --- doc/examples/plot_matching.py | 143 ++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 doc/examples/plot_matching.py diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py new file mode 100644 index 00000000..abed0ee9 --- /dev/null +++ b/doc/examples/plot_matching.py @@ -0,0 +1,143 @@ +""" +============================ +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 measure. 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, shear=0, + 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) + + +# 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 i in range(len(src)): + if i in inliers: + ax_idx = 0 + color = 'g' + else: + ax_idx = 1 + color = 'r' + ax[ax_idx].plot((src[i, 1], dst[i, 1] + 200), (src[i, 0], dst[i, 0]), '-', + color=color) + ax[ax_idx].plot(src[i, 1], src[i, 0], '.', markersize=10, color=color) + ax[ax_idx].plot(dst[i, 1] + 200, dst[i, 0], '.', markersize=10, color=color) + +plt.show() From 15e079a8c2743eeaa303a347065a3e1e8f9cd44e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 4 May 2013 11:46:07 +0200 Subject: [PATCH 39/54] Remove unused shear argument --- doc/examples/plot_matching.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index abed0ee9..0546d6f5 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -41,8 +41,7 @@ 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, shear=0, - translation=(20, -10)) +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) From 43b4e9efb193ee666303169e212207187cf540ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 07:02:37 +0200 Subject: [PATCH 40/54] Remove extra space --- doc/examples/plot_matching.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index 0546d6f5..d5839fe7 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -4,7 +4,7 @@ Robust matching using RANSAC ============================ In this simplified example we first generate two synthetic images as if they -were taken from different view points. +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 measure. Note, From ef7b82e116e53a182b90956342f8c36feee45837 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 07:04:17 +0200 Subject: [PATCH 41/54] Clarify SSD description --- doc/examples/plot_matching.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index d5839fe7..57919080 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -7,9 +7,10 @@ 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 measure. Note, -that this measure is only robust towards linear radiometric and not geometric -distortions and is thus only usable with slight view point changes. +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 From d6c49ea6b81425b714923b89d17bc474c8571828 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 07:08:50 +0200 Subject: [PATCH 42/54] Add comment to explain commented out code --- skimage/measure/fit.py | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 3f0636c8..621eca98 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -191,6 +191,7 @@ class CircleModel(BaseModel): 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 From 5656feb61be869e5e78da56dce65ef81887b20fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 07:10:05 +0200 Subject: [PATCH 43/54] Use double quotes for code in doc string --- skimage/measure/fit.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 621eca98..9b9d4526 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -43,7 +43,7 @@ class LineModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. """ @@ -76,7 +76,7 @@ class LineModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. Returns ------- @@ -166,7 +166,7 @@ class CircleModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. """ @@ -211,7 +211,7 @@ class CircleModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. Returns ------- @@ -291,7 +291,7 @@ class EllipseModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. """ @@ -361,7 +361,7 @@ class EllipseModel(BaseModel): Parameters ---------- data : (N, 2) array - N points with `(x, y)` coordinates, respectively. + N points with ``(x, y)`` coordinates, respectively. Returns ------- From 9e7fb07d67b58382cea197f1597fff1ba921e933 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 11:36:38 +0200 Subject: [PATCH 44/54] Add missing double colons for equations --- skimage/measure/fit.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 9b9d4526..cd739316 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -18,18 +18,18 @@ class LineModel(BaseModel): """Total least squares estimator for 2D lines. - Lines are parameterized using polar coordinates as functional model: + 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: + 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: + The `_params` attribute contains the parameters in the following order:: dist, theta @@ -143,16 +143,16 @@ class CircleModel(BaseModel): """Total least squares estimator for 2D circles. - The functional model of the circle is: + 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: + circle:: min{ sum((r - sqrt((x_i - xc)**2 + (y_i - yc)**2))**2) } - The `_params` attribute contains the parameters in the following order: + The `_params` attribute contains the parameters in the following order:: xc, yc, r @@ -260,7 +260,7 @@ class EllipseModel(BaseModel): """Total least squares estimator for 2D ellipses. - The functional model of the ellipse is: + 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) @@ -270,14 +270,14 @@ class EllipseModel(BaseModel): shortest distance from the point to the ellipse. This estimator minimizes the squared distances from all points to the - ellipse: + 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: + The `_params` attribute contains the parameters in the following order:: xc, yc, a, b, theta @@ -513,7 +513,7 @@ def ransac(data, model_class, min_samples, residual_threshold, ------- model : object Best model with largest consensus set. - inliers : (N,) array + inliers : (N, ) array Indices of inliers. References From 252c889427ccbf02d099b9b13e10780262b881a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 11:38:46 +0200 Subject: [PATCH 45/54] Add missing double quotes for equations and code --- skimage/measure/fit.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index cd739316..ff2c6420 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -23,13 +23,13 @@ class LineModel(BaseModel): 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`. + 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:: + The ``_params`` attribute contains the parameters in the following order:: dist, theta @@ -152,7 +152,7 @@ class CircleModel(BaseModel): min{ sum((r - sqrt((x_i - xc)**2 + (y_i - yc)**2))**2) } - The `_params` attribute contains the parameters in the following order:: + The ``_params`` attribute contains the parameters in the following order:: xc, yc, r @@ -266,18 +266,18 @@ class EllipseModel(BaseModel): 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. + 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`. + 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:: + The ``_params`` attribute contains the parameters in the following order:: xc, yc, a, b, theta From 63577f275c6d06d0b5395f6f7601c4e9f449abd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 11:41:08 +0200 Subject: [PATCH 46/54] Update doc string of ransac with new validation functions --- skimage/measure/fit.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index ff2c6420..966009da 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -483,14 +483,15 @@ def ransac(data, model_class, min_samples, residual_threshold, 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)` and - `is_degenerate(*data)` must all take each data array as separate - arguments. + 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)` + * ``estimate(*data)`` + * ``residuals(*data)`` min_samples : int The minimum number of data points to fit a model. From f18eef6bc0d2faa0167420394548bfe63bbe604e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 11:43:52 +0200 Subject: [PATCH 47/54] Add title of wikipedia reference --- skimage/measure/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 966009da..5f2d328f 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -494,7 +494,7 @@ def ransac(data, model_class, min_samples, residual_threshold, * ``residuals(*data)`` min_samples : int - The minimum number of data points to fit a model. + 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 @@ -519,7 +519,7 @@ def ransac(data, model_class, min_samples, residual_threshold, References ---------- - .. [1] http://en.wikipedia.org/wiki/RANSAC + .. [1] "RANSAC", Wikipedia, http://en.wikipedia.org/wiki/RANSAC Examples -------- From 434eb48620a0da75d773cbb2cb376023dfb26202 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 11:51:20 +0200 Subject: [PATCH 48/54] Add test case for underdetermined LineModel estimation --- skimage/measure/tests/test_fit.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 7faa3c0b..ad45d02f 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -36,6 +36,11 @@ def test_line_model_estimate(): assert_almost_equal(model0._params, model_est._params, 1) +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))) From 231596b35148242282a2b2eab6c5eb51c9076985 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 12:08:12 +0200 Subject: [PATCH 49/54] Add test cases for residuals --- skimage/measure/tests/test_fit.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index ad45d02f..f0401028 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -36,6 +36,17 @@ def test_line_model_estimate(): 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) @@ -74,6 +85,15 @@ def test_circle_model_estimate(): 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))) @@ -107,6 +127,15 @@ def test_ellipse_model_estimate(): 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(): # generate original data without noise model0 = CircleModel() From 785e602aba9a9556bc71f99ceb5f06f52f42caf2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 15:32:33 +0200 Subject: [PATCH 50/54] Fix last remaining PEP8 errors --- skimage/measure/fit.py | 32 ++++++++++++++++--------------- skimage/measure/tests/test_fit.py | 2 +- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index 5f2d328f..fa78fb82 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -25,7 +25,8 @@ class LineModel(BaseModel): 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:: + This estimator minimizes the squared distances from all points to the + line:: min{ sum((dist - x_i * cos(theta) + y_i * sin(theta))**2) } @@ -51,14 +52,15 @@ class LineModel(BaseModel): 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 + 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 + else: # under-determined raise ValueError('At least 2 input points needed.') # angle perpendicular to line angle @@ -266,8 +268,8 @@ class EllipseModel(BaseModel): 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. + 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:: @@ -451,19 +453,19 @@ class EllipseModel(BaseModel): 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): - """Fits a model to data with the RANSAC (random sample consensus) algorithm. + """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: + 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 + 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 @@ -481,9 +483,9 @@ def ransac(data, model_class, min_samples, residual_threshold, 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)``, + 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. diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index f0401028..0f21b49b 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -162,7 +162,7 @@ def test_ransac_geometric(): # generate original data without noise src = 100 * np.random.random((50, 2)) model0 = AffineTransform(scale=(0.5, 0.3), rotation=1, - translation=(10, 20)) + translation=(10, 20)) dst = model0(src) # add some faulty data From 01124f5bcc75e07c3f7e13b0aefc8942ada2f76d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 18:08:21 +0200 Subject: [PATCH 51/54] Use boolean mask for inlier return value of RANSAC --- doc/examples/plot_matching.py | 16 ++++++---------- skimage/measure/fit.py | 9 +++------ skimage/measure/tests/test_fit.py | 11 +++++------ 3 files changed, 14 insertions(+), 22 deletions(-) diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index 57919080..d01231fe 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -105,6 +105,7 @@ 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 @@ -128,16 +129,11 @@ ax[1].axis('off') ax[1].axis((0, 400, 200, 0)) ax[1].set_title('Faulty correspondencies') -for i in range(len(src)): - if i in inliers: - ax_idx = 0 - color = 'g' - else: - ax_idx = 1 - color = 'r' - ax[ax_idx].plot((src[i, 1], dst[i, 1] + 200), (src[i, 0], dst[i, 0]), '-', + +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[i, 1], src[i, 0], '.', markersize=10, color=color) - ax[ax_idx].plot(dst[i, 1] + 200, dst[i, 0], '.', markersize=10, 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() diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py index fa78fb82..87febf70 100644 --- a/skimage/measure/fit.py +++ b/skimage/measure/fit.py @@ -517,7 +517,7 @@ def ransac(data, model_class, min_samples, residual_threshold, model : object Best model with largest consensus set. inliers : (N, ) array - Indices of inliers. + Boolean mask of inliers classified as ``True``. References ---------- @@ -598,8 +598,6 @@ def ransac(data, model_class, min_samples, residual_threshold, # number of samples N = data[0].shape[0] - data_idxs = np.arange(N) - for _ in range(max_trials): # choose random sample set @@ -623,12 +621,11 @@ def ransac(data, model_class, min_samples, residual_threshold, sample_model_residuals = np.abs(sample_model.residuals(*data)) # consensus set / inliers - sample_model_inliers = data_idxs[sample_model_residuals - < residual_threshold] + 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 = sample_model_inliers.shape[0] + sample_inlier_num = np.sum(sample_model_inliers) if ( # more inliers sample_inlier_num > best_inlier_num diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index 0f21b49b..a2a40f58 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -167,17 +167,16 @@ def test_ransac_geometric(): # add some faulty data outliers = (0, 5, 20) - dst[0] = (10000, 10000) - dst[1] = (-100, 100) - dst[2] = (50, 50) + 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, 10) + model_est, inliers = ransac((src, dst), AffineTransform, 2, 20) # test whether estimated parameters equal original parameters assert_almost_equal(model0._matrix, model_est._matrix) - for outlier in outliers: - assert outlier not in inliers + assert np.all(np.nonzero(inliers == False)[0] == outliers) def test_ransac_is_data_valid(): From 54e3757cb0d2931a3a7969be703b180811eb0f35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 18:08:46 +0200 Subject: [PATCH 52/54] Wrap lines --- doc/examples/plot_matching.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py index d01231fe..6d3ea72f 100644 --- a/doc/examples/plot_matching.py +++ b/doc/examples/plot_matching.py @@ -8,9 +8,9 @@ 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. +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 @@ -134,6 +134,7 @@ 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) + ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10, + color=color) plt.show() From 13bfd603f391c2e573200e9765354e91eb8f52de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 18:25:15 +0200 Subject: [PATCH 53/54] Add simple RANSAC example script --- doc/examples/plot_ransac.py | 55 +++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 doc/examples/plot_ransac.py diff --git a/doc/examples/plot_ransac.py b/doc/examples/plot_ransac.py new file mode 100644 index 00000000..26fd78c0 --- /dev/null +++ b/doc/examples/plot_ransac.py @@ -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() From 219d6217250ff176937dd220013f96d928bcb52a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 6 May 2013 18:34:34 +0200 Subject: [PATCH 54/54] Set random seed to avoid random test failures --- skimage/measure/tests/test_fit.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py index a2a40f58..bde1fcc7 100644 --- a/skimage/measure/tests/test_fit.py +++ b/skimage/measure/tests/test_fit.py @@ -137,6 +137,8 @@ def test_line_model_residuals(): def test_ransac_shape(): + np.random.seed(1) + # generate original data without noise model0 = CircleModel() model0._params = (10, 12, 3) @@ -159,6 +161,8 @@ def test_ransac_shape(): 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, @@ -180,6 +184,8 @@ def test_ransac_geometric(): 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) @@ -188,6 +194,8 @@ def test_ransac_is_data_valid(): 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,