Delphi實現直線和圓的最小二承法擬合
阿新 • • 發佈:2019-02-03
在工程計算中,經常會遇到資料的擬合處理,擬合處理用的最多的是最小二承法。我曾經做過一個數據處理的軟體,其中用到了直線和圓的最小二承法擬合算法,是用delphi實現的,現將原始碼貼出來。
一直線擬合
- {
- 直線的方程為形式:y=k*x+b ,x=c;
- X,y為需要處理的資料
- }
- function FitLine(x,y: arrayofdouble;len: integer; var k,b,c: double):boolean;
- var
- i: integer;
- sx,sy,s2x,sxy:double;
- begin
- result:=true;
- if len<
- begin
- result:=false;
- exit;
- end;
- sx:=0;
- sy:=0;
- s2x:=0;
- sxy:=0;
- for i:=0to len-1do
- 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=0then
- c:=x[0]
- else
- begin
- c:=infinity;
- k:=(sy*sx-sxy*len)/(sx*sx-len*s2x);
- b:=(sy-k*sx)/len;
- end;
- end;
二、圓的擬合<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />
- {
- 圓的方程為:(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;