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)
% 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)