Delphi实现直线和圆的最小二承法拟合

在工程计算中，经常会遇到数据的拟合处理，拟合处理用的最多的是最小二承法。我曾经做过一个数据处理的软件，其中用到了直线和圆的最小二承法拟合算法，是用delphi实现的，现将源代码贴出来。

一 直线拟合

{

直线的方程为形式：y=k*x+b ,x=c;

X,y为需要处理的数据

}

function FitLine(x,y: array of double;len: integer; var k,b,c: double):boolean;

var

i: integer;

sx,sy,s2x,sxy:double;

begin

result:=true;

if len<2 then //少于两点无法拟合

begin

result:=false;

exit;

end;

sx:=0;

sy:=0;

s2x:=0;

sxy:=0;

for i:=0 to len-1 do

begin

sx:=sx+x[i];

sy:=sy+y[i];

s2x:=s2x+x[i]*x[i];

sxy:=sxy+y[i]*x[i];

end;

if sx*sx-len*s2x=0 then

c:=x[0]

else

begin

c:=infinity;

k:=(sy*sx-sxy*len)/(sx*sx-len*s2x);

b:=(sy-k*sx)/len;

end;

end;

二、圆的拟合

{

圆的方程为：(x-xr)* (x-xr) +(y-yr)* (y-yr)=rad*rad

}

function FitCircle(x,y:array of double; len:integer; var xr,yr,rad: double):boolean;

function D2(x,y:array of double; len:integer):double;

var

i:integer;

sa,sb,sab:double;

begin

sa:=0;

sb:=0;

sab:=0;

for i:=0 to len-1 do

begin

sa:=sa+x[i];

sb:=sb+y[i];

sab:=sab+x[i]*y[i];

end;

result:=sab-(sa*sb)/len;

end;

function D3(x,y:array of double; len:integer):double;

var

i:integer;

sa,sb2,sab2:single;

begin

sa:=0;

sb2:=0;

sab2:=0;

for i:=0 to len-1 do

begin

sa:=sa+X[i];

sb2:=sb2+y[i]*y[i];

sab2:=sab2+X[i]*y[i]*y[i];

end;

result:=sab2-(sa*sb2)/len;

end;

var

i,n:integer;

t1,t2,t3,t4:single;

sum:single;

begin

result:=true;

if len<3 then

begin

result:=false;//少于三点,拟合没有意义退出

exit;

end;

t1:=(D3(x,x,len)+D3(x,y,len))*D2(y,y,len);

t2:=(D3(y,x,len)+D3(y,y,len))*D2(x,y,len);

t3:=2*(D2(x,x,len)*D2(y,y,len)-D2(x,y,len)*D2(x,y,len));

if t3=0 then

begin

result:=false;

exit;

end

else

xr:=(t1-t2)/t3;

t1:=(D3(y,y,len)+D3(y,x,len))*D2(x,x,len);

t2:=(D3(x,y,len)+D3(x,x,len))*D2(x,y,len);

t3:=2*(D2(x,x,len)*D2(y,y,len)-D2(x,y,len)*D2(x,y,len));

if t3=0 then

begin

result:=false;

exit;

end

else

yr:=(t1-t2)/t3;

sum:=0;

for i:=0 to len-1 do

sum:=sum+sqr(x[i]-xr)+sqr(y[i]-yr);

rad:=sqrt(sum/len);

end;