function means = kmeans(data,n)
% function means = kmeans(data,n)
%
% returns a set of n means (each a row of the means matrix) estimated
% by running the k-means algorithm on the data

mins = min(data,[],1);
maxs = max(data,[],1);
npts = size(data,1);
d = size(data,2);
delta = maxs-mins;
% distribute the means uniformly across the axis-aligned bounding box for
% the data
means = rand(n,d).*repmat(delta,n,1) + repmat(mins,n,1);
oldmeans = means - 1;
% itterate until the means stop changing
while(sum(sum(oldmeans ~= means))),
	oldmeans = means;

% uncomment the following to see the code run ... its neat!
	if (d==2),
		plot(data(:,1),data(:,2),'gx');
		hold on;
		plot(means(:,1),means(:,2),'ro');
		hold off;
		pause;
	end;

	s = nearest(data,means,eye(n,n));

	% this part is slightly subtle... to insure that we have n means,
	% we have to replace any center that doesn't capture any point
	% to help get this center restarted, we place it right on top of
	% a random point.  In case two replaced ones end up on the same
	% random point, we must repeat until all centers have captured a
	% point
	okay = sum(s)>0;
	while(sum(okay)<n), % not all of the means have captured a point
		for i=1:n,
			if (~okay(i)),
				means(i,:) = data(ceil(rand(1)*npts),:);
			end;
		end;
		s = nearest(data,means,eye(n,n));
		okay = sum(s)>0;
	end;

	for i=1:n,
		means(i,:) = mean(data(s(:,i)>0,:));
	end;
end;
