function [clusterSignals,Lowest1StMomentSig,TF_sel,signalsAll,signalsAllsel,idxSelClu] = DoKmeansClustering(data4D,TF,timeAxis,pass,clusters,evalFM,visu)
%DoKmeansClustering(data4D,TF_ROU,timeAxis,pass,clusters)
disp(['k-means běh ' num2str(pass) ' se ' num2str(clusters) ' shluky začal']);
%% 
[nRows,nCols,nSlices,nSamples]=size(data4D);
TF_4D=repmat(TF,[1 1 1 nSamples]);
signals=data4D(TF_4D);
rSig=reshape(signals,[],nSamples);
nSignals=size(rSig,1);
rSigNorm=rSig./repmat(sum(rSig,2),[1 nSamples]);

[corR,corC,corS] = ind2sub(size(TF),find(TF == 1));

opts    = statset('Display','off');    %show just results
[idx,~] = kmeans(rSigNorm,clusters,'Distance','sqeuclidean','Replicates',15,'Options',opts);
disp(['k-means běh ' num2str(pass) ' dokončen']);

% 'sqeuclidean' (default) | 'cityblock' | 'cosine' | 'correlation' | 'hamming'
colorMap=uint8(zeros(nRows,nCols,nSlices));

for i=1:nSignals
    colorMap(corR(i),corC(i),corS(i))=idx(i);
end

% dataPeak=squeeze(data4D(:,:,:,10));
dataPeak=mean(data4D(:,:,:,3:25),4);

% cmap=hsv(clusters);
% cmap=[0 0 0;cmap];


clusterSignals=zeros(clusters,nSamples);
clusterSignalsNorm=zeros(clusters,nSamples);
nSignalsInCluster=zeros(clusters,1);

signalsAll=cell(clusters,2);

FirstMoment=zeros(clusters,1);
Leg=cell(clusters,1);

for i=1:clusters
    signalsTmp=rSig(idx==i,:);% get signal from original data (with no normalization)
    clusterSignals(i,:)=mean(signalsTmp,1); %mean cluster signal
    signalsTmpNorm=rSigNorm(idx==i,:); % get signal from AUC normalized data
    
    signalsAll{i,1}=signalsTmp;
    signalsAll{i,2}=signalsTmpNorm;
    
    nSignalsInCluster(i,1)=size(signalsTmpNorm,1); %number of voxels included in the cluster
    clusterSignalsNorm(i,:)=mean(signalsTmpNorm,1);    %mean normalized data
    clusterSignalsNorm(i,:)=clusterSignalsNorm(i,:)/sum(clusterSignalsNorm(i,:)); %normalize again
    Leg{i,1}=['Cluster ' num2str(i)];
    if nSignalsInCluster(i,1)<=3
        %         if there are just two signals in cluster -> disqualify from analysis
        if evalFM
            FirstMoment(i,1)=Inf;
        else
            FirstMoment(i,1)=nSamples;
        end
    else
        if evalFM
            FirstMoment(i,1)=sum(clusterSignalsNorm(i,:).*timeAxis);%get 1st moment
        else
            [~,FirstMoment(i,1)]=max(clusterSignalsNorm(i,:)); %get position of maximum (peak)
        end
    end
    %     [~,FirstMoment(i,1)]=max(clusterSignalsNorm(i,:)); %get position of maximum (peak)
end

[~,idxSelClu]=min(FirstMoment);
disp(['Shluk s maximální pozicí je č. ' num2str(idxSelClu)]);

Lowest1StMomentSig=clusterSignals(idxSelClu,:); %Extract signal cluster with lowest 1st moment

signalsAllsel=signalsAll(idxSelClu,:);
%New reduced mask of cluster with the lowest 1st moment
TF_sel=false(size(TF));
corRnew=corR(idx==idxSelClu);
corCnew=corC(idx==idxSelClu);
corSnew=corS(idx==idxSelClu);
nOfNewVoxels=size(corSnew,1);
for i=1:nOfNewVoxels
    TF_sel(corRnew(i),corCnew(i),corSnew(i))=true;
end

colorMapSel=zeros(nRows,nCols,nSlices);
for i=1:nOfNewVoxels
    colorMapSel(corRnew(i),corCnew(i),corSnew(i))=idxSelClu;
end
if visu
    showClusterMap2(dataPeak,colorMap,TF, TF_sel,pass,clusters,nSignalsInCluster,idxSelClu,clusterSignalsNorm,Lowest1StMomentSig);
end

end