function L = cclabel(BW,mode)
%CCLABEL Label connected components.
%	CCLABEL(BW), where BW is a matrix containing 0's and 1's,
%	returns a matrix the same size containing labels for the
%	8-connected components in BW.
%
%	CCLABEL(BW,4) returns a label matrix for the 4-connected
%	components in B.
%
%	Example:
%	>> I = [1 1 0 1 1 0 1
%	        1 1 0 1 1 0 1
%	        0 0 1 0 0 0 1];
%	>> cclabel(I)
%	ans =
%	     1     1     0     1     1     0     2
%	     1     1     0     1     1     0     2
%	     0     0     1     0     0     0     2
%	>> cclabel(I,4)
%	ans =
%	     1     1     0     3     3     0     4
%	     1     1     0     3     3     0     4
%	     0     0     2     0     0     0     4

%	Steven L. Eddins, August 1995
%	Copyright (c) 1993-1996 by The MathWorks, Inc.
%	$Revision: 1.2 $  $Date: 1995/09/06 15:06:13 $

if (nargin < 2)
  mode = 8;
end

[M,N] = size(BW);

%
% Compute run-length encoding
%
rowZeros = zeros(1,size(BW,2));
BW = [rowZeros ; BW ; rowZeros];
d = diff(BW);
[sr,sc] = find(d == 1);  % start row and column indices
[er,ec] = find(d == -1); % end row and column indices
numRuns = length(sr);
if (numRuns == 0)
  runs = [];
else
  runs = [sc sr (er-1)];
end

% 
% First labeling pass
%
labels = zeros(numRuns,1);
equivTable = [];
nextEquivTableIdx = 1;
currentLabel = 1;
for k = 1:numRuns
  column = runs(k,1);
  rowStart = runs(k,2);
  rowEnd = runs(k,3);
  idx = find(runs(:,1) == column - 1); % find runs on previous column
  if (isempty(idx))
    labels(k) = currentLabel;
    currentLabel = currentLabel + 1;
  else
    if (mode == 8)
      % Find eight-connected objects
      overlaps = find(((rowEnd >= (runs(idx,2)-1)) & ...
	  (rowStart <= (runs(idx,3)+1))));
    else
      overlaps = find(((rowEnd >= runs(idx,2)) & ...
	  (rowStart <= runs(idx,3))));
    end
    if (isempty(overlaps))
      labels(k) = currentLabel;
      currentLabel = currentLabel + 1;
    elseif (length(overlaps) == 1)
      labels(k) = labels(idx(overlaps));
    else
      labels(k) = labels(idx(overlaps(1)));
      for n = 2:length(overlaps)
	if (labels(idx(overlaps(1))) ~= labels(idx(overlaps(n))))
	  if (nextEquivTableIdx > (size(equivTable,1)+1))
	    equivTable = [equivTable ; zeros(size(equivTable))];
	  end
	  equivTable(nextEquivTableIdx:nextEquivTableIdx+1,:) =  ...
	      [labels(idx(overlaps(1))) labels(idx(overlaps(n))) ;
	       labels(idx(overlaps(n))) labels(idx(overlaps(1))) ]; 
	  nextEquivTableIdx = nextEquivTableIdx + 2;
	end
      end
    end
  end
end
numLabels = currentLabel - 1;
equivTable = [equivTable(1:nextEquivTableIdx-1,:);
  repmat((1:numLabels)', 1, 2)];


%
% Equivalence class resolution
%
% This next section is borrowed from components.m.
% John Gilbert, Xerox PARC, 8 June 1992.
% Copyright (c) 1990-1995 by Xerox Corporation.  All rights reserved.

A = sparse(equivTable(:,1), equivTable(:,2), ones(size(equivTable,1),1), ...
    numLabels, numLabels);

[p,p,r,r] = dmperm(A);

sizes = diff(r);                        % Sizes of components, in vertices.
numObjs = length(sizes);      % Number of components.
 
% Now compute an array "blocks" that maps vertices of A to components;
% First, it will map vertices of A(p,p) to components...
 
blocks = zeros(1,numLabels);
blocks(r(1:numObjs)) = ones(1,numObjs);
blocks = cumsum(blocks);
 
% Second, permute it so it maps vertices of A to components.
 
blocks(p) = blocks;

%
% Finally!  Fill in each run in the output image with a label.
%

L = zeros(M,N);
for k = 1:numRuns
  column = runs(k,1);
  rowStart = runs(k,2);
  rowEnd = runs(k,3);
  value = blocks(labels(k));
  L(rowStart:rowEnd,column) = value(ones(rowEnd-rowStart+1,1),1);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1997-1999. The Regents of the University
% of California (Regents). All Rights Reserved.
% 
% Permission to use, copy, modify, and distribute this software and its
% documentation for educational, research, and not-for-profit purposes,
% without fee and without a signed licensing agreement, is hereby
% granted, provided that the above copyright notice, this paragraph and
% the following two paragraphs appear in all copies, modifications, and
% distributions. Contact The Office of Technology Licensing, UC Berkeley,
% 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, (510)
% 643-7201, for commercial licensing opportunities. Created by
% Chad Carson, Serge Belongie, and Jitendra Malik, Department of Electrical
% Engineering and Computer Sciences, University of California, Berkeley.
% For technical information, contact carson@eecs.berkeley.edu.
% 
% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF
% ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION
% TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

