Grayscale image segmentation: 4 clusters
Important: possible input JPG or BMP image
13-3-2022;  29-1-2023;  7-2-2023

Modules

Needs["ComputationalGeometry`"]

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

Princip

princip[x_,cen_]:=cen[[Position[Table[(x-cen[[i]])^2,{i,Length[cen]}],Min[Table[(x-cen[[i]])^2,{i,Length[cen]}]]][[1,1]]]]

WF

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

Indexes

(* Klasteri trebaju biti poredani od najmanjeg (najsvjetlije) prema najvećem (najtamnijem) centru *)
MD1[PI_,p_,Ind_]:=Module[{n=1,k=Length[PI],z,q,edge,dist},
z=q=Table[0,{j,k}];
Do[
z[[j]]=Mean[PI[[j,All,2]]];
q[[j]]=Quantile[Abs[PI[[j,All,2]]-z[[j]]],p] ;
,{j,Length[PI]}];
If[Ind≠0,Print["z=",z,"\nq=",Round[q,-.01],";  p=",p]];
edge={Min[z[[1]]+q[[1]],Max[PI[[1,All,2]]]],Max[z[[1]]-q[[1]],Min[PI[[1,All,2]]]]};
Do[
edge=Append[edge,Min[z[[j]]+q[[j]],Max[PI[[j,All,2]]]]];
edge=Append[edge,Max[z[[j]]-q[[j]],Min[PI[[j,All,2]]]]];
,{j,2,k}];
(* Sortiranje rubova od najmanjeg prema najvećem *)
edge=Sort[edge];
If[Ind≠0,Print["edge=",Round[edge,-.001]]];
dist={edge[[3]]-edge[[2]]}; 
If[Length[edge]>4,
Do[
dist=Append[dist,edge[[j+1]]-edge[[j]]];
,{j,4,Length[edge]-2,2}]
	]; 
If[Ind≠0,Print["dist=",Round[dist,-.001]]];
{(m-kk)Min[dist]/WF[PI],q,edge,dist}
]

(* Calinski-Harabasz index *)
Clear[VCH]
VCH[PI_]:=Module[{k=Length[PI],m,FF,mj,F,G},
FF=Table[0,{j,Length[PI]}];
m=Sum[Length[PI[[j]]],{j,Length[PI]}];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
set=Flatten[PI,1];
cc=Mean[WeightedData[set[[All,2]],set[[All,1]]]];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-z[[j]] ]^2,{s,mj}];
   ,{j,Length[PI]}];
F=Total[FF];
G=Sum[Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}]*Norm[z[[j]]-cc]^2,{j,Length[PI]}];
(m-Length[PI])*G/((Length[PI]-1)F)
                                     ]

(* Davies-Bouldin *)
Clear[VDB]
VDB[PI_]:=Module[{FF,VV,mj,sum=0},
FF=VV=Table[0,{j,Length[PI]}];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,i,1]]*Norm[PI[[j,i,2]]-z[[j]]]^2,{i,mj}];
VV[[j]]=Sqrt[FF[[j]]/Sum[PI[[j,i,1]],{i,mj}]]
   ,{j,Length[PI]}];
Do[
max={};
Do[
If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/Norm[z[[i]]-z[[j]]]]]
,{i,Length[PI]}];
sum=sum+Max[max]
,{j,Length[PI]}];
sum/Length[PI]
            ]

(* Simplified Silhouette With Criterion *)
Clear[SSWC]
SSWC[PI_]:=Module[{k=Length[PI],a,b,dd,w},
set=Flatten[PI,1];m=Length[set];
a=b=Table[0,{i,m}];
w=set[[All,1]];A=set[[All,2]];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Do[
dd=Sort[Table[w[[i]] Norm[A[[i]]-z[[j]]]^2,{j,k}]];
a[[i]]=dd[[1]];b[[i]]=dd[[2]];
,{i,m}];
Sum[(b[[i]]-a[[i]])/Max[{a[[i]],b[[i]]}],{i,m}]/m
]

(* Dunn index *)
Clear[VD]
VD[PI_]:=Module[{z,DD,diam},
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
DD=Table[If[i≠j,Min[Table[Norm[PI[[i,p,2]]-PI[[j,q,2]]],{p,Length[PI[[i]]]},{q,Length[PI[[j]]]}]],Infinity],{i,Length[PI]},{j,Length[PI]}];
diam=Table[Max[Table[Norm[PI[[i,p,2]]-PI[[i,q,2]]],{p,Length[PI[[i]]]},{q,Length[PI[[i]]]}]],{i,Length[PI]}];
Min[DD/Max[diam]]
]

Weighted k - means   for the data with additional information

WKMeansDop[w_,B_,z_,Ind_:0]:=Module[{m=Length[B],k=Length[z],n=Length[B[[1]]],c0,Pod,imin,PI,It=0},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,m}]; 
c0=z; 
PI=Table[{},{i,k}]; 
			(* Minimal distance principle *)
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
(*Print["Initial centers: ",z,"; 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"]];
Chop[Norm[c-c0]]≠0&&It<5,c0=c;It=It+1;
				(* New partition *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
               ];
{PI,c0,WF[PI],It}]

Important !!!
Image consists of clusters from WKMeans in GrayLevel of its centers

SelectedGrayLevel[sol_]:=Module[{AAA,podDop0,sl0,pp0},
		(* The whole image put in GrayLevel [[1]] *)
Do[
AAA=Table[Table[{1,sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[1]]]}];
		(* jj-th cluster put in color sol[[2,jj]] *)
Which[
jj==2,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==3,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==4,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==5,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==6,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==7,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==8,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}];
	];
		(* Image construction *)
podDop0=Flatten[AAA,1];
sl0=SortBy[podDop0,Last];
pp0=sl0[[All,1]];
slikaN0[jj]=Table[Table[pp0[[i+k bb]],{i,bb}],{k,0,aa-1}];
(* Print[Image[slikaN0[jj],ImageSize->Small]]; *)
		(* Color of center sol[[jj]] *)
(* Print[Graphics[{GrayLevel[sol[[2,jj]]],Disk[]},ImageSize->20]];*)
,{jj,2,Length[sol[[2]]]}]; 
						(* All clusters together *)
		(* The whole image put in GrayLevel [[1]] *)
AAA=Table[Table[{1,sol[[1,j,s,3]]},{s,Length[sol[[1,j]]]}],{j,Length[sol[[1]]]}];
		(* For all jj, jj-th cluster put in color sol[[2,jj]] *)
 Do[
Which[
jj==2,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==3,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==4,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==5,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==6,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==7,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}],
jj==8,
AAA[[jj]]=Table[{sol[[2,jj]],sol[[1,jj,s,3]]},
{s,Length[sol[[1,jj]]]}];
	];
,{jj,2,kk}];
			(* Image construction *)
podDop0=Flatten[AAA,1];
sl0=SortBy[podDop0,Last];
pp0=sl0[[All,1]];
slikaNAll=Table[Table[pp0[[i+k bb]],{i,bb}],{k,0,aa-1}];
Print["Each cluster in GrayLevel of its center:"];
Print[Image[slikaNAll,ImageSize->Small]]
]

Image processing (input is possiblein BMP or JPG format)

Processing[Fig_,Color_]:=Module[{},
or=ImageResize[Fig,150];
fig=ImageData[or];
pod=Join@@fig;  
				(* BMP image *)
If[Length[Dimensions[fig]]==3,
{aa,bb,cc}=Dimensions[fig];
			(* Additional information *)
slikaDop=Table[{fig[[i,j,1]],{i,j}},{i,aa},{j,bb}];
	];
			(* JPG image *)
If[Length[Dimensions[fig]]==2,
{aa,bb}=Dimensions[fig];
		(* Additional information *)
slikaDop=Table[{fig[[i,j]],{i,j}},{i,aa},{j,bb}];
	];
				(* Continue *)
If[Length[Dimensions[fig]]==2,format="JPG",format="BMP"];
podDop=Flatten[slikaDop,1];
w=Table[1,{i,Length[podDop]}];
				(* Image construction *)
pp=podDop[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}];
Print["\nOriginal: ",format," Dimension image: ",{aa,bb,cc}=Dimensions[slikaDop],";  |pod| = ",Length[podDop]];
(* sl0=Print[Image[slikaN,ImageSize->Small]]; *)
		(* Defining of objective function and domenas + SymDIRECT *)
AA=podDop;Nf=0;m=Length[AA]; 
IndMD1=IndSSWC=IndCH=IndDB={};
Do[
kk=it;
If[kk==2,
FFCompile=Compile[{{c1,_Real},{c2,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{kk}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==3,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{kk}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==4,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real},{c4,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2,(c4-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{kk}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==5,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real},{c4,_Real},{c5,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2,(c4-AA[[i,1]])^2,(c5-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{5}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==6,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real},{c4,_Real},{c5,_Real},{c6,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2,(c4-AA[[i,1]])^2,(c5-AA[[i,1]])^2,(c6-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{6}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==7,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real},{c4,_Real},{c5,_Real},{c6,_Real},{c7,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2,(c4-AA[[i,1]])^2,(c5-AA[[i,1]])^2,(c6-AA[[i,1]])^2,(c7-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{7}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];,
If[kk==8,
FFCompile=Compile[{{c1,_Real},{c2,_Real},{c3,_Real},{c4,_Real},{c5,_Real},{c6,_Real},{c7,_Real},{c8,_Real}},Evaluate[Sum[Min[{(c1-AA[[i,1]])^2,(c2-AA[[i,1]])^2,(c3-AA[[i,1]])^2,(c4-AA[[i,1]])^2,(c5-AA[[i,1]])^2,(c6-AA[[i,1]])^2,(c7-AA[[i,1]])^2,(c8-AA[[i,1]])^2}],{i,m}]]];
dom=Table[{0,1},{8}];
Sym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->50,Eps->10^(-4)];
]]]]]]];
Print["\n",Style["SymDIRECT ",16,Blue],Style[kk,16,Blue],Style[" centers: ",16,Blue],Round[ Sym[[1]],-.001],"; F=",Sym[[2]],"; Iter = ",iter];
Print[Table[Graphics[{GrayLevel[Sym[[1,j]]],Disk[]},ImageSize->20],{j,Length[Sym[[1]]]}]];
image=Image[Map[(princip[#,Sym[[1]]])&,fig,{2}],ImageSize->Small];
							(* WKMeans *) 
sol=WKMeansDop[w,AA,Sym[[1]],0];
Print["WKMeans correctiona: ", Round[sol[[2]],-.001],"; F=",sol[[3]],"; Iter = ",sol[[4]]];
Print[Table[Graphics[{GrayLevel[sol[[2,j]]],Disk[]},ImageSize->20],{j,Length[sol[[2]]]}]];
Print["Clusters size: ",Table[Length[sol[[1,i]]],{i,Length[sol[[2]]]}]];
							(* Indexes *)
vmd1=MD1[sol[[1]],p,0];  IndMD1=Append[IndMD1,vmd1[[1]]];
vss=SSWC[sol[[1]]];IndSSWC=Append[IndSSWC,vss];
vch=VCH[sol[[1]]];IndCH=Append[IndCH,vch];
vdb=VDB[sol[[1]]];IndDB=Append[IndDB,vdb];
Print["Gaps between p-clusters from left to right: ",Round[vmd1[[4]],-.001],"\nThe smallest ",Min[Round[vmd1[[4]],-.001]]," is attain between ",aaa=Ordering[vmd1[[4]],1][[1]]," and ",aaa+1," cluster from the left"];
Print["kk=",kk,"; MD1=",vmd1[[1]],"; CH=",vch,"; DB=",vdb,"; SSWC=",vss];
qq=vmd1[[2]];
If[Color==1,
				(* Color figure *)
data0=Table[Table[{sol[[1,jj,s,2]],0},{s,Length[sol[[1,jj]]]}],{jj,Length[sol[[1]]]}];
sldata0= Show[Table[ListPlot[data0[[j]],PlotStyle->{boja[[j]],PointSize[.007],Opacity[.8]}],{j,Length[sol[[1]]]}]];
slcen0=ListPlot[Table[{sol[[2,jj]],0},{jj,Length[sol[[1]]]}],PlotStyle->{Black,Opacity[.8],PointSize[.02]},PlotMarkers->{"✶",20}];
sl1=Show[sldata0,slcen0];
Print[Show[sl1,Table[Graphics[{boja[[s]],  Opacity[0.25],Rectangle[{vmd1[[3,2s-1]],-0.025},{vmd1[[3,2s]],0.025}]}],{s,kk}],AxesOrigin->{0,0},Axes->{True,False},PlotRange->{{0,1},{-.1,.1}},AspectRatio->.1,ImageSize->500]],
				(* Grayscale figure *)
data0=Table[Table[{sol[[1,jj,s,2]],0},{s,Length[sol[[1,jj]]]}],{jj,Length[sol[[1]]]}];
sldata0= Show[Table[ListPlot[data0[[j]],PlotStyle->{GrayLevel[.05],PointSize[.005],Opacity[.8]}],{j,Length[sol[[1]]]}]];
slcen0=ListPlot[Table[{sol[[2,jj]],0},{jj,Length[sol[[1]]]}],PlotStyle->{Orange,Opacity[.8],PointSize[.02]},PlotMarkers->{"✶",15}];
sl1=Show[sldata0,slcen0];
Print[Show[Table[Graphics[{GrayLevel[sol[[2,s]]], Rectangle[{vmd1[[3,2s-1]],-0.05},{vmd1[[3,2s]],0.05}]}],{s,kk}],sl1,AxesOrigin->{0,0},Axes->{True,False},PlotRange->{{0,1},{-.1,.1}},AspectRatio->.1,ImageSize->200]];
]
				(* Clusters in k-partition *) 
SelectedGrayLevel[sol];
,{it,2,KK}];
Print["\nMD1: ",IndMD1,"; MAPart=",Ordering[IndMD1,-1][[1]]+1,"\nIndCH: ",IndCH,"; MAPart=",Ordering[IndCH,-1][[1]]+1, 
"\nIndDB: ",IndDB,"; MAPart=",Ordering[IndDB,1][[1]]+1, 
"\nIndSSWC: ",IndSSWC,"; MAPart=",Ordering[IndSSWC,-1][[1]]+1];
Print[GraphicsGrid[{{ListPlot[Table[{i+1,IndMD1[[i]]},{i,Length[IndMD1]}],Joined->True,GridLines->Automatic,AxesOrigin->{0,.5},PlotLabel->"MD1"], ListPlot[Table[{i+1,IndCH[[i]]},{i,Length[IndCH]}],Joined->True,GridLines->Automatic,AxesOrigin->{0,40000},PlotLabel->"CH"],ListPlot[Table[{i+1,IndDB[[i]]},{i,Length[IndDB]}],Joined->True,GridLines->Automatic,AxesOrigin->{0,.5},PlotLabel->"DB"],ListPlot[Table[{i+1,IndSSWC[[i]]},{i,Length[IndSSWC]}],Joined->True,GridLines->Automatic,AxesOrigin->{0,.8},PlotLabel->"IndSSWC"]}}]];
]

Image segmentation

An example : letter A

Fig=Image[Import["Desktop/Priv/Taler1.jpg"]]
(* kk=3;   Number of clusters *)
p=.7; KK=6;
boja={Red,Green,Blue,Orange,Purple,Magenta,Gray, Brown};
color=0;
Processing[Fig,color] (* If color ≠ 0, figures of clusters are in color; Otherwise figures of clusters are grayscale *)

tab={1,2,3,4,5,6};
ListPlot[Table[{i+1,IndMD1[[i]]},{i,Length[IndMD1]}],Joined->True,GridLines->{tab,{20,30,40}},
AxesOrigin->{0,10},Ticks->{Automatic,{20,30,40}},ImageSize->150]
ListPlot[Table[{i+1,IndCH[[i]]},{i,Length[IndCH]}],Joined->True,GridLines->{tab,{80000,100000,120000}},AxesOrigin->{0,70000},Ticks->{Automatic,{80000,100000,120000}},ImageSize->150]
ListPlot[Table[{i+1,IndDB[[i]]},{i,Length[IndDB]}],Joined->True,GridLines->{tab,{.45,.50,.55}},AxesOrigin->{0,.40},Ticks->{Automatic,{.45,.50,.55}},ImageSize->150]
ListPlot[Table[{i+1,IndSSWC[[i]]},{i,Length[IndSSWC]}],Joined->True,GridLines->{tab,{.89,.90,.91}},AxesOrigin->{0,.88},Ticks->{Automatic,{.89,.90,.91}},ImageSize->150]

Angiography
Import image.pre and call Processing
Import image.post and call Processing

Fig=Image[Import["Desktop/Priv/BelcicPre-13-10-2017.jpg"]]
p=.7; KK=8; (* Number of clusters: 2,...,8 *)
boja={Red,Green,Blue,Orange,Purple,Magenta,Yellow, Brown};
Processing[Fig,0] (* If color ≠ 0, figures of clusters are in color; Otherwise figures of clusters are grayscale *)