function  [R,sizes] = lyngby_img_regident(V)

% lyngby_img_regindent - Region identification
%
%       function  [R, sizes] = lyngby_img_regident(V)
%
%       Input:  V       Binary image with values in {0,1}
%      
%       Output: R       Labeled image
%               sizes   Number of elements in region
%  
%       This function computes a list of connected regions (clusters)
%       by labelling all pixels by a region number, ie., so-called
%       region identification (coloring or labelling).
%
%       The regions in R is sorted so that the largest cluster has the
%       label=1, the second largest label=2, ... The Non-labeled region
%       is returned with zero.
%
%       The number of regions is max(R) or length(sizes).

% cvs : $Id: lyngby_img_regident.m,v 1.1 1999/05/27 20:41:17 fnielsen Exp $
%       $Revision: 1.1 $


    [rV, cV] = size(V);
    V(1,:) = zeros(size(V(1,:)));

    R = zeros(size(V));     
 
    MA = sum(sum(V)) + rV;
    list = ones(MA,1);        
    hlist = ones(MA,1);     % Equivalence list
    sizes = zeros(MA,1);

    % First pass
    label = 1;
    for j = 2:rV
      if V(j,1) == 1
	% First row
	R(j,1) = label;
      end
      for i = 2:cV
	if V(j,i) == 1 
	  if V(j,i-1) == 1
	    % Same region as previous
	    R(j,i) = label;
	  else
	    % New region
	    label = label + 1;
	    R(j,i) = label;
	    list(label) = label;
	    hlist(label) = label;
	  end
	  if V(j-1,i) == 1
	    % Label collision
	    hlist(label) = hlist(list(R(j-1,i)));
	  end,
	end
      end
      label = label + 1;
    end 
    
    % Second pass
    for j = 2:(label-1),
      hlist(j) = hlist(min(max(1,hlist(j)),MA));
    end
    for j = 1:label
      list(j) = hlist(list(j));
    end
    for j = 1:rV
      for i = 1:cV
	if V(j,i) == 1
	  R(j,i) = list(R(j,i));
	  sizes(R(j,i)) = sizes(R(j,i))+1;
	end
      end
    end

    % Collapse labels
    [sortSizes,iSizes] = sort(-sizes);
    sortSizes = -sortSizes;
    nRegions = sum(sizes~=0);
    for n = 1:nRegions
      R(find(R==iSizes(n))) = n;
    end
    sizes = sortSizes(1:nRegions);









