function l = label4(a)
%
% label4: assigns unique sequential labels to 4-connected components in array
% On input:
%     a: 2D array
% On output:
%     l: labeled image (unique, consecutive labels on 4-connected sets)
% Call:
%     l = label4(a)  
%

[m n] = size(a);
l = zeros(m,n);

indexl = 0;
for i1 = 1:m
  for j1 = 1:n
    if a(i1,j1) > 0
      iup = i1-1; lup = 0;
      jleft = j1-1; lleft = 0;
      if iup>0      % if up pixel labeled
        lup = l(iup,j1);
      end
      if jleft>0    % if left pixel labeled
        lleft = l(i1,jleft);
      end
      if lup==0 & lleft==0    % no labeled neighbor
        indexl = indexl+1;    % add new label
%        eq(indexl,indexl) = 1;
        l(i1,j1) = indexl;
      elseif lup>0 & lleft==0   % only up neighbor
        l(i1,j1) = lup;
      elseif lup==0 & lleft>0   % only left neighbor
        l(i1,j1) = lleft;
      elseif lup==lleft         % both same label
        l(i1,j1) = lup;
      elseif lup~=lleft         % labels differ
        lmin = min(lup,lleft);
        lmax = max(lup,lleft);
        l(i1,j1) = lmin;
%        eq(lleft,lup) = 1;
%        eq(lup,lleft) = 1;
%        eq(lmax,lmax) = lup;
      end
    end
  end       % Initial labeling: non-unique, not compacted
end        % % Builds the equivalence table, eq

lmax = 0;
for i1 = 1:m
  for j1 = 1:n
    l1 = l(i1,j1); lup = 0; lleft = 0; ldown = 0; lright = 0;
    if i1>1
      lup = l(i1-1,j1);
    end
    if j1>1
      lleft = l(i1,j1-1);
    end
    if i1<m
      ldown = l(i1+1,j1);
    end
    if j1<n
      lright = l(i1,j1+1);
    end
    lnew = l1;
    if lup>0 & lup<lnew
      lnew = lup;
    end
    if lleft>0 & lleft<lnew
      lnew = lleft;
    end
    if ldown>0 & ldown<lnew
      lnew = ldown;
    end
    if lright>0 & lright<lnew
      lnew = lright;
    end
    if lnew ~= l1
      i2 = i1;    % do same row
      for j2 = 1:n
        if l(i2,j2)==l1
          l(i2,j2) = lnew;
        end
      end
      i2 = i1-1; count = 1;  % do previous rows
      while i2>0 & count>0
        count = 0;
        for j2 = 1:n
          if l(i2,j2)==l1
            l(i2,j2) = lnew;
            count = count+1;
          end
        end
        i2 = i2-1;
      end
      i2 = i1+1; count = 1;  % do following rows
      while i2<m+1 & count>0
        count = 0;
        for j2 = 1:n
          if l(i2,j2)==l1
            l(i2,j2) = lnew;
            count = count+1;
          end
        end
        i2 = i2+1;
      end
    end
    if lnew>lmax 
      lmax = lnew;
    end
  end
end

logl = zeros(lmax+1,1);
for i1 = 1:lmax       % index label table
  logl(i1) = i1;
end

lhist = zeros(lmax+1,1);     % calculate histogram for labels
for i1 = 1:m
  for j1 = 1:n
    lval = l(i1,j1);
    if (lval>0) & (lval<=lmax)
        lhist(lval)= lhist(lval)+1;
    end
  end
end

for i1 = 1:lmax
  while (lhist(i1) == 0) & (sum(lhist(i1:lmax)) ~= 0)
    for i2 = i1:lmax
      lhist(i2) = lhist(i2+1);
      logl(i2) = logl(i2+1);
    end
  end
end

for i1 = 1:lmax
  lval = logl(i1); 
  if lval>0
    logil(lval) = i1;
  end
end

for i1 = 1:m
  for j1 = 1:n
    lval = l(i1,j1);
    if lval>0
      l(i1,j1) = logil(lval);
    end
  end
end
