function h = TGhillsl(f,K)
%
% TGhillsl returns non-zero at each pixel with surrounding gradient vectors
%         all within 90 degrees of each other
% On input:
%     f: 2D image function
%     K: window dimension
% On output:
%     h: hill labels (value is max angle dif in window)
% Call:
%     h = TGhillsl(f,K)
%

[m,n] = size(f);
k = (K+1)/2;
h = zeros(m,n);

[fc,fr] = gradient(f);
theta = atan2(fr,fc);
ptheta = posori(theta);

for r = k:m-k+1
  for c = k:n-k+1
    win = ptheta(r-k+1:r+k-1,c-k+1:c+k-1);
    tmax = max(max(win));
    tmin = min(min(win));
    tdif = tmax-tmin;
    if tdif<pi/2
      h(r,c) = tdif;
    end;
  end;
end;
