function [c1,c2] = ctopo4(f)
%
% ctopo4: find ravines by computing the curl of the 90 degree rotated 
%         gradient vectors 
% On input:
%   f:  2D image function
% On output:
%   c1:  curl of the 90 degree rotated gradient (equivalent to the negative
%       Laplacian) in one principal curvature direction
%   c2:  curl of the 90 degree rotated gradient (equivalent to the negative
%       Laplacian) in one principal curvature direction
% Call:
%     [c1,c2] = ctopo4(f)
%
tic
[fc,fr] = gradient(f);
[fcc,fcr] = gradient(fc);
[frc,frr] = gradient(fr);

c = -fcc-frr;
cmask2 = c*0;
cmask1 = c*0;
angle = posori(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
	  cmask1(i1,j1) = 0;
	  cmask2(i1,j1) = 0;
	  if locmaxdir(c,i1,j1,angle(i1,j1))
	    cmask1(i1,j1) = 1;
	  end;
	  if locmaxdir(c,i1,j1,rem((angle(i1,j1)+pi/2),2*pi))
	    cmask2(i1,j1) = 1;
	  end;
	end
  end
end;

c1 = c.*cmask1;
c2 = c.*cmask2;
toc
