Correct implementation of Gaussian blur via FFT

Question

I have read a lot of questions on SO about Gaussian blur and FFT, but there aren't answer how to implement steps of it (but there are comments like "it's your homework"). I want to know, how to properly pad the kernel and use FFT and IFFT on kernel and image. Can you provide some pseudocode or implementation in any language like Java, Python, etc. how to do this, or at least some good tutorial how to understand it:

1. FFT the image
2. FFT the kernel, padded to the size of the image
3. multiply the two in the frequency domain (equivalent to convolution in the spatial domain)
4. IFFT (inverse FFT) the result

Steps copied from Gaussian blur and FFT

Wyrick

Answer 1

A Matlab Example. It should be a good place to start for you.

Load the picture:

%Blur Demo

%Import image in matlab default image set.
origimage = imread('cameraman.tif');

%Plot image
figure, imagesc(origimage)
axis square
colormap gray
title('Original Image')
set(gca, 'XTick', [], 'YTick', [])

The whole process:

%Blur Kernel
ksize = 31;
kernel = zeros(ksize);

%Gaussian Blur
s = 3;
m = ksize/2;
[X, Y] = meshgrid(1:ksize);
kernel = (1/(2*pi*s^2))*exp(-((X-m).^2 + (Y-m).^2)/(2*s^2));

%Display Kernel
figure, imagesc(kernel)
axis square
title('Blur Kernel')
colormap gray

%Embed kernel in image that is size of original image
[h, w] = size(origimage);
kernelimage = zeros(h,w);
kernelimage(1:ksize, 1:ksize) = kernel;

%Perform 2D FFTs
fftimage = fft2(double(origimage));
fftkernel = fft2(kernelimage);

%Set all zero values to minimum value
fftkernel(abs(fftkernel)<1e-6) = 1e-6;

%Multiply FFTs
fftblurimage = fftimage.*fftkernel;

%Perform Inverse 2D FFT
blurimage = ifft2(fftblurimage);

%Display Blurred Image
figure, imagesc(blurimage)
axis square
title('Blurred Image')
colormap gray
set(gca, 'XTick', [], 'YTick', [])

Before Image:

After Image:

Note, because the zero filling isn't placing the kernel centered, you get an offset. This answer explains the wrapping issue. gaussian blur with FFT

Inflation

Great example code. I do have a few suggestions make this even better. (1) Show an example result using an actual image we can test with. You're using imread with a local image which we don't have access to. (2) Don't simply use the red channel to convert to grayscale. Use rgb2gray to do that for you. (3) If you have time, extend to colour images where you would perform the filtering separately per channel. (4) Setting all zero values to a minimum value with find is a bit inefficient. Use logical indexing instead: fftkernel(fftkernel == 0) = 1e-6; You don't need find here. – Polite