Dopasowanie elispy - pomoc w implementacji

Dopasowanie elispy - pomoc w implementacji
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/rbf/CVonline/LOCAL_COPIES/PILU1/demo.html - applet w javie
http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/PILU1/ElliFit.java - kod źrodłowy.

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

Kopiuj
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.

Azarien
  • Rejestracja:około 21 lat
  • Ostatnio:2 minuty
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”.

Kliknij, aby dodać treść...

Pomoc 1.18.8

Typografia

Edytor obsługuje składnie Markdown, w której pojedynczy akcent *kursywa* oraz _kursywa_ to pochylenie. Z kolei podwójny akcent **pogrubienie** oraz __pogrubienie__ to pogrubienie. Dodanie znaczników ~~strike~~ to przekreślenie.

Możesz dodać formatowanie komendami , , oraz .

Ponieważ dekoracja podkreślenia jest przeznaczona na linki, markdown nie zawiera specjalnej składni dla podkreślenia. Dlatego by dodać podkreślenie, użyj <u>underline</u>.

Komendy formatujące reagują na skróty klawiszowe: Ctrl+B, Ctrl+I, Ctrl+U oraz Ctrl+S.

Linki

By dodać link w edytorze użyj komendy lub użyj składni [title](link). URL umieszczony w linku lub nawet URL umieszczony bezpośrednio w tekście będzie aktywny i klikalny.

Jeżeli chcesz, możesz samodzielnie dodać link: <a href="link">title</a>.

Wewnętrzne odnośniki

Możesz umieścić odnośnik do wewnętrznej podstrony, używając następującej składni: [[Delphi/Kompendium]] lub [[Delphi/Kompendium|kliknij, aby przejść do kompendium]]. Odnośniki mogą prowadzić do Forum 4programmers.net lub np. do Kompendium.

Wspomnienia użytkowników

By wspomnieć użytkownika forum, wpisz w formularzu znak @. Zobaczysz okienko samouzupełniające nazwy użytkowników. Samouzupełnienie dobierze odpowiedni format wspomnienia, zależnie od tego czy w nazwie użytkownika znajduje się spacja.

Znaczniki HTML

Dozwolone jest używanie niektórych znaczników HTML: <a>, <b>, <i>, <kbd>, <del>, <strong>, <dfn>, <pre>, <blockquote>, <hr/>, <sub>, <sup> oraz <img/>.

Skróty klawiszowe

Dodaj kombinację klawiszy komendą notacji klawiszy lub skrótem klawiszowym Alt+K.

Reprezentuj kombinacje klawiszowe używając taga <kbd>. Oddziel od siebie klawisze znakiem plus, np <kbd>Alt+Tab</kbd>.

Indeks górny oraz dolny

Przykład: wpisując H<sub>2</sub>O i m<sup>2</sup> otrzymasz: H2O i m2.

Składnia Tex

By precyzyjnie wyrazić działanie matematyczne, użyj składni Tex.

<tex>arcctg(x) = argtan(\frac{1}{x}) = arcsin(\frac{1}{\sqrt{1+x^2}})</tex>

Kod źródłowy

Krótkie fragmenty kodu

Wszelkie jednolinijkowe instrukcje języka programowania powinny być zawarte pomiędzy obróconymi apostrofami: `kod instrukcji` lub ``console.log(`string`);``.

Kod wielolinijkowy

Dodaj fragment kodu komendą . Fragmenty kodu zajmujące całą lub więcej linijek powinny być umieszczone w wielolinijkowym fragmencie kodu. Znaczniki ``` lub ~~~ umożliwiają kolorowanie różnych języków programowania. Możemy nadać nazwę języka programowania używając auto-uzupełnienia, kod został pokolorowany używając konkretnych ustawień kolorowania składni:

```javascript
document.write('Hello World');
```

Możesz zaznaczyć również już wklejony kod w edytorze, i użyć komendy  by zamienić go w kod. Użyj kombinacji Ctrl+`, by dodać fragment kodu bez oznaczników języka.

Tabelki

Dodaj przykładową tabelkę używając komendy . Przykładowa tabelka składa się z dwóch kolumn, nagłówka i jednego wiersza.

Wygeneruj tabelkę na podstawie szablonu. Oddziel komórki separatorem ; lub |, a następnie zaznacz szablonu.

nazwisko;dziedzina;odkrycie
Pitagoras;mathematics;Pythagorean Theorem
Albert Einstein;physics;General Relativity
Marie Curie, Pierre Curie;chemistry;Radium, Polonium

Użyj komendy by zamienić zaznaczony szablon na tabelkę Markdown.

Lista uporządkowana i nieuporządkowana

Możliwe jest tworzenie listy numerowanych oraz wypunktowanych. Wystarczy, że pierwszym znakiem linii będzie * lub - dla listy nieuporządkowanej oraz 1. dla listy uporządkowanej.

Użyj komendy by dodać listę uporządkowaną.

1. Lista numerowana
2. Lista numerowana

Użyj komendy by dodać listę nieuporządkowaną.

* Lista wypunktowana
* Lista wypunktowana
** Lista wypunktowana (drugi poziom)

Składnia Markdown

Edytor obsługuje składnię Markdown, która składa się ze znaków specjalnych. Dostępne komendy, jak formatowanie , dodanie tabelki lub fragmentu kodu są w pewnym sensie świadome otaczającej jej składni, i postarają się unikać uszkodzenia jej.

Dla przykładu, używając tylko dostępnych komend, nie możemy dodać formatowania pogrubienia do kodu wielolinijkowego, albo dodać listy do tabelki - mogłoby to doprowadzić do uszkodzenia składni.

W pewnych odosobnionych przypadkach brak nowej linii przed elementami markdown również mógłby uszkodzić składnie, dlatego edytor dodaje brakujące nowe linie. Dla przykładu, dodanie formatowania pochylenia zaraz po tabelce, mogłoby zostać błędne zinterpretowane, więc edytor doda oddzielającą nową linię pomiędzy tabelką, a pochyleniem.

Skróty klawiszowe

Skróty formatujące, kiedy w edytorze znajduje się pojedynczy kursor, wstawiają sformatowany tekst przykładowy. Jeśli w edytorze znajduje się zaznaczenie (słowo, linijka, paragraf), wtedy zaznaczenie zostaje sformatowane.

  • Ctrl+B - dodaj pogrubienie lub pogrub zaznaczenie
  • Ctrl+I - dodaj pochylenie lub pochyl zaznaczenie
  • Ctrl+U - dodaj podkreślenie lub podkreśl zaznaczenie
  • Ctrl+S - dodaj przekreślenie lub przekreśl zaznaczenie

Notacja Klawiszy

  • Alt+K - dodaj notację klawiszy

Fragment kodu bez oznacznika

  • Alt+C - dodaj pusty fragment kodu

Skróty operujące na kodzie i linijkach:

  • Alt+L - zaznaczenie całej linii
  • Alt+, Alt+ - przeniesienie linijki w której znajduje się kursor w górę/dół.
  • Tab/⌘+] - dodaj wcięcie (wcięcie w prawo)
  • Shit+Tab/⌘+[ - usunięcie wcięcia (wycięcie w lewo)

Dodawanie postów:

  • Ctrl+Enter - dodaj post
  • ⌘+Enter - dodaj post (MacOS)