Color image segmentation - Mahalanobis Hard Clustering
All indexes are calculated on the basis of p-clusters

Moduli

Needs["ComputationalGeometry`"]
Needs["MultivariateStatistics`"]

Squares Mahalanobis distance

dm[x_,y_,S_]:=(x-y).Inverse[S].(x-y);
dM[x_,y_,S_]:=(Det[S]^(1/n))(x-y).Inverse[S].(x-y);

Modul IntersectTest

IntersectTest[c1_,cov1_,c2_,cov2_]:=Module[{FF,max},
FF[lambda_]:=lambda (1-lambda) (c2-c1).Inverse[lambda cov2+(1-lambda) cov1].(c2-c1);
max=First[FindMaximum[{FF[x],0<x<1},x]];
max>1
]

Discrete Hausdorff Distance

ClearAll[DiscreteHausdorffDistance];
Options[DiscreteHausdorffDistance] = Options[Nearest];
DiscreteHausdorffDistance[s1_, s2_, opts : OptionsPattern[]] := Max[
	Nearest[s1 -> "Distance", s2, 1, opts],
	Nearest[s2 -> "Distance", s1, 1, opts]
];

Position of the Minimal and the Maximal element of the upper triangular matrix

MinEl[R_]:=Module[{n=Length[R[[1]]],min=Infinity,red,ord,imin},
Do[
red=Drop[R[[i]],{1,i}];
ord=Ordering[red,1][[1]];
If[red[[ord]]<min,min=red[[ord]];imin={i,ord+i}];
(*Print[i,red,ord,":  ",imin];*)
,{i,n-1}];
{min,imin}]

MaxEl[R_]:=Module[{n=Length[R[[1]]],max=-Infinity,red,ord,imax},
Do[
red=Drop[R[[i]],{1,i}];
ord=Ordering[red,-1][[1]];
If[red[[ord]]>max,max=red[[ord]];imax={i,ord+i}];
(* Print[i,red,ord,":  ",imax]; *)
,{i,n-1}];
{max,imax}]

Minimal and the Maximal element of the upper triangular matrix

MinR[R_]:=Module[{n=Length[R[[1]]],min=Infinity},
Do[
red=Drop[R[[i]],{1,i}];
minn=Min[red];
If[minn<min,min=minn];
,{i,n-1}];
min]

MaxR[R_]:=Module[{n=Length[R[[1]]],max=-Infinity},
Do[
red=Drop[R[[i]],{1,i}];
maxn=Max[red];
If[maxn>max,max=maxn];
,{i,n-1}];
max]

Important !!!
Image composed of selected clusters from WKMeans

(* JJ: indeksi izabranih klastera *)
SelectedColors[sol_,JJ_]:=Module[{AA,podDop0,sl0,pp0,slikaN0,zC},
AA=Table[Table[{{1,1,1},sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[1]]]}];
Do[
AA[[jj]]=Table[{sol[[1,jj,s,2]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}];
,{jj,JJ}];
podDop0=Flatten[AA,1];
sl0=SortBy[podDop0,Last];
pp0=sl0[[All,1]];
slikaN0=Table[Table[pp0[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print["All centers: "];
Print[Table[Graphics[{RGBColor[sol[[2,j]]],Disk[]},ImageSize->20],{j,Length[sol[[2]]]}]];
Print["Choosen centers: "];
zC=Table[sol[[2,j]],{j,JJ}];
Print[Table[Graphics[{RGBColor[zC[[i]]],Disk[]},ImageSize->20],{i,Length[zC]}]];
Print[Image[slikaN0,ImageSize->Small]]
]

Objective functions

MWF[PI_,c_,Scov_]:=Module[{m=Length[Flatten[PI,1]],n=Length[Scov]},
Return[Sum[PI[[j,s,1]]*dM[PI[[j,s,2]],c[[j]],Scov[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

WF[PI_]:=Module[{m=Length[Flatten[PI,1]],c},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Return[Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-c[[j]]]^2,{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

Mahalanobis Weighted k - means  for the data with additional information

(* Koristi se dM *)
MWKMeansAdd[w_,B_,z_,Ind_:0]:=Module[{m=Length[w],k=Length[z],W=Total[w],PI0,PIpriv,tab,min,imin,Pod,c,c0,n=Length[z[[1]]],Scov,pom,corr,It=0,IT=20,epsilon=.000000001,eps=.008},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,m}]; 
AG=Table[{0,1},{i,3}];
c0=z;
						(* Računanje particije *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];pom=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
                                                           (* Računanje centara *)
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
F0=WF[PI];
If[Ind≠ 0,Print[It,": F0 = ",F0]];
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition"]];
							(* LOOP *)
While[
c0=c;It=It+1;
				(* Racunanje matrica Scov *)
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,All,2]][[s]]-c[[j]]),(PI[[j,All,2]][[s]]-c[[j]])],{s,Length[PI[[j]]]}]/mj[[j]]
,{j,k}];
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-tog centra (",c[[j]],"): ",Scov[[j]]," je singularna"];Print["|PI[j]|: ",Length[PI[[j]]],", Eigen_min=",Min[Eigenvalues[Scov[[j]]]]];Scov[[j]]=IdentityMatrix[n] ]
,{j,k}];
PIpriv=PI;
PI=Table[{},{i,k}];
corr=Table[Log[m^2 Det[Scov[[jj]]]/mj[[jj]]^2],{jj,k}];
Do[
Do[
tab=Table[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
PI=DeleteCases[PI,{}];
k=Length[PI];
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
F1= MWF[PI,c,Scov];delF=Abs[F0-F1]/F0;absc=Norm[c-c0,Infinity];
If[Ind≠ 0,Print[It,": F1 = ",F1,"; |c-c0| = ",absc,"; delF = ",delF]];
absc>eps&&It<IT, F0=F1; 
]; 
				(* Slucaj smanjenja broja klastera u particiji *)
If[Length[c]<Length[c0],
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[B[[i]]-c[[s]]]^2,{s,Length[c]}],Norm[B[[i]]-{u,v}]^2]],{i,m}];
dom=AG;
								(* Novi DIRECT centar*)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];Print["DIRECT: ",min];
c=Append[c,min[[1]]];
Print["Novi optimirani centar: ",c];
k=Length[c]; Print["Korigirani k: ",k];
								(* Nova particija *)
PI=Table[{},{j,k}];
Do[
tab=Table[Norm[Pod[[i]][[2]]-c[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It,F1}
]

Mahalanobis Incremental algorithm   for the data with additional information

(* Ako je poznata originalna particija PI0, možemo je usporediti s particijom dobivenom u brk-koraku *)
MIncAdd[w_,B_,cen_,KK_,Ind_:0]:=Module[{centri,centriN,p,dom,Scov,min,PI1},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,m}]; 
AG=Table[{0,1},{i,3}];
(* Initialization *)
centri=cen;
IndMDB=IndMHV=IndMVolume=IndMMD={};
Scov=Table[IdentityMatrix[n],{i,Length[cen]}]; 
sol=MWKMeansAdd[w,B,cen]; 
centri=sol[[2]];
brk0=Length[centri];
					(* Iterativni postupak *)
Do[
If[Ind≠0,Print[Style["\nIteracija ",Red,Bold,12],Style[brk,Red,Bold,14]," Na postojecih ",Length[centri]," centara dodaje se 1 novi"]];
k=Length[centri];
						(* Compile *)
Clear[u1,u2,u3,F2];
FCompile=Compile[{{u1,_Real},{u2,_Real},{u3,_Real}},
Evaluate[Sum[w[[i]]Min[Append[Table[Norm[Pod[[i,2]]-centri[[s]]]^2,{s,Length[centri]}],Norm[Pod[[i,2]]-{u1,u2,u3}]^2]],{i,m}] ]];	
F2[u1_,u2_,u3_]:=Apply[FCompile,{u1,u2,u3}];	
dom=AG;
					(* DIRECT *)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)]//Timing; min=minp[[2]];If[Ind≠0,Print["DIRECT: ",min]];
centriN=Append[centri,min[[1]]];
If[Ind≠0,Print["Novi optimirani centri: ",centriN]];
(* Ovdje se vrijednost funkcije cilja ne moze racunati kao MWFmin jer jos nemam Scov[[k+1]] *)
k=Length[centriN];
					(* Poziv Mahalanobis k-means *)
t1=Timing[sol=MWKMeansAdd[w,B,centriN]];
tab=Table[{sol[[1,j]],sol[[2,j]],sol[[3,j]],Norm[sol[[2,j]]],j},{j,Length[sol[[2]]]}];
sort=Sort[tab,#1[[4]]<#2[[4]]&];
sol[[1]]=sort[[All,1]];sol[[2]]=sort[[All,2]];sol[[3]]=sort[[All,3]];
If[Ind≠0,Print["Novi MWKMeans profinjeni centri: ",sol[[2]]]];
Scov=sol[[3]];
If[Ind≠0,Print["MWKMeans iteracija = ",sol[[4]],"; MWF = ",FF=MWF[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",mjx=Table[Length[sol[[1,j]]],{j,Length[PI]}]]];
If[Ind≠0,Print["Time (DIRECT+WKM) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]]];
Print[brk+1,"-particija:\n",{brk+1,Niter,Round[min[[2]],-.001],Round[minp[[1]],-.01],sol[[4]],Round[FF,-.01],Round[t1[[1]],-.01],mjx}];
				(* p-clusters, p-centers, p-covariance matrices *)
PI0=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];  (* Clusters bez koordinata pixela *)
ell=Table[Ellipsoid[sol[[2,j]], kappa sol[[3,j]]],{j,Length[sol[[2]]]}];
pPI0=pclusters[PI0,sol[[2]],sol[[3]],kappa];
pcen=Table[Mean[pPI0[[j]]],{j,Length[sol[[1]]]}];
pScov=Table[Covariance[pPI0[[j]]],{j,Length[sol[[1]]]}];				
						(* p-INDEXES *)
pvdb=pMDB[pPI0,pcen,pScov];IndMDB=Append[IndMDB,pvdb];
phv=pMHV[pPI0,pcen,pScov];IndMHV=Append[IndMHV,phv];
pvolume=pMVolume[pPI0,pcen,pScov];IndMVolume=Append[IndMVolume,pvolume];
pmmd=pMMD[pPI0,pcen,pScov,Ind2];
IndMMD=Append[IndMMD,pmmd];
If[Ind≠0,Print["\nMDB=",Round[pvdb,-.0001],"\nMHV=",Round[phv,-.000001],"\nMVolume=",pvolume"\nMMD=",Round[pmmd,-.001]]];
			(* FIGURE: Ellipsoids + Nova segmentirana slika *)
pell=Table[Ellipsoid[pcen[[j]], kappa pScov[[j]]],{j,Length[pcen]}];
slpell=Show[Table[Graphics3D[{RGBColor[pcen[[j]]],Opacity[0.4],ell[[j]]}],{j,Length[sol[[2]]]}],Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"}];
slptoc=Table[ListPointPlot3D[PI0[[j]],PlotStyle->{RGBColor[pcen[[j]]],PointSize[.004],Opacity[.5]}],{j,Length[pcen]}];
Print[Show[slpell,slptoc, ImageSize->Small]];
time=Timing[PI=Flatten[Table[Table[{sol[[2,j]],sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[2]]]}],1];
sl=SortBy[PI,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
];
Print["Segmentirana slika: br.clust.=",Length[sol[[2]]],  ";  CPU-time = ", time[[1]],"\nOptimalni centri: ",sol[[2]],"; F = ",MWF[sol[[1]],sol[[2]],sol[[3]]]];
Print[slN=Image[slikaN,ImageSize->Small]];
Print[Table[Graphics[{RGBColor[sol[[2,j]]],Disk[]},ImageSize->20],{j,Length[sol[[2]]]}]];
centri=sol[[2]];     
,{brk,KK}];
                                               	(* Graphs of Indexes *)
Print["\nIndMDB: ",IndMDB,"; MAPart=",Ordering[IndMDB,1][[1]]+1,"\nIndMMD: ",IndMMD,"; MAPart=",Ordering[IndMMD,-1][[1]]+1,"\nIndMHV: ",IndMHV,"; MAPart=",Ordering[IndMHV,1][[1]]+1,"\nIndMVolume: ",IndMVolume,"; MAPart=",Ordering[IndMVolume,1][[1]]+1];
Print[GraphicsGrid[{{ListPlot[Table[{i+1,IndMDB[[i]]},{i,Length[IndMDB]}],Joined->True,GridLines->Automatic,PlotLabel->"MDB",AxesOrigin->{1,0}],
ListPlot[Table[{i+1,IndMMD[[i]]},{i,Length[IndMMD]}],Joined->True,GridLines->Automatic,PlotLabel->"MMD",AxesOrigin->{1,0}]},
{ListPlot[Table[{i+1,IndMHV[[i]]},{i,Length[IndMHV]}],Joined->True,GridLines->Automatic,PlotLabel->"MHV",AxesOrigin->{1,0}],
ListPlot[Table[{i+1,IndMVolume[[i]]},{i,Length[IndMVolume]}],Joined->True,GridLines->Automatic,PlotLabel->"MVolume",AxesOrigin->{1,0}]}},ImageSize->Medium]];
{PI,sol[[2]],sol[[3]]}
]

p - clusters

pclusters[partition_,centers_,covariances_,kappa_]:=Module[{invCovs=Inverse/@covariances,k=Length[centers]},
Table[
Module[{clusterData=partition[[j]],diff,dists},
If[clusterData==={},{},
diff=clusterData-ConstantArray[centers[[j]],Length[clusterData]];
dists=Total[(diff.invCovs[[j]])*diff,{2}];
Pick[clusterData,UnitStep[kappa-dists],1]]
],{j,k}
]
];

Indexes - Mahalanobis

Clear[MDB]  
pMDB[PI_,z_,Scov_]:=Module[{k=Length[z],FF,VV,mj,tab,sum=0},
VV=Table[0,{j,k}];
Do[
VV[[j]]=Sum[dM[z[[j]],PI[[j,s]],Scov[[j]]],{s,Length[PI[[j]]]}]/Length[PI[[j]]]
,{j,k}];
Do[
max={};
Do[
If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/dM[z[[i]],z[[j]],Scov[[i]]+Scov[[j]]] ]]
,{i,k}];
sum=sum+Max[max]
,{j,k}];
sum/k
            ]

(*  Hypervolume index - čista determinanta (ne dijeli se s brojem točaka) *)
Clear[MHV]
pMHV[PI_,z_,Scov_]:=Module[{},
 Sum[Sqrt[Det[Scov[[j]]]],{j,Length[PI]}]
]

(*  MVolume index - determinanta se dijeli s brojem točaka u klasteru *)
Clear[pMVolume]
pMVolume[PI_,z_,Scov_]:=Module[{},
 Sum[Sqrt[Det[Scov[[j]]]/Length[PI[[j]]]],{j,Length[PI]}]
]

Clear[pMMD]
pMMD[PI_,z_,Scov_,Ind2_:0]:=Module[{k=Length[PI],min,dist,P1,P2},
(* x, y ne smiju biti u popisu internih varijabli *)
DDD=Table[0,{r,k},{s,k}];
				(* Distances between p-clusters *)
Do[
Do[
If[Ind2==0,DDD[[r,s]]=Min[Flatten[Table[EuclideanDistance[a,b],{a,PI[[r]]},{b,PI[[s]]}],1]],
dist=DiscreteHausdorffDistance[pPI0[[r]],pPI0[[s]]];
DDD[[r,s]]=dist[[1]];
{P1,P2}={pPI0[[r,dist[[2,1]]]], pPI0[[s,dist[[2,2]]]]};
If[n==2,
Print[Show[
Graphics[{boja[[r]],EdgeForm[Directive[Thickness[.004],boja[[r]]]],Opacity[0.075],ell[[r]]}],
Graphics[{boja[[s]],EdgeForm[Directive[Thickness[.004],boja[[s]]]],Opacity[0.075],ell[[s]]}],
ListPlot[PI0,PlotStyle->Table[{PointSize[.008],boja[[j]]},{j,Length[sol[[2]]]}]], 
ListPlot[{P1,P2}, PlotStyle->{Red,PointSize[.025]}], 
Graphics[{Red,Thin,Line[{P1,P2}]}],
Axes->True, AxesOrigin->{0,0},AxesStyle->{Thickness[.005],Thickness[.005]},GridLines->{tabx,taby},PlotRange->AG,AspectRatio->Automatic,ImageSize->150]],
Print[Show[
Graphics3D[{RGBColor[z[[r]]],Opacity[0.5],ell[[r]]}],Graphics3D[{RGBColor[z[[s]]],Opacity[0.5],ell[[s]]}],
ListPointPlot3D[PI0,PlotStyle->Table[{PointSize[.001],Opacity[.5],RGBColor[z[[j]]]},{j,Length[sol[[2]]]}]],
ListPointPlot3D[{P1,P2}, PlotStyle->{Red,PointSize[Large]}],
Graphics3D[{Red,Thin,Line[{P1,P2}]}],
Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},ImageSize->150]]
		];
			];
,{s,r+1,k}]
,{r,k}];
Print["DDD=",MatrixForm[DDD]," - distance matrix between clusters"];
pm=Sum[Length[pPI0[[j]]],{j,k}];
(pm-Length[PI])MinR[DDD]/MWF[pPI0,pcen,pScov]
]

Broj piksela izabrane boje

NrC[slcolor_]:=Module[{slika,data,nrcolor},
slika=ImageData[slcolor];nrcolor=0;
data=Flatten[slika,1];
Do[
If[data[[i]]≠{1.,1.,1.},nrcolor=nrcolor+1]
,{i,Length[data]}];
nrcolor
]

DIRECT Modul

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="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
(* skaliranje funkcije na hiperkocku duljine stranice 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* početni centar i promjeri *)
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]]};
(* GLAVNA PETLJA *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* KORAK 1: traženje potencijalno optimalnih hiperkvadara *)
dmaxrazni=Union[dmax];
(* svi podaci grupirani po d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* podaci s najmanjim f grupirani po d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* provjera ostalih uvjeta *)
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];
(* KORAK 2: dijeljenje potencijalno optimalnih *)
Do[
(* smjerovi po kojima se dijeli *)
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}];
Niter=iter;
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* slika za 2D slučaj *)
(*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]}]
(* vraćanje u početnu domenu *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* ispis slika za 2D slučaj *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Hard Mahalanobis Clustering

Input slike (Insert -> Picture -> From File)

Fig=Image[]; 
dim=ImageDimensions[Fig];
Print["Originalna slika: ImageDimensions = ",dim," --- ",dim[[1]] dim[[2]] ]
If[Max[dim]>128,
resize=128; Fig=ImageResize[Fig,resize];
dim=ImageDimensions[Fig];
Print["Resize slika: ImageDimensions = ",dim," --- ",dim[[1]] dim[[2]] ]
]

Added information

(* CPU Time za MWKMeansDop i MIncAdd algorithm ovisi o brju piksela na izabranoj slici or *)
(* or=ImageResize[Fig,100]; ImageDimensions[or] *)
slika=ImageData[Fig];
Print["Originalna slika"]
slor=Image[slika,ImageSize->Small]
(* p and kappa *)
n=3;  p=.6;  kappa=Quantile[ChiSquareDistribution[n],p];
Print["n = ",n,";  p = ",p,";  kappa = ",kappa];
(* Podaci *)
pod=Join@@slika;   
pod=Flatten[slika,1];
Print["Dim slike: ",{aa,bb,cc}=Dimensions[slika],";  |pod| = ",Length[pod]]
Print["Prijedolg  centara:"]
SeedRandom[3]
boja=Table[pod[[RandomInteger[{1,aa bb}]]],{i,10}];
Table[Graphics[{RGBColor[boja[[i]]],Disk[]},ImageSize->20],{i,Length[boja]}]

Dodavanje dopunskih informacija i kontrola (crtanje slike)

slikaDop=Table[{slika[[i,j]],{i,j}},{i,aa},{j,bb}];
podDop=Flatten[slikaDop,1];
m=Length[podDop];
w=Table[1,{i,Length[podDop]}];
Print["Dim slike: ",{aa,bb,cc}=Dimensions[slikaDop],";  |pod| = ",Length[podDop]]
(*
pp=podDop[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Image[slikaN,ImageSize->Small]
*)
Print["Prijedolg  centara:"]
SeedRandom[3]
boja=Table[pod[[RandomInteger[{1,aa bb}]]],{i,10}]
Table[Graphics[{RGBColor[boja[[i]]],Disk[]},ImageSize->20],{i,Length[boja]}]

Segmentacija slike Mahalanobis WKMeans algorithm
Pokrenuti samo po potrebi

(* iz Prijedloga centara izabrati nekoliko boja i unijeti ih u z;
Pozicija ellipsoida određena je centrom sol[[2,j]]; Veličina i smjer elipsoida određena je sa kappa sol[[3,j]], gdje je kappa=3.5 *)
z={boja[[3]],boja[[2]]}; k=Length[z];
Print["Initial centers:"]
Table[Graphics[{RGBColor[z[[i]]],Disk[]},ImageSize->20],{i,Length[z]}]
time=Timing[
sol=MWKMeansAdd[w,podDop,z,0];
                                 (* Definiranje nove segmentirane slike *)
PI=Flatten[Table[Table[{sol[[2,j]],sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,k}],1];
sl=SortBy[PI,Last];
];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
slN=Image[slikaN,ImageSize->Small]
Print["Segmentirana slika: br.clust.=",k, "; MSE = ",ImageDistance[slor,slN,DistanceFunction->"MeanSquaredEuclideanDistance"]";  CPU-time = ", time[[1]],"\nOptimalni centri: ",sol[[2]],"; F = ",sol[[5]]]
						 (* Ellipsoids *)
kappa=3.5;
ell=Table[0,{j,Length[sol[[2]]]}];
slEll=Show[Table[Graphics3D[{RGBColor[sol[[2,j]]],Opacity[0.4],ell[[j]]=Ellipsoid[sol[[2,j]], kappa sol[[3,j]]]}],{j,Length[sol[[2]]]}],Axes->True];
sltoc=Table[ListPointPlot3D[sol[[1,j,All,2]],PlotStyle->{RGBColor[sol[[2,j]]],PointSize[.008],Opacity[.5]}],{j,Length[sol[[1]]]}];
Print["Ellipsoidal clusters:"];
Print[Show[slEll,sltoc,ImageSize->250]];
Print["Izabrane boje:"];
Print[Table[Graphics[{RGBColor[sol[[2,j]]],Disk[]},ImageSize->20],{j,Length[sol[[2]]]}]];
Do[
Print["Ellipse ",j,": Center = ",sol[[2,j]],"; M-rad = ",Surd[Det[sol[[3,j]]],4] "; |pi| = ",Length[sol[[1,j]]]];
nell=0;
Do[
If[Element[sol[[1,j,All,2]][[s]],ell[[j]]]==True,nell=nell+1]
,{s,Length[sol[[1,j]]]}];
Print["    Elemens {in, all}: ",br={nell,Length[sol[[1,j]]]},"; % in = ",100. br[[1]]/br[[2]]]
,{j,Length[sol[[2]]]}];

Mahalanobis Incremental algorithm

Print["Initial center(s):"];
z0={boja[[3]]}; k=Length[z0];
KK=2; Ind=1; Ind2=0;
Table[Graphics[{RGBColor[z0[[i]]],Disk[]},ImageSize->20],{i,Length[z0]}]
(* w - weights; podDop - podaci dopunjenis koordinatama; cen - center(s); KK - broj iteracija; Ind - ako je "0", ne ispisuju se rezultati međukoraka *)
inc=MIncAdd[w,podDop,z0,KK,1];

For the text

Print[TableForm[Table[{i+1,Round[IndMDB[[i]],-.001],Round[IndMHV[[i]],-.000001],Round[IndMVolume[[i]],-.00000001],Round[IndMMD[[i]],-.00001]},{i,Length[IndMMD]}],TableHeadings->{None,{"k", "MDB","MHV","MVolume","MMD"}}]]

ListPlot[Table[{i+1,IndMDB[[i]]},{i,Length[IndMDB]}],Joined->True,GridLines->Automatic,AxesOrigin->{1,0},ImageSize->Small]
ListPlot[Table[{i+1,IndMVolume[[i]]},{i,Length[IndMVolume]}],Joined->True,GridLines->Automatic,AxesOrigin->{1,.0000005},ImageSize->Small]
ListPlot[Table[{i+1,IndMMD[[i]]},{i,Length[IndMMD]}],Joined->True,GridLines->Automatic,AxesOrigin->{1,0},ImageSize->Small]

Redoslijed klastera

sol2={{0.13435403951523475`,0.0709367082288051`,0.0623107549830894`},{0.27800386813057576`,0.2055350056011913`,0.11595567092945559`},{0.4506574043707705`,0.31988213430329115`,0.15546843167069252`},{0.7062071870437586`,0.4642214868834843`,0.16073362532706606`},{0.6734824537291568`,0.5952259330523464`,0.2862458768554153`},{0.9227967299772649`,0.725561497326201`,0.28692359702501746`}};
tab=Table[{sol2[[j]],Norm[sol2[[j]]],j},{j,Length[sol2]}];
sort=Sort[tab,#1[[2]]<#2[[2]]&]
z=sort[[All,1]];
Table[Graphics[{RGBColor[z[[i]]],Disk[]},ImageSize->20],{i,Length[z]}]

MWKMeans algorithm for the chosen centers

(* Za MWKMeansDop biram 5-partition jer je to particija s maksimalnim brojem dobro razdvojenih klastera *)
Ind2=0;
z={{0.16419859436103243`,0.1020742746191925`,0.07587396650896751`},{0.44160679991464435`,0.30419897100315474`,0.14174431562226014`},{0.6817272509324401`,0.558500858751998`,0.24665715554613188`},{0.9005896647845358`,0.7097955828747264`,0.2884080242060367`},{0.9325934147243795`,0.9294857565667773`,0.9120976692563814`}}; k=Length[z];
Print["Centri:"]
Table[Graphics[{RGBColor[z[[i]]],Disk[]},ImageSize->20],{i,Length[z]}]
time=Timing[
sol=MWKMeansAdd[w,podDop,z];
(* Definiranje nove segmentirane slike *)
PI=Flatten[Table[Table[{sol[[2,j]],sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,k}],1];
sl=SortBy[PI,Last];
];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
slN=Image[slikaN,ImageSize->Small]
Print["Centri Novi:"]
Print[Table[Graphics[{RGBColor[sol[[2,i]]],Disk[]},ImageSize->20],{i,Length[sol[[2]]]}]]
Print["Segmentirana slika: br.clust.=",k,  "; MSE = ",ImageDistance[slor,slN,DistanceFunction->"MeanSquaredEuclideanDistance"],";  CPU-time = ", time[[1]],"\nOptimalni centri: ",sol[[2]],"; F = ",sol[[5]]]
Do[
Print["Ellipsoid ",j,": Center = ",sol[[2,j]],"; M-rad = ",Surd[Det[sol[[3,j]]],4] "; |pi| = ",Length[sol[[1,j]]],"; Poluosi: ",Eigenvalues[sol[[3,j]]],"; Var: ",Variance[Table[sol[[2,j]]-sol[[1,j,s,2]],{s,Length[sol[[1,j]]]}]]];
,{j,Length[sol[[2]]]}];

Analiza slike po bojama iz Mahalanobis WKMeans algorithma

podDopCl=Table[0,{j,Length[sol[[1]]]}];
Do[
AA=Table[Table[{{1,1,1},sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[1]]]}];
AA[[jj]]=Table[{sol[[1,jj,s,2]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}];
podDopCl[[jj]]=AA[[jj]];
(* Slika satavljena od pi_jj (ostalo je prazno) *)
podDopjj=Flatten[AA,1];
sljj=SortBy[podDopjj,Last];
ppjj=sljj[[All,1]];
slikaNjj=Table[Table[ppjj[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print["\nSlika klastera s centrom: ",Graphics[{RGBColor[sol[[2,jj]]],Disk[]},ImageSize->20]];
Print[Image[slikaNjj,ImageSize->Small]]; 
	,{jj,Length[sol[[1]]]}]
(* Slika sastavljena od svih klastera *)
AA=Table[Table[{sol[[1,j,s,2]],sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[1]]]}];
podDop0=Flatten[AA,1];
sl0=SortBy[podDop0,Last];
pp0=sl0[[All,1]];
slikaN0=Table[Table[pp0[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print["\nSlika svih klastera: "]
Print[Image[slikaN0,ImageSize->Small]]

sol[[2]]

Important !!!
Image composed of selected clusters from WKMeans

JJ={1,2,3};
SelectedColors[sol,JJ]