%im = phantom(256);
%im = double(imread('../../StandardTestImages/house.png','PNG'));
im = double(imread('../../StandardTestImages/Cameraman256.png','PNG'));

% add noise
%im = im.*(1 + 0.1*randn(size(im)));

[gx,gy] = gradient(im);
[gxx,gxy] = gradient(gx);
[gyx,gyy] = gradient(gy);
gg = sqrt(gx.^2 + gy.^2);

figure(1); clf; hold on;
subplot(2,3,1);colormap(gray); imagesc(im);title('original');
subplot(2,3,2);colormap(gray);  imagesc(gx);title('f_x');
subplot(2,3,3);colormap(gray);  imagesc(gy);title('f_y');
subplot(2,3,4);colormap(gray); imagesc(gxx);title('f_{xx}');
subplot(2,3,5);colormap(gray);  imagesc(gxy);title('f_{xy}');
subplot(2,3,6);colormap(gray);  imagesc(gyy);title('f_{yy}');


%%
for t = 1:10:360
    sn = sin(t*2*pi/360);
    cs = cos(t*2*pi/360);
    
figure(2); clf;
subplot(2,3,1);colormap(gray); imagesc(im);title('original');
subplot(2,3,2);colormap(gray);  imagesc(gx);title('f_x');
subplot(2,3,3);colormap(gray);  imagesc(gy);title('f_x');

subplot(2,3,4);colormap(gray);  imagesc(gg);title('| grad f |');
subplot(2,3,4);colormap(gray);  imagesc(cs.*gx + sn.*gy);title(['f_{\theta}, \theta = ',num2str(t),' degrees']);
subplot(2,3,6);colormap(gray);  imagesc(cs^2.*gxx + 2*cs.*sn.*gxy + sn^2.*gyy);title(['f_{\theta,\theta}, \theta = ',num2str(t),' degrees']);
pause(0.5);
end

%%
eps = 1e-3;
gc = gx./(gg+eps);
gs = gy./(gg+eps);
gnn = (gc.^2).*gxx + 2*gc.*gs.*gxy + (gs.^2).*gyy;
unn = (gx.^2).*gxx + 2*gx.*gy.*gxy + (gy.^2).*gyy;
unn2 = sqrt(abs(unn)).*sign(unn);
figure(3); clf;
subplot(2,3,1);colormap(gray); imagesc(im);title('original');
subplot(2,3,2);colormap(gray); imagesc(gnn);title('f_{nn}');
subplot(2,3,3);colormap(gray); imagesc(unn);title('|grad f|^2 f_{nn}');
subplot(2,3,4);colormap(gray); imagesc(gc);title('f_x/|grad f|');
subplot(2,3,5);colormap(gray); imagesc(gs);title('f_y/|grad f|');
subplot(2,3,6);colormap(gray); imagesc(unn2);title(' sqrt(||grad f|^2 f_{nn}|');