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