Ellipse-detection
7-1-2023

Moduli

Mahalanobis distance-like function

dM[xpod_,ypod_,S_]:=Sqrt[Det[S]](xpod-ypod).Inverse[S].(xpod-ypod);
dm[xpod_,ypod_,S_]:=(xpod-ypod).Inverse[S].(xpod-ypod);

Distance from a point to the circle

(* LAD distance *)
Dist[T_,c_,r_]:=Module[{},Return[Abs[Norm[T-c]-r] ]
                                             ];ind1=1;
			(*  TLS quasi distance *)
Dist[T_,c_,r_]:=Module[{},Return[(Norm[T-c]-r)^2 ]
                                             ];ind1=2;
			(* Algebraic fit criterion  *)
Dist[T_,c_,r_]:=Module[{},Return[(Norm[T-c]^2-r^2)^2]
                                             ];ind1=3;

Distance from a point to the M-circle

(* LAD distance *)
indikator=1;
DistM[T_,Circ_]:=Module[{},Return[Abs[Sqrt[dM[T,Circ[[1]],Circ[[3]]]]-Circ[[2]]]]
                                             ];
(* TLS quasi distance *)
indikator=2;
DistMTLS[T_,Circ_]:=Module[{},Return[Abs[Sqrt[dM[T,Circ[[1]],Circ[[3]]]]-Circ[[2]]]^2  ]
                                             ];
(* Algebraic fit criterion *)
indikator=3;
DistM[T_,Circ_]:=Module[{},Return[Abs[dM[T,Circ[[1]],Circ[[3]]]-Circ[[2]]^2] ^2]
                                             ];

Objective functions for circles and M-circles

WFLS[PI_,Circ_]:=Module[{},
Return[
Sum[Sum[PI[[j,s,1]]Dist[PI[[j,s,2]],Circ[[j,1]],Circ[[j,2]]],{s,Length[PI[[j]]]}]
,{j,Length[PI]}]
]]

Fmin[w_,A_,Circ_]:=Module[{k,m},
m=Length[A];k=Length[Circ];
Sum[w[[i]]Min[Table[Dist[A[[i]],Circ[[j,1]],Circ[[j,2]]],{j,k}]],{i,m}]
]

WFLSM[PI_,Circ_]:=Module[{},
Return[Sum[PI[[j,s,1]]*DistM[PI[[j,s,2]],Circ[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

Weighted k - closest circles algorithm

WKC[w_,A_,Circ_,Iter_]:=Module[{It=0,m=Length[w],k=Length[Circ],Kr,PI,tab,fc,min,imin,Pod,F0,F1,p,q,r,dom,min1,fCompile,FF},
Kr=Table[0,{j,k}];
Pod=Table[{w[[i]],A[[i]]},{i,m}];
						(* Odredjivanje particije *)
PI=Table[{},{i,k}];
Do[
tab=Table[Dist[Pod[[i,2]],Circ[[j,1]],Circ[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F1=WFLS[PI,Circ];
                                              (* Računanje kružnica-centara *)
Do[
 {p0,q0}=Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]];
r0=Sqrt[Sum[Norm[{p0,q0}-PI[[j,s,2]]]^2,{s,Length[PI[[j]]]}]/Length[PI[[j]]]];
fc[p_,q_,r_]:=Sum[PI[[j,s,1]]*Dist[PI[[j,s,2]],{p,q},r],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,r],{p,p0},{q,q0},{r,r0}];
Kr[[j]]={{p/.min[[2]],q/.min[[2]]},r/.min[[2]]};
,{j,k}];
F0=WFLS[PI,Kr]; 
While[
Chop[F0-F1]≠0&&It<Iter,
F0=F1;It=It+1;
PI=Table[{},{i,k}];
Do[
tab=Table[Dist[Pod[[i,2]],Kr[[j,1]],Kr[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Kr=Table[0,{j,k}];
Do[
{p0,q0}=Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]];
r0=Sqrt[Sum[Norm[{p0,q0}-PI[[j,s,2]]]^2,{s,Length[PI[[j]]]}]/Length[PI[[j]]]];
fc[p_,q_,r_]:=Sum[PI[[j,s,1]]*Dist[PI[[j,s,2]],{p,q},r],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,r],{p,p0},{q,q0},{r,r0}];
Kr[[j]]={{p/.min[[2]],q/.min[[2]]},r/.min[[2]]};
,{j,k}];
F1=WFLS[PI,Kr];
];
{PI,Kr,F1,It}
]

Weighted k - closest Mahalanobis circles algorithm

WKMC[w_,A_,Circ_,Iter_]:=Module[{m=Length[w],k=Length[Circ],n=Length[A[[1]]],PI,tab,centroid,p,q,r,imin,Pod,F0,F1,IT=0,r0},
Clear[dist,p,q,r,r0,al1,al2,al3,Kr,Cov];
Kr=Cov=Table[{},{j,k}]; 
Pod=Table[{w[[i]],A[[i]]},{i,m}];
						(* Odredjivanje particije *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[DistM[Pod[[i,2]],Circ[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}]; 
k=Length[PI];
F1=WFLSM[PI,Circ];  
dist[jj_,p_,q_,r_,al1_,al2_,al3_]:=Sum[PI[[jj,s,1]]*DistM[PI[[jj,s,2]],{{p,q},r,{{al1,al2},{al2,al3}} }],{s,Length[PI[[jj]]]}];
                           (* Računanje centara i radijusa Mahalanobis kruznica *)
centroid=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}]; 
Do[Cov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,s,2]]-centroid[[j]]),(PI[[j,s,2]]-centroid[[j]])],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
,{j,k}];
Do[
If[indikator==1,r0=Median[Sqrt[Table[dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]],{s,Length[PI[[j]]]}]]],
If[indikator==2,r0=Sum[PI[[j,s,1]]*Sqrt[dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],r0=Sqrt[Sum[PI[[j,s,1]]*dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}]]
]];
min=FindMinimum[dist[j,p,q,r,al1,al2,al3],{{p,centroid[[j,1]]},{q,centroid[[j,2]]},{r,r0},{al1,Cov[[j,1,1]]},{al2,Cov[[j,1,2]]},{al3,Cov[[j,2,2]]}},Method->"PrincipalAxis"]; 
Kr[[j]]={{p/.min[[2]],q/.min[[2]]},r/.min[[2]],{{al1/.min[[2]],al2/.min[[2]]},{al2/.min[[2]],al3/.min[[2]]}}};
,{j,k}];
F0=WFLSM[PI,Kr];
(*Print["{IT,F0}:",{IT,F0}];*)
(*
slcentroid=ListPlot[Table[Kr[[j,1]],{j,k}],PlotStyle->{Red,PointSize[.03],Opacity[.5]}];
Do[
sle[j]=ContourPlot[dM[Kr[[j,1]],{x,y},Kr[[j,3]]]==Kr[[j,2]]^2,
{x,ax,bx},{y,ay,by},ContourStyle->{Red}],
{j,k}];
Print[Show[slpod,slc,slcentroid,Table[sle[j],{j,k}],elipse0, AxesOrigin->{0,0},AspectRatio->Automatic,ImageSize->200]];
*)
						(*  *)
While[
Chop[F0-F1]≠0&&IT<Iter,IT=IT+1;
F0=F1;
PI=Table[{},{i,k}];
Do[
tab=Table[DistM[Pod[[i,2]],Kr[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
centroid=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}];
Do[Cov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,s,2]]-centroid[[j]]),(PI[[j,s,2]]-centroid[[j]])],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
,{j,k}];
Do[
If[indikator==1,r0=Median[Sqrt[Table[dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]],{s,Length[PI[[j]]]}]]],
If[indikator==2,r0=Sum[PI[[j,s,1]]*Sqrt[dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],r0=Sqrt[Sum[PI[[j,s,1]]*dM[centroid[[j]],PI[[j,s,2]],Cov[[j]]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}]]
]];
min=FindMinimum[dist[j,p,q,r,al1,al2,al3],{{p,centroid[[j,1]]},{q,centroid[[j,2]]},{r,r0},{al1,Cov[[j,1,1]]},{al2,Cov[[j,1,2]]},{al3,Cov[[j,2,2]]}},Method->"PrincipalAxis"]; 
Kr[[j]]={{p/.min[[2]],q/.min[[2]]},r/.min[[2]],{{al1/.min[[2]],al2/.min[[2]]},{al2/.min[[2]],al3/.min[[2]]}}};
,{j,k}];
F1=WFLSM[PI,Kr];
(*
slcentroid=ListPlot[Table[Kr[[j,1]],{j,k}],PlotStyle->{Red,PointSize[.03],Opacity[.5]}];
Do[
sle[j]=ContourPlot[dM[Kr[[j,1]],{x,y},Kr[[j,3]]]==Kr[[j,2]]^2,
{x,ax,bx},{y,ay,by},ContourStyle->{Red}],
{j,k}];
Print[Show[slpod,slc,slcentroid,Table[sle[j],{j,k}],elipse0, AxesOrigin->{0,0},AspectRatio->Automatic,ImageSize->200]];
*)
];
Krug=Kr;(*Do[Krug[[j,1]]=KK[[j,1]],{j,k}];*)
(*Print["F0:",WFLSM[PI,Krug]];*)
{PI,Kr,F1,IT}
]

Theta

Theta[A_]:=
Module[{a=A[[1,1]],b=A[[2,2]],c=A[[1,2]],d=A[[2,1]],k}, 
Which[c≠d, Print["Matrica nije simetrična"],
a<0 || a b-c^2<0, Print["Matrica nije pozitivno definitna"],
a==b && c==0, Print["Matrici ne možemo odrediti kut"],
a>b && c≥ 0, k=(1/2) ArcTan[((2c)/(a-b))]//N,
a≤ b && c>0,  k=(1/2) ArcCos[((a-b)/(Sqrt[(((a-b))^2)+4 (c^2)]))]//N,
a<b && c≤0, k=(π/2)+(1/2) ArcTan[((2c)/(a-b))]//N, 
a≥ b && c<0, k=π+(1/2) ArcSin[ ((2 c)/(Sqrt[(((a-b))^2)+4 (c^2)]))]//N
];k];

Transformation : M - circle -> Ellipse (parametrical)

EllParam[Mcircle_]:=Module[{wwN,ewwN,thN},
wwN=Mcircle[[3]]/Sqrt[Det[Mcircle[[3]]]];
ewwN=Eigensystem[wwN];
(* theta New *)
 thN=Theta[Mcircle[[3]]];
{Mcircle[[1]],Mcircle[[2]]*Sqrt[ewwN[[1]]],thN}
]

Similarity

(* EllS: Original ellipses zadana parametarski; EllN: calculated ellipses zadana kao M-kružnica*)

Scn2[EllS_,EllN_,Ind_]:=Module[{sim,sim1,komp,komp1,EllSM,EllS0,EllS1,sa,sb,sp,eigS,eigN,Nr=0},
EllSM=Table[{EllS[[j,1]],Sqrt[EllS[[j,2,1]] EllS[[j,2,2]]], Q[EllS[[j,3]]]*DiagonalMatrix[EllS[[j,2]]^2]*Transpose[Q[EllS[[j,3]]]]},{j,Length[EllS]}];
EllS1=WKMC[w,A,EllSM,10];ELLS0=EllS1[[2]];
Do[
 eigN=Eigensystem[EllN[[s,3]]/Sqrt[Det[EllN[[s,3]]]]];
sax=EllN[[s,2]]*Sqrt[eigN[[1]]]; 
sim=0;komp=0;
Do[
eigS=Eigensystem[ELLS0[[j,3]]];
sa=Norm[ELLS0[[j,1]]-EllN[[s,1]]];
sb=Abs[ELLS0[[j,2]]-EllN[[s,2]]];
sp=Abs[Dot[eigS[[2,1]],eigN[[2,1]]]];
komp1={sa,sb,(sax[[1]]/sax[[2]]-1)(1-sp)};sim1=Exp[-(sa+sb+(sax[[1]]/sax[[2]]-1)(1-sp))];
If[sim1≥ sim,sim=sim1;komp=komp1;par={s,j}];
,{j,Length[ELLS0]}];
(*Print[par,": komp=",komp,"; sim=",sim];*)
If[sim>.9,Nr=Nr+1;
Print[par[[1]],". calculated: ",EllParam[EllN[[par[[1]]]]]," coincides with\n", par[[2]],". Original    ",EllParam[ELLS0[[par[[2]]]]]]]
,{s,Length[EllN]}];
Nr
]

NEW Method II for Ellipses

(* 
itKM - maximal number of iteration WKM;
itDIR - maximal number of iteration DIRECT;
Ind ≠ 0 - show results DIRECT and WKM;
Sl ≠ 0 - show figures (original, DIRECT, WKC, Comparison);
Primjedba-1: U DIRECT se ne računaju radijusi, već se stavlja RR=1;
Primjedba-2: U WKMC modulu koristim FindMinimum uz dobru početnu aproksimaciju, a u NewCirc modulu koristim DIRECT;
Primjedba-3: u popis internih varijabli ne smije doci ime modula koji se koristi;
*)

NewEllipse[k_,w_,A_,itKM_:50,itDIR_:10,Ind_:0,Sl_:0]:=Module[{m=Length[A],FCompile,f,min,B,kCl=k,z,zz,centri,rad,cenDIR,sol,time,vdb,rand},
aaa=Floor[Min[A[[All,1]]]]; bbb=Ceiling[Max[A[[All,1]]]];
ccc=Floor[Min[A[[All,2]]]]; ddd=Ceiling[Max[A[[All,2]]]];
rrr=1;RR=1;
tab0={bbb-aaa,ddd-ccc};
tab1={aaa,ccc}; 
II={0,1}; 
dom=Table[II,{j,2k}];	
AA1=DiagonalMatrix[Flatten[Table[tab0,{j,k}],1]];
AA=DiagonalMatrix[Flatten[Table[1/tab0,{j,k}],1]];
fu=Flatten[Table[tab1,{j,k}],1];
TF[x_]:=AA.(x-fu);
TF1[x_]:=AA1.x+fu;
If[k==1,
FCompile=Compile[{{z1,_Real},{z2,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR]}],{i,m}]]];	
f[z1_,z2_]:=Apply[FCompile,TF1[{z1,z2}]],
If[k==2,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_]:=Apply[FCompile,TF1[{z1,z2,z4,z5}]],
If[k==3,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real},{z7,_Real},{z8,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR],Dist[A[[i]],{z7,z8},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_,z7_,z8_]:=Apply[FCompile,TF1[{z1,z2,z4,z5,z7,z8}]],
If[k==4,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real},{z7,_Real},{z8,_Real},{z10,_Real},{z11,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR],Dist[A[[i]],{z7,z8},RR],Dist[A[[i]],{z10,z11},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_,z7_,z8_,z10_,z11_]:=Apply[FCompile,TF1[{z1,z2,z4,z5,z7,z8,z10,z11}]],
If[k==5,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real},{z7,_Real},{z8,_Real},{z10,_Real},{z11,_Real},{z13,_Real},{z14,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR],Dist[A[[i]],{z7,z8},RR],Dist[A[[i]],{z10,z11},RR],Dist[A[[i]],{z13,z14},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_,z7_,z8_,z10_,z11_,z13_,z14_]:=Apply[FCompile,TF1[{z1,z2,z4,z5,z7,z8,z10,z11,z13,z14}]],
If[k==6,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real},{z7,_Real},{z8,_Real},{z10,_Real},{z11,_Real},{z13,_Real},{z14,_Real},{z16,_Real},{z17,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR],Dist[A[[i]],{z7,z8},RR],Dist[A[[i]],{z10,z11},RR],Dist[A[[i]],{z13,z14},RR],Dist[A[[i]],{z16,z17},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_,z7_,z8_,z10_,z11_,z13_,z14_,z16_,z17_]:=Apply[FCompile,TF1[{z1,z2,z4,z5,z7,z8,z10,z11,z13,z14,z16,z17}]],
If[k==7,
FCompile=Compile[{{z1,_Real},{z2,_Real},{z4,_Real},{z5,_Real},{z7,_Real},{z8,_Real},{z10,_Real},{z11,_Real},{z13,_Real},{z14,_Real},{z16,_Real},{z17,_Real},{z19,_Real},{z20,_Real}},
Evaluate[Sum[Min[{Dist[A[[i]],{z1,z2},RR],Dist[A[[i]],{z4,z5},RR],Dist[A[[i]],{z7,z8},RR],Dist[A[[i]],{z10,z11},RR],Dist[A[[i]],{z13,z14},RR],Dist[A[[i]],{z16,z17},RR],Dist[A[[i]],{z19,z20},RR]}],{i,m}]]];	
f[z1_,z2_,z4_,z5_,z7_,z8_,z10_,z11_,z13_,z14_,z16_,z17_,z19_,z20_]:=Apply[FCompile,TF1[{z1,z2,z4,z5,z7,z8,z10,z11,z13,z14,z16,z17,z19,z20}]]  ]]]]]]];
						(* DIRECT *)
tDIR=Timing[min=DIRECT[f,dom,DMin->10.^(-4),MaxIter->itDIR, Eps->10^(-4)]];
If[Ind≠0,Print["DIRECT: F=",min[[2]],"; Nr.ITER=",len=Length[ITER],"; Del_F2=",Abs[ITER[[len,2]]-ITER[[len-1,2]]]]];
					(* Defining of cenDIR *)
zz=TF1[min[[1]]];
centri=rad=Table[0,{j,k}];
Do[
kk=2j-1;
centri[[j]]={zz[[kk]],zz[[kk+1]]};
,{j,k}];
cenDIR=Table[{centri[[j]],RR},{j,k}]; 
							(* WKC *)
tWKC=Timing[solC=WKC[w,A,cenDIR,itKM];];
cenWKC=Table[{solC[[2,j,1]],solC[[2,j,2]],IdentityMatrix[2]},{j,k}];
							(* WKMC *)
tWKMC=Timing[sol=WKMC[w,A,cenWKC,itKM];];
time=tDIR[[1]]+tWKC[[1]]+tWKMC[[1]];
(*B=Table[sol[[1,j,All,2]],{j,Length[sol[[2]]]}]; 
If[Length[B]<Length[pod],rand=Rand[B,pod],rand=Rand[pod,B]]; *)
If[Ind≠0,Print["WKMC: F=",sol[[3]], "; Nr.It(WKC,WKMC)=",{solC[[4]],sol[[4]]},"; TotalTime=",{tDIR[[1]],tWKC[[1]],tWKMC[[1]]},"=",time(*,"\n{Rand, Jaccard} =",rand[[2]]*)]];
If[Length[sol[[1]]]<kCl,Print["Došlo je do smanjenja broja klastera: od ",kCl," na ", Length[sol[[1]]]]];
cen=sol[[2]]; 
If[Sl≠0, 
AG={{aaa,bbb},{ccc,ddd}};
gx=Table[i,{i,AG[[1,1]],AG[[1,2]]}]; 
gy=Table[i,{i,AG[[2,1]],AG[[2,2]]}]; 
(* SLIKA DIRECT kružnice*)
slc=ListPlot[cenDIR[[All,1]],PlotStyle->{Black,PointSize[.04],Opacity[.5]}];
kruznice=Show[Table[Graphics[Circle[cenDIR[[j,1]],cenDIR[[j,2]]]],{j,Length[cenDIR]}]];
sl1=Show[slpod,slc,kruznice,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"DIRECT"];
(*  For the paper *)
Print[Show[slpod,slc,kruznice,AspectRatio->Automatic,Ticks->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]]; 
(* FIGURE of WKC kružnice*)
slc=ListPlot[cenWKC[[All,1]],PlotStyle->{PointSize[.05],Opacity[.5]}];
kruznice=Show[Table[Graphics[Circle[cenWKC[[j,1]],cenWKC[[j,2]]]],{j,Length[cenWKC]}]];
slC=Show[kruznice,slc,slpod,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"KCC"];
(* FIGURE of WKMC kružnice*)
slc=ListPlot[cen[[All,1]],PlotStyle->{PointSize[.05],Opacity[.5]}];
(*kruznice=Show[Table[Graphics[{Red,Circle[cen[[j,1]],cen[[j,2]]]}],{j,Length[cen]}]];*)
kruznice=Table[ContourPlot[dM[sol[[2,j,1]],{x,y},sol[[2,j,3]]]==sol[[2,j,2]]^2,
{x,AG[[1,1]],AG[[1,2]]},{y,AG[[2,1]],AG[[2,2]]},ContourStyle->{Red}],{j,Length[cen]}];
(* ALl figures *)
sl2=Show[kruznice,slc,slpod,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"KMCC"];
(*    ZA TEXT RADA *)
Print[Show[slpod,kruznice,slc,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},AxesLabel->{x,y},ImageSize->Small]];
pts=Show[slE,kruznice,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"Komparacija elipsi"];
Print[GraphicsGrid[{{(*sl0,*)sl1,slC,sl2,pts}},ImageSize->Large]]
];
ell=Table[EllParam[sol[[2,j]]],{j,Length[sol[[2]]]}];
vdb=VDB[sol[[1]],sol[[2]]];  vgo=GOindex[sol[[1]],sol[[2]]];  
{sol[[1]],sol[[2]],sol[[3]],sol[[4]],time,vdb,vgo,(*rand[[2]]*),{itDIR,solC[[4]],sol[[4]]},{tDIR[[1]],tWKC[[1]],tWKMC[[1]]}}
]

Davies - Bouldin Index with Hausforff distance betwwen two ellipses

(* Slučajno se izabere nh=50 točaka po obje elipse pa se računa Haudorffova udaljenost tih skupova točaka *)

Clear[VDB]
VDB[PI_,z_]:=Module[{k=Length[z],FF=Table[0,{j,k}],VV=Table[0,{j,k}],mj,sum=0},
nh=50;tran=RandomReal[{0,2Pi},nh];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,s,1]]*Sqrt[DistM[PI[[j,s,2]],z[[j]]]],{s,mj}];
VV[[j]]=FF[[j]]/Sum[PI[[j,s,1]],{s,mj}];
   ,{j,k}];
Do[
max={};
Do[
(*If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/Hdist[z[[i]],z[[j]]]]]*)
If[s≠j,
ells=TMCIR[z[[s]]];ellj=TMCIR[z[[j]]];
tocS=Table[ells[[1]]+Q[ells[[3]]].{ells[[2,1]]Cos[t],ells[[2,2]]Sin[t]},{t,tran}];
tocJ=Table[ellj[[1]]+Q[ellj[[3]]].{ellj[[2,1]]Cos[t],ellj[[2,2]]Sin[t]},{t,tran}];max=Append[max,(VV[[s]]+VV[[j]])/(Hd[tocS,tocJ]) ]];
,{s,k}]; 
sum=sum+Max[max]
,{j,k}];
sum/k
            ]

Davies - Bouldin Index with PRL distance betwwen two ellipses

(* Korištena je izvedena Hausdorffova udaljenost
 dH(C1,C2)=dM(C1,C2,Sigma1+Sigma2)+|r2-r1|^2 kod uporabe algebarske udaljenosti D,a u indeksu da je Sqrt[D] i dH^2 *)

Clear[VDB]
VDB[PI_,z_]:=Module[{k1,FF,VV,mj,sum},
k1=Length[z];FF=Table[0,{j,k1}];VV=Table[0,{j,k1}];sum=0;
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,s,1]]*DistM[PI[[j,s,2]],z[[j]]],{s,mj}];
VV[[j]]=FF[[j]]/Sum[PI[[j,s,1]],{s,mj}];
   ,{j,k1}];
Do[
max={};
Do[
(*If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/Hdist[z[[i]],z[[j]]]]]*)
If[s≠j,max=Append[max,(VV[[s]]+VV[[j]])/(dM[z[[s,1]],z[[j,1]],z[[s,3]]+z[[j,3]]]+ Abs[z[[s,2]]-z[[j,2]]]^2)]];
,{s,k1}]; 
sum=sum+Max[max]
,{j,k1}]; 
sum/k1
            ]

Geometrical - Index

(* An ellipse as M-circle *)

GOindex[PI_,Ellipses_]:=Module[{k1,h,rho={},msd={},rj,GOM},
k1=Length[Ellipses]; 
Do[
mj=Length[PI[[j]]];
ell=EllParam[Ellipses[[j]]];
(*dEl=(4NIntegrate[Sqrt[ell[[2,1]]^2Cos[t]^2+ell[[2,2]]^2Sin[t]^2],{t,0,Pi/2}]);*)
h=(ell[[2,1]]-ell[[2,2]])^2/(ell[[2,1]]+ell[[2,2]])^2;dEl=Pi*(ell[[2,1]]+ell[[2,2]])*(1+3h/(10+Sqrt[4-3h]));
ud=Sum[DistM[PI[[j,s,2]],Ellipses[[j]]],{s,mj}];
msd=Append[msd,dEl*ud/(mj^2)]; 
,{j,Length[Ellipses]}]; 
Print["M-circles: ",Ellipses];
Print["|pi|: ",LPI=Table[Length[PI[[j]]],{j,k1}]];
Print[k1,": MSD=",msd]; 
GOM=k1 Total[msd]
];

DIRECT Modul

Mathematica codes written by I.Vazler

Options[DIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
DIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,ind,novic,novid,novifc,red,TMin,iter=0,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="Maximal number of iterations is reached `1`.";
n=Length[dom];
(* Scaling the function f on the hypercube with side length 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* Initial center and diameter *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};
TMin={c[[1]],fc[[1]]};ITER={};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: searching for the potentially optimal hyperrectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with smallest f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other conditions *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* STEP 2: dividing potentially optimal hyperrectangles *)
Do[
(* Dividing directions *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
If[novifc[[2red[[1]]-1]]<TMin[[2]],TMin=N[{novic[[2red[[1]]-1]],novifc[[2red[[1]]-1]]},10]];
If[novifc[[2red[[1]]]]<TMin[[2]],TMin=N[{novic[[2red[[1]]]],novifc[[2red[[1]]]]},10]];
,{pot,pot2}];
ITER=AppendTo[ITER,{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];
(*Print[iter,": ",ITER];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure in 2D case *)
(*If[n==2,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}]
(* Back in initial domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Images printing in 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Loading of one data set

SetDirectory[NotebookDirectory[]]; 
DATAELL=ReadList["DATAELL-1.txt"];
Q[th_]:={{Cos[th],-Sin[th]},{Sin[th],Cos[th]}}
ii=49;  (* DATAELL-1: 39, 3, 49, 17  separated *)
		(* DATAELL-2: 39, 4, 27, 1  intersecting *)
Print["\nFrom the set DATAELL of ",Length[DATAELL]," examples, ",Style["example:",16,Blue],Style[ii,18,Blue]," is choosen"];
k=Length[DATAELL[[ii,2]]];
Print["Centers, Semiaxes, Theta:"]
Ell=DATAELL[[ii,2]]
A=DATAELL[[ii,1]]; 
m=Length[A]; w=Table[1,{i,m}];
(* Borders + Figures *)
aaa=Floor[Min[A[[All,1]]]]; bbb=Ceiling[Max[A[[All,1]]]];
ccc=Floor[Min[A[[All,2]]]]; ddd=Ceiling[Max[A[[All,2]]]];
AG={{aaa,bbb},{ccc,ddd}};
gx=Table[i,{i,AG[[1,1]],AG[[1,2]]}]; 
gy=Table[i,{i,AG[[2,1]],AG[[2,2]]}];
slE=Table[ParametricPlot[Ell[[j,1]]+Q[Ell[[j,3]]].{Ell[[j,2,1]]*Cos[tt],Ell[[j,2,2]]*Sin[tt]},{tt,0,2Pi},PlotStyle->{Blue}],{j,Length[Ell]}];
slcenters=ListPlot[Table[Ell[[j,1]],{j,Length[Ell]}],PlotStyle->{Black,PointSize[.03],Opacity[.5]}];
slpod=ListPlot[A,PlotStyle->{Blue,PointSize[0.01]},AxesLabel->{x,y}];
sl0=Show[slpod,slE,slcenters,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},(* PlotLabel->"ORIGINAL", *)ImageSize->Small]
(* For the paper 
Show[slE,slcenters,slpod,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->Small]*)

Searching for an optimal k - partition

p=2; k=4;
Print["Traženje optimalne ",k," - particije:"];
new0=NewEllipse[k,w,A,50,30,1,1];
Print["vdb(",k,") = ",new0[[6]],"; vgo(",k,") = ",new0[[7]]]
sim=Scn2[PODELL[[ii,2]],new0[[2]],1];
Print["Nr. of detected ellipses: ",sim]

Reconstruction of one data set with DB-index

(* 
1. DB-index nema svojstvo monotonosti pa zato moramo racunati sve particije(Primjer ii=3 pokazuje da DB index ukazuje na k=3 jer je DB(4)>DB(3) iako je rjesenje k=5; Slično, Primjer ii=49 pokazuje da DB index ukazuje na k=2 jer je DB(3)>DB(2) iako je rjesenje k=4 
2. Kod WKMC ispisuje se broj iteracija za WKC i WKMC i vrijeme za DIRECT, WKC, WKMC;
*)

(* Racunam sve particije pa treba staviti nk=broj originalnih elipsi + 1 *)
p=2; nk=5;
new={};s=1;
Do[
Print["Searchin for a optimal ",s," - partition:"];
k=s;new=Append[new,NewEllipse[s,w,A,50,15,1,1]];
,{s,nk}];
new[[1,6]]=new[[1,3]]/m;
Print["DB=",vdb=Table[new[[s,6]],{s,nk}]];
kk=Ordering[vdb,1][[1]];
new0=new[[kk]];
(* Results *)
(*If[vdb1<vdb0,new0=new1;kk=k,kk=k-1];*)
EllN=Table[0,{j,kk}];
Print["Partition with the most appropriate number of ellipses has ",kk," ellipses:\nCenters, Radia, Covariance Matrices:",new0[[2]]]
Print["Centers, Semiaxes, Theta: "]
Do[
(* Eigensystem[ScovN]*)
wwN=new0[[2,j,3]]/Sqrt[Det[new0[[2,j,3]]]];
ewwN=Eigensystem[wwN];
EllN[[j]]={new0[[2,j,1]],new0[[2,j,2]]*Sqrt[ewwN[[1]]],Theta[new0[[2,j,3]]]}
,{j,kk}]
EllN
						(* FIGURES *)
slEN=Table[ParametricPlot[EllN[[j,1]]+Q[EllN[[j,3]]].{EllN[[j,2,1]]*Cos[tt],EllN[[j,2,2]]*Sin[tt]},{tt,0,2Pi},PlotStyle->{Black}],{j,Length[EllN]}];
slcenters=ListPlot[Table[EllN[[j,1]],{j,Length[EllN]}],PlotStyle->{Black,PointSize[.03],Opacity[.5]}];
slRez=Show[slEN,slcenters,slpod,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"Rekonstrukcija",ImageSize->Small]
                                    (* Rezultati detekcije i komparacija *)
Print["\nComparison originala and calculated ellipses:"];
sim=Scn2[PODELL[[ii,2]],new0[[2]],0];
Print["Nr. of detected ellipses: ",sim]

Reconstruction of one data set with GO-index

(* 
1. GO-index ima gotovo uvijek svojstvo monotonosti pa ne moramo racunati sve particije. Zaustavimo se kada GO-index počinje rasti  
2. Kod WKMC ispisuje se broj iteracija za WKC i WKMC i vrijeme za DIRECT, WKC, WKMC;
*)

p=2; k=1;
Print[ii,": Searchin for a optimal ",k," - partition:"];
new0=NewEllipse[k,w,A,50,15,1,1];
vgo0=new0[[7]]; goAll={vgo0};
Print["vgo(",k,") = ",vgo0,"\n"]
Fwkmc={new0[[3]]};
While[
k=k+1; 
new1=NewEllipse[k,w,A,50,15,1,1];
Fwkmc=Append[Fwkmc,new1[[3]]];
vgo1=new1[[7]];goAll=Append[goAll,vgo1];
Print["vgo(",k,") = ",vgo1,"\n"];
vgo1<vgo0,vgo0=vgo1; new0=new1
];
					(* RESULTS *)
kk=k-1;
Print["GO-All=",goAll];
Print["Fwkmc=",Fwkmc];
EllN=Table[0,{j,kk}];
Print["Partition with the most appropriate number of ellipses has ",kk," ellipses:\nCenters, Radia, Covariance Matrices:",new0[[2]]]
Print["Centers, Semiaxes, Theta: "]
Do[
(* Eigensystem[ScovN]*)
wwN=new0[[2,j,3]]/Sqrt[Det[new0[[2,j,3]]]];
ewwN=Eigensystem[wwN];
thN=Theta[new0[[2,j,3]]];
EllN[[j]]={new0[[2,j,1]],new0[[2,j,2]]*Sqrt[ewwN[[1]]],thN}
,{j,kk}]
EllN
						(* FIGURES *)
slEN=Table[ParametricPlot[EllN[[j,1]]+Q[EllN[[j,3]]].{EllN[[j,2,1]]*Cos[tt],EllN[[j,2,2]]*Sin[tt]},{tt,0,2Pi},PlotStyle->{Black}],{j,Length[EllN]}];
slcenters=ListPlot[Table[EllN[[j,1]],{j,Length[EllN]}],PlotStyle->{Black,PointSize[.03],Opacity[.5]}];
slRez=Show[slEN,slcenters,slpod,AspectRatio->Automatic,GridLines->{gx,gy},PlotRange->AG,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotLabel->"Rekonstrukcija",ImageSize->Small]
                                    (* Results andcomparison *)
Print["\nComparison originala and calculated ellipses:"];
sim=Scn2[PODELL[[ii,2]],new0[[2]],0];
Print["Nr. of detected ellipses: ",sim]