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