function labels = TGridges3(f)
%
% TGridges3: Tom's Gradient Topographic Sketch
%    labels = TGridges3(f)

tablef = [1,2,4,6,7,8,5,3];        % forward pixel types
tableb = [1,3,5,9,10,11,4,2];      % backward pixel types
tablel = [12,13,4,14,15,16,5,17];  % left pixel types
tabler = [12,17,5,16,15,14,4,18];  % right pixel types

dr = [0,-1,-1,-1,0,1,1,1];    % offsets to pixel in directions
dc = [1,1,0,-1,-1,-1,0,1];    % 0,1,2,3,4,5,6,7  neighbors

[m,n] = size(f);
rmin = 3; rmax = m-2; cmin = 3; cmax = n-2;  % assumes 3x3 neighborhood
[gc,gr] = gradient(f);
gdir = eightdir(gc,gr);   % convert gradient vector into 8 neighbor dirs
theta = atan2(gr,gc);
labels = zeros(m,n);
for r = rmin:rmax
  for c = cmin:cmax
    label = 0;
    drc = gdir(r,c);     %  pixel under consideration
    angrc = theta(r,c);
    tfr = r + dr(drc+1);  % forward pixel
    tfc = c + dc(drc+1);
    dirf = gdir(tfr,tfc);
%       tf = tablef(dir8dif(drc,dirf)+1);
        angf = theta(tfr,tfc);
        tlr = r + dr(rem(drc+2,8)+1);  % left pixel
        tlc = c + dc(rem(drc+2,8)+1);
        dirl = gdir(tlr,tlc);
        angl = theta(tlr,tlc);
%       tl = tablel(dir8dif(drc,dirl)+1);
        tbr = r + dr(rem(drc+4,8)+1);  % back pixel
        tbc = c + dc(rem(drc+4,8)+1);
        dirb = gdir(tbr,tbc);
        angb = theta(tbr,tbc);
%       tb = tableb(dir8dif(drc,dirb)+1);
        trr = r + dr(rem(drc+6,8)+1);  % right pixel
        trc = c + dc(rem(drc+6,8)+1);
        dirr = gdir(trr,trc);
        angr = theta(trr,trc);
%       tr = tabler(dir8dif(drc,dirr)+1);
        tlbr = tlr + dr(rem(dirl+4,8)+1); % left then back pixel
        tlbc = tlc + dc(rem(dirl+4,8)+1);
        dirlb = gdir(tlbr,tlbc);
        anglb = theta(tlbr,tlbc);
%       tlb = tableb(dir8dif(dirl,dirlb)+1);
        tlfr = tlr + dr(dirl+1);  % left then forward pixel
        tlfc = tlc + dc(dirl+1);
        dirlf = gdir(tlfr,tlfc);
        anglf = theta(tlfr,tlfc);
%       tlf = tablef(dir8dif(dirl,dirlf)+1);
        trbr = trr + dr(rem(dirr+4,8)+1);  % right then back pixel
        trbc = trc + dc(rem(dirr+4,8)+1);
        dirrb = gdir(trbr,trbc);
        angrb = theta(trbr,trbc);
%       trb = tableb(dir8dif(dirr,dirrb)+1);
        trfr = trr + dr(dirr+1);  % right then forward pixel
        trfc = trc + dc(dirr+1);
        dirrf = gdir(trfr,trfc);
        angrf = theta(trfr,trfc);
%       trf = tablef(dir8dif(dirr,dirrf)+1);

        tf3 = 0;
        tf2r = tfr + dr(drc+1);
        tf2c = tfc + dc(drc+1);
        tf3r = tf2r + dr(drc+1);
        tf3c = tf2c + dc(drc+1);
        if tf3r>0 & tf3r<=m & tf3c>0 & tf3c<=n
          dirf3 = gdir(tf3r,tf3c);
          tf3rr = tf3r + dr(rem(dirf3+6,8)+1);
          tf3rc = tf3c + dc(rem(dirf3+6,8)+1);
          if tf3rr>0 & tf3rr<=m & tf3rc>0 & tf3rc<=n
            tf3lr = tf3r + dr(rem(dirf3+2,8)+1);
                tf3lc = tf3c + dc(rem(dirf3+2,8)+1);
                if tf3lr>0 & tf3lr<=m & tf3lc>0 & tf3lc<=n
                  angf3 = theta(tf3r,tf3c);
                  angf3r = theta(tf3rr,tf3rc);
                  angf3l = theta(tf3lr,tf3lc);
                  tf3 = 1;
                end;
          end;
        end;

        tf4 = 0;
        tf4r = tf3r + dr(drc+1);
        tf4c = tf3c + dc(drc+1);
        if tf4r>0 & tf4r<=m & tf4c>0 & tf4c<=n
          dirf4 = gdir(tf4r,tf4c);
          tf4rr = tf4r + dr(rem(dirf4+6,8)+1);
          tf4rc = tf4c + dc(rem(dirf4+6,8)+1);
          if tf4rr>0 & tf4rr<=m & tf4rc>0 & tf4rc<=n
            tf4lr = tf4r + dr(rem(dirf4+2,8)+1);
                tf4lc = tf4c + dc(rem(dirf4+2,8)+1);
                if tf4lr>0 & tf4lr<=m & tf4lc>0 & tf4lc<=n
                  angf4 = theta(tf4r,tf4c);
                  angf4r = theta(tf4rr,tf4rc);
                  angf4l = theta(tf4lr,tf4lc);
                  tf4 = 1;
                end;
          end;
        end;

% label topographic type

% Type 1:          -> -> <-
%                  -> -> <-
%                  -> -> <-

    if (dftype(r,c,angrc,tfr,tfc,angf)==3) ...
         & (dftype(tlr,tlc,angl,tlfr,tlfc,anglf)==3) ...
     & (dftype(trr,trc,angr,trfr,trfc,angrf)==3) ...
         & (dftype(tbr,tbc,angb,r,c,angrc)==1) ...
         & (dftype(tlbr,tlbc,anglb,tlr,tlc,angl)==1) ...
         & (dftype(trbr,trbc,angrb,trr,trc,angr)==1) ...
         & (dstype(r,c,angrc,tlr,tlc,angl)==1) ...
         & (dstype(r,c,angrc,trr,trc,angr)==1)
           label = 1;
        end;

% Type 2:          -> . <-
%                  -> . <-
%                  -> . <-
    if (dftype(tlbr,tlbc,anglb,tlfr,tlfc,anglf)==3) ...
         & (dftype(tbr,tbc,angb,tfr,tfc,angf)==3) ...
     & (dftype(trbr,trbc,angrb,trfr,trfc,angrf)==3) ...
         & (dstype(tbr,tbc,angb,tlbr,tlbc,anglb)==1) ...
         & (dstype(tbr,tbc,angb,trbr,trbc,angrb)==1) ...
         & (dstype(tfr,tfc,angf,tlfr,tlfc,anglf)==1) ...
         & (dstype(tfr,tfc,angf,trfr,trfc,angrf)==1)
           label = 2;
        end;

        offctr = 0;
% Typ 3:          -> . . <-
%                 -> . . <-
%                 -> . . <-
        if tf3 > 0   %  pixels exist 3 forward
          if (dftype(r,c,angrc,tf3r,tf3c,angf3)==3) ...
           & (dftype(tlr,tlc,angl,tf3rr,tf3rc,angf3r)==3) ...
           & (dftype(trr,trc,angr,tf3lr,tf3lc,angf3l)==3) ...
           & (dstype(r,c,angrc,tlr,tlc,angl)==1) ...
           & (dstype(r,c,angrc,trr,trc,angr)==1) ...
           & (dstype(tf3r,tf3c,angf3,tf3rr,tf3rc,angf3r)==1) ...
           & (dstype(tf3r,tf3c,angf3,tf3lr,tf3lc,angf3l)==1)
             label = 3;
             offctr = 1;
          end;
        end;

% Typ 4:          -> . . . <-
%                 -> . . . <-
%                 -> . . . <-
        if tf4 > 0 & offctr==0   %  pixels exist 4 forward
          if (dftype(r,c,angrc,tf4r,tf4c,angf4)==3) ...
           & (dftype(tlr,tlc,angl,tf4rr,tf4rc,angf4r)==3) ...
           & (dftype(trr,trc,angr,tf4lr,tf4lc,angf4l)==3) ...
           & (dstype(r,c,angrc,tlr,tlc,angl)==1) ...
           & (dstype(r,c,angrc,trr,trc,angr)==1) ...
           & (dstype(tf4r,tf4c,angf4,tf4rr,tf4rc,angf4r)==1) ...
           & (dstype(tf4r,tf4c,angf4,tf4lr,tf4lc,angf4l)==1)
             label = 4;
             offctr = 1;
          end;
        end;

        if offctr>0
          labels(tf2r,tf2c) = label;
        else
          labels(r,c) = label;
        end;
  end
end
