1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ For this size of a kernel for a single image you better, performance wise, apply the filtration in the Spatial Domain. $\endgroup$ – Royi Oct 18 at 4:45
1
$\begingroup$

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. 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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.