function [c1,c2] = cridges(f)
%
% cridges: find ridges in a 2D image (using curl of rotated gradient)
% On input:
%   f:  2D image function
% On output:
%   c:  non-hillside, local image maximum of positive curl of the 90 degree 
%       locally rotated gradient (equivalent to the negative Laplacian), with
%     c1: maxima in first principal direction
%     c2: maxima in second principal direction
% Call:
%     [c1,c2] = cridges(f)
%

tic
[fc,fr] = gradient(f);
[fcc,fcr] = gradient(fc);
[frc,frr] = gradient(fr);

c = -fcc-frr;
cmask2 = c*0;
cmask1 = c*0;
angle = posori2(atan2(fcr,fcc-frr)/2);
cmask = c>0;
[m,n] = size(c);
for i1 = 2:m-1
  for j1 = 2:n-1
    if cmask(i1,j1)>0
      if locmaxdir(f,i1,j1,angle(i1,j1))
        cmask1(i1,j1) = 1;
      end;
      if locmaxdir(f,i1,j1,rem((angle(i1,j1)+pi/2),2*pi))
        cmask2(i1,j1) = 1;
      end;
    end
  end
end;
h = TGhillsl(fc,fr,3);
c1 = c.*cmask1.*(h<=0);
c2 = c.*cmask2.*(h<=0);
toc
