精华区文章阅读

发信人: zjliu (秋天的萝卜), 信区: Matlab
标  题: 圆的拟合的方法zz
发信站: 哈工大紫丁香 (Wed Dec 31 09:16:35 2003), 站内信件

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

Peter Boettcher wrote:
>

> "Steven" <s_f_lu@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

> boettcher@ll.mit.edu

> (781) 981-5275


--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝

※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 202.118.229.162]

Matlab 版 (精华区)