function [lb, lanes] = lanes2(bands)
% lanes = lanes2(bands)
%   bands(n1,19) has the number of the band above band n1
[m n] = size(bands);
lb = bands;
minrow = min(bands(:,15));  % should exclude bands(:,20) ~= :
maxrow = max(bands(:,15));
nottop = zeros(m,1); below = zeros(m,1);
for n1 = 1:m   % find tops of lanes (no band above)
  ab = bands(n1,19);   % and add pointers to band below
  if ab>0 & bands(n1,20)>0
    nottop(n1) = 1;
    below(ab) = n1;
    lb(ab,18) = n1;
  end
end
numlanes = 0; tops = [];
for n1 = 1:m
  if nottop(n1)==0 & bands(n1,20)>0
    tops = [tops,n1];    % list of tops of lanes
    numlanes = numlanes+1;
  end
end
lanes = zeros(numlanes,8);
for l1 = 1:numlanes    % get lane info
  n1 = tops(l1); count = 1;    % first band in lane
  lanes(l1,1) = bands(n1,15); lanes(l1,3) = bands(n1,15);
  lanes(l1,2) = bands(n1,16); lanes(l1,4) = bands(n1,16);
  glsum = bands(n1,6);
  n1 = below(n1);   % follow pointers down lane
  while n1>0
    lanes(l1,3) = bands(n1,15);  % set end of lane row and col
    lanes(l1,4) = bands(n1,16);
    glsum = glsum + bands(n1,6);
    count = count+1;
    n1 = below(n1);
  end
  lanes(l1,9) = count;
  lanes(l1,5) = glsum/count;  % avg gray level of bands in lane
end
for l1 = 1:numlanes   % compute rho, theta of line through lane
  x1 = lanes(l1,4); y1 = -lanes(l1,3); x2 = lanes(l1,2); y2 = -lanes(l1,1);
  lanes(l1,7) = atan2(x1-x2,y2-y1);  % translate to origin and rotate 90 deg
  if (x1~=x2) | (y1~=y2)  % min distance from origin to line
    v1 = sqrt(x1^2 + y1^2);
    v2 = sqrt((x2-x1)^2 + (y2-y1)^2);
    phi = acos((-x1*(x2-x1)-y1*(y2-y1))/(v1*v2));
    lanes(l1,6) = v1*sin(phi);
  else
    lanes(l1,6) = x1;  % should set to average of nearest lanes
  end
end
