function b = bands1(pl,rl,f,fa,fb,fc,fd,fe,cx,ch,cr,rx,rh,rr,r1);
% b = bands1(pl,rl,f,fa,fb,fc,fd,fe,cx,ch,cr,rx,rh,rr,r1);
%
% Creates a band matrix; each row has info on a band
%
% On input:
%           pl: 2D thinned parabola labeled image
%           rl: 1D thinned row parabola labeled image
%           f:  original gray level image
%           fa: a coefficient from 2D parabola fit (on x^2)
%           fb: b coefficient from 2D parabola fit (on y^2)
%           fc: c coefficient from 2D parabola fit (on x)
%           fd: d coefficient from 2D parabola fit (on x)
%           fe: e coefficient from 2D parabola fit (constant)
%           cx: 1D parab fit exists image
%           ch: 1D parab fit height image
%           cx: 1D parab fit radius image
%           cx: 1D parab fit exists image
%           cx: 1D parab fit height image
%           cx: 1D parab fit radius image
%           cx: ridge labels image (1: valley, 2: flat, 4: rising to right,
%                                   8: rising to left, 16 peak)
%
% On output:
%            b: Has one row of 20 elements for each band where:
%
%       col  1: left column terminator
%       col  2: right column terminator
%       col  3: count of left support (number of ridges in row)
%       col  4: count of right support (number of ridges in row)
%       col  5: column extent
%       col  6: gray level average
%       col  7: 2D parab height avg (e coefficient)
%       col  8: avg gray level to right
%       col  9: avg gray level up
%       col 10: avg gray level to left
%       col 11: avg gray level down
%       col 12: avg gray level of surround:  (c8+c9+c10+c11)/4
%       col 13: relative height: c6-c12
%       col 14: 1D ridge label at peak location
%       col 15: 2D parab peak row value
%       col 16: 2D parab peak col value
%       col 17: avg col parab radius value
%       col 18: index of band below this one
%       col 19: index of band above this one
%       col 20: equivalence to another band

[m n] = size(f);
numpks = max(max(pl));
b = zeros(numpks,20);

bindex = 0;
for i1 = 1:m
  for j1 = 1:n
    if pl(i1,j1)>0    % for every pixel labeled in the 2D parabola set
      bindex = bindex+1;
      j2 = j1-1; more = 1; lct = 0; lcol = j1;  % assign it the connected
      rlabel = rl(i1,j1); b(bindex,14) = rlabel;  % ridge label 
      b(bindex,15) = i1; b(bindex,16) = j1;    % set this as its center
      while j2>0 & more>0       % scan left
        more = 0;
        if r1(i1,j2)>0 & rlabel == rl(i1,j2)
          more = 1;
          rtype = r1(i1,j2);                    % set up a first round left
          if rtype==2 | rtype == 4              % and right column boundary
            lct = lct+1;                        % for the band; based on
            lcol = j2;                          % ridge responses in the row
          end
          j2 = j2-1;
        end
      end
      j2 = j1+1; more = 1; rct = 0; rcol = j1;   % scan right
      while j2<n+1 & more>0
        more = 0;
        if r1(i1,j2)>0 & rlabel==rl(i1,j2)
          more = 1;
          rtype = r1(i1,j2);
          if rtype==2 | rtype == 8
            rct = rct+1;
            rcol = j2;
          end
          j2 = j2+1;
        end
      end
      b(bindex,1) = lcol;     % left col of band
      b(bindex,2) = rcol;     % right col of band
      dist = 100;
      for j2 = lcol:rcol    % if there's a 1D row parabola peak, make it the
        d12 = abs(j2-j1);   % center and make it 11 pixels wide
        if rx(i1,j2)>0 & d12<dist
          lcol = max(1,j2-min(5,round(rr(i1,j2))));
          rcol = min(j2+min(5,round(rr(i1,j2))),n);
          b(bindex,1) = lcol;
          b(bindex,2) = rcol;
          dist = d12;
        end
      end
% Should insert code here to calculate a row parabola if no rx exists in band
%-------   close to good code 
% [a,b,c] = qf1d(f(i1,lcol:rcol),rcol-lcol+1);
% jcenter = (rcol-lcol)/2;
% if c<0
%   jpeak = jcenter- (b/(2*c));
%-------   end of close to good code 
% If that doesn't fit well, and no rx in band, look for 16 (peak) in r1
% If that doesn't work, look for 4 to 8 transition in r1
      b(bindex,3) = lct;    % left support of band
      b(bindex,4) = rct;    % right support of band
      b(bindex,5) = rcol-lcol+1;
      b(bindex,16) = (lcol+rcol)/2;
      glsum = 0; htsum =0; radsum = 0;
      for j2 = lcol:rcol
        glsum = glsum+f(i1,j2);
        htsum = htsum+fe(i1,j2);
        radsum = radsum+cr(i1,j2);
      end
      b(bindex,6) = glsum/b(bindex,5);   % gray level avg
      b(bindex,7) = htsum/b(bindex,5);   % row height avg
      b(bindex,17) = radsum/b(bindex,5);   % col parab radius avg
      imin = max(1,i1-1); imax = min(m,i1+1);
      jmin = min(n,rcol+1); jmax = min(n,rcol+2); sum = 0;
      for i2 = imin:imax    % gray level avg to right
        for j2 = jmin:jmax
          sum = sum + f(i2,j2);
        end
      end
      b(bindex,8) = sum/((imax-imin+1)*(jmax-jmin+1));
      imin = max(1,i1-2); imax = max(1,i1-1);
      jmin = lcol; jmax = rcol; sum = 0;
      for i2 = imin:imax    % gray level avg uo
        for j2 = jmin:jmax
          sum = sum + f(i2,j2);
        end
      end
      b(bindex,9) = sum/((imax-imin+1)*(jmax-jmin+1));
      imin = max(1,i1-1); imax = min(m,i1+1);
      jmin = max(1,lcol-2); jmax = max(1,lcol-1); sum = 0;
      for i2 = imin:imax    % gray level avg to left
        for j2 = jmin:jmax
          sum = sum + f(i2,j2);
        end
      end
      b(bindex,10) = sum/((imax-imin+1)*(jmax-jmin+1));
      imin = min(m,i1+1); imax = min(m,i1+2);
      jmin = lcol; jmax = rcol; sum = 0;
      for i2 = imin:imax    % gray level avg down
        for j2 = jmin:jmax
          sum = sum + f(i2,j2);
        end
      end
      b(bindex,11) = sum/((imax-imin+1)*(jmax-jmin+1));  
      b(bindex,12) = (b(bindex,8)+b(bindex,9)+ ...
         b(bindex,10)+b(bindex,11))/4;  % gray level avg around
      b(bindex,13) = b(bindex,6) - b(bindex,12);  % relative height
    end
  end
end
