% Matlab code to denoise an image by the Rudin-Osher-Fatemi model
% Min TV(u) + lambda*(L^2[(f-u)])^2
% lambda = coefficient of the fidelity term
% lambda needs to be tunned to obtain a sharp denoised image
% h = space discretization step
% The stationary problem (no time-regularization) is used here
% A fixed point iteration is used
% The numerical scheme is implicit in the max norm and unconditionally stable
% A stopping criteria can be included, or the "IterMax" must be given.

% LAST MODIFIED: 2005-04-19

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT u0 = noisy image

% OUTPUT: u(i,j) = denoised image

% read the noisy data u0 in BMP format

u0=imread('fic.noise.bmp');
u0=double(u0);
[M N]=size(u0);

% visualize the image u0 in Matlab (rescaled)
imagesc(u0); axis image; axis off; colormap(gray);

%---------------------------------------------------------------------------
% initialize u by u0 (not necessarily) or by a constant, or by a random image
%---------------------------------------------------------------------------

u=u0;
[M,N]=size(u);

%-----------------------------------------------
% PARAMETERS
%-----------------------------------------------

% coefficient of the TV norm (needs to be adapted for each image)
lambda=0.028;

% space discretization
h=1.;

% number of iterations (depends on the image)
IterMax=100;

% needed to regularize TV at the origin
eps=0.000001;

%-----------------------------------------------------------
% BEGIN ITERATIONS IN ITER
%-----------------------------------------------------------

for Iter=1:IterMax,

Iter

for i=2:M-1,
for j=2:N-1,

%-----------computation of coefficients co1,co2,co3,co4---------

ux=(u(i+1,j)-u(i,j))/h;
uy=(u(i,j+1)-u(i,j))/h;
Gradu=sqrt(eps*eps+ux*ux+uy*uy);
co1=1./Gradu;

ux=(u(i,j)-u(i-1,j))/h;
uy=(u(i-1,j+1)-u(i-1,j))/h;
Gradu=sqrt(eps*eps+ux*ux+uy*uy);
co2=1./Gradu;

ux=(u(i+1,j)-u(i,j))/h;
uy=(u(i,j+1)-u(i,j))/h;
Gradu=sqrt(eps*eps+ux*ux+uy*uy);
co3=1./Gradu;

ux=(u(i+1,j-1)-u(i,j-1))/h;
uy=(u(i,j)-u(i,j-1))/h;
Gradu=sqrt(eps*eps+ux*ux+uy*uy);
co4=1./Gradu;

co=1.+(1/(2*lambda*h*h))*(co1+co2+co3+co4);

div=co1*u(i+1,j)+co2*u(i-1,j)+co3*u(i,j+1)+co4*u(i,j-1);

u(i,j)=(1./co)*(u0(i,j)+(1/(2*lambda*h*h))*div);

end
end

%------------ FREE BOUNDARY CONDITIONS IN u -------------------
for i=2:M-1,
u(i,1)=u(i,2);
u(i,N)=u(i,N-1);
end

for j=2:N-1,
u(1,j)=u(2,j);
u(M,j)=u(M-1,j);
end

u(1,1)=u(2,2);
u(1,N)=u(2,N-1);
u(M,1)=u(M-1,2);
u(M,N)=u(M-1,N-1);
%----------------------------------------------------------------------

%---------------- END ITERATIONS IN Iter ------------------------------

%%% Compute the discrete energy at each iteration
en=0.0;
for i=2:M-1,
for j=2:N-1,
ux=(u(i+1,j)-u(i,j))/h;
uy=(u(i,j+1)-u(i,j))/h;
fid=(u0(i,j)-u(i,j))*(u0(i,j)-u(i,j));
en=en+sqrt(eps*eps+ux*ux+uy*uy)+lambda*fid;
end
end
%%% END computation of energy

Energy(Iter)=en;

end

% visualize the image in Matlab (re-scaled)
figure
imagesc(u); axis image; axis off; colormap(gray);

% use export to save the image in ps or eps format

% Plot the Energy versus iterations
figure
plot(Energy);legend('Energy/Iterations');