matlab高斯2D模糊的函式
% 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