Matlab code: gaussian blur in domain
I am going through a code that does blur images. However, I am having a hard time understanding the code, and I was wondering if anyone could help me understand what the code is about.
Here the variable "Iref" is an image.
Imin = min(Iref(:));
Iref_fft = Iref-Imin;
Iref_fft = fftshift(Iref_fft,1);
Iref_fft = fftshift(Iref_fft,2);
Iref_fft = fft(Iref_fft,[],1);
Iref_fft = fft(Iref_fft,[],2);
Iref_fft = fftshift(Iref_fft,1);
Iref_fft = fftshift(Iref_fft,2);
Here I am already confused as to what it means to apply fftshift to an image that is not already in the fourier domain. So, I can tell that it performs Fourier transform along each of the axes, but why does it do fftshift before and after?
The code follows like this:
Nx_r = 32;
Ny_r = 32;
sigma = 1.5;
wx = reshape(gausswin(Nx_r,sigma), [1 Nx_r]);
wy = reshape(gausswin(Ny_r,sigma), [Ny_r 1]);
wx_rep = repmat(wx, [Ny_r 1]);
wy_rep = repmat(wy, [1 Nx_r]);
Window = wx_rep .* wy_rep;
xIndices = floor((Nx-Nx_r)/2)+1 : floor((Nx-Nx_r)/2)+Nx_r;
yIndices = floor((Ny-Ny_r)/2)+1 : floor((Ny-Ny_r)/2)+Ny_r;
Iref_blurred = zeros(Ny,Nx);
Iref_blurred(yIndices,xIndices,:) = Iref_fft(yIndices,xIndices) .* Window;
Iref_blurred = fftshift( ifft2( fftshift(Iref_blurred) ) );
Iref_blurred = abs(Iref_blurred)+Imin;
In the next steps I think we are doing a Gaussian blur. However, I thought the kernel should be generated in the fourier domain in the same way that we could multiply them, as in the line:
Iref_blurred(yIndices,xIndices,:) = Iref_fft(yIndices,xIndices) .* Window;
I'm not sure if Window is the Fourier transform of a Gaussian convolution kernel or at least can't tell it from the code.
So, I'm a little confused as to how this is achieved with gaussian blur. Any help understanding this code would be appreciated.
source to share
You are correct that this code is missing a Gaussian FFT, but the thing to remember (or learn) is that the Fourier space representation of a Gaussian is also a Gaussian
, with only the reciprocal standard deviation. Whoever wrote this code probably knew about it, or just forgot and got lucky.
See the gausswin
docs section Gaussian Window and Fourier Transform . A compressed version of the gausswin
example in the documentation:
N = 64; n = -(N-1)/2:(N-1)/2; alpha = 8; w = gausswin(N,alpha); nfft = 4*N; freq = -pi:2*pi/nfft:pi-pi/nfft; wdft = fftshift(fft(w,nfft)); plot(n,w) hold on plot(freq/pi,abs(wdft) / 10,'r') title('Gaussian Window and FFT') legend({'win = gausswin(64,8)','0.1 * abs(FFT(win))'})
So, interpreting the output gausswin
as Fourier space at once, without execution and FFT, equates to a Gaussian window with a much larger sigma in the spatial domain.
source to share