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

% Image filter based on total variation minimalization.
% This function uses ROF filtration model, which is solved by Elueler-Lagrange method. 
% 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 [u, PSNR, SNR] = EL_TV(u)
%%
e = noise_est(u);                                               % Noise estimation 
lambda = e*3;                                                   % Computing of filtration parameter
dt = 0.2;                                                       % Differences step in each iteration
iter = 10;                                                      % Number of itertions 
f = u;
%% Filtration iterations 
[x,y] = size(u);
for i = 1:dt:iter

%%% First and second differences in directions x,y and xy.
ux = ( u(:,[2:y,y])-u(:,[1,1:y-1]))/2;
uy = (u([2:x,x],:)-u([1,1:x-1],:))/2;
uxy = (u([2:x,x],[2:y,y])+u([1,1:x-1],[1,1:y-1])-u([1,1:x-1],[2:y,y])-u([2:x,x],[1,1:y-1]))/4;
uxx = u(:,[2:y,y])-2*u+u(:,[1,1:y-1]);
uyy = u([2:x,x],:)-2*u+u([1,1:x-1],:);
%%%

Nom =((uxx.*((uy).^2))-(2*(ux.*uy.*uxy))+(uyy.*((ux).^2)));      % Nominator of filtration equation
Den =((((ux.^2)+(uy.^2)).^(3/2))+0.01);                          % Deominator of filtration equation
uNewPom = (Nom./Den);                                            
u = u+dt*(uNewPom)-2*(1/lambda)*(u-f);                           % Filtration in one iteration


end
%% Objective rating parameters
SNR = abs(10*log10((sum(u(:)))/(sum(sum(f-u)))));                           

MSE1=sum(sum((u-f).^2))/(x*y);                                             
PSNR= 10*log10(1/MSE1);
end
