Dopasowanie elispy - pomoc w implementacji

Odpowiedz Nowy wątek
2011-10-20 19:25
Progrmaistyczny Nowi
0

Cześć,
Próbuje bezskutecznie zaimplementować algorytm napisany w javie, służący do dopasowywania elipsy do zbioru zadanych punktów.

http://homepages.inf.ed.ac.uk[...]/LOCAL_COPIES/PILU1/demo.html - applet w javie
http://homepages.inf.ed.ac.uk[...]CAL_COPIES/PILU1/ElliFit.java - kod źrodłowy.

Poniżej wklejam moją własną implementację, która z jakichś powodów nie działa:

unit EllipseFit;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Math, StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Memo1: TMemo;
    Label7: TLabel;
    Label8: TLabel;
    procedure Button1Click(Sender: TObject);

  private
    { Private declarations }
  public
    { Public declarations }
  end;
  TDoubleM = array of array of Double;
  TDoubleV = array  of Double;
  TPointsV = record
    x,y:integer;
  end;
var
  Form1: TForm1;

  ControlPoints: array of TPointsV;

implementation

{$R *.dfm}

procedure multMatrix(m,g,mg:TDoubleM);
var i,j,k:integer;
begin
  for i:=0 to 3 do
    for j:=0 to 1 do
      mg[i,j]:=0;

  for i:=0 to 3 do
    for j:=0 to 1 do
      for k:=0 to 3 do
        mg[i,j]:=mg[i,j] + (m[i,k]*g[k,j]);
end;

procedure rotate(a:TDoubleM; i,j,k,l:integer; tau,s:double);
var  g,h:double;
begin
  g:=a[i,j];
  h:=a[k,l];
  a[i,j]:=g-s*(h+g*tau);
  a[k,l]:=h+s*(g-h*tau);
end;

procedure jacobi(a:TDoubleM; n:integer; d:TDoubleV; v:TDoubleM; nrot:integer);
var j,iq,ip,i,licz1,licz2:integer;
    tresh,theta,tau,t,sm,s,h,g,c:double;
    b,z:TDoubleV;
begin
  SetLength(b,n+1);
  SetLength(z,n+1);

  for ip:=1 to n do
  begin
    for iq:=1 to n do v[ip,iq]:=0;
    v[ip,ip]:=1;
  end;

  for ip:=1 to n do
  begin
    b[ip]:=d[ip];
    d[ip]:=a[ip,ip];
    z[ip]:=0;
  end;
  nrot:=0; licz1:=0; licz2:=0;
  for i:=1 to 50 do
  begin
    sm:=0;
    for ip:=1 to n-1 do
    begin
      for iq:=ip+1 to n do sm:=sm+abs(a[ip,iq]);
    end;
    if sm=0 then
    begin
      licz1:=licz1+1;
      Form1.Label7.Caption:='N '+inttostr(licz1);
    end else begin
      licz2:=licz2+1;
      Form1.Label8.Caption:='T '+inttostr(licz2);
    end;
    if (i<4) then tresh:=0.2*sm/(n*n) else tresh:=0;
    for ip:=1 to n-1 do
    begin
      for iq:=ip+1 to n do
      begin
        g:=100*abs(a[ip,iq]);
        if (i>4)and(abs(d[ip])+g = abs(d[ip]))and(abs(d[iq])+g = abs(d[iq])) then
          a[ip,iq]:=0
        else if abs(a[ip,iq]) > tresh then
        begin
          h:=d[iq]-d[ip];
          if (abs(h)+g = abs(h)) then t:=(a[ip,iq])/h
          else begin
            theta:=0.5*h/(a[ip,iq]);
            t:=1/(abs(theta)+sqrt(1+sqr(theta)));
            if (theta < 0) then t:= -t;
          end;
          c:=1/sqrt(1+sqr(t));
          s:=t*c;
          tau:=s/(1+c);
          h:=t*a[ip,iq];
          z[ip]:=z[ip] - h;
          z[iq]:=z[iq] + h;
          d[ip]:=d[ip] - h;
          d[iq]:=d[iq] + h;
          a[ip,iq]:=0;
          for j:=1 to ip-1 do rotate(a,j,ip,j,iq,tau,s);
          for j:=ip+1 to iq-1 do rotate(a,ip,j,j,iq,tau,s);
          for j:=iq+1 to n do rotate(a,ip,j,iq,j,tau,s);
          for j:=1 to n do rotate(v,j,ip,j,iq,tau,s);
          nrot:=nrot+1;
        end;
      end;
    end;
    for ip:=1 to n do
    begin
      b[ip]:=b[ip] + z[ip];
      d[ip]:=b[ip];
      z[ip]:=0;
    end;
  end;

  b:=Nil;
  z:=Nil;
end;

procedure choldc(a:TDoubleM; n:integer; l:TDoubleM);
var i,j,k:integer;
    sum:double;
    p:TDoubleV;
begin
  SetLength(p,n+1);

  for i:=1 to n do
  begin
    for j:=i  to n do
    begin
      sum:=a[i,j];
      k:=i-1;
      while k >= 1 do
      begin
        sum:=sum - a[i,k]*a[j,k];
        k:=k-1;
      end;
      if i=j then
      begin
        if sum <= 0 then
        else p[i]:=sqrt(sum);
      end else
        a[j,i]:=sum/p[i];
    end;
  end;
  for i:=1 to n do
    for j:=1 to n do
      if i=j then l[i,i]:=p[i]
      else begin
        l[j,i]:=a[j,i];
        l[i,j]:=0;
      end;

  p:=Nil;
end;

function inverse(TB,InvB:TDoubleM; N:integer):integer;
var k,i,j,p,q:integer;
    mult,D,temp,maxpivot,eps:double;
    npivot:integer;

    B,A,C:TDoubleM;
begin
  eps:=10e-20;
  SetLength(B,N+1,N+2);
  SetLength(A,N+1,2*N+2);
  SetLength(C,N+1,N+1);

  for k:=1 to N do
    for j:=1 to N do B[k,j]:=TB[k,j];

  for k:=1 to N do
  begin
    for j:=1 to N+1 do A[k,j]:=B[k,j];
    for j:=N+2 to 2*N+1 do A[k,j]:=0;
    A[k,k-1+N+2]:=1;
  end;
  for k:=1 to N do
  begin
    maxpivot:=abs(A[k,k]);
    npivot:=k;
    for i:=k to N do
      if (maxpivot < abs(A[i,k])) then
      begin
        maxpivot:=abs(A[i,k]);
        npivot:=i;
      end;
      if maxpivot > eps then
      begin
        if npivot = k then
          for j:=k to 2*N+1 do
          begin
            temp:=A[npivot,j];
            A[npivot,j]:=A[k,j];
            A[k,j]:=temp;
          end;
        D:=A[k,k];
        for j:=2*N+1 downto k do A[k,j]:=A[k,j]/D;
        for i:=1 to N do
        begin
          if i<>k then
          begin
            mult:=A[i,k];
            for j:=2*N+1 downto k do A[i,j]:=A[i,j]-mult*A[k,j];
          end;
        end;
      end else Result:=-1;
  end;

  p:=1;
  for k:=1 to N do
  begin
    q:=1;
    for j:=N+2 to 2*N+1 do
    begin
      InvB[p,q]:=A[k,j];
      q:=q+1;
    end;
    p:=p+1;
  end;
  Result:=0;

  A:=Nil; B:=Nil; C:=Nil;
end;

procedure AperB(_A,_B,_res:TDoubleM; _righA,_colA,_righB,_colB:integer);
var p,q,l:integer;
begin
  for p:=1 to _righA do
    for q:=1 to _colB do
    begin
      _res[p,q]:=0;
      for l:=1 to _colA do _res[p,q]:=_res[p,q]+_A[p,l]*_B[l,q];
    end;
end;

procedure A_TperB(_A,_B,_res:TDoubleM; _righA,_colA,_righB,_colB:integer);
var p,q,l:integer;
begin
  for p:=1 to _colA do
    for q:=1 to _colB do
    begin
      _res[p,q]:=0;
      for l:=1 to _righA do _res[p,q]:=_res[p,q]+_A[l,p]*_B[l,q];
    end;
end;

procedure AperB_T(_A,_B,_res:TDoubleM; _righA,_colA,_righB,_colB:integer);
var p,q,l:integer;
begin
  for p:=1 to _colA do
    for q:=1 to _colB do
    begin
      _res[p,q]:=0;
      for l:=1 to _righA do _res[p,q]:=_res[p,q]+_A[p,l]*_B[q,l];
    end;
end;

procedure draw_conic(pvec:TDoubleV; nptsk:integer; points:TDoubleM);
var npts:integer;
    u,Aiu,L,B,Xpos,Xneg,ss1,ss2:TDoubleM;
    lambda:TDoubleV;
    uAiu,A,Ai,Aib,bs,r1:TDoubleM;

    pi,theta,kk:double;
    i,j:integer;

    Ao,Ax,Ay,Axx,Ayy,Axy:double;
begin
  SetLength(u,3,npts+1);
  SetLength(Aiu,3,npts+1);
  SetLength(L,3,npts+1);
  SetLength(B,3,npts+1);
  SetLength(Xpos,3,npts+1);
  SetLength(Xneg,3,npts+1);
  SetLength(ss1,3,npts+1);
  SetLength(ss2,3,npts+1);
  SetLength(lambda,npts+1);
  SetLength(uAiu,3,npts+1);
  SetLength(A,3,3);
  SetLength(Ai,3,3);
  SetLength(Aib,3,2);
  SetLength(bs,3,2);
  SetLength(r1,2,2);

  pi:=3.14781;
  npts:=round(nptsk/2);

  Ao:=pvec[6];
  Ax:=pvec[4];
  Ay:=pvec[5];
  Axx:=pvec[1];
  Ayy:=pvec[3];
  Axy:=pvec[2];

  A[1,1]:=Axx; A[1,2]:=Axy/2;
  A[2,1]:=Axy/2; A[2,2]:=Ayy;
  b[1,1]:=Ax; b[2,1]:=Ay;

  theta:=0;
  for i:=1 to npts do
  begin
    u[1,i]:=cos(theta);
    u[2,i]:=sin(theta);
    theta:=theta + pi/npts;
  end;

  inverse(A,Ai,2);

  AperB(Ai,bs,Aib,2,2,2,1);
  A_TperB(bs,Aib,r1,2,1,2,1);
  r1[1,1]:=r1[1,1] - 4*Ao;

  AperB(Ai,u,Aiu,2,2,2,npts);
  for i:=1 to 2 do
    for j:=1 to npts do uAiu[i,j]:=u[i,j]*Aiu[i,j];

  for j:=1 to npts do
  begin
    kk:=r1[1,1]/(uAiu[1,j]+uAiu[2,j]);
    if (kk >= 0) then
      lambda[j]:=sqrt(kk)
    else
      lambda[j]:= -1;
  end;

  for j:=1 to npts do
  begin
    L[1,j]:=L[2,j];
    L[2,j]:=lambda[j];
  end;
  for j:=1 to npts do
  begin
    B[1,j]:=bs[1,1];
    B[2,j]:=bs[2,1];
  end;

  for j:=1 to npts do
  begin
    ss1[1,j]:=0.5*(L[1,j]*u[1,j]-B[1,j]);
    ss1[2,j]:=0.5*(L[2,j]*u[2,j]-B[2,j]);
    ss2[1,j]:=0.5*(-L[1,j]*u[1,j]-B[1,j]);
    ss2[2,j]:=0.5*(-L[2,j]*u[2,j]-B[2,j]);
  end;

  AperB(Ai,ss1,Xpos,2,2,2,npts);
  AperB(Ai,ss2,Xneg,2,2,2,npts);

  for j:=1 to npts do
  begin
    if lambda[j]=-1 then
    begin
      points[1,j]:=-1;
      points[2,j]:=-1;
      points[1,j+npts]:=-1;
      points[2,j+npts]:=-1
    end else begin
      points[1,j]:= Xpos[1,j];
      points[2,j]:= Xpos[2,j];
      points[1,j+npts]:= Xneg[1,j];
      points[2,j+npts]:= Xneg[2,j];
    end;
  end;

  u:=Nil; Aiu:=Nil; L:=Nil; B:=Nil; Xpos:=Nil; Xneg:=Nil; ss1:=Nil; ss2:=Nil;
  lambda:=Nil; uAiu:=Nil; A:=Nil; Ai:=Nil; Aib:=Nil; bs:=Nil; r1:=Nil;
end;

procedure TForm1.Button1Click(Sender: TObject);
var np,i,j,k:integer;
    D,S,Constr,temp,L,C,invL,V,sol,XY:TDoubleM;
    ds,pvec:TDoubleV;
    tx,ty,mod1,zero,minev:double;
    nrot,npts,solind:integer;
begin
  np:=10;
  nrot:=0;
  npts:=50;

  SetLength(D,np+1,7);
  SetLength(S,7,7);
  SetLength(Constr,7,7);
  SetLength(temp,7,7);
  SetLength(L,7,7);
  SetLength(C,7,7);
  SetLength(invL,7,7);
  SetLength(ds,7);
  SetLength(V,7,7);
  SetLength(sol,7,7);
  SetLength(XY,3,npts+1);
  SetLength(pvec,7);

  Constr[1,3]:=-2;
  Constr[2,2]:=1;
  Constr[3,1]:=-2;

  randomize;
  for i:=1 to np do
  begin
    tx:=random(10)-random(10);//ControlPoints[i-1].x;
    ty:=random(10)-random(10);//ControlPoints[i-1].y;
    D[i,1]:=sqr(tx);
    D[i,2]:=tx*ty;
    D[i,3]:=sqr(ty);
    D[i,4]:=tx;
    D[i,5]:=ty;
    D[i,6]:=1;
  end;

  A_TperB(D,D,S,np,6,np,6);
  choldc(S,6,L);

  Form1.Memo1.Clear; k:=1;
  for i:=0 to 6 do
    for j:=0 to 6 do
    begin
      Form1.Memo1.Lines.Add('');
      Form1.Memo1.Lines[k]:=floattostr(Constr[j,i]);
      k:=k+1;
    end;

  inverse(L,invL,6);
  AperB_T(Constr,invL,temp,6,6,6,6);
  AperB(invL,temp,C,6,6,6,6);
  jacobi(C,6,ds,V,nrot);
  A_TperB(invL,V,sol,6,6,6,6);

  for j:=1 to 6 do
  begin
    mod1:=0;
    for i:=1 to 6 do mod1:=mod1+sol[i,j]*sol[i,j];
    for i:=1 to 6 do sol[i,j]:=sqrt(mod1);
  end;

  zero:=10e-20;
  minev:=10e+20;
  solind:=0;
  for i:=1 to 6 do
    if (ds[i]<0)and(abs(ds[i])>zero) then solind:=1;

  for j:=1 to 6 do
    pvec[j]:=sol[j,solind];

  Form1.Label1.Caption:=FloattoStr(pvec[1]);
  Form1.Label2.Caption:=FloattoStr(pvec[2]);
  Form1.Label3.Caption:=FloattoStr(pvec[3]);
  Form1.Label4.Caption:=FloattoStr(pvec[4]);
  Form1.Label5.Caption:=FloattoStr(pvec[5]);
  Form1.Label6.Caption:=FloattoStr(pvec[6]);

  D:=Nil; S:=Nil; Constr:=Nil; temp:=Nil; L:=Nil; C:=Nil; invL:=Nil;
  ds:=Nil; V:=Nil; sol:=Nil; XY:=Nil; pvec:=Nil;
end;

end.

Generalnie dużo tego nie ma, i jeśli ktoś chciałby samemu spróbować napisać ten algorytm na podstawie podanych powyżej źródeł to byłbym bardzo wdzięczny. Myślę, że taki algorytm przyda się wielu osobom, więc dobrze by było aby w sieci taki kod był dostępny :-)

Serdecznie pozdrawiam.

Pozostało 580 znaków

2011-10-21 10:14
0

sprawdź czy poszczególne funkcje (jacobi, inverse, multMatrix) zwracają te same wartości dla tych samych danych w javie i delphi.

zlokalizuj problem, nie poprzestawaj na jedynie „nie działa”.

Pozostało 580 znaków

Odpowiedz
Liczba odpowiedzi na stronę

1 użytkowników online, w tym zalogowanych: 0, gości: 1, botów: 0