function cs = kmeans(ps, k)
minpt = min(ps');
maxpt = max(ps');
cs = (minpt'*ones(1,k)) + (maxpt-minpt)'*ones(1,k).*rand(size(ps,1),k);
for j = 1:10
        figure(107);
    plot(ps(1,:), ps(2,:), 'ko');
    hold on
    plot(cs(1,:), cs(2,:), 'ro', 'LineWidth',4);
    hold off
    keyboard
    D = [];
    for i = 1:k
        c = cs(:,i);
        Di = sum((c*ones(1,length(ps))-ps).^2);
        D = [D; Di];
    end
    cs = [];
    for i = 1:k
        mems = D(i,:) < min(D([1:(i-1),(i+1):k], :));
        cs = [cs, mean((ps(:,find(mems)))')'];
    end
end
    figure(107);
    plot(ps(1,:), ps(2,:), 'ko');
    hold on
    plot(cs(1,:), cs(2,:), 'ro', 'LineWidth',4);
    hold off
    keyboard