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

 

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

 

一 直线拟合

1. {
2. 直线的方程为形式：y=k*x+b ,x=c;
3. X,y为需要处理的数据
4. }
5. function  FitLine(x,y: array of double;len: integer; var k,b,c: double):boolean;
6. var
7.   i:            integer;
8.   sx,sy,s2x,sxy:double;
9. begin
10.   result:=true;
11.   if len<2 then    //少于两点无法拟合
12.   begin
13.     result:=false;
14.     exit;
15.   end;
16.   sx:=0;
17.   sy:=0;
18.   s2x:=0;
19.   sxy:=0;
20.   for i:=0 to len-1 do
21.   begin
22.     sx:=sx+x[i];
23.     sy:=sy+y[i];
24.     s2x:=s2x+x[i]*x[i];
25.     sxy:=sxy+y[i]*x[i];
26.   end;
27.   if sx*sx-len*s2x=0 then
28.     c:=x[0]
29.   else
30.   begin
31.     c:=infinity;
32.     k:=(sy*sx-sxy*len)/(sx*sx-len*s2x);
33.     b:=(sy-k*sx)/len;
34.   end;
35. end;

二、圆的拟合

1. {
2.  圆的方程为：(x-xr)* (x-xr) +(y-yr)* (y-yr)=rad*rad 
3. }
4. function FitCircle(x,y:array of double; len:integer; var xr,yr,rad: double):boolean;
5.   function D2(x,y:array of double; len:integer):double;
6.   var
7.     i:integer;
8.     sa,sb,sab:double;
9.   begin
10.     sa:=0;
11.     sb:=0;
12.     sab:=0;
13.     for i:=0 to len-1 do
14.     begin
15.       sa:=sa+x[i];
16.       sb:=sb+y[i];
17.       sab:=sab+x[i]*y[i];
18.     end;
19.     result:=sab-(sa*sb)/len;
20.   end;
21.   function D3(x,y:array of double; len:integer):double;
22.   var
23.     i:integer;
24.     sa,sb2,sab2:single;
25.   begin
26.     sa:=0;
27.     sb2:=0;
28.     sab2:=0;
29.     for i:=0 to len-1 do
30.     begin
31.       sa:=sa+X[i];
32.       sb2:=sb2+y[i]*y[i];
33.       sab2:=sab2+X[i]*y[i]*y[i];
34.     end;
35.     result:=sab2-(sa*sb2)/len;
36.   end;
37. var
38.   i,n:integer;
39.   t1,t2,t3,t4:single;
40.   sum:single;
41. begin
42.   result:=true;
43.   if len<3 then
44.   begin
45.     result:=false;//少于三点,拟合没有意义退出
46.     exit;
47.   end;
48.   t1:=(D3(x,x,len)+D3(x,y,len))*D2(y,y,len);
49.   t2:=(D3(y,x,len)+D3(y,y,len))*D2(x,y,len);
50.   t3:=2*(D2(x,x,len)*D2(y,y,len)-D2(x,y,len)*D2(x,y,len));
51.   if t3=0 then
52.   begin
53.     result:=false;
54.     exit;
55.   end
56.   else
57.     xr:=(t1-t2)/t3;
58.   t1:=(D3(y,y,len)+D3(y,x,len))*D2(x,x,len);
59.   t2:=(D3(x,y,len)+D3(x,x,len))*D2(x,y,len);
60.   t3:=2*(D2(x,x,len)*D2(y,y,len)-D2(x,y,len)*D2(x,y,len));
61.   if t3=0 then
62.   begin
63.     result:=false;
64.     exit;
65.   end
66.   else
67.     yr:=(t1-t2)/t3;
68.   sum:=0;
69.   for i:=0 to len-1 do
70.     sum:=sum+sqr(x[i]-xr)+sqr(y[i]-yr);
71.   rad:=sqrt(sum/len);
72. end;