function result = hk(input) 
% function result = hk(input)
%
% Apply the Hoshen-Kopelman Algorithm for Cluster Identification
%
% Based on http://splorg.org/people/tobin/projects/hoshenkopelman/
%
% Note: there's probably a better way to do this in MATLAB, since this code
% does not utilize vectorization at all.
%
% Note: I had some trouble with this as I am assuming that all variables are
% integers.  I'm not sure how to tell Matlab to actually use integers instead 
% of floating point numbers. 
%
% 2002-06-16 Tobin Fricke <tobin@splorg.org> Weizmann Institute of Science
% $Id: hk.m,v 1.1 2002/06/16 17:10:33 tobin Exp tobin $

matrix = round(input);

[m,n] = size(matrix);

% labels = zeros(m*n/2,1);  % allocate space for the labels

n_labels = 0;             % number of labels used so far

labels=[];

%  /* scan the matrix */
  
for i=1:m,
  for j=1:n
    if matrix(i,j) == 1
	
      % check the neighbors of this cell, up and to the left

      up = 0;
      left = 0;
    
      if (i > 1), 
        up = matrix(i-1,j);
      end
    
      if (j > 1)
        left = matrix(i,j-1);
      end
    
      switch ((up>0)  + (left>0))
        
      case 0 %  new cluster 
        n_labels = n_labels + 1;           % make the new label number
	    labels( n_labels ) = n_labels;     % no known alias yet (label[x]=x)
	    matrix(i,j) = n_labels;            % apply this label to this position 

	  case 1 %  same cluster 
	    matrix(i,j) = max([up,left]);      % whichever is nonzero is already labelled; use that label

	  case 2 % two clusters 
       
        while labels(up) < up,      % this scheme could probably be improved
          % dereference the label until it points to itself
          up = labels(up);
        end
      
	    while labels(left) < left,
          left = labels(left);
        end
        
	    label_min = min([up,left]);  % this is the label that we will use
	    label_max = max([up,left]);  % note that this is an alias for label_min  

	    matrix(i,j) = label_min;
	    labels(label_max) = label_min;
      otherwise
        error('error!')
      end % switch(up+left)
    end % if
  end % for
end % for

%  Renumber the labels so that they're continuous, and eliminate cluster aliases

disp('relabeling'),
j = 1;
for i=1:n_labels
    if (labels(i) == i)
	    labels(i) = j;
        j = j+1;
    else
	    labels(i) = labels(labels(i));
    end
end

%  apply the relabeling to the matrix 

for i=1:m
    for j=1:n
        if not (matrix(i,j) == 0)
          matrix(i,j) = labels(matrix(i,j)+0);
        end
    end
end

result = matrix;


