Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB

Question

I have created a function that filters an image (250x250) with a gaussian blur kernel (5x5) using FFT and IFFT. I am trying to get my filtered image to equal exactly the filtered image created by the 'imfilter' function using the 'replicate' option. There are small differences in the two filtered images. They look identical visually but upon inspecting the pixel values, they are slightly different.

This is the general logic my code follows:

Zero pad the borders of the input image to size 254x254

Zero pad the kernel to size 254x254

'clamp to edge' the edge pixels of the input image. The outermost input image matrix rows and columns are extended to the borders, replacing the zeros. I am doing this to emulate the 'replicate' option in 'imfilter'.

Use 'fft2' to convert both the input image and kernel to frequency domain

Element wise multiply the two together

Convert the result to the time domain using 'ifft2'.

Crop away the border pixels, converting the filtered image from 254x254 ---> 250x250

What am I doing wrong here? Thanks.

Answer 1

The trick is to properly compensate for the fact that Frequency Domain multiplication applies a convolution with the circular boundary conditions in the spatial domain.

You may use the following code:

clear('all');
close('all');

gaussianKernelStd       = 2;
gaussianKernelRadius    = ceil(5 * gaussianKernelStd);

mI = im2double(imread('cameraman.tif'));
mI = mI(:, :, 1);

numRows = size(mI, 1);
numCols = size(mI, 2);

vX = [-gaussianKernelRadius:gaussianKernelRadius].';
vK = exp(-(vX .* vX) ./ (2 * gaussianKernelStd * gaussianKernelStd));
mK = vK * vK.';
mK = mK ./ sum(mK(:)); %<! The Gaussian Kernel

mIPad = padarray(mI, [gaussianKernelRadius, gaussianKernelRadius], "replicate", 'both'); %<! For the replicate array

mKC = CircularExtension2D(mK, size(mIPad, 1), size(mIPad, 2)); %<! Circular extension for the 2D Kernel

startIdx = gaussianKernelRadius + 1;

mIFiltered      = ifft2(fft2(mIPad) .* fft2(mKC), 'symmetric');
mIFiltered      = mIFiltered(startIdx:(startIdx + numRows - 1), startIdx:(startIdx + numCols - 1)); %<! Removing the padding
mIFilteredRef   = imfilter(mI, mK, 'replicate', 'same', 'conv'); %<! Reference

figure(); imshow(mIFiltered);
figure(); imshow(mIFilteredRef);

max(abs(mIFilteredRef(:) - mIFiltered(:))) %<! Should be very very low

The function CircularExtension2D() is given In my StackExchange Signal Processing Q38542 GitHub Repository (Look at the SignalProcessing\Q38542 folder). It was taken from my answer to Applying Image Filtering (Circular Convolution) in Frequency Domain.

The steps the code implements are as following:

1. Pad the image in order to have Replicate boundary condition convolution.
2. Convert the spatial domain kernel into a form which matches the image in frequency domain. We assume top left of the image is (0, 0) in spatial domain. So we need the (0, 0) of the kernel to also be in the top left corner.
3. Apply circular convolution using frequency domain.

As you can see, the result is prefect.

In my answer to How Much Zero Padding Do We Need to Perform Filtering in the Fourier Domain? I implemented a MATLAB Function, ImageFilteringFrequencyDomain(), to apply Frequency Domain Convolution with the border conditions supported in imfilter().

Remark
In your case, the kernel is 5x5 which is very small.
For small kernels and a single image it is better to apply the convolution in the spatial domain. It will be much faster.