matlab code: Gaussian blurring in fourier domain

Question

I am looking at some code which performs blurring of images. However, I am having trouble understanding the code and I was wondering if someone could help me understand what the code is doing roughly.

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 the fftshift on an image which is not already in the fourier domain. So, I can tell it is doing a fourier transform along each of the axes but why does it do a fftshift before and after?

The code follows as:

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 subsequent steps, I think we are doing a Gaussian blurring. However, I thought the kernel has to be generated in the fourier domain as well before we could multiply them like the line:

Iref_blurred(yIndices,xIndices,:) = Iref_fft(yIndices,xIndices) .* Window;

I am not sure if the Window is the fourier transform of the Gaussian convolution kernel or at least not being able to tell it from the code.

So, I am a bit confused as to how this is achieving the Gaussian blurring. Any help in understanding this code would be appreciated.

Answer 1

You are correct that there is no FFT of the Gaussian going on in this code, but the thing to remember (or learn) is that the Fourier space representation of a Gaussian is also a Gaussian, just with the reciprocal standard deviation. Whoever wrote this code probably knew this, or they just forgot and got lucky.

See the section of gausswin docs called Gaussian Window and the Fourier Transform. Condensed version of gausswin example in 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 of gausswin as the Fourier space right away, without performing and FFT, equates to a Gaussian window in the with a much larger sigma in the spacial domain.