function labels = TGravines(f)
%
% TGravines: Tom's Gradient Topographic Sketch
%    labels = TGravines(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
r
  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);
    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);
    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);
    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);
    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);
    tlfr = tlr + dr(dirl+1);  % left then forward pixel
    tlfc = tlc + dc(dirl+1);
    dirlf = gdir(tlfr,tlfc);
    anglf = theta(tlfr,tlfc);
    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);
    trfr = trr + dr(dirr+1);  % right then forward pixel
    trfc = trc + dc(dirr+1);
    dirrf = gdir(trfr,trfc);
    angrf = theta(trfr,trfc);
    tbbr = tbr + dr(rem(drc+4,8)+1);  % back then back pixel
    tbbc = tbc + dc(rem(drc+4,8)+1);
    dirbb = gdir(tbbr,tbbc);
    angbb = theta(tbbr,tbbc);
    trbbr = trbr + dr(rem(dirr+4,8)+1);  % right then back then back pixel
    trbbc = trbc + dc(rem(dirr+4,8)+1);
	trbb =0;
	if trbbr>0 & trbbr<=m & trbbc>0 & trbbc<=n
      dirrbb = gdir(trbbr,trbbc);
      angrbb = theta(trbbr,trbbc);
	  trbb = 1;
	end;
    tlbbr = tlbr + dr(rem(dirl+4,8)+1); % left then back then back pixel
    tlbbc = tlbc + dc(rem(dirl+4,8)+1);
	tlbb = 0;
	if tlbbr>0 & tlbbr<=m & tlbbc>0 & tlbbc<=n
       dirlbb = gdir(tlbbr,tlbbc);
       anglbb = theta(tlbbr,tlbbc);
	   tlbb = 1;
	end;

    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)==1) ...
         & (dftype(tlr,tlc,angl,tlfr,tlfc,anglf)==1) ...
     & (dftype(trr,trc,angr,trfr,trfc,angrf)==1) ...
         & (dftype(tbr,tbc,angb,r,c,angrc)==2) ...
         & (dftype(tlbr,tlbc,anglb,tlr,tlc,angl)==2) ...
         & (dftype(trbr,trbc,angrb,trr,trc,angr)==2) ...
         & (dstype(r,c,angrc,tlr,tlc,angl)==1) ...
         & (dstype(r,c,angrc,trr,trc,angr)==1)
           label = 1;
        end;

% Type 2:          <- . ->
%                  <- . ->
%                  <- . ->
    if tlbb>0 & trbb>0
	  if (dftype(r,c,angrc,tbbr,tbbc,angbb)==2) ...
	   & (dftype(tlr,tlc,angl,tlbbr,tlbbc,anglbb)==2) ...
	   & (dftype(trr,trc,angr,trbbr,trbbc,angrbb)==2) ...
	   & (dstype(r,c,angrc,tlr,tlc,angl)==1) ...
	   & (dstype(r,c,angrc,trr,trc,angr)==1) ...
	   & (dstype(tbbr,tbbc,angbb,tlbbr,tlbbc,anglbb)==1) ...
	   & (dstype(tbbr,tbbc,angbb,trbbr,trbbc,angrbb)==1)
	    label = 2;
	  end;
	end;
% Type 3:          <- . . ->
%                 <- . . ->
%                 <- . . ->
% Typ 4:          -> . . . <-
%                 -> . . . <-
%                 -> . . . <-
    labels(r,c) = label;
  end
end
