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

% Image filter based on total variation minimalization.
% This function uses Poisson noise model, which replaced Gaussian model
% in ROF equation 
% 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] = TVPoiss(u)
% Filtration parameters settings
nest = noise_est(u);                                            % Noise estimation
nest = abs(1-((nest-0.7)/90));                                  % Computing of fist filtration parameter h1
if nest>0.7
    nest = 0.85;
end

f = u;
iter = 50;                                                      % Number of iterations                                               
pom1 = nest;
pom2 = pom1*0.05;                                               % Computing of second filtration parameter h2
u=f;
by = zeros(size(u)); bx = by;
b2 = by; 
%% Individually iterations 
for it = 1:iter
    Lap = LAP(u);                                               % Lamplacian
    dx = sign(diffx(u)+bx).*max(abs(diffx(u)+bx)-1/(nest),0);   % X direction D subproblem
    dy = sign(diffy(u)+by).*max(abs(diffy(u)+by)-1/(nest),0);   % y direction D subproblem
 
 
    z = (((u-(pom1)+b2)/2)+(sqrt(((u-(pom1)+b2).^2)+(pom1).*4.*f))./2); % Subproblem Z solving
    u = (z-b2)-(pom2).*((divX(dx-bx)+divY(dy-by))-Lap);         % Subproblem Z solving
    by=(diffy(u))+by-dy;                                        % error parameter of U and D sub.probl. in Y direction
    bx=(diffx(u))+bx-dx;                                        % error parameter of U and D sub.probl. in X direction
    b2 = b2+u-z;                                                % error parameter of Z subproblem
       
end
filtU = u;
%%Objective rating parameters
SNR = 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

