%Doug Van Nort  MUMT 611  Final Project: Timbre Recognition via HMM's


function [c,mu] = kmeans(x,K,mu)
% kmeans : K-means clustering
% [c,mu] = kmeans(x,K)
%	x  - d*n samples
%	K  - number of clusters wanted
%	mu	- optional d*K initial guess for cluster centroids
% returns:
%	c  - 1*n calculated membership vector where c(j) \in 1..K
%	mu	- d*K cluster centroids
%
%This algorithm was originally taken from Sing-Tze Bow, 'Pattern Recognition'
% and is based on an implementation by Frank Dellaert Copyright (c) 1995-2001
% All rights Reserved

%------------------------------------------------------------------------
% Check arguments
%------------------------------------------------------------------------
if ndims(x)~=2, error('x must be a matrix'); end
[d,n] = size(x);
if ~isequal(size(K),[1 1]), error('K must be a scalar');end
K=round(K);
if (K<1 | K>d), error(sprintf('K must be at least one and at most n=%d',d)); end

%------------------------------------------------------------------------
% Initialize variables
%------------------------------------------------------------------------
oldmu = Inf*ones(K,d);
c = zeros(1,d);
D = zeros(K,n);

%------------------------------------------------------------------------
% step 1: Arbitrarily choose K samples as the initial cluster centers
%------------------------------------------------------------------------
if nargin<3
   p=randperm(d);
   mu = x(p(1:K),:);
end

while(1),
   %------------------------------------------------------------------------
   % step 2: distribute the samples x to the chosen cluster domains
   %         based on which cluster center is nearest
   %------------------------------------------------------------------------
   
   %changed the distance calculations and matrix handling - doug van nort 4/04
   
   for j=1:K,                        % for every cluster        
      center = mu(j,:);              % get cluster center
      size(center);
      if ~isequal(center,oldmu(j,:)) % has it moved ? 
          
          for l=1:d
         D(j,l) = dist2(x(l,:),center); % calculate sqrDist from x to center
                    
            end
      end
   end
   oldmu = mu;
   
   % find minimum
   [Dmin,index] = min(D);
   moved = sum(index~=c);
   %fprintf(2,'moved = %d\n',moved);
   c = index;
   
   %------------------------------------------------------------------------
   % step 3: Update the cluster centers
   %------------------------------------------------------------------------
   for i=1:K
      ci=find(c==i);
      mu(i,:)=mean(x(ci,:),1);
   end
   
   %------------------------------------------------------------------------
   % step 4: Check convergence
   %------------------------------------------------------------------------
   if (moved==0), break, end
   
   %------------------------------------------------------------------------
end % while(1)


