WKMeans  algorithm  for the data with additional information
19-12-2025, 4-3-2026

(* Koristi se artificijelnim setovima podataka iz [0,1]^2 i na podacima color image melanoma;
Na setovima artificijelnih podataka iz [0,1]^2 opcije WKMeansAdd koristi klasični k-means algorithm a, opcija WKTriAdd koristi Lemmu; 
Pozivom FigureKM iscrtavaju se tri slike Podaci Initial approximation Results; 
Color image melanoma od slike konstruira set s elementima {{RGB-trojke},{Pixels}};
Prelazom na jednobojnu sliku korištenjem SymDIRECT određuje se i.a. za kk-partition;
Pozivom "WKMeansAdd" s 1 iter određuju se sivi klasteri u sivoj kk-partition;
Korištenjem datacolor i odgovarajućih centroida u "WKMeansAdd" ili "WKTriAdd" određuje se kk-partition sa Spherical color p-Clusters;
 Korištenjem datacolor i istih centroida u "MWKMeansDop" određuje se kk-partition s Ellipsoidal color Dropp-Clusters; 
Iscrtava se siva kk-part i color WKMeans kk-part i color MWKMeans kk-part s p-clusters; *)

Modules

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

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);

Objective functions

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]]*EuclideanDistance[PI[[j,s,2]],c[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

WG[PI_]:=Module[{k=Length[PI], A=Flatten[PI,1],c},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Return[Sum[Length[PI[[j]]]*EuclideanDistance[c[[j]],Mean[WeightedData[A[[All,2]],A[[All,1]]] ]],{j,k}]]
]

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]]]}] ]
];

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]
];

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}
]
];

Weighted k - means  for data with additional information

(* Elements of data set B are of the form {{},{},{}};
For example: {0.5654,{0.870588, 0.819608,0.862745},{1,1}} *)
WKMeansAdd[w_,B_,z_,IT_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1,1]]],c0,Pod,q,imin,PI,It=0,eps=.00005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,m}]; 
c0=z; 
PI=Table[{},{i,k}]; 
If[IT!=1,Print[Style["WKMeans ",16,Blue],k,Style["-partition with Spherical p-clusters:",16,Blue]]];
			(* Minimal distance principle *)
Do[
tab=Table[EuclideanDistance[Pod[[i,2]],c0[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
						(* New centers *)
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition
				in the case n=2 - can be extended to an arbitrary n *)
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition ","\nTry with a new initial centers"];
	];
absc=ChessboardDistance[c,c0];
absc>eps&&It<IT,c0=c;It=It+1;
				(* New partition *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[EuclideanDistance[Pod[[i,2]],c0[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
If[Ind≠ 0,Print[It,":  absc=",absc,";   F = ",WF[PI]]];
               ];
PIP=PI; rrest={};
q=Table[{},{j,k}]; 
Do[
q[[j]]=.5EuclideanDistance[c0[[j]],Nearest[Drop[c0,{j}],c0[[j]]][[1]]];
,{j,k}]; 
{PI,c0,WF[PI],It,q}]

Separability for spherical 2-partition

(* Elements of data set B are of the form {{},{},{}};
Ako je (z2-z1,a-pol)<=0, kut između z2-z1 i a-pol je tupi: element a pripada PI[[1]] *)
Clear[q];WKSep2[w_,B_,z_,IT_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1,1]]],cen,imin,It=1,eps=.0005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,m}]; 
cen=z;  
PI=q=Table[{},{i,k}]; 
Print[Style["WKSep ",16,Blue],2,Style["-partition with Spherical p-clusters - Separability ",16,Blue]];
					(* Grouping based on Lemma *)
pol=Mean[cen];
zz=cen[[2]]-cen[[1]];
PI[[1]]=Select[Pod,zz.(#[[2]]-pol)<=0&];
(* PI[[2]]=Select[Pod,zz.(#[[2]]-pol)>0&]; *)
PI[[2]]=Complement[Pod,PI[[1]]];
q[[1]]=q[[2]]=Table[.5EuclideanDistance[cen[[1]],cen[[2]]],{j,k}];
If[Ind≠ 0, 
Print[0,". |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  F=",WF[PI]]; ];
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition *)
If[Length[c]<Length[cen],
Print["There was a case of reducing the number of clusters in the partition ","\nTry with a new initial centers"];
	];
absc=ChessboardDistance[c,cen];
absc>eps&&It<IT,cen=c;It=It+1;
					(* New partition *)
PI=Table[{},{j,k}]; 
pol=Mean[cen];
zz=cen[[2]]-cen[[1]];
PI[[1]]=Select[Pod,zz.(#[[2]]-pol)<=0&];
(*PI[[2]]=Select[Pod,zz.(#[[2]]-pol)>0&]; *)
PI[[2]]=Complement[Pod,PI[[1]]];
q[[1]]=q[[2]]=Table[.5EuclideanDistance[cen[[1]],cen[[2]]],{j,k}];
If[Ind≠ 0,Print[It,": |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  absc=",absc,";   F = ",WF[PI],"\ncen=",cen]];
];
PIP=PI; rrest={};
{PI,c,WF[PI],It,q}
		]

WKVorAdd

(* Elements of data set B are of the form {{},{},{}}; *)
Clear[q];WKVorAdd[w_,B_,z_,IT_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1,1]]],cen,imin,It=0,eps=.0005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,m}]; 
cen=z;  
PI=q=Table[{},{i,k}]; 
Print[Style["WKVorAdd ",16,Blue],k,Style["-partition with Spherical clusters-q-Nearest ",16,Blue]];
					(* Grouping based on Lemma *)
rest=Pod;
Do[
q[[j]]=.5EuclideanDistance[cen[[j]],Nearest[Drop[cen,{j}],cen[[j]]][[1]]];
PI[[j]]=Select[rest,EuclideanDistance[#[[2]],cen[[j]]]<=q[[j]]&];
rest=Complement[rest,PI[[j]]]; 
,{j,k}]; 
				(* MDP *)
Do[
tab=Table[EuclideanDistance[rest[[i,2]],cen[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],rest[[i]]],
{i,Length[rest]}];
				(* p-clusters *)
If[Ind≠ 0,part=Table[PI[[j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,cen,cov,kappa];
FigureVor[PI,cen,q,pclu];
Print[0,". |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  F=",WF[PI]]; ];
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition *)
If[Length[c]<Length[cen],
Print["There was a case of reducing the number of clusters in the partition ","\nTry with a new initial centers"];
	];
absc=ChessboardDistance[c,cen];
absc>eps&&It<IT,cen=c;It=It+1;
					(* New partition *)
PI=Table[{},{j,k}]; 
rest=Pod;
Do[
q[[j]]=.5EuclideanDistance[cen[[j]],Nearest[Drop[cen,{j}],cen[[j]]][[1]]];
PI[[j]]=Select[rest,EuclideanDistance[#[[2]],cen[[j]]]<=q[[j]]&];
rest=Complement[rest,PI[[j]]];  
,{j,k}]; 
rrest=rest; PIP=PI;
			(* k-means on the Rest *)
Do[
tab=Table[EuclideanDistance[rrest[[i,2]],cen[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],rrest[[i]]],
{i,Length[rrest]}];
					(* p-clusters *)
If[Ind≠ 0, part=Table[PI[[j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,cen,cov,kappa];
FigureVor[PI,cen,q,pclu];
Print[It,". |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  F=",WF[PI]]; ];
];
{PI,c,WF[PI],It,q}]

WKTriAdd

(* Elements of data set B are of the form {{},{},{}};
For example: {0.5654,{0.870588, 0.819608,0.862745},{1,1}} *)
WKTriAdd[w_,B_,z_,IT_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1,1]]],d,cen,Pod,imin,PI,It=0,eps=.0005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,m}]; 
cen=z; 
PI=q=Table[{},{i,k}]; 
Print[Style["WKTri ",16,Blue],k,Style["-partition with Spherical p-clusters: q-Nearest ",16,Blue]];
					(* Grouping based on Lemma *)
rest=Pod;
Do[
q[[j]]=.5EuclideanDistance[cen[[j]],Nearest[Drop[cen,{j}],cen[[j]]][[1]]];
PI[[j]]=Select[rest,EuclideanDistance[#[[2]],cen[[j]]]<=q[[j]]&];
rest=Complement[rest,PI[[j]]];  
,{j,k}]; 
If[Ind≠ 0, 
Print[0,". |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  F=",WF[PI]]; ];
			(* The rest are grouped by MDP *)
Do[
tab=Table[EuclideanDistance[rest[[i,2]],cen[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],rest[[i]]],
{i,Length[rest]}];
If[Ind!=0,Print[" |rest| = ",Length[rest]," podijelio se prema MDP"]];
						(* New centers *)
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition
				in the case n=2 - can be extended to an arbitrary n *)
If[Length[c]<Length[cen],
Print["There was a case of reducing the number of clusters in the partition ","\nTry with a new initial centers"];
	];
absc=ChessboardDistance[c,cen];
absc>eps&&It<IT,cen=c;It=It+1;
If[Ind≠ 0,Print[It,": F = ",WF[PI],"\ncen=",cen]];
					(* New partition *)
PI=Table[{},{i,k}]; 
			(* Grouping based on Lemma *)
rest=Pod;
Do[
q[[j]]=.5EuclideanDistance[cen[[j]],Nearest[Drop[cen,{j}],cen[[j]]][[1]]];
PI[[j]]=Select[rest,EuclideanDistance[#[[2]],cen[[j]]]<=q[[j]]&];
rest=Complement[rest,PI[[j]]]; 
,{j,k}]; 
If[Ind≠ 0,Print[It,": |clu|:",el=Table[Length[PI[[j]]],{j,k}],"=",Total[el],";  |Rest|=",Length[rest],";  absc=",absc,";   F = ",WF[PI],"\ncen=",cen]];
rrest=rest; PIP=PI;
			(* The rest is grouped by MDP *)
Do[
tab=Table[EuclideanDistance[rest[[i,2]],cen[[j]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],rest[[i]]],
{i,Length[rest]}];
If[Ind!=0,Print[It," |rest| = ",Length[rest]," podijelio se prema MDP"]];
If[Ind≠ 0,Print[It,":  absc=",absc,";   F = ",WF[PI]]];
               ];
{PI,cen,WF[PI],It,q}]

Mahalanobis Weighted k-Sep-2   za podatke s dopunskim informacijama

MWKMeans2[w_,B_,z_,Ind_:0]:=Module[{m=Length[B],k=Length[z],W=Total[w],tab,min,imin,PI,c0,c,Scov,PIpriv,Pod,n=Length[z[[1]]],pom,It=0,epsilon=.000000001,eps=.0025},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],i},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
c0=z;
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];mj=Table[0,{j,k}]; 
							(* k-triangle *)
rest=Pod;
pol=Mean[c0];
zz=c0[[2]]-c0[[1]];
PI[[1]]=Select[rest,zz.(#[[2]]-pol)<=0&];
(*PI[[2]]=Select[rest,zz.(#[[2]]-pol)>0&];*)
PI[[2]]=Complement[rest,PI[[1]]];
					 (* 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," |rest| = ",Length[rest],": F0(WF) = ",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]]=Covariance[PI[[j,All,2]]]
,{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}];
Do[
Do[
tab=Table[Sqrt[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
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=ChessboardDistance[c,c0];
If[Ind≠ 0,Print[It,": F1(MWF) = ",F1,"; |c-c0| = ",absc,"; delF = ",delF]];
absc>eps&&It<IT, F0=F1; 
]; 
{PI,c,Scov,It,F1}
]

Mahalanobis Weighted k - means   for data with additional information

MWKMeansAdd[w_,B_,z_,IT_,Ind_:0]:=Module[{m=Length[B],k=Length[z],W=Total[w],PIpriv,tab,min,imin,PI,Scov,Pod,c,c0,n=Length[z[[1]]],pom,corr,It=0,epsilon=.000000001,eps=.0005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
c0=z;
						(* Partition determining by MDP *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[EuclideanDistance[Pod[[i,2]],c0[[j]]],{j,k}]//N;
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];API=PI;
                                                           (* Centers determining *)
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;
				(* Determination of covariance matrices Scov *)
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Covariance[PI[[j,All,2]]]
,{j,k}]; 
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-th center (",c[[j]],"): ",Scov[[j]]," is singular"];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[Sqrt[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]]],{j,k}];
imin=PositionSmallest[tab][[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=ChessboardDistance[c,c0];
If[Ind≠ 0,Print[It,": F1 = ",F1,"; |c-c0| = ",absc,"; delF = ",delF]];
absc>eps&&It<IT, F0=F1; 
]; 
			(* The case of reducing the number of clusters in a partition *)
If[Length[c]<Length[c0],
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[SquaredEuclideanDistance[B[[i]],c[[s]]],{s,Length[c]}],Norm[B[[i]]-{u,v}]^2]],{i,m}];
dom=AG;
						(* New DIRECT center *)
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];
						(* New partition *)
PI=Table[{},{j,k}];
Do[
tab=Table[EuclideanDistance[Pod[[i]][[2]],c[[j]]],{j,k}]//N;
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It,F1}
]

GS k - partition

GSkpart[dataGS_,kk_,Ind_]:=Module[{},
(* Defining of domena and Compile Function FCompile((Subscript[x, 1]),...,(Subscript[x, k])) *)
dom=Table[{0,1},{kk}];
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(s-1),1+s-1}]-dataGS[[i,1]]]^2,{s, 1,kk}]],{i,Length[dataGS]}]]];
FCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,kk}]]; 
						(* Call SymDIRECT *)
Nf=0;
tSym=Timing[
minSym=SymDIRECT[FCompile,dom,DMin->10^(-4),MaxIter->100,Eps->10^(-4)]
];
cenSym=Sort[minSym[[1]]];
If[Ind≠ 0,Print[Style["\nSymDIRECT: ",16,Blue],"centers: ",cenSym,Table[Graphics[{GrayLevel[cenSym[[j]]],Disk[]},ImageSize->20],{j,Length[cenSym]}],"\nF = ",minSym[[2]],"; Nf = ",Nf,"; It = ",iter,"; Time = ",tSym[[1]]]];
		(* Call WKMeans: correction + clusters determining *)
                  (* dataGS: 1.komp-graylevel; 2.komp-pixel; 3.komp-color *)
time=Timing[solGS=WKMeansAdd[w,dataGS,cenSym,1,0]];
Print[Style["WKMeans-1-iter: GS-centers+clusters: ",16,Blue],solGS[[2]],"  ",Table[Graphics[{GrayLevel[solGS[[2,j]]],Disk[]},ImageSize->20],{j,Length[solGS[[2]]]}],"\nF = ",solGS[[3]],"; Iteracija=",solGS[[4]],"; Time = ",time[[1]]];
		(* Centroids calculation by using cluster data of graylevel image !!! *)
centroids0=Table[Mean[WeightedData[solGS[[1,j,All,4]],solGS[[1,j,All,1]]]],{j,Length[solGS[[1]]]}];
dataColor=Table[{podDop[[i,1]],podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}]; 
(* dataColor is defined formal, because WKMeansDop calls data with 3 components *)
                                            (* Call WKMeans - Spherical Clusters *)
timeSph=Timing[solColor=WKVorAdd[w,dataColor,centroids0,5,0]];
							(* p-clusters *)
part=Table[solColor[[1,j,All,2]],{j,kk}]; 
cov=Table[Covariance[part[[j]]],{j,kk}];
pclu=pclusters[part,solColor[[2]],cov,kappa];
							(* Results *)
Print[Style["Extended results:",16,Blue],"\nCentroids: ",solColor[[2]],"  ",Table[Graphics[{RGBColor[solColor[[2,j]]],Disk[]},ImageSize->20],{j,Length[solColor[[2]]]}],"\nCLU: ",el=Table[Length[solColor[[1,j]]],{j,kk}],"~",Total[el],"\nF = ",solColor[[3]],"; Iteracija=",solColor[[4]],"; Time = ",timeSph[[1]]];   
FigureVor[solColor[[1]],solColor[[2]],solColor[[5]],pclu];
			      (* Call MWKMeans - Ellipsoidal Clusters  *)                      
time=Timing[mwkm=MWKMeansAdd[w,dataColor,centroids0,IT,0] ];Print[Style["MWKMeans ",16,Blue],kk,Style["-partition with Ellipsoidal p-clusters:\n",16,Blue],"centroids: ",mwkm[[2]],"  ",Table[Graphics[{RGBColor[mwkm[[2,j]]],Disk[]},ImageSize->20],{j,Length[mwkm[[2]]]}],"\nCLU: ",el=Table[Length[mwkm[[1,j]]],{j,kk}],"--",Total[el],"\nF = ",mwkm[[5]],"; Iteracija=",mwkm[[4]],"; Time = ",time[[1]]];
(*  For graylevel image testing; 
Print["MWKM-B-",kk,mwkm[[2]],"; F=",mwkm[[5]],"; it=",mwkm[[4]], "; CPU=", time[[1]]]; *)
{solGS,centroids0,solColor,mwkm}
]

SymDIRECT

(* For a given center and length of the side of the hyperrectangle the function test returns:
  1 If the whole hyperrectangle is located in the region \Delta
0 If a hyperrectangle lies in the region \Delta only partially
-1 If a hyperrectangle appears outside the region \Delta (or if a hyperrectangle has only one point in \Delta)*)
test[c_,dd_]:=Module[{n=Length[c],S,d=dd/2,unutra=1},
S=Table[1,{n}];
Do[
If[c[[i]]-d[[i]]<c[[i+1]]+d[[i+1]],
unutra=0;Break[];
]
,{i,1,n-1}];
If[unutra==0,
Do[
If[c[[i]]+S[[i]]*d[[i]]<c[[i+1]]+d[[i+1]],
If[c[[i]]+S[[i]]*d[[i]]>c[[i+1]]-d[[i+1]],
S[[i+1]]=-1,
unutra=-1;Break[];
]
]
,{i,1,n-1}];
];
Return[unutra];
]

Options[SymDIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
SymDIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,temp2,ind,novic,novid,novifc,noviunutra,red,TMin,dtest,anim={},unutra,brisi},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
SymDIRECT::maxiter="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
iter=0;
(* scaling the function on the hypercube with length of the side 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
unutra={test[c[[1]],d[[1]]]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: searching for 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 *)
brisi={};
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]}];Nf=Nf+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]}];
(* Eliminating those that are not inside *)
If[unutra[[pot]]==0,
(* Central hyperrectangle can go out or stay completely inside *)
If[test[c[[pot]],d[[pot]]]==-1,AppendTo[brisi,{pot}],unutra[[pot]]=test[c[[pot]],d[[pot]]]];
temp={};noviunutra={};
Do[
temp2=test[novic[[i]],novid[[i]]];
If[temp2==-1,AppendTo[temp,{i}],AppendTo[noviunutra,temp2]];
,{i,Length[novic]}];
novic=Delete[novic,temp];
novid=Delete[novid,temp];
novifc=Delete[novifc,temp];
,
(* If the divided hyperrectangle lies completely inside, all subhyperrectangles remain inside too *)
noviunutra=Table[1,{Length[novic]}];
];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
unutra=Join[unutra,noviunutra];
(* Update of actual minimum *)
temp=Position[novifc,Min[novifc]][[1,1]];
If[novifc[[temp]]<TMin[[2]],TMin=N[{novic[[temp]],novifc[[temp]]},10]];
,{pot,pot2}];
(* Deletion of those hyperrectangles which went outside *)
c=Delete[c,brisi];
d=Delete[d,brisi];
dmax=Delete[dmax,brisi];
fc=Delete[fc,brisi];
unutra=Delete[unutra,brisi];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[SymDIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure in 2D case *)
(*If[n==2 && iter<100,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Tooltip[Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10],unutra[[i]]]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* Back to the 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];
]

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];
]

Figures

(* Illustration for k-means algorithm *)
FigureKM[dol_,sol_,PIP_,rest_]:=Module[{k=Length[sol[[2]]]},
				(* Input: given data *)
AG={{0,1.01},{0,1}};
slpart0=ListPlot[Table[dol[[1,j,All,2]],{j,Length[dol[[1]]]}],PlotStyle->{PointSize[.014]}, GridLines->Automatic,PlotRange->AG];
slcen0=ListPlot[dol[[2]],PlotStyle->{Orange,PointSize[.03],Opacity[.8]},PlotRange->AG];
sl0=Show[slpart0,slcen0,DiagramPlot[cen],AspectRatio->Automatic,PlotLabel->"Data"<>" ("<>ToString[k mc]<>")"];
				(* Prema Voronoi *)
slPI=ListPlot[Table[PIP[[j,All,2]],{j,k}],PlotStyle->{PointSize[.014]}, GridLines->Automatic,PlotRange->AG];
slrest=ListPlot[rrest[[All,2]],PlotStyle->{Purple,PointSize[.016]}, GridLines->Automatic,PlotRange->AG];
slDisk=Table[Graphics[{LightBlue,Opacity[.2],Disk[sol[[2,j]],sol[[5,j]]]}],{j,k}];slCir=Table[Graphics[{Thin,Circle[sol[[2,j]],sol[[5,j]]]}],{j,k}];
slcen1=ListPlot[sol[[2]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.03],Opacity[.8]}];
sl1=Show[slPI,slDisk,slCir,slrest,slcen1,DiagramPlot[sol[[2]]],AspectRatio->Automatic,PlotLabel->"Voronoi + rest "<>" (|rest|="<>ToString[Length[rrest]]<>")"];
			           (* Output: Results *)
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}],PlotStyle->{PointSize[.014]}, GridLines->Automatic,PlotRange->AG];
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.03],Opacity[.8]},PlotRange->AG];
sl2=Show[slPI,slc,DiagramPlot[sol[[2]]],AspectRatio->Automatic,PlotLabel->"Result"];
Print[GraphicsGrid[{{sl0,sl1,sl2}},ImageSize->Large]];
If[Length[rrest]!=0,mm=Length[Flatten[PIP,1]]+Length[rest];
Print[Style["Lemma~2 obuhvaća ",13,Blue],Style[Round[100. (1-Length[rrest]/mm)],14,Red],Style["% podataka. Na ostatak (purple točkice) primijenjen je MDP; CPU=",13,Blue],Style[time[[1]],14,Red]]];
]

(* Illustration for k-Voronoi method *)
FigureVor[PI_,cen_,q_,pclu_]:=Module[{k=Length[cen]},
A=Flatten[Table[PI[[j,All,2]],{j,k}],1];
If[n==2,
AG={{0,1.01},{0,1}};
slVcirc=Show[
ListPlot[A,PlotStyle->{Gray,PointSize[.012]}],DiagramPlot[cen], 
Table[Graphics[{LightBlue,Disk[cen[[j]],q[[j]]]}],{j,k}],Table[Graphics[{Thin,Circle[cen[[j]],q[[j]]]}],{j,k}],
ListPlot[cen,PlotMarkers->{"✶",15},PlotStyle->{Black,PointSize[.025],Opacity[.6]}],
ListPlot[Table[PI[[j,All,2]],{j,Length[PI]}],PlotStyle->{PointSize[.014]}],
GridLines->Automatic,AspectRatio->Automatic,PlotRange->AG,Frame->True,PlotLabel->"Voronoi circles + all data"];
slpclu=Show[
Show[Table[ConvexHullMesh[pclu[[j]]],{j,k}]],DiagramPlot[cen], 
ListPlot[pclu,PlotStyle->{PointSize[.015]}],
ListPlot[cen,PlotMarkers->{"✶",13},PlotStyle->{Black,PointSize[.025],Opacity[.6]}],
GridLines->Automatic,AspectRatio->Automatic,PlotRange->AG,Frame->True,PlotLabel->"p-clusters"<>" (p="<>ToString[p]<>")"];
Print[GraphicsGrid[{{slVcirc,slpclu}},ImageSize->Medium]],
							(* For n=3 *)
AG={{0,1},{0,1},{0,1}};
sldata=Show[
ListPointPlot3D[cen,PlotStyle->{Black,PointSize[.04],Opacity[.6]}],
ListPointPlot3D[A,PlotStyle->{PointSize[.013],Opacity[.6]}],PlotRange->AG,BoxRatios->{1,1,1},Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},PlotLabel->"Data"<>" ("<>ToString[Length[A]]<>")"];
slsol=Show[
Graphics3D[{Opacity[.5],Table[Sphere[solColor[[2,j]],1.1q[[j]]],{j,Length[solColor[[2]]]}]}],
ListPointPlot3D[solColor[[2]],PlotStyle->{Black,PointSize[.04]},PlotRange->AG], 
Table[ListPointPlot3D[PIP[[j,All,2]],PlotStyle->{RGBColor[solColor[[2,j]]],PointSize[.013],Opacity[.7]}],{j,Length[PIP]}],
ListPointPlot3D[rest[[All,2]],PlotStyle->{Red,PointSize[.013],Opacity[.7]}],
PlotRange->AG,BoxRatios->{1,1,1},Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},PlotLabel->"Lema~2 + rest"];
slpclu=Show[
slSph=Graphics3D[{Opacity[.6],Table[Sphere[solColor[[2,j]],1.1q[[j]]],{j,Length[solColor[[2]]]}]}],ListPointPlot3D[solColor[[2]],PlotStyle->{Black,PointSize[.05],Opacity[.5]},PlotRange->AG],
Table[ListPointPlot3D[pclu[[j]],PlotStyle->{RGBColor[solColor[[2,j]]],PointSize[.015],Opacity[.7]}],{j,Length[solColor[[1]]]}],
PlotRange->AG,BoxRatios->{1,1,1},Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},PlotLabel->"p-clusters"<>" (p="<>ToString[p]<>")"];Print[GraphicsGrid[{{sldata,slsol,slpclu}},ImageSize->Large]];
If[Length[rrest]!=0,Print["Lemma~2 obuhvaća ",Style[Round[100. (1-Length[rrest]/Length[A])],14,Red],"% podataka. Na ostatak (crvene točkice) je primijenjen MDP"]];
]
]

Examples with artificial data sets

Example - 2  clusters

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd, and k - separability method (only for k = 2) WKSep2

(* dol - initial data;  sol - solution *)
(* cpu:  WKMeansAdd (0.22);  WKTriAdd (0.125);  WKVorAdd (0.17); WKSep2 (0.10) *)
SeedRandom[13];mc=500;
cen={{.45,.35},{.6,.75}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",Table[Length[dol[[1,j]]],{j,k}]];
					(* MethodDop *)
n=2;  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
z0={{.3,.4},{.7,.6}};
time=Timing[sol=WKSep2[w,B,z0,10,0]];
				(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results Extended *)
Print[Style["Extended results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; Total-CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],";   d(z0,cen1)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol,cen1)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\ncen1: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 4  clusters

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 -solution *)
SeedRandom[13];mc=100;
cen={{.4,.3},{.25,.55},{.6,.75},{.75,.3}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",tab=Table[Length[dol[[1,j]]],{j,k}]," = ", Total[tab]];
					(* WKMeansDop *)
z0={{.3,.3},{.4,.5},{.5,.8},{.6,.3}};
n=Length[cen[[1]]];  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
(* Parameter Ind određuje hoće li se printati slike i podaci međurezultata *)
time=Timing[sol=WKVorAdd[w,B,z0,10,0]];
					(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results *)
Print[Style["Results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; |Rest|=",Length[rrest],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],";   d(z0,z^*)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol[[2]],z^*)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\nz^*: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 6  clusters

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 - solution *)
SeedRandom[13];mc=500;
cen={{.8,.3},{.65,.5},{.6,.7},{.3,.8},{.4,.2},{.2,.4}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",el=Table[Length[dol[[1,j]]],{j,k}]," = ",Total[el]];
					 (* MethodAdd *)
z0={{.7,.3},{.6,.4},{.5,.7},{.2,.8},{.4,.3},{.3,.5}};
n=2;  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
time=Timing[sol=WKVorAdd[w,B,z0,20,0]];
						(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results Extended *)
Print[Style["Results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; |Rest|=",Length[rrest],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],"\nd(z0,z^*)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol[[2]],z^*)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\nz^*: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 8  clusters

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 -  solution *)
SeedRandom[13];mc=500;
cen={{.7,.2},{.6,.4},{.5,.6},{.25,.75},{.35,.2},{.18,.4},{.83,.6},{.6,.8}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",el=Table[Length[dol[[1,j]]],{j,k}]," = ",Total[el]];
					(* MethodAdd *)
z0={{.8,.2},{.7,.5},{.4,.5},{.2,.7},{.2,.2},{.25,.5},{.9,.7},{.5,.9}};
n=2;  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
time=Timing[sol=WKTriAdd[w,B,z0,5,0]];
						(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results *)
Print[Style["Results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; |Rest|=",Length[rrest],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],";   d(z0,z^*)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol[[2]],z^*)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\nz^*: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 10  clusters

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 -solution *)
SeedRandom[13];mc=500;
cen={{.78,.82},{.83,.55},{.55,.75},{.23,.82},{.5,.5},{.75,.2},{.6,.37},{.3,.6},{.35,.2},{.18,.35}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",el=Table[Length[dol[[1,j]]],{j,k}]," = ",Total[el]];
					(* MethodAdd *)
z0={{.7,.85},{.8,.5},{.55,.75},{.3,.9},{.5,.6},{.8,.3},{.6,.3},{.3,.5},{.2,.1},{.18,.35}}; 
n=2;  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
time=Timing[sol=WKTriAdd[w,B,z0,5,0]];
					(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results *)
Print[Style["Results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; |Rest|=",Length[rrest],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],";   d(z0,z^*)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol[[2]],z^*)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\nz^*: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 5 big   clusters in R^2

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 -solution *)
SeedRandom[13];mc=1000;
cen={{.4,.25},{.25,.55},{.8,.7},{.75,.3},{.5,.8}}; k=Length[cen]; n=Length[cen[[1]]];
kov=Table[IdentityMatrix[2].005,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],kov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",Table[Length[dol[[1,j]]],{j,k}]];
					(* WKMeansDop *)
z0={{.3,.3},{.4,.5},{.5,.8},{.7,.7},{.6,.3}};
n=2;  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
time=Timing[sol=WKVorAdd[w,B,z0,6,0]];
						(* p-clusters *)
part=Table[sol[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,sol[[2]],cov,kappa];
					(* Results *)
Print[Style["Results:",16,Blue],"\nF= ",sol[[3]], ";  absc=",absc,"; It = ",sol[[4]],"; |Rest|=",Length[rrest],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[sol[[1,j]]],{j,k}],"=",Total[el],";   d(z0,z^*)]=",DiscreteHausdorffDistance[sol[[2]],z0],";  d(dol[[2]],z^*)=",DiscreteHausdorffDistance[dol[[2]],sol[[2]]],"\nz^*: ",sol[[2]]  ];
	                                            (* Figures *)                                     
FigureVor[sol[[1]],sol[[2]],sol[[5]],pclu]; 
FigureKM[dol,sol,PIP,rrest]

Example - 3 big   clusters in R^3

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd

(* dol - initial data;  sol - cen1 -solution *)
SeedRandom[13];mc=2000;
cen={{.3,.3,.2},{.7,.4,.5},{.5,.8,.7}}; k=Length[cen]; n=Length[cen[[1]]];
cov=Table[IdentityMatrix[3].015,{j,k}];
data=Table[RandomReal[MultinormalDistribution[cen[[j]],cov[[j]]],mc],{j,k}];
A=Flatten[data,1];
B=Table[{A[[i]],i,i},{i,Length[A]}];
w=Table[1,{i,Length[A]}];
dol=WKMeansAdd[w,B,cen,1,0];
Print["Given centroids(A) cen = ",cen=dol[[2]]];
Print["|given clu|: ",Table[Length[dol[[1,j]]],{j,k}]];
					(* WKMeansDop *)
z0={{.4,.2,.3},{.6,.5,.4},{.6,.7,.6}};
n=Length[cen[[1]]];  p=.95;  kappa=Quantile[ChiSquareDistribution[n],p];
time=Timing[solColor=WKMeansAdd[w,B,z0,5,0]];
							(* p-clusters *)
part=Table[solColor[[1,j,All,2]],{j,k}]; 
cov=Table[Covariance[part[[j]]],{j,k}];
pclu=pclusters[part,solColor[[2]],cov,kappa];
					(* Results *)
Print[Style["Results:",16,Blue],"\nF= ",solColor[[3]], ";  absc=",absc,"; It = ",solColor[[4]],"; CPU=",time[[1]]];
Print["|clu|:",el=Table[Length[solColor[[1,j]]],{j,k}],"=",Total[el],";   d(z0,z^*)]=",DiscreteHausdorffDistance[solColor[[2]],z0],";  d(dol,z^*)=",DiscreteHausdorffDistance[dol[[2]],solColor[[2]]],"\nz^*: ",solColor[[2]] ];
FigureVor[solColor[[1]],solColor[[2]],solColor[[5]],pclu];

Color images

For the calculation you can choose : ordinary k - means algorithm WKMeansAdd, k - triangle method WKTriAdd, k - Voronoi method WKVorAdd, and k - separability method (only for k = 2) WKSep2

(*ISIC-0000140*)      (* ISIC-0000117 *)     (* Lena:*)

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]] ]
]		
                                                               (* Data *)
podaci=ImageData[Fig];
pod=Flatten[podaci,1];
Print["Image Dimension: ",{aa,bb,cc}=Dimensions[podaci],";  |pod| = ",aa bb]
			(* Adding additional informations (pixels coordinates) *)
slikaDop=Table[{podaci[[i,j]],{i,j}},{i,aa},{j,bb}];
podDop=Flatten[slikaDop,1];
w=Table[1,{i,Length[podDop]}];
m=Length[podDop]; n=3;  p=.95;  IT=10;
(* Ako je Ind2!=0, different images are shown *)
Ind2=1;  kappa=Quantile[ChiSquareDistribution[n],p];
Print["m = ",m,";  n = ",n,";  p = ",p,";  kappa = ",kappa];
					 (* Transition to graylevel image *)
GL={0.2989,0.5870,0.1140}; (* parameter za prijelaz na sivu sliku *)
dataGS=Table[{GL.podDop[[i,1]],podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}];
(* For R - G - B color image:
dataGS=Table[{GL.podDop[[i,1]],podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}];
dataGS=Table[{Norm[ReplacePart[podDop[[i,1]],1->3*podDop[[i,1,1]]]]/4,podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}]; (* R *)
dataGS=Table[{Norm[ReplacePart[podDop[[i,1]],2->3*podDop[[i,1,2]]]]/4,podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}];(* G *)
dataGS=Table[{Norm[ReplacePart[podDop[[i,1]],3->3*podDop[[i,1,3]]]]/4,podDop[[i,2]],podDop[[i,1]]},{i,Length[podDop]}]; (* B *)
*)
				(* Construction of graylevel image *)
If[Ind2!=0,
Print["Grayscale image"];
pp=dataGS[[All,1]];
slikaGS=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print[Image[slikaGS,ImageSize->100]]
  ];
							  (* Call GSkpart *)
            (* Based on GS-centers(SymDIR) an i.a. for color centroids is determined *)
kk=5;gs=GSkpart[dataGS,kk,1];
				(* Plot kk-partition of graylevel image *)
If[Ind2!=0,
Print["\n",kk,Style["-partition of grayscale image: \n",16,Blue]];
PI=Flatten[Table[Table[{gs[[1,2,j]],gs[[1,1,j,s,3]]},{s,Length[gs[[1,1,j]]]}],{j,Length[gs[[1,2]]]}],1];
sl=SortBy[PI,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print[slN=Image[slikaN,ImageSize->Small]];
  ];
			(* Crtanje spherical kk-partitopn of color image *)	
If[Ind2!=0,
Print["\nSpherical ",kk,Style["-partition of color image: \n",16,Blue]];
PI=Flatten[Table[Table[{gs[[3,2,j]],gs[[3,1,j,s,3]]},{s,Length[gs[[3,1,j]]]}],{j,Length[gs[[3,2]]]}],1];
sl=SortBy[PI,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print[slN=Image[slikaN,ImageSize->Small]];
  ];	
				(* Plot Mahalanobis kk-particije color image  *)
If[Ind2!=0,
Print["\nMahalanobis ",kk,Style["-partition of color image: \n",16,Blue]];
PI=Flatten[Table[Table[{gs[[4,2,j]],gs[[4,1,j,s,3]]},{s,Length[gs[[4,1,j]]]}],{j,Length[gs[[4,2]]]}],1];
sl=SortBy[PI,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print[slN=Image[slikaN,ImageSize->Small]];
  ];