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

来源：互联网发布：小米云服务擦除数据编辑：程序博客网时间：2024/05/21 21:34

在工程计算中，经常会遇到数据的拟合处理，拟合处理用的最多的是最小二承法。我曾经做过一个数据处理的软件，其中用到了直线和圆的最小二承法拟合算法，是用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;