function labels = TGTS(f)
%
% TGTS: Tom's Gradient Topographic Sketch
%    labels = TGTS(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
lab = zeros(3,3);
labels = zeros(m,n);
for r = rmin:rmax
  for c = cmin:cmax
%  r
%  c
    label = 0;
    drc = gdir(r,c);
    tfr = r + dr(drc+1);
    tfc = c + dc(drc+1);
    dirf = gdir(tfr,tfc);
    tf = tablef(dir8dif(drc,dirf)+1);
    tlr = r + dr(rem(drc+2,8)+1);
    tlc = c + dc(rem(drc+2,8)+1);
    dirl = gdir(tlr,tlc);
    tl = tablel(dir8dif(drc,dirl)+1);
    tbr = r + dr(rem(drc+4,8)+1);
    tbc = c + dc(rem(drc+4,8)+1);
    dirb = gdir(tbr,tbc);
    tb = tableb(dir8dif(drc,dirb)+1);
    trr = r + dr(rem(drc+6,8)+1);
    trc = c + dc(rem(drc+6,8)+1);
    dirr = gdir(trr,trc);
    tr = tabler(dir8dif(drc,dirr)+1);
    tlbr = tlr + dr(rem(dirl+4,8)+1);
    tlbc = tlc + dc(rem(dirl+4,8)+1);
    dirlb = gdir(tlbr,tlbc);
    tlb = tableb(dir8dif(dirl,dirlb)+1);
    tlfr = tlr + dr(dirl+1);
    tlfc = tlc + dc(dirl+1);
    dirlf = gdir(tlfr,tlfc);
    tlf = tablef(dir8dif(dirl,dirlf)+1);
    trbr = trr + dr(rem(dirr+4,8)+1);
    trbc = trc + dc(rem(dirr+4,8)+1);
    dirrb = gdir(trbr,trbc);
    trb = tableb(dir8dif(dirr,dirrb)+1);
    trfr = trr + dr(dirr+1);
    trfc = trc + dc(dirr+1);
    dirrf = gdir(trfr,trfc);
    trf = tablef(dir8dif(dirr,dirrf)+1);
% alternative version of ridge
    tblr = tbr + dr(rem(dirb+2,8)+1);
    tblc = tbc + dc(rem(dirb+2,8)+1);
    dirbl = gdir(tblr,tblc);
    tbl = tablel(dir8dif(dirb,dirbl)+1);
    tbrr = tbr + dr(rem(dirb+6,8)+1);
    tbrc = tbc + dc(rem(dirb+6,8)+1);
    dirbr = gdir(tbrr,tbrc);
    tbr1 = tabler(dir8dif(dirb,dirbr)+1);
	
    tflr = tfr + dr(rem(dirf+2,8)+1);
    tflc = tfc + dc(rem(dirf+2,8)+1);
    dirfl = gdir(tflr,tflc);
    tfl = tablel(dir8dif(dirf,dirfl)+1);
    tfrr = tfr + dr(rem(dirf+6,8)+1);
    tfrc = tfc + dc(rem(dirf+6,8)+1);
    dirfr = gdir(tfrr,tfrc);
    tfr1 = tabler(dir8dif(dirf,dirfr)+1);
%    lab(1,1) = tlb; lab(1,2) = tl; lab(1,3) = tlf;
%    lab(2,1) = tb; lab(2,2) = 0; lab(2,3) = tf;
%    lab(3,1) = trb; lab(3,2) = tl; lab(3,3) = trf;
%	lab

% label topographic type

%  hillside
    if ((tlb==1)|(tlb==2)|(tlb==3)) & ((tl==12)|(tl==13)|(tl==17)) ...
	 & ((tlf==1)|(tlf==2)|(tlf==3)) & ((tb==1)|(tb==2)|(tb==3))    ...
	 & ((tf==1)|(tf==2)|(tf==3))    & ((trb==1)|(trb==2)|(trb==3)) ...
	 & ((tr==12)|(tr==18)|(tr==17)) & ((trf==1)|(trf==2)|(trf==3)) 
	  label = 1;     % hillside
	end
%  peak
    if ((tlb==1)|(tlb==2)|(tlb==4)) & ((tl==5)|(tl==16)|(tl==17)) ...
	 & ((tlf==7)|(tlf==8)|(tlf==5)) & ((tb==1)|(tb==2)|(tb==3))    ...
	 & ((tf==7)|(tf==6)|(tf==8))    & ((trb==1)|(trb==5)|(trb==3)) ...
	 & ((tr==5)|(tr==16)|(tr==17)) & ((trf==7)|(trf==6)|(trf==4)) 
	  label = label + 2;     % peak
	end
%  pit
    if ((tlb==5)|(tlb==9)|(tlb==10)|(tlb==11)) & ((tl==4)|(tl==13)|(tl==14)) ...
	 & ((tlf==1)|(tlf==2)|(tlf==3)|(tlf==4)) & ((tb==9)|(tb==10)|(tb==11)) ...
	 & ((tf==1)|(tf==2)|(tf==3)) & ((trb==9)|(trb==10)|(trb==11)|(trb==4)) ...
	 & ((tr==4)|(tr==14)|(tr==18)) & ((trf==1)|(trf==2)|(trf==3)|(trf==5)) 
	  label = label + 4;     % pit
	end
%  ridge 1
    if ((tlb==1)|(tlb==2)|(tlb==3)) & ((tl==12)|(tl==13)|(tl==17)) ...
	 & ((tlf==7)|(tlf==8)|(tlf==6)) & ((tb==1)|(tb==2)|(tb==3))    ...
	 & ((tf==7)|(tf==6)|(tf==8))    & ((trb==1)|(trb==2)|(trb==3)) ...
	 & ((tr==12)|(tr==18)|(tr==17)) & ((trf==7)|(trf==6)|(trf==8)) 
	  label = label + 8;     % ridge
	end
%  ridge 2
    if ((tbl==5)|(tbl==16)|(tbl==17)) & ((tl==5)|(tl==16)|(tl==17)) ...
	 & ((tfl==5)|(tfl==16)|(tfl==17)) & ((tb==1)|(tb==2)|(tb==3))    ...
	 & ((tf==1)|(tf==2)|(tf==3))    & ((tbr1==5)|(tbr1==16)|(tbr1==17)) ...
	 & ((tr==5)|(tr==16)|(tr==17)) & ((tfr1==5)|(tfr1==16)|(tfr1==17)) 
	  label = label + 32;     % ridge
	end
%  ravine
    if ((tlb==10)|(tlb==9)|(tlb==11)) & ((tl==12)|(tl==13)|(tl==17)) ...
	 & ((tlf==1)|(tlf==2)|(tlf==3)) & ((tb==9)|(tb==10)|(tb==11))    ...
	 & ((tf==1)|(tf==2)|(tf==3))    & ((trb==9)|(trb==10)|(trb==11)) ...
	 & ((tr==12)|(tr==18)|(tr==17)) & ((trf==1)|(trf==2)|(trf==3)) 
	  label = label + 16;     % ravine
	end
	labels(r,c) = label;
  end
end
