1. 程式人生 > >Delphi實現直線和圓的最小二承法擬合

Delphi實現直線和圓的最小二承法擬合

 

在工程計算中,經常會遇到資料的擬合處理,擬合處理用的最多的是最小二承法。我曾經做過一個數據處理的軟體,其中用到了直線和圓的最小二承法擬合算法,是用delphi實現的,現將原始碼貼出來。

直線擬合

  1. {
  2. 直線的方程為形式:y=k*x+b ,x=c;
  3. X,y為需要處理的資料
  4. }
  5. function  FitLine(x,y: arrayofdouble;len: integervar k,b,c: double):boolean;
  6. var
  7.   i:            integer;
  8.   sx,sy,s2x,sxy:double;
  9. begin
  10.   result:=true;
  11. if len<
    2then//少於兩點無法擬合
  12. begin
  13.     result:=false;
  14.     exit;
  15. end;
  16.   sx:=0;
  17.   sy:=0;
  18.   s2x:=0;
  19.   sxy:=0;
  20. for i:=0to len-1do
  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=0then
  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;

二、圓的擬合<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />

  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;