function [U,T]=ROF(Im,Init,lambda,tau,tol)

%ROF    An implementation of the Rudin-Osher-Fatemi (ROF) denoising model
%       using the numerical procedure presented in Eq. (11) of A. Chambolle
%       (2005). Observe that this version is implemented using periodic boundary 
%       conditions (essentially turning the rectangular image domain into a torus!).
%
%Output variables:
%U      - Denoised and detextured image (also the primal variable)
%T      - Texture residual
%
%Input variables:
%Im     - Noisy input image (grayscale)
%Init   - Initial guess for U
%lambda - Weight of the TV-regularizing term
%tau    - Steplength in the Chambolle algorithm (0< tau <=1/4; default=2.49)
%tol    - Tolerance for determining the stop criterion (default=10e-3)
%
%Program variables:
%U      - The primal field (also the denoised image)
%Px     - X-component of the dual field
%Py     - Y-cpmponent of the dual field
%
%Copyright Niels Chr Overgaard, Lund Institure of Techology and Malmoe
%University.
%Latest update: 10 November 2008


%---Initialization---------

[N,M] = size(Im); %size of noisy image.
U = Init; %the primal variable.
Px = U; %zeros(N,M); %x-component to the dual field.
Py = U; %zeros(N,M); %y-component of the dual field.
Err = 1; 
iter = 0;


%---Main iteration

while Err > tol
    
    %Save old iteration
    Uold = U;
    
    %Gradient of primal variable
    LxU = [U(2:N,:); U(1,:)]; %Left translation w.r.t. the x-direction
    LyU = [U(:,2:M) U(:,1)]; %Left translation w.r.t. the y-direction
    GradUx = LxU-U; %x-component of U's gradient
    GradUy = LyU-U; %y-component of U's gradient
    
    %First we update the dual varible
    PxNew = Px + (tau/lambda)*GradUx; %Non-normalized update of x-component (dual)
    PyNew = Py + (tau/lambda)*GradUy; %Non-normalized update of y-component (dual)
    NormNew = max(1,sqrt(PxNew.^2+PyNew.^2));
    Px = PxNew./NormNew; %Update of x-component (dual)
    Py = PyNew./NormNew; %Update of y-component (dual)
    
    %Then we update the primal variable
    RxPx = [Px(N,:) ; Px(1:(N-1),:)]; %Right x-translation of x-component
    RyPy = [Py(:,M) Py(:,1:(M-1))]; %Right y-translation of y-component
    DivP = (Px-RxPx)+(Py-RyPy); %Divergence of the dual field.
    U = Im + lambda*DivP; %Update of the primal variable
    
    %Update of error-measure
    Err = norm(U-Uold,'fro')/sqrt(N*M);
    iter = iter +1;
    
end

%The texture residual
T = Im - U;

%Displays number of iterations used
ITER = int2str(iter);
%disp(['ROF iterations =',ITER])