function [segs,clusters] = CV_kmeans(im,f_dist_func,k)
%
% CV_kmeans - K-means clustering
% On input:
%     im (mxnxp array): input image
%     f_dist_func (string): name of distance function
%     k (int): number of clusters
% On output:
%     segs (mxn array): segmented image (distinct integer for each class)
%     clusters (kxp array): center value (mean) for each cluster
% Call:
%    [s,c] = CV_kmeans(house1,CV_Euclidean,5);
% Author:
%  Tom Henderson
%  UU
%  Fall 2004
%  Modified: Fall 2012
%

MAX_D = 2;

[num_rows,num_cols,num_feat] = size(im);
clusters = zeros(k,num_feat);
pts = [floor(num_rows*rand(k,1))+1,floor(num_cols*rand(k,1))+1];
for p = 1:k
    clusters(p,:) = im(pts(p,1),pts(p,2),:);
end
unchanged = 0;
counts = zeros(k,1);
sums = zeros(k,num_feat);

while (unchanged==0)
    for r = 1:num_rows
        for c = 1:num_cols
            dist = zeros(k,1);
            for cl = 1:k
                im_vec = im(r,c,:);
                im_vec = reshape(im_vec,1,num_feat);
                dist(cl) = norm(clusters(cl,:)-im_vec);
            end
            [val,index] = min(dist);
            sums(index,:) = sums(index,:) + im_vec;
            counts(index) = counts(index)+1;
        end
    end
    for cl = 1:k
        if (counts(cl)==0)
            counts(cl) = 1;
            sums(cl,:) = im(floor(num_rows*rand)+1,...
                floor(num_cols*rand)+1,:);
        end
    end
    for cl = 1:k
        new_clusters(cl,:) = sums(cl,:)/counts(cl);
    end
    d = sum(norm(new_clusters-clusters));
    if (d<MAX_D)
        unchanged = 1;
    else
        clusters = new_clusters;
    end
end

segs = zeros(num_rows,num_cols);
for r = 1:num_rows
    for c = 1:num_cols
        dist = zeros(k,1);
        for cl = 1:k
            im_vec = im(r,c,:);
            im_vec = reshape(im_vec,1,num_feat);
            dist(cl) = norm(clusters(cl,:)-im_vec);
        end
        [val,index] = min(dist);
        segs(r,c) = index;
    end
end
