%
% K-means clustering algorithm
%
%  [ xhat, alloc ] = kmeans( data, kmean [, means]  )
%
%   xhat  = cluster centres
%   alloc = allocation of data to clusters
%   data  = data
%   kmean = number of clusters
%   means = optional starting means
%
% bob newell, april 1997
%
%--------------------------------------------------------
function [ xhat, alloc ] = kmeans( data, kmean, means )
%-----------------------
% initialise
hh = kmean;
[ kk, n ] = size( data );
%----------------------
% step 1 - given or random means
if nargin > 2,
  xhat = means;
else
  notuniq = 1;
  while notuniq,
    sel = round( rand( 1, hh ) * ( kk - 1 ) ) + 1;
    notuniq = 0;
    for i = 2:hh,
      for j = 1:i-1,
        if sel(i) == sel(j), notuniq = 1; end;
      end;
    end;
  end;
  xhat = data( sel, : );
end;
%-----------------------
% convergence check
crit = 1;
while abs( crit ) > 0,
%-----------------------
% step 2 - allocate data
  for i = 1:kk,
    dmin = 1.0e30; jmin = 0;
    for j = 1:hh,
      d = norm( data(i,:) - xhat(j,:) );
      if d < dmin,
        dmin = d; jmin = j;
      end;
    end; 
    alloc(i) = jmin;
  end;
%-----------------------------
% step 3 - calculate new means
  crit = 0;
  dtot = 0;
  for i = 1:hh,
    sel = find( alloc == i );
    for j = 1:length( sel ),
      dtot = dtot + norm( data(sel(j),:) - xhat(i,:) );
    end;
    crit = crit + norm( xhat(i,:) - mean( data( sel, : ) ) );
    xhat(i,:) = mean( data( sel, : ) );
  end;
%  [ crit dtot ]
end;
%----------------------------
for i = 1:hh,
  sel = find( alloc == i );
%  disp( [ 'cluster ' num2str(i) ' = ' num2str( length( sel ) ) ] );
%  disp( xhat(i,:) );
end;
