least squares circle

[Steven]

Hello,

Could anyone tell me how to fit a circle to a set of data points? Thanks.

Steven

[Peter Boettcher]

"Steven" <s_f...@21cn.com> writes:

> Could anyone tell me how to fit a circle to a set of data points? Thanks.

Here is a copy of the entry I just added to the FAQ
(http://www.mit.edu/~pwb/cssm/):

Q5.5: How can I fit a circle to a set of XY data?

An elegant chunk of code to perform least-squares circle fitting was
written by Bucher Izhak and has been floating around the newgroup for
some time. The first reference to it that I can find is in
msgid:<3A13371D.A732886D%40skm.com.au>:

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R. A is an optional
% output describing the circle's equation:
%
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0

% by Bucher izhak 25/oct/1991

n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));

--
Peter Boettcher
MIT Lincoln Laboratory
boet...@ll.mit.edu
(781) 981-5275

[Tom Davis]

Total least squares (TLS) curve fitting minimizes the sum of squared
orthogonal distances from the data to the curve. The result is a
maximum likelihood estimator when both the x and y data are subject to
errors that are independently and identically distributed with zero mean
and common variance. The TLS circle fit cannot be accomplished in
closed form.

function [xc,yc,r,fval,exitflag,output] = tlscirc(x,y)
% TLS circle fit
options=optimset('TolFun',1e-10);
[x0,y0]=circfit(x,y); % center estimate
[z,fval,exitflag,output]=fminsearch(@circobj,[x0,y0],options,x,y);
[f,r]=circobj(z,x,y); xc=z(1); yc=z(2);

function [f,r] = circobj(z,x,y)
% TLS circle fit objective f and radius r
n=length(x);
X=x-z(1); Y=y-z(2); X2=X'*X; Y2=Y'*Y;
r=sum(sqrt(X.^2+Y.^2))/n; f=X2+Y2-n*r^2;

The circfit function is used above as an initial estimator. The
objective values of tlscirc and circfit differ by nearly 14% in the
example below.

x=[1;2;3;5;7;9]; y=[7;6;7;8;7;5];
plot(x,y,'bo'), hold on

[xc,yc,r,f]=tlscirc(x,y)
% [xc,yc,r,f]=[4.7398, 2.9835, 4.7142, 1.2276]
rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','b')

[xc,yc,r] = circfit(x,y)
X=x-xc; Y=y-yc; f=sum((sqrt(X.^2+Y.^2)-r).^2)
% [xc,yc,r,f]=[4.7423, 3.8351, 4.1088, 1.3983]
rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','r')

hold off, axis equal

[Tom Davis]

Let X=x-xc and Y=y-yc. Then CIRCFIT minimizes sum((X.^2+Y.^2-r^2).^2).
This problem has been considered by many investigators, e.g. Crawford
(1983), Nucl. Instrum. Methods, 211, 223-225; Moura and Kitney (1991),
Comput. Phys. Comm., 64, 57-63; Coope (1993), J. Optim. Theory Appl.,
76(2), 381-388; and Nievergelt (1994), Comput. Phys. Comm., 81, 342-350.
The geometric interpretation of this minimization is unclear.

In contrast, the total least squares (TLS), or orthogonal distance
regression (ODR), objective is min sum((sqrt(X.^2+Y.^2)-r).^2), i.e.
minimize the sum of squared orthogonal distances from the data to the
circle.

CIRCFIT is an excellent approximation for small residual problems with
well distributed data, a good approximation for large residual problems
with well distributed data, and a poor approximation for problems with
data clustered in a segement of the circle, but nonetheless, a good
initial estimate to start TLSCIRC.

% small residual problem
x=[0.7;3.3;5.6;7.5;6.4;4.4;0.3;-1.1];
y=[4;4.7;4;1.3;-1.1;-3;-2.5;1.3];
figure(1)

plot(x,y,'bo'), hold on

[xc,yc,r,f]=tlscirc(x,y)
% [xc,yc,r,f]=[3.0432, 0.7457, 4.1059, 0.2959]

rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','b')

[xc,yc,r] = circfit(x,y)
X=x-xc; Y=y-yc; f=sum((sqrt(X.^2+Y.^2)-r).^2)

% [xc,yc,r,f]=[3.0603, 0.7436, 4.1091, 0.2973]

rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','r')
hold off, axis equal

% large residual problem
x=[5;10;7;9;2;-1;0;-1];
y=[-3;3;5;9;8;7;3;-1];
figure(2)

plot(x,y,'bo'), hold on

[xc,yc,r,f]=tlscirc(x,y)
% [xc,yc,r,f]=[4.1714, 2.8790, 5.7389, 12.9080]

rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','b')

[xc,yc,r] = circfit(x,y)
X=x-xc; Y=y-yc; f=sum((sqrt(X.^2+Y.^2)-r).^2)

% [xc,yc,r,f]=[4.3160, 3.2446, 5.8362, 13.6179]

rectangle('Position', [xc-r,yc-r,2*r,2*r],...
'Curvature', [1,1],'EdgeColor','r')
hold off, axis equal

[Steven]

Hello,

Thanks for your help.

Steven