function u = TV_L2(u0,lambda, IterMax)
% 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 

u0=double(u0);
[M N]=size(u0);

% visualize the image u0 in Matlab (rescaled) 
figure(1); subplot(1,3,1);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  
%-----------------------------------------------
if nargin < 2
% coefficient of the TV norm (needs to be adapted for each image) 
lambda=0.028; 

%  number of iterations (depends on the image) 
IterMax=100;
end
%  needed to regularize TV at the origin 
eps=0.000001;

%  space discretization 
h=1.;

%-----------------------------------------------------------
%     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 -------------------
u(:,1)=u(:,2);
u(:,N)=u(:,N-1);
u(1,:)=u(2,:);
u(M,:)=u(M-1,:);
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(1); subplot(1,3,2);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(1);subplot(1,3,3);plot(Energy);legend('Energy/Iterations');