function [f,p_new] = rof( f0, beta, tau, h, theta, maxNumIter )
    [m,n] = size(f0)

    f = f0(:);

    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;

    for i=1:maxNumIter
        disp(i);
        p = p_new;

        tmp = div_x*p(:,1) + div_y*p(:,2) - 2*f/beta;
        grad = [grad_x*tmp grad_y*tmp];

        denom = 1 + tau * sqrt( grad(:,1).^2 + grad(:,2).^2 );
        p_new = ( p + tau * grad ) ./ ( denom * ones (1, 2) ); 
        
        %% stopping criterion
        if ( max(max(abs(p_new-p))) < theta )
            break;
        end
    end
 
    %% compute denoised image
    f = -(beta/2) * (div_x*p(:,1) + div_y*p(:,2)) + f;
  
    f = reshape(f,m,n);

end


%% computes matrix representation of the x-derivative
function dx = xDerivative( m, n, h )
    
    rows = m*n;
    nzs = 2*m*(n-1);
    row = zeros(nzs, 1);
    col = zeros(nzs, 1);
    val = zeros(nzs, 1);
    count = 0;
    for i=1:(n-1)
        for j=1:m
            count = count+1;
            row(count) = (i-1)*m+j;
            col(count) = (i-1)*m+j;
            val(count) = -1/h;
            count = count+1;
            row(count) = (i-1)*m+j;
            col(count) = i*m+j;
            val(count) = 1/h;
        end
    end
    dx = sparse(row, col, val, rows, rows);
end

%% computes matrix representation of the y-derivative
function dy = yDerivative( m, n, h )
    
    rows = m*n;
    nzs = 2*(m-1)*n;
    row = zeros(nzs, 1);
    col = zeros(nzs, 1);
    val = zeros(nzs, 1);
    count = 0;
    for i=1:n
        for j=1:(m-1)
            count = count+1;
            row(count) = (i-1)*m+j;
            col(count) = (i-1)*m+j;
            val(count) = -1/h;
            count = count+1;
            row(count) = (i-1)*m+j;
            col(count) = (i-1)*m+j+1;
            val(count) = 1/h;
        end
    end
    dy = sparse(row, col, val, rows, rows);
end

