%%  set parameters

impath = '../../StandardTestImages/';

BLURDATA = 0.01; % blur image to give less severe structure
%im = phantom(256);
%im = ones(256,256);im(50:100,100:150) = 2; im(150:200,100:150) = 0;
%pim = imread('pebbles.png','png'); im = imresize(sum(pim,3),[256,256]);
%im = double(imread([impath 'Cameraman256.png'],'PNG'));
im = double(imread([impath 'house.png'],'PNG'));
%im = double(imread([impath 'ot.pgm'],'pgm'));
%im = double(imread([impath 'Lena512.png'],'PNG'));

nlevel =  0;%0.01; % noiselevel
sigma = 0; %0.001; %0.002; % blur kernel in fraction of image width
imax = max(max(im));
im = im./imax; % normalise image
beta = 1e-2; % additional for divide by zero;
if BLURDATA > 0
     im = imblur(im,BLURDATA);
end
im = im + nlevel*randn(size(im));

%% calculate curvatures
[imLvvUn,g,skap,imstr] = imcurvature(im,beta,sigma,1);
[imLvv,g,skap,imstr] = imcurvature(im,beta,sigma,2);
[imdivnorm,g,skap,imstr] = imcurvature(im,beta,sigma,3);
[imdivtang,g,skap,imstr] = imcurvature(im,beta,sigma,4);
kmax = max(imdivnorm(:));
kmin = min(imdivnorm(:));

%% display some system parameters
figure(1); clf; colormap(hot);
subplot(3,4,1); imagesc(g); title('gradient');
subplot(3,4,2); imagesc(imstr.mask); title('gradient mask');
subplot(3,4,3); imagesc(imstr.mask1); title('mask dilated 1');
subplot(3,4,4); imagesc(imstr.mask2); title('mask dilated 2');
subplot(3,4,5); imagesc(imstr.simx); title('f_x');
subplot(3,4,6); imagesc(imstr.simy); title('f_y');
subplot(3,4,7); imagesc(imstr.simxx); title('f_{xx}');
subplot(3,4,8); imagesc(imstr.simyy); title('f_{yy}');

subplot(3,4,9); imagesc(imstr.sx); title('n_x');
subplot(3,4,10); imagesc(imstr.sy); title('n_y');
subplot(3,4,11); imagesc(imstr.sxx); title('n_{xx}');
subplot(3,4,12); imagesc(imstr.syy); title('n_{yy}');

%% display
figure(2); clf; colormap(hot); %colormap(bone);
subplot(2,4,1); imagesc(im); colorbar; % (gray);
subplot(2,4,2); imagesc(g);colorbar;  title('gradient');
subplot(2,4,3); imagesc(skap);colorbar;  title('reciprocal gradient');% colormap(gray);
subplot(2,4,5); imagesc(imLvvUn);colorbar;  title('Unnormalised L_{vv} image curvature'); %colormap(gray);
subplot(2,4,6); imagesc(imLvv,[kmin,kmax]); colorbar; title('Normalised L_{vv}/(L_n^3) image curvature');%colormap(gray);
%subplot(2,4,6); imagesc(sign(imLvv).*log(abs(imLvv))); colorbar; title('Normalised L_{vv}/(L_n^3) image curvature');
subplot(2,4,7); imagesc(imdivnorm); colorbar; title('Divergence of normal');% \nabla\cdot\hat{N}'); %colormap(gray);
subplot(2,4,8); imagesc(imdivtang); colorbar; title('Divergence of tangent');% \nabla\cdot\hat{N}'); %colormap(gray);
subplot(2,4,4); hold on; plot(imLvv(:,150),'k'); plot(imdivnorm(:,150),'--r');
%plot(imdivtang(:,150),'--go')