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

% Image filter based on weighted averaging. 
% This function using floating window and in defined area is compared similarity
% between defined window and window aroud central pixel.
% Default the window size is set on 51 pixels and searching area is set on 3, but it
% can be changed by user.
% 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 [uout,PSNR, SNR] = NLM_filt(u,win,area)
if nargin<2
    win=51;
end
if nargin<3
    area = 3;
end
e = noise_est(u);                                               % Noise estimation
h=60+(1/e)*100;                                                 % Decay parameter
area= floor(area/2);                                            % Half size of area                                         % Half size of window
win=floor(win/2);                                               % Half size of window  
h=h*h;
u2 = padarray(u,[win,win],'symmetric');
x = size(u,1);
y = size(u,2);
uout=zeros(x,y);


K = make_kernel(win);                                           % Kernel with Gaussian distribution,which makes stronger
                                                                % centrals  pixels of window stronger  
                                                                
%% Chosing of central pixels for filtrations                                                              
parfor k=1:x
for l=1:y 
        kpom = k+win;
        lpom = l+win;
        
        win_W = u2(kpom-win:kpom+win,lpom-win:lpom+win);        % Chosen window

        %%% Position of central pixel on image edges
        kmin = max(kpom-area,win+1);
        kmax = min(kpom+area,x+win);
        lmin = max(lpom-area,win+1);
        lmax = min(lpom+area,y+win);
        AVG = 0;
        Z = 0;
        %%%
        %%% Compring between window and odher windows in defined area.
        for r = kmin:1:kmax
             for s = lmin:1:lmax

                    win_W2 = u2(r-win:r+win,s-win:s+win);       % Chosen other window area
                    E = (K.*(win_W-win_W2)).^2;                 % Euclidean distance 
                    E = (sum(E(:)));
                    w = exp(-E*h);                              % Weight computation with Decay parameter
                    Z = Z+w;
                    AVG = AVG + w*u2(r,s);                      % Weighted average
            end
       end
           if Z>0
                 uout(k,l) = AVG/(Z);
           else
               uout(k,l) = u(k,l);
           end
end
end
%% Objective rating parameters.
SNR = abs(10*log10((sum(u(:)))/(sum(sum(u-uout)))));
MSE = abs(sum(sum((uout-u).^2))/(size(uout,1)*size(uout,2)));    
PSNR = 10*log10(1/MSE);
 
 end
 