% Bc. Zbynek Holub, Diploma thesis: Denoising of Images from Electron Microscope

% Image filter based on total variation minimalization.
% This function uses ROF model, computed by split Bregman algorithm.
% On input is only image and filtration parameter is set automatically by
% standart deviation, which is calculated by noise estimation.. For estimation is used function from
% https://www5.cs.fau.de/our-team/schwemmer-chris/software/. Noise est. function
% is decribed in artickle http://electronicimaging.spiedigitallibrary.org/article.aspx?doi=10.1117/1.3476329
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [filtU, PSNR, SNR] = TVROF(u)
% noise estimation 
nest = noise_est(u);
% Filtration parameters settings 
lambda = (60/nest)*5;
if lambda<12
    lambda = 12;
end

gama =  1- 60/(nest*100);

if gama > 0.8
    gama = 0.8;
end
%Number of iterations 
iter = 50; 

f = u;
by = (zeros(size(u))); bx = by;

% total variation iterations - split Bregman
for it = 1:iter
    Lap = LAP(u);                                                           % Laplacian computing
    dx = sign(diffx(u)+bx).*max(abs(diffx(u)+bx)-1/gama,0);                 % X direction D subproblem
    dy = sign(diffy(u)+by).*max(abs(diffy(u)+by)-1/gama,0);                 % Y direction D subproblem

    u = (lambda*f -( (divX(dx-bx))+(divY(dy-by)))+gama*(Lap))./(lambda);    % U subproblem computing, result is
                                                                            % filtrated image in current iteration.

    by=(diffy(u))+by-dy;                                                    % error parameter in Y direction
    bx=(diffx(u))+bx-dx;                                                    % error parameter in X direction
       
end
filtU = u;                                                                  % Denoised image
% Objective rating parameters
SNR = abs(10*log10((sum(u(:)))/(sum(sum(f-u)))));
MSE = abs(sum(sum((u-f).^2))/(size(u,1)*size(u,2)));                                             
PSNR = 10*log10(1/MSE);
end

