function labels = TGrav(f,fh,d)
%
% TGrav: Tom's Gradient Topographic Sketch for ridges
% On input:
%     f:   image function (2D)
%    fh:   labeled hill function (2D): 1 if hillside pixel
%     d:   scale of ridge
% On output:
%     labels:   0  unknown
%               1  ridge
% Calling sequence:
%    labels = TGrav(f,fh,d)

[m,n] = size(f);
d1p5 = 1.5*d;
rmin = ceil(d1p5+1); rmax = floor(m-d1p5); cmin = rmin; cmax = floor(n-d1p5);
[gc,gr] = gradient(f);
theta = atan2(gr,gc);
labels = zeros(m,n);
r45 = ceil(d1p5*sin(pi/4));;
step = (pi/4)/r45;
for r = rmin:rmax
r
  for c = cmin:cmax
    for angle = 0:step:pi
	  delr = ceil(d1p5*sin(angle));  delc = ceil(d1p5*cos(angle));
	  r1 = r-delr;  c1 = c+delc;
	  r2 = r+delr;  c2 = c-delc;
	  if (r1>0) & (r1<=m) & (c1>0) & (c1<=n) ...
	   & (r2>0) & (r2<=m) & (c2>0) & (c2<=n)
	     if dftype(r1,c1,theta(r1,c1),r2,c2,theta(r2,c2))==2 ...
		  & fh(r1,c1)>0 & fh(r2,c2)>0
		   labels(r,c) = 1;
		   break;
		 end;
	  end;
	end;
  end;
end;
