Code covered by the BSD License

# fitcircle.m

### Richard Brown (view profile)

20 May 2007 (Updated )

Fits circles to 2D data using nonlinear least squares to minimise geometric error

fitcircle(x, varargin)
```function [z, r, residual] = fitcircle(x, varargin)
%FITCIRCLE    least squares circle fit
%
%   [Z, R] = FITCIRCLE(X) fits a circle to the N points in X minimising
%   geometric error (sum of squared distances from the points to the fitted
%   circle) using nonlinear least squares (Gauss Newton)
%       Input
%           X : 2xN array of N 2D points, with N >= 3
%       Output
%           Z : center of the fitted circle
%           R : radius of the fitted circle
%
%   [Z, R] = FITCIRCLE(X, 'linear') fits a circle using linear least
%   squares minimising the algebraic error (residual from fitting system
%   of the form ax'x + b'x + c = 0)
%
%   [Z, R] = FITCIRCLE(X, Property, Value, ...) allows parameters to be
%   passed to the internal Gauss Newton method. Property names can be
%   supplied as any unambiguous contraction of the property name and are
%   case insensitive, e.g. FITCIRCLE(X, 't', 1e-4) is equivalent to
%   FITCIRCLE(X, 'tol', 1e-4). Valid properties are:
%
%       Property:                 Value:
%       --------------------------------
%       maxits                    positive integer, default 100
%           Sets the maximum number of iterations of the Gauss Newton
%           method
%
%       tol                       positive constant, default 1e-5
%           Gauss Newton converges when the relative change in the solution
%           is less than tol
%
%   [X, R, RES] = fitcircle(...) returns the 2 norm of the residual from
%   the least squares fit
%
%   Example:
%       x = [1 2 5 7 9 3; 7 6 8 7 5 7];
%       % Get linear least squares fit
%       [zl, rl] = fitcircle(x, 'linear')
%       % Get true best fit
%       [z, r] = fitcircle(x)
%
%   Reference: "Least-squares fitting of circles and ellipses", W. Gander,
%   G. Golub, R. Strebel - BIT Numerical Mathematics, 1994, Springer

% This implementation copyright Richard Brown, 2007, but is freely
% available to copy, use, or modify as long as this line is maintained

error(nargchk(1, 5, nargin, 'struct'))

% Default parameters for Gauss Newton minimisation
params.maxits = 100;
params.tol    = 1e-5;

% Check x and get user supplied parameters
[x, fNonlinear, params] = parseinputs(x, params, varargin{:});

% Convenience variables
m  = size(x, 2);
x1 = x(1, :)';
x2 = x(2, :)';

% 1) Compute best fit w.r.t. algebraic error using linear least squares
%
% Circle is represented as a matrix quadratic form
%     ax'x + b'x + c = 0
% Linear least squares estimate found by minimising Bu = 0 s.t. norm(u) = 1
%     where u = [a; b; c]

% Form the coefficient matrix
B = [x1.^2 + x2.^2, x1, x2, ones(m, 1)];

% Least squares estimate is right singular vector corresp. to smallest
% singular value of B
[U, S, V] = svd(B);
u = V(:, 4);

% For clarity, set the quadratic form variables
a = u(1);
b = u(2:3);
c = u(4);

z = -b / (2*a);
r = sqrt((norm(b)/(2*a))^2 - c/a);

% 2) Nonlinear refinement to miminise geometric error, and compute residual
if fNonlinear
[z, r, residual] = fitcircle_geometric(x, z, r);
else
residual = norm(B * u);
end

% END MAIN FUNCTION BODY

% NESTED FUNCTIONS
function [z, r, residual] = fitcircle_geometric(x, z0, r0)
% Use a simple Gauss Newton method to minimize the geometric error
fConverged = false;

% Set initial u
u     = [z0; r0];

% Delta is the norm of current step, scaled by the norm of u
delta = inf;
nIts  = 0;

for nIts = 1:params.maxits
% Find the function and Jacobian
[f, J] = sys(u);

% Solve for the step and update u
h = -J \ f;
u = u + h;

% Check for convergence
delta = norm(h, inf) / norm(u, inf);
if delta < params.tol
fConverged = true;
break
end
end

if ~fConverged
warning('fitcircle:FailureToConverge', ...
'Gauss Newton iteration failed to converge');
end
z = u(1:2);
r = u(3);
f = sys(u);
residual = norm(f);

function [f, J] = sys(u)
%SYS   Nonlinear system to be minimised - the objective
%function is the distance to each point from the fitted circle
%contained in u

% Objective function
f = (sqrt(sum((repmat(u(1:2), 1, m) - x).^2)) - u(3))';

% Jacobian
denom = sqrt( (u(1) - x1).^2 + (u(2) - x2).^2 );
J = [(u(1) - x1) ./ denom, (u(2) - x2) ./ denom, repmat(-1, m, 1)];
end % sys

end % fitcircle_geometric

% END NESTED FUNCTIONS

end % fitcircle

function [x, fNonlinear, params] = parseinputs(x, params, varargin)
% Make sure x is 2xN where N > 3
if size(x, 2) == 2
x = x';
end

if size(x, 1) ~= 2
error('fitcircle:InvalidDimension', ...
'Input matrix must be two dimensional')
end

if size(x, 2) < 3
error('fitcircle:InsufficientPoints', ...
'At least 3 points required to compute fit')
end

% determine whether we are measuring geometric error (nonlinear), or
% algebraic error (linear)
fNonlinear = true;
switch length(varargin)
% No arguments means a nonlinear least squares with defaul parameters
case 0
return

% One argument can only be 'linear', specifying linear least squares
case 1
if strncmpi(varargin{1}, 'linear', length(varargin{1}))
fNonlinear = false;
return
else
error('fitcircle:UnknownOption', 'Unknown Option')
end

% Otherwise we're left with user supplied parameters for Gauss Newton
otherwise
if rem(length(varargin), 2) ~= 0
error('fitcircle:propertyValueNotPair', ...
'Additional arguments must take the form of Property/Value pairs');
end

% Cell array of valid property names
properties = {'maxits', 'tol'};

while length(varargin) ~= 0
property = varargin{1};
value    = varargin{2};

% If the property has been supplied in a shortened form, lengthen it
iProperty = find(strncmpi(property, properties, length(property)));
if isempty(iProperty)
error('fitcircle:UnkownProperty', 'Unknown Property');
elseif length(iProperty) > 1
error('fitcircle:AmbiguousProperty', ...
'Supplied shortened property name is ambiguous');
end

% Expand property to its full name
property = properties{iProperty};

switch property
case 'maxits'
if value <= 0
error('fitcircle:InvalidMaxits', ...
'maxits must be an integer greater than 0')
end
params.maxits = value;
case 'tol'
if value <= 0
error('fitcircle:InvalidTol', ...
'tol must be a positive real number')
end
params.tol = value;
end % switch property
varargin(1:2) = [];
end % while

end % switch length(varargin)

end %parseinputs```