Melanoma Color  Image: Mahalanobis  Incremental Algorithm  -  02 - 12 - 2025, 18-1-2025

eps = 0.0005;  MWKMeans2 : only  2-partition;  MWKMeansDop : k-partition  for  k > 2
i. a. for  2-partition  using  a  grayscale  image  and  SymDIRECT

(* Najprije se identificira područje melanoma na slici; i.a. određuje modulom IASymDIR[];
Nakon toga, pozivom MIncAdd algoritma određuje se sljedeći centar DIRECT optimizacijom; Nakon toga poziva se MWKMeans algorithm s do 10 iteracija;
Izabere se samo klaster melanoma (tamnije boje) i na njemu se provede MIncAdd algorithm s KK=4 iteracije;
Tada se određuju eigenvalues cijelih ili p-clusters.
 *)

Modules

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

Squares Euclidean and 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);

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

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

Initial  approximation - 2-partition

IASymDIR[kk_,w_,dataGS_,Ind_:0]:=Module[{z0},
dom=Table[{0,1},{kk}];
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[EuclideanDistance[Table[(Subscript[x, j]),{j,1+(s-1),1+s-1}],dataGS[[i,1]]],{s, 1,kk}]],{i,Length[dataGS]}]]];
FCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,kk}]]; 
						(* Call SymDIRECT *)
Nf=0;
minSym=SymDIRECT[FCompile,dom,DMin->10^(-4),MaxIter->70,Eps->10^(-4)];
cenSym=Sort[minSym[[1]]];
sol2=WKMeansAdd[w,dataGS,cenSym,1,0]; (* Determining clusters size *)Print[Style["IA-SymDIRECT: ",16,Blue],"centers: ",cenSym,Table[Graphics[{GrayLevel[cenSym[[j]]],Disk[]},ImageSize->20],{j,Length[cenSym]}],"\nF(WF) = ",minSym[[2]],"; Nf = ",Nf,"; It = ",iter,"; {|pi_1|,|pi_2|}:",{Length[sol2[[1,1]]],Length[sol2[[1,2]]]}];
z0=Table[Mean[WeightedData[sol2[[1,j,All,4]],sol2[[1,j,All,1]]]],{j,Length[sol2[[1]]]}];
{z0}]

Mahalanobis Weighted k-2 means  for the data with additional information

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}]; 
							(* 2-separable *)
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&];
					 (* 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," |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;
				(* 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}];
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 the data with additional information

(* It is dM used *)
MWKMeansAdd[w_,B_,z_,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=.0025},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{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(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;
				(* 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(MWF) = ",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_,om_]:=Sum[w[[i]]Min[Append[Table[SquaredEuclideanDistance[B[[i,1]],c[[s]]],{s,Length[c]}],Norm[B[[i,1]]-{u,v,om}]^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["New optimized center: ",c];
k=Length[c]; Print["Corrected 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}
]

Weighted k - means  for the 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,imin,PI,It=0,eps=.0005},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],B[[i,3]]},{i,m}]; 
c0=z; 
PI=Table[{},{i,k}]; 
			(* 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]]];
               ];
{PI,c0,WF[PI],It}]

MIncAdd - clustering

(* Mahalanobis Incremental Algorithm for the data with additiona informations *)
MIncAdd[w_,B_,cen_,KK_,Ind_:0]:=Module[{centri,dom,Scov,min},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
				(* Initialization *)
centri=cen;   
Do[
If[brk>1,
If[Ind≠0,Print[Style["Iteration ",Red,Bold,12],Style[brk,Red,Bold,14]," To the existing ",Length[centri]," centers 1 new center is added"]];
k=Length[centri];  
					(* DIRECT Compile *)
time1=Timing[
Clear[u1,u2,u3,F2];
FCompile=Compile[{{u1,_Real},{u2,_Real},{u3,_Real}},
Evaluate[Sum[w[[i]]Min[Append[Table[EuclideanDistance[Pod[[i,2]],centri[[s]]],{s,Length[centri]}],EuclideanDistance[Pod[[i,2]],{u1,u2,u3}]]],{i,Length[B]}] ]];	
F2[u1_,u2_,u3_]:=Apply[FCompile,{u1,u2,u3}];	
dom=AG;					
min=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)];
];
centriN=Append[centri,min[[1]]],
centriN=cen];
							(*  *)
If[Ind≠0,Print["New optimized centers: ",centriN]];
k=Length[centriN];
			(* Call Mahalanobis k-means *)
time2=Timing[ 
If[k==2,sol=MWKMeans2[w,B,centriN,Ind],
 sol=MWKMeansAdd[w,B,centriN,Ind]]
(* sol=MWKMeansAdd[w,B,centriN,Ind]  *)
                    ]; 
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["New MWKMeans refined centers ",sol[[2]]]];
Scov=sol[[3]];
If[Ind≠0,Print["MWKMeans: iter=",sol[[4]],"; MWF = ",FF=MWF[sol[[1]],sol[[2]],sol[[3]]],"; |pi|: ",mjx=Table[Length[sol[[1,j]]],{j,Length[sol[[1]]]}]]];
                                          (* Construction of MWKM image *)
If[Ind2≠0,
dol=sol;
Do[dol[[1,2,s,2]]={1,1,1},{s,Length[dol[[1,2]]]}];PIP=Flatten[Table[Table[{dol[[1,j,s,2]],dol[[1,j,s,3]]},{s,Length[dol[[1,j]]]}],{j,Length[dol[[2]]]}],1];
sl=SortBy[PIP,Last];
pp=sl[[All,1]];
slikaN=Table[Table[pp[[i+k bb]],{i,bb}],{k,0,aa-1}] ;
Print["MWKM",Length[sol[[2]]],"-image.",  " Optimized centers: ",sol[[2]],"; F(MWF) = ",MWF[sol[[1]],sol[[2]],sol[[3]]],"\np-clusters size: ",{Length[sol[[1,1]]],Length[sol[[1,2]]]},"; iter=",sol[[4]]];
Print[slN=Image[slikaN,ImageSize->Small]];
        ];
                                                           (* pPI0, cov, Eigenvalues *)
If[Ind1≠0,
time3=Timing[
pPI0=Table[{},{j,Length[sol[[1]]]}];
PI0=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];  (* Clusters bez koordinata pixela *)
cc=Table[Mean[PI0[[j]]],{j,Length[sol[[1]]]}];
cov=Table[Covariance[PI0[[j]]],{j,Length[sol[[1]]]}];
(*ell=Table[Ellipsoid[cc[[j]], 11.3449  cov[[j]]],{j,Length[PI0]}];*)(* p=.99 *)
pell=Table[Ellipsoid[cc[[j]], kappa  cov[[j]]],{j,Length[sol[[1]]]}];
Do[
Do[
If[PI0[[j,s]]∈pell[[j]] ,pPI0[[j]]=Append[pPI0[[j]],PI0[[j,s]]]  ];
,{s,Length[PI0[[j]]]}]
,{j,Length[sol[[1]]]}];
cov=Table[Covariance[pPI0[[j]]],{j,Length[sol[[1]]]}];
ei=Table[Eigenvalues[cov[[j]]],{j,Length[sol[[1]]]}];
                                                    (* Ellipsoids *)
If[Ind3!=0,
slell=Show[Table[Graphics3D[{RGBColor[cc[[j]]],Opacity[0.4],pell[[j]]}],{j,Length[sol[[2]]]}],Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"}];
sltoc=Table[ListPointPlot3D[pPI0[[j]],PlotStyle->{RGBColor[cc[[j]]],PointSize[.008],Opacity[.8]}],{j,Length[cc]}];
Print[Show[slell,sltoc, ImageSize->150]];
  ];
];
];
					(* Ellipsoids + Eigenvalues *)
If[Ind1==0,
time3=Timing[
PI0=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];                  (* Clusters without pixels coordinates *)                                    
cc=Table[Mean[PI0[[j]]],{j,Length[sol[[1]]]}];
cov=Table[Covariance[PI0[[j]]],{j,Length[sol[[1]]]}];
ei=Table[Eigenvalues[cov[[j]]],{j,Length[sol[[1]]]}];
ell=Table[Ellipsoid[cc[[j]], 11.3449  cov[[j]]],{j,Length[sol[[1]]]}];
pell=Table[Ellipsoid[cc[[j]], kappa  cov[[j]]],{j,Length[sol[[1]]]}];
                                              (* Ellipsoids *)
If[Ind3!= 0,slell=Show[Table[Graphics3D[{RGBColor[cc[[j]]],Opacity[0.4],pell[[j]]}],{j,Length[sol[[1]]]}],Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"}];
sltoc=Table[ListPointPlot3D[PI0[[j]],PlotStyle->{RGBColor[cc[[j]]],PointSize[.008],Opacity[.8]}],{j,Length[cc]}];
Print[Show[slell,sltoc, ImageSize->150]];
   ];
];
];	
centri=sol[[2]];
If[Ind5≠0,
Print[brk+1,Style["-partition of the set data2:",13,Blue]];
Print["Centroids: ",centri, "; MWF= ",sol[[5]]];
Print[Style["Eigenvalues: ",13,Blue],ei];
EIG[[brk]]=Append[EIG[[brk]],ei];
];
,{brk,KK}];
{sol[[1]],sol[[2]],sol[[3]],ei}]

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

The main program (IT of MWKM=5,  TotalTime,  Sqrt[dM])

(*ISIC_0000117*)  (* ISIC_0000140*)

Fig=Image[];     
dim=ImageDimensions[Fig];
Print["Original image: ImageDimensions = ",dim," --- ",dim[[1]] dim[[2]] ]
If[Max[dim]>128,
resize=128; Fig=ImageResize[Fig,resize];
dim=ImageDimensions[Fig];
]							
                                                                      (* Data *)
podaci=ImageData[Fig];
pod=Flatten[podaci,1];
Print["Resize Image Dimension: ",{aa,bb,cc}=Dimensions[podaci],";  |pod| = ",dim[[1]] dim[[2]]]
			(* 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;  kappa=Quantile[ChiSquareDistribution[n],p];
Print["m = ",m,";  n = ",n,";  p = ",p,";  kappa = ",kappa];
						(*  *)
(* definiran je poseban postupak izbora i.a. za MIncAdd; 
IT određuje maksimalni broj iteracija u MWKM i time bitno utječe na CPU;
Ako je Ind≠0, ispisuju se detalji izračunavanja;
						(*  *)
Ako je Ind1=0, određuju se cijeli klasteri, u suprotnom određuju se p-klasteri;
tj. određuju se Eigenvalues cijelih ell, a u protivnom Eig p-ellipsoida;
U oba slučaja crataju se p-clusters jer bi cijeli clu zauzeli veliki dio slike; 
						(*  *)
Ako je Ind2≠0, crta se slika klastera melanoma, u protivnom, ne crta se; 
Ako je Ind3≠0, iscrtavaju se ellipsoids, u protivnom, ne iscrtavaju se; 
						(*  *)
Ako je Ind5=0, obrađuje se bazična slika, u protivnom obrađuje se klaster melanoma 
*)
IT=10; KK=1;  (* KK: number of clusters in the partition *)
Ind=0; Ind1=0; Ind2=1; Ind3=1;Ind5=0;
EV2=EV3=EV4={};
ispis1={};
						(* Grayscale 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]}];
				(* Construction of grayscale 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]]
  ];  *)
				(* Initial approximation *)
time0=Timing[z0=IASymDIR[2,w,dataGS,Ind][[1]]  ]; 
Print["Initial approximation: ",z0]
					(* Call MIncAdd *)
Print[Style["Searching for the melanoma cluster",Red,15]];
minc2=MIncAdd[w,podDop,z0,KK,Ind]; 
			(* Melanoma cluster: data2 + Eigenvalues *)
time4=Timing[
data2=Table[{minc2[[1,1,i,2]],minc2[[1,1,i,3]],minc2[[1,1,i,2]]},{i,Length[minc2[[1,1]]]}];
data2GS=Table[{GL.minc2[[1,1,i,2]],minc2[[1,1,i,3]],minc2[[1,1,i,2]]},{i,Length[minc2[[1,1]]]}];
w=Table[1,{i,Length[minc2[[1,1]]]}];
];
If[Ind1==0,Print[Style["Eigenvalues of the melanoma cluster:  ",13,Blue],minc2[[4,1]]],
Print[Style["Eigenvalues of the melanoma p-cluster:  ",13,Blue],minc2[[4,1]]];
];
			(* |clu|  *)
Print[Style["{|clu|, eps, Iter, F(MWF)}: ",13,Blue],Table[{Length[sol[[1,j]]]} ,{j,Length[sol[[1]]]}],";  eps=",absc,"; Iter=",sol[[4]],"; F(MWF)=",sol[[5]]]
ispis1=Append[ispis1,minc2[[4,1]]];      
Print[Style["Centroids: ",13,Blue],minc2[[2]]] 
Print[Style["CPU {i.a, MWKM2, PI+ei, data2}: ",13,Blue],TT={time0[[1]],time2[[1]],time3[[1]],time4[[1]]},"; Total = ",TTT=Total[TT]];  
(*  *)
					(* Clustering Melanoma data *)
(*  *)
IT=5; KK=3;
Ind=1; Ind1=0; Ind2=0; Ind3=1;Ind5=1;
EIG=Table[{},{j,KK}];					
Print[Style["\nSearching for the melanoma partitions",Red,15]];
				(* Initial approximation *)
time00=Timing[z0=IASymDIR[2,w,data2GS,Ind][[1]]  ];
Print[Style["data2: Initial approximation: ",13,Blue],z0];
mincKK=MIncAdd[w,data2,z0,KK,Ind];

Print[Style["CPU {i.a., DIR, MWKM, PI+ei}: ",13,Blue],TT={time00[[1]],time1[[1]],time2[[1]],time3[[1]]},"; Total = ",TTT=TTT+Total[TT]];  
EV2=Append[EV2,EIG[[1]]];
EV3=Append[EV3,EIG[[2]]];
EV4=Append[EV4,EIG[[3]]];
Print[Style["Eigenvalues of melanoma cluster (the first cluster of podDop) ",16,Red],ispis1];
Print[Style["Eigenvalues of partitions(2,3,4) clusters ",16,Red]];
Print["EV2: ",X2=EV2[[All,1]]]
Print["EV3: ",X3=EV3[[All,1]]]
Print["EV4: ",X4=EV3[[All,1]]]