function ROF

    function f = rof( f0, lambda, h, maxNumIter, stopEps )
        [m,n] = size(f0);
        u = f0(:);
        
        tau = 0.125;
        sigma = tau;
        theta = 1.0;
 
        grad_x = xDerivative( m, n, h );
        grad_y = yDerivative( m, n, h );
        
        div_x = -grad_x';
        div_y = -grad_y';
        
        p = zeros(n*m, 2);
        p_new = p;
        u_new = u;
        
        for i=1:maxNumIter
            if mod(i,10) == 0
                fprintf('%d steps. ||u^(n+1)-u^n || = %f \n',[i,threshold]);
            end
            p = p_new;
            u = u_new;
            
            %% proximal operator of H^*
            tmp = p + sigma*[grad_x*u grad_y*u];
            denom = max ( 1, sqrt ( tmp(:,1).^2 + tmp(:,2).^2 ) );
            p_new = [tmp(:,1)./denom tmp(:,2)./denom];
            
            %% proximal operator of G
            u_new = ( u + tau * ( div_x*p_new(:,1) + div_y*p_new(:,2) + f0(:)/lambda) ) / ( 1 + tau/lambda);
            
            %% acceleration step
            u_new = u_new + theta * ( u_new - u );

            %% stopping criterion
            threshold = max(abs(u_new(:)-u(:)));
            if ( threshold < stopEps )
                break;
            end
        end
        
        f = reshape(u_new,m,n);
        
    end

    %% computes matrix representation of the x-derivative
    function dx = xDerivative( m, n, h )
        d = spdiags(ones(m,1)*[-1,1],[0,1],m,m)/h;
        dx = sparse(kron(speye(n),d));
    end

    %% computes matrix representation of the y-derivative
    function dy = yDerivative( m, n, h )
        d = spdiags(ones(n,1)*[-1,1],[0,1],n,n)/h;
        dy = sparse(kron(d,speye(m)));
    end

clear
close all

%% Parameter
lambda=0.1;
h = 1;
maxIter = 1000;
stopEps = 1e-3;

%% Read input image
f = double(imread('mandril.png'))/255;

%% Rudin-Osher-Fatemi-denoising
fDenoised = rof(f,lambda,h,maxIter,stopEps);

subplot(1,2,1)
imshow(f)
title('Input image')

subplot(1,2,2)
imshow(fDenoised)
title('ROF result')

end