function J = KMeansBW(I, k, partnum)
I = double(I);
switch partnum
    case 1
        % Here k will actually be a vector of cluster centers.
        J = Part1(I,k);
    case 2
        % Also k will be the cluster centers.
        as = Part1(I,k);
        J = result_image(k, as);
        figure(1);
        colormap(gray);
        imshow(uint8(round(J)));
    case 3
        % In this case k will contain the assignments that we computed for
        % each cluster.  J, the output will contain the new cluster
        % centers.
        J = compute_centers(k,I);
    case 4
        J = Part4(I,k);
end

function J = Part1(I, Cs)
% for part one, the input will be an image and a representation of the clusters centers.
% We'll return a representation of which pixel belongs to which cluster.
D = distance_to_centers(Cs,I);
J = assign_centers(D);

function J = Part4(I,k)
Cs = random_centers(k);
cont = 1;
while cont
    D = distance_to_centers(Cs, I);
    as = assign_centers(D);
    newCs = compute_centers(as, I);
    cont = not_converged(Cs, newCs);
    Cs = newCs
end
J = result_image(Cs, as);

function Cs = random_centers(k)
Cs = 1+ floor(255*rand(1,k));

function D = distance_to_centers(Cs,I)
% Cs is kx3, I is nx3.  We'll return an nxk matrix with distances to
% centers.
D = [];
for c = Cs
    D = cat(3, D, (I-c).^2);
end

function as = assign_centers(D)
k = size(D,3);
as = [];
for i = 1:k
    d = D(:,:,i);
    mems = ones(size(D(:,:,1)));
    for j = 1:k
        if j ~= i
            mems = mems & (d < D(:,:,j));
        end
    end
    as = cat(3, as, mems);
end

function Cs = compute_centers(as, I)
Cs = [];
for i = 1:size(as,3)
    Cs = [Cs, mean(I(logical(as(:,:,i))))];
end

function o = not_converged(Cs1, Cs2)
o = any(abs(Cs1-Cs2) > .0001);

function J = result_image(Cs, as)
J = zeros(size(as(:,:,1)));
for i = 1:size(as,3)
    J = J + as(:,:,i)*Cs(i);
end