function [ f_blurring ] = Func_GaussianBlurring2D( I,FWHM, Method,Pixel_Size,Pad)


%   Pixel_Size = 1e-6;     %   Pixel Size [m]
%   FWHM = 5;              %   [pixels]
%   Method = 1;


sigma = FWHM/2.0/sqrt(2.0*log(2.0));    % in pixel
SIGMA = FWHM*Pixel_Size/2.0/sqrt(2.0*log(2.0));  % in physical dimension


%% Sptial Convolution
% Method 1 is less optimized than method 2, beacuse the discretized Gassian
% kernal in spatial domain may not be accurate enough


if Method == 1  


    H = fspecial('gaussian',2^nextpow2(sigma*64),sigma);
    f_blurring = imfilter(I, H, 'replicate');
    
    return 
        
end
    
%% Frequency Domain    




[rows,cols]=size(I);       % rows and cols must be even
f_pad = padarray(I, [rows/2, cols/2], Pad);     % Pad the input wavefield matrix


delta_x = Pixel_Size;


%% The FT of the Gaussian kernel is known


if Method == 2
 
F_max = 0.5*(1/delta_x);
    delta_k = (2*F_max)/(rows*2);
    [Kx, Ky] = meshgrid((-rows:1:rows-1)*delta_k,(-cols:1:cols-1)*delta_k);
    F_gauss_kernel = exp( - (2*pi^2)*SIGMA^2*(Kx.^2+Ky.^2)) ; 
            %   The FT expression of the standard normalized Gassian distribution
    F_gauss_kernel = fftshift(F_gauss_kernel);
            %   In order to match the arrangement of F_pad
    F_pad = fft2(f_pad);
    
    f_pad_blurring = ifft2(F_pad.*F_gauss_kernel);


    f_blurring = abs(f_pad_blurring(rows/2+1:end-rows/2,cols/2+1:end-cols/2));   
    
    
    return
end


%% A more precise description of matlab FFT operations
if Method == 2.1


F_max = 0.5*(1/delta_x);
    delta_k = (2*F_max)/(rows*2);
    [Kx, Ky] = meshgrid((-rows:1:rows-1)*delta_k,(-cols:1:cols-1)*delta_k);    
    F_gauss_kernel = exp( - (2*pi^2)*SIGMA^2*(Kx.^2+Ky.^2)); 
            %   The FT expression of the standard normalized Gassian distribution
            %   arraging from -Fs_max to Fs_max
        
    F_pad = ifftshift(fft2(f_pad)*delta_x^2);
            %   The true fourier transform from a spatial domain function 
            %   fftshift is to rearrange it from -Fs_max to Fs_max
    
    f_pad_blurring = ifft2(ifftshift(F_pad.*F_gauss_kernel))*(rows*2*cols*2)*(delta_k)^2;
            %   The true inverse Fourier transform from a frequency domain function
            %   ifftshift is to rearrange it from 0 to 2*Fs_max
            %   The result starts from index 0 in the spatial domain
            
    f_blurring = abs(f_pad_blurring(rows/2+1:end-rows/2,cols/2+1:end-cols/2));   


end
%% The spatial domain expression of the Gaussian kernel is known 
% Method 3 is less optimized than method 2, beacuse in some cases where
% the sigma is too small, the spatial sampling is not enough to evaluate
% the Gaussian kernel


if Method == 3


    [x, y] = meshgrid((-rows:1:rows-1)*delta_x,(-cols:1:cols-1)*delta_x);
    gauss_kernel = 1/(2*pi*SIGMA^2) * exp(-(x.^2+y.^2)/(2*SIGMA^2));
            %   the standard normalized Gassian distribution    
    gauss_kernel = fftshift(gauss_kernel);
            %   rearrange the gassian kernel profile to start with the peak
            %   otherwise the kernel will introduce pixels delay
    F_gauss_kernel = fft2(gauss_kernel) * delta_x^2;    
            %  FFT needs multiplying delta_x^2 to perform ture FT magnitude
    
    F_pad = fft2(f_pad);
    
    f_pad_blurring = ifft2(F_pad.*F_gauss_kernel);    


    f_blurring = abs(f_pad_blurring(rows/2+1:end-rows/2,cols/2+1:end-cols/2));   
    
    return;
    
end
  


%% A more precise description of matlab FFT operations


if Method == 3.1
    
F_max = 0.5*(1/delta_x);
    delta_k = (2*F_max)/(rows*2);    
    [x, y] = meshgrid((-rows:1:rows-1)*delta_x,(-cols:1:cols-1)*delta_x);
    gauss_kernel = 1/(2*pi*SIGMA^2) * exp(-(x.^2+y.^2)/(2*SIGMA^2));
            %   the standard normalized Gassian distribution    
    gauss_kernel = fftshift(gauss_kernel);
            %   rearrange the gassian kernel profile to start with the peak
            %   otherwise the kernel will introduce pixels delay
    F_gauss_kernel = fft2(gauss_kernel) * delta_x^2;    
            %  FFT needs multiplying delta_x^2 to perform ture FT magnitude
            %  Here it starts with the frequency k=0
    F_gauss_kernel = fftshift(F_gauss_kernel);  
            %  Rearrange the result to start with -Fs_max
    
    F_pad = fft2(f_pad) * delta_x^2;
    F_pad = fftshift(F_pad);
    
    F_pad_blurring = F_pad.*F_gauss_kernel;
    
    f_pad_blurring = ifft2(ifftshift(F_pad_blurring))*(rows*2*cols*2)*(delta_k)^2;


    f_blurring = abs(f_pad_blurring(rows/2+1:end-rows/2,cols/2+1:end-cols/2));   
    
    return;
    
end


end