最小二乘法拟合圆

最小二乘圆拟合matlab程序:
t=0:0.01:pi;
a=20;%设定圆心X轴数值
b=30;%设定圆心Y轴数值
r=5;%设定圆半径数值
x=a+r*cos(t)+randn(1,315);
y=b+r*sin(t)+randn(1,315);
plot(x,y);
hold on;
x=x(:);
y=y(:);
m=[x y ones(size(x))]\[-(x.^2+y.^2)];
xc = -.5*m(1)%拟合圆心X轴数值
yc = -.5*m(2)%拟合圆心Y轴数值
R = sqrt((m(1)^2+m(2)^2)/4-m(3))%拟合半径数值
plot(xc,yc,'r-x',(xc+R*cos(t)),(yc+R*sin(t)),'r-');
axis equal;



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

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));