%--------------------------------------------------------------------------
% Created: 04/24/2014 Simon Nguyen
%
% Revision History:
%
% Purpose: This function takes an input image and integer k that produces k
% regions on the image.
%
% Variables:
%   Inputs:
%       im = input image
%       k  = integer value for number of regions
%   Outputs:
%       img = output image with boundaries drawn
%
% function img = mykmeans(im, k)
%--------------------------------------------------------------------------

function img = mykmeans(im, k)
% Choose k random points from within the image
[r, c, ~] = size(im);
pts = [randi(r, 1, k); randi(c, 1, k)];

% Find the colors at the chosen points
m = zeros(size(pts)); m = [m; m(1,:)];
for kk = 1 : k
    pt = pts(:,kk);
    m(:,kk) = [im(pt(1), pt(2), 1);
               im(pt(1), pt(2), 2);
               im(pt(1), pt(2), 3)];
end

% Create a temporary image containing the clusters.
kcluster = uint8(zeros(r,c));

% Limit iterations to 100
ic = 1;
while ic < 100

    % After each iteration, values within the regions are averaged and
    % chosen as the next set of m-values
    if ic > 1
        prevkm = m;

        for kk = 1 : k
            curkc = uint8(kcluster == kk);
            curr = curkc.*im(:,:,1);
            curg = curkc.*im(:,:,2);
            curb = curkc.*im(:,:,3);

            nvals = sum(kcluster(:) == kk);
            m(:,kk) = [sum(curr(:))/nvals;
                       sum(curg(:))/nvals;
                       sum(curb(:))/nvals];
        end
    end

    % Iterate through each pixel of the image and find its color distance
    % with respect to each of the m-values, assigning the pixel to the
    % corresponding group.
    for cc = 1 : c
        for rr = 1 : r
            codiff = zeros(1, length(m(1,:)));

            for mm = 1 : length(m(1,:))
                km = m(:,mm);
                curPixel = double([im(rr, cc, 1);
                            im(rr, cc, 2);
                            im(rr, cc, 3)]);

                codiff(:,mm) = sqrt( (km(1) - curPixel(1))^2 + ...
                                    (km(2) - curPixel(2))^2 + ...
                                    (km(3) - curPixel(3))^2);
            end

            codmin = find(codiff == min(codiff), 1);
            kcluster(rr, cc) = codmin;

        end
    end

    % After each iteration, check if the difference between the previous
    % clustering is similar enough to the current iteration.
    if ic > 1
        kmDiff = abs((prevkm-m)./m);

        if sqrt(kmDiff(1,:).^2 + kmDiff(2,:).^2 + kmDiff(3,:).^2) < .05
            break;
        end
    end

    ic = ic + 1;

end

% Create the boundaries for each region
kcvals = unique(kcluster);
newim = im;
tempim = label2rgb(kcluster);

for kk = 1 : length(kcvals)

    curkc = kcluster == kk;
    nkc = (curkc - imerode(curkc, strel('square', 2)));
    [r, c] = find(nkc == 1);

    for rr = 1 : length(r)
        aa = r(rr);
        bb = c(rr);
        newim(aa, bb, 1) = tempim(aa, bb, 1);
        newim(aa, bb, 2) = tempim(aa, bb, 2);
        newim(aa, bb, 3) = tempim(aa, bb, 3);
    end
end

img = newim;
end