% denoise a noisy step function v by minimizing the functional
% J(u) = \int (u-v)^2 + mu \int sqrt(|u'|^2+ep^2)

% start fresh
clear all;
clc;

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

% grid
N  = 128;
x  = (0:N)/N;
h  = 1/N;

% regularization
mu = 1.0e-2;
% rounding in cost
ep = 1.0e-3;
% noise level
ns = 0.05;
% max number of iterations
imax = 100;
% termination criterion for iterations
tol = 1.0e-5;

% theoretical convergence criterion
if ((4*mu/(ep*h^2)) > 1)
    warning(strcat('should be less than one:',num2str(4*mu/(ep*h^2))));
end

% derived parameters
N1      = N+1;
muh2    = mu/h^2;
muh2ep2 = muh2*ep^2;

% exact step function
ux = ones(N1,1);
ux(1:round(N/2)) = 0;

% data
randn('state',0);
v  = ux + ns*randn(N1,1);

% initialize solution
% u  = v;
u  = ux;

% fixed-point iteration
for i=1:imax
    Q2 = ((u(2:N1) - u(1:N))/h).^2;
    E  = 1./sqrt(Q2 + ep^2);
    D  = diag([E;0]+[0;E]) - diag(E,1) - diag(E,-1);
% evaluation of fixed-point function
    G  = v - muh2*D*u;
% solution update
    du = G - u;
    u  = u + du;
    if (norm(du) < tol*norm(u))
        break;
    end
end

% report performance
if (i < imax)
    disp(strcat('iterations to convergence: ',num2str(i)));
else
    warning('maximum iterations reached');
end

% show result
plot(x,ux,'b',x,v,'r',x,u,'g');
legend('exact','noisy','denoised');
title('denoising');
xlabel('spatial dimension');
ylabel('intensity');
