% visualize the gradient and newton fields of the quadratic cost functional,
%   J(x,y) = (x-0.1).^2/2 + (y-0.5).^2 + mu*(x-y).^2/2
% and the non-quadratic cost functional
%   J(x,y) = (x-0.1).^2/2 + (y-0.5).^2 + mu*sqrt((x-y).^2 + ep^2);
% which are 2D model functions for the image processing counterparts.  Take
%   X = (x,y),    F(X) = DJ(X),  DF(X) = D^2 J(X)

% start fresh
clear all;
clc;

% figure
h1 = figure(1);
close(h1);
h1 = figure(1);
set(h1,'Position',[10 10 900 600]);

% grid
N = 50;
x = (0:N)/N;
y = x';
xx = kron(x,ones(size(y)));
yy = kron(ones(size(x)),y);

% regularization
mu = 1.0e0;
ep = 1.0e-2;

% number of contours
nc=101;
% vector scaling
vs=2;

% functions
f = 0.5*(xx-0.1).^2 + 0.5*(yy-0.5).^2 + mu*sqrt((xx-yy).^2 + ep^2);
g = 0.5*(xx-0.1).^2 + 0.5*(yy-0.5).^2 + 0.5*mu*(xx-yy).^2;

% gradient of f
fx = (xx-0.1) + mu*(xx-yy)./sqrt((xx-yy).^2 + ep^2);
fy = (yy-0.5) + mu*(yy-xx)./sqrt((yy-xx).^2 + ep^2);

% gradient of g
gx = (xx-0.1) + mu*(xx-yy);
gy = (yy-0.5) + mu*(yy-xx);

% hessian of f
fxx = 1 + mu*ep^2./sqrt((xx-yy).^2 + ep^2).^3;
fxy =   - mu*ep^2./sqrt((xx-yy).^2 + ep^2).^3;
fyy = 1 + mu*ep^2./sqrt((xx-yy).^2 + ep^2).^3;

% modified hessian of f
fmxx = 1 + mu./sqrt((xx-yy).^2 + ep^2);
fmxy =   - mu./sqrt((xx-yy).^2 + ep^2);
fmyy = 1 + mu./sqrt((yy-xx).^2 + ep^2);

% hessian of g
gxx = (1+mu)*ones(size(xx));
gxy =   -mu *ones(size(xx));
gyy = (1+mu)*ones(size(xx));

% Newton field for f
dd=fxx.*fyy-fxy.^2;
dd=(dd ~= 0)./(dd + (dd == 0));
uf=(fy.*fxy-fx.*fyy).*dd;
vf=(fx.*fxy-fy.*fxx).*dd;

% Newton field for fm
dd=fmxx.*fmyy-fmxy.^2;
dd=(dd ~= 0)./(dd + (dd == 0));
ufm=(fy.*fmxy-fx.*fmyy).*dd;
vfm=(fx.*fmxy-fy.*fmxx).*dd;

% Newton field for g
dd=gxx.*gyy-gxy.^2;
dd=(dd ~= 0)./(dd + (dd == 0));
ug=(gy.*gxy-gx.*gyy).*dd;
vg=(gx.*gxy-gy.*gxx).*dd;

% show gradient field for f
subplot(2,3,1);
contour(xx,yy,f,nc)
title('-DF for |DX|_ep');
hold on;

quiver(xx,yy,-fx,-fy,vs,'k');
hold off;
axis equal;
axis fill;
axis tight;
axis off;

% show gradient field for g
subplot(2,3,4);
contour(xx,yy,g,nc)
title('-DF for |DX|^2');
hold on;

quiver(xx,yy,-gx,-gy,vs,'k');
hold off;
axis equal;
axis fill;
axis tight;
axis off;

% show Newton field for f
subplot(2,3,2)
contour(xx,yy,f,nc)
title('-DF\\F for |DX|_ep');
hold on;

m=1.0e0;
quiver(xx,yy,min(max(uf,-m),m),min(max(vf,-m),m),vs,'k');
hold off;
axis equal;
axis fill;
axis tight;
axis off;

% show modified Newton field for f
subplot(2,3,3)
contour(xx,yy,f,nc)
title('-DFm\\F for |DX|_ep');
hold on;

m=1.0e0;
quiver(xx,yy,min(max(ufm,-m),m),min(max(vfm,-m),m),vs,'k');
hold off;
axis equal;
axis fill;
axis tight;
axis off;

% show Newton field for g
subplot(2,3,5)
contour(xx,yy,g,nc)
title('-DF\\F for |DX|^2');
hold on;

quiver(xx,yy,ug,vg,vs,'k');
hold off;
axis equal;
axis fill;
axis tight;
axis off;
