function [imkap,sg,skap,imstr] = imcurvature(im,beta,sd,CFlag)
% 
% return Image Curvature, together with gradient, and normalising factor
%
% Cflag is a switch
%    Cflag == 1 : This is : 'Lvv curvature'
%    Cflag == 2 ; This is "normalised Lvv curvature" Lvv/Ln^3
%    Cflag == 3 : This is "divergence of normal"
%    Cflag == 4 : This is "divergence of tangent"
%
n = size(im);
[imstr.imx,imstr.imy] = gradient(im);
g = sqrt(imstr.imx.^2 + imstr.imy.^2);
if(sd == 0)
    imstr.sim = im;
    imstr.simx = imstr.imx;
    imstr.simy = imstr.imy;
else
    % smoothed version. Use Fourier smoothing
    imstr.sim= imblur(im,sd);
    [imstr.simx,imstr.simy] = gradient(imstr.sim);
end
% if beta==0
%     sg = sqrt((simx).^2 + (simy).^2);
%     skap = 1./sg;
%     gzind = find(sg==0);
%     skap(gzind) = 0;
% else
% skap = beta./sqrt((simx).^2 + (simy).^2 + beta^2);
%end

sg = sqrt(imstr.simx.^2 + imstr.simy.^2);
skap = 1./sg;
gzind = find(sg <= beta);
skap(gzind) = 0;
% make a mask of regions where gradient is non-zero
imstr.mask = ones(size(im)); imstr.mask(gzind) = 0;
% now expand by one pixel
[mx,my] = gradient(imstr.mask); gz1 = find(mx ~= 0 | my ~=0);
gz1 = union(gzind,gz1);
imstr.mask1 = imstr.mask; imstr.mask1(gz1) = 0;
% and again...
[mx,my] = gradient(imstr.mask1); gz2 = find(mx ~= 0 | my ~=0);
gz2 = union(gz1,gz2);
imstr.mask2 = imstr.mask1; imstr.mask2(gz2) = 0;

imstr.sx = skap.*imstr.simx; 
imstr.sy = skap.*imstr.simy; % "sine and cosine"
[imstr.simxx,imstr.simxy] = gradient(imstr.simx);
[imstr.simyx,imstr.simyy] = gradient(imstr.simy);
[imstr.sxx,imstr.sxy] = gradient(imstr.sx);
[imstr.syx,imstr.syy] = gradient(imstr.sy);

imstr.sx(gz1) = 0;imstr.sy(gz1) = 0;
imstr.sxx(gz2) = 0;imstr.sxy(gz2) = 0;
imstr.syx(gz2) = 0;imstr.syy(gz2) = 0;
imstr.simx(gz1) = 0;imstr.simy(gz1) = 0;
imstr.simxx(gz2) = 0;imstr.simxy(gz2) = 0;
imstr.simyx(gz2) = 0;imstr.simyy(gz2) = 0;

switch CFlag
    case 1
        imkap = imstr.simy.^2 .* imstr.simxx + imstr.simx.^2 .* imstr.simyy - imstr.simx.*imstr.simy.*(imstr.simxy+imstr.simyx);
    case 2
%         imkap = simy.^2 .* simxx + simx.^2 .* simyy - simx.*simy.*(simxy+simyx);
%         imkap = imkap.*(skap.^3);
        imkap = imstr.sx.^2 .*imstr.sxx + imstr.sy.^2 .*imstr.syy - imstr.sx.*imstr.sy.*(imstr.sxy+imstr.syx);
    case 3  % divergence of normal
        %        imkap = divergence(sx,sy);
        imkap = (imstr.sxx + imstr.syy);
    case 4  % divergence of tangent
        %        imkap = divergence(sx,sy);
        imkap = imstr.sxy - imstr.syx;
    otherwise
        disp(['Unknown method. Clflag=',num2str(CFlag)]);
end
if beta==0
    imkap(gzind) = 0;
end