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

% if direction is taken into account, then need to compute:
%     alpha = atan2(0.5,1);
%     beta = theta - alpha
%     gamma = theta + alpha
%     phi = theta + pi - alpha
%     sigma = theta + pi + alpha
%     q = p*sqrt(5)/2    where p is the distance in pixels to look forward
% where theta is direction of gradient at pixel, and p=d/2
% Then the four corners of the rectangle aimed in the theta direction:
%    r1 = q*sin(beta);   c1 = q*cos(beta);
%    r2 = q*sin(gamma);  c2 = q*cos(gamma);
%    r3 = q*sin(phi);    c3 = q*cos(phi);
%    r4 = q*sin(sigma);  c4 = q*cos(sigma);

[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);
h = zeros(m,n);
for r = rmin:rmax
  for c = cmin:cmax
    grc = [gr(r,c),gc(r,c)];
	count = 0;
	for r2 = r-dhalf:r+dhalf
	  for c2 = c-dhalf:c+dhalf
	    if dot([gr(r2,c2),gc(r2,c2)],grc)>del
		  count = count + 1;
		end;
	  end;
	end;
	h(r,c) = count;
  end;
end;
h = h/(d*d);
