function [u,s_tot] = AnisoFilteringPCG(f)

% convert f to double
f = double(f);

%% some parameters
% size of the kernel used for calculating the smoothed edge indicator
kernel_size = 31;
kernel_var = 3;
filter_ker = fspecial('gaussian',kernel_size,kernel_var);

% 'steepness parameter' for the anisotropy
% smaller parameters lead to sharper edges, but also makes the problem more
% difficult to solve
epsilon_aniso = 1;

% regularisation parameter
lambda = 1500;

% maximal number of CG steps (ideally, the iteration should stop earlier)
n_max = 1000;

% size of the (squared) residual used for stopping the iteration
s_stop = 1;

%% computation of anisotropy
% smooth the image and compute the size of its gradient
f_smooth = imfilter(f,filter_ker,'symmetric','same','conv');
[f_x,f_y] = num_grad(f_smooth);
grad_size = f_x.^2 + f_y.^2;
tot_grad = sum(grad_size,3);

% actual anisotropy
aniso = 1./(tot_grad + epsilon_aniso^2);

%% computation of diagonal
aniso_diag = 2*aniso;
aniso_diag(2:end,:) = aniso_diag(2:end,:) + aniso(1:end-1,:);
aniso_diag(:,2:end) = aniso_diag(:,2:end) + aniso(:,1:end-1);
aniso_tot = zeros(size(f));
aniso_tot(:,:,1) = aniso_diag;
aniso_tot(:,:,2) = aniso_diag;
aniso_tot(:,:,3) = aniso_diag;

Adiag = ones(size(f)) + lambda*aniso_tot;
Adiag_inv = 1./Adiag;


%% Initialisation
% initial guess = 0
% it would be also possible to use f as initial guess
u = zeros(size(f));

% computation of initial residual
% note that the evaluation of div(a grad u) is performed in several steps
% also, the actual matrix of the problem is never computed
% (if u = 0, one could shorten this to the simple declaration r = f)
r = f - u;
[u_x,u_y] = num_grad(u);
au_x = aniso.*u_x;
au_y = aniso.*u_y;
dau = num_divergence(au_x,au_y);
r = r + lambda*dau;

% multiply the residual with the inverse of the preconditioner
z = Adiag_inv.*r;

% initialisation of search direction
p = z;

% initialise the size of the residual
s = sum(sum(sum(z.*r)));

% just for evaluation of the method: all the residual sizes are stored in
% some vector of increasing size - note that s contains the residuals for
% the preconditioned system!
s_tot = sum(sum(sum(r.*r)));

for k=1:n_max
    % evaluate the operator at p
    % again, computation of div(a grad p) takes several steps
    w = p;
    [p_x,p_y] = num_grad(p);
    ap_x = aniso.*p_x;
    ap_y = aniso.*p_y;
    dap = num_divergence(ap_x,ap_y);
    w = w - lambda*dap;

    % computation of step length
    alpha = s/sum(sum(sum(w.*p)));
    
    % updates of u and the residual
    u = u + alpha*p;
    r = r - alpha*w;

    % multiply the residual with the inverse of the preconditioner
    z = Adiag_inv.*r;
    
    % computation of new residual size (for the preconditioned system!)
    s_old = s;
    s = sum(sum(sum(z.*r)));
    
    % for evaluation of the method: attach the current residual size to
    % s_tot (again, we need to compute the residual for the actual equation
    % here)
    size_res = sum(sum(sum(r.*r)));
    s_tot = [s_tot,size_res];
    
    if(size_res < 1)
        k
        break;
    end
    
    % compute the next update direction
    beta = s/s_old;
    p = z + beta*p;
end

% plot the result
figure(3);imagesc(uint8(u));