function [L,C,iters] = kmeans3(X,k,C,distance,convTh,maxIter,tol,spamitmax)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% K-means, K-medians, K-spatialmedians clustering for given C
% Joonas H, 26.6.2017
% Modified, 27.6.2017 (Tommi K.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

iters = 0;
converged = 0;

if(strcmp(distance,'seuclidean'))
    [~,L] = min(pdist2(X,C),[],2);
    % K-means
    while (~converged && (iters <maxIter))
        Cref = C;
        for i = 1:k
            labels = L==i;
            if (isempty(labels))
                disp('Empty cluster. Returning. ');
                iters = -1;
                return
            elseif (sum(labels) == 1)
                C(i,:) = X(labels,:);
            else
                C(i,:) = sum(X(labels,:),1)/sum(labels);
            end
        end
        [~,L] = min(pdist2(X,C),[],2);
        iters = iters + 1;
        maxCmove = max(min(pdist2(C,Cref).^2,[],2));
        if(maxCmove <= convTh || k==1)
            converged = 1;
        end
    end

elseif(strcmp(distance,'cityblock'))
    [~,L] = min(pdist2(X,C,distance),[],2);
    % K-medians
    while (~converged && (iters <maxIter))
        Cref = C;
        for i = 1:k
            labels = L==i;
            if (isempty(labels))
                disp('Empty cluster. Returning. ');
                iters = -1;
                return
            elseif (sum(labels) == 1)
                C(i,:) = X(labels,:);
            else
                C(i,:) = median(X(labels,:));
            end
        end
        [~,L] = min(pdist2(X,C,distance),[],2);
        iters = iters + 1;
        maxCmove = max(min(pdist2(C,Cref,distance),[],2));
        if(maxCmove <= convTh || k==1)
            converged = 1;
        end
    end
    
elseif(strcmp(distance,'euclidean'))
    [~,L] = min(pdist2(X,C,distance),[],2);
    % K-spatialmedians
    while (~converged && (iters <maxIter))
        Cref = C;
        for i = 1:k
            labels = find(L==i);
            if (isempty(labels))
                disp('Empty cluster. Returning. ');
                iters = -1;
                return
            elseif (sum(labels) == 1)
                C(i,:) = X(labels,:);
            else
                [C(i,:),~] = SpatMed3(X(labels,:),tol,spamitmax);
            end
        end
        [~,L] = min(pdist2(X,C,distance),[],2);
        iters = iters + 1;
        maxCmove = max(min(pdist2(C,Cref,distance),[],2));
        if(maxCmove <= convTh || k==1)
            converged = 1;
        end
    end
end
