function [l,cpit,cravine,cridge,cpeak] = churlabel(f,t1,t2,k)
%
% churlabel: uses the curl of the 90 degree rotated gradient vectors
%            to determine if there is a ridge (positive curl) or a
%            ravine (negative curl)
% On input:
%   f:  2D image function
%   t1: threshold for positive curl
%   t2: threshold for negative curl
%   k:  window size for non-maxima suppression for peaks and pits
% On output:
%   l:  all labels in a bit code
%       0: no label
%       1: pit
%       2: ravine
%       4: ridge
%       8: peak
%   cpit: pit image
%   crav: ravine image
%   crid: ridge image
%   cpeak: peak image
%

[fc,fr] = gradient(f);
[fcc,fcr] = gradient(fc);
[frc,frr] = gradient(fr);
[m,n] = size(f);
l = zeros(m,n);
cridge = l; cravine = l; cpeak = l; cpit = l;

curl = -fcc-frr;

cridge  = 4*(curl>t1);
cravine = 2*(curl<t2);
cpeak = 8*(nmaxs2(f.*cridge,k)>0);
cpit = 1*(nmins2(f.*cravine,k)>0);

l = cridge+cravine+cpeak+cpit;
