function [aclus,centers,clustercounts,iters] = ...
           mycluster(asdata,centers,tol,maxiters)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [aclus,centers,clustercounts,iters] = ...     %
%            mycluster(asdata,centers,tol,maxiters)      %
% Cluster data using the k-means algorithm               %
%                                                        %
% Input:                                                 %
% asdata = mxpxq array of values                         %
% centers = qxk array of initial guesses at centers      %
%           of the k clusters                            %
% tol = tolerance on movement of centers to define       %
%       convergence of the iteration                     %
% maxiters = maximum number of iterations allowed        %
% Output:                                                %
% aclus = mxp array of the cluster number for each value %
% centers = qxk array of the centers of the k clusters   %
% clustercounts = 1xk array of number of points in each  %
%                 cluster                                %
% iters = number of iterations                           %
% DPO 03/03                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[m,p,q] = size(asdata);
[q,k] = size(centers);

oldcenters = centers - 5*tol;
iters = 0;

while (norm(oldcenters-centers,'fro') > tol) & (iters < maxiters)

   iters = iters + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Assign each point to the nearest cluster.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   [aclus,clustercounts] = map_to_cluster(asdata,centers);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Update the cluster centers.                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   ssum = zeros(q,k);

   for i=1:m,
   for j=1:p,
     kk = aclus(i,j);
     ssum(:,kk) = ssum(:,kk) + reshape(asdata(i,j,:),q,1);
   end
   end
   
   oldcenters = centers;
   for kk=1:k,
      if (clustercounts(kk) > 0)
        centers(:,kk) = ssum(:,kk)/ clustercounts(kk);
      else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Reinitialize a lost center to a data point1.  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ii = min(floor(rand*m+1),m);
        jj = min(floor(rand*p+1),p);
        centers(:,kk) = reshape(asdata(ii,jj,:),q,1);
      end
   end

end %while