function h = TGhills2(f,d,del)
%
% TGhills: Tom's Gradient Topographic Sketch hill finder
% On input:
%     f:   image function (2D)
%     d:   analyze d by d area
%     del: same direction if within del degrees 
% On output:
%     h:   hill likelihood in interval [0,d*d]    
% Calling sequence:
%    h = TGhills2(f,d,del)

r2d = 180/pi;
[m,n] = size(f);
dhalf = floor((d-1)/2);
rmin = dhalf+1; rmax = m-dhalf; cmin = dhalf+1; cmax = n-dhalf;
[gc,gr] = gradient(f);
theta = round(atan2(gr,gc)*r2d);
h = zeros(m,n);
for r = rmin:rmax
  for c = cmin:cmax
    t1 = theta(r,c); tup = t1+del; tdn = t1-del;
	count = 0;
	for r2 = r-dhalf:r+dhalf
	  for c2 = c-dhalf:c+dhalf
	    t2 = theta(r2,c2);
	    if (t2<tup) | (t2>tdn)
	       count = count + 1;
	    end;
	  end;
	end;
	h(r,c) = count;
  end;
end;
h = h/(d*d);
