Issue with lowpass image filter

Question

I need to apply a simple low-pass filter on an image, however this needs to be done in the frequency domain. However the final result is a very noisy image, suggesting that my process is incorrect somewhere.

My track of thought has been

1) Center the image (Testpatt) and apply a 2D FFT

tpa1 = size(Testpatt, 1); tpa2 = size(Testpatt, 2);
dTestpatt = double(Testpatt);


for i = 1:tpa1
    for j = 1:tpa2

        dTestpatt(i,j) = (dTestpatt(i,j))*(-1)^(i + j);
    end
end

fftdTestpatt = fft2(dTestpatt);

2) Generate the low-pass filter and multiply it with the Fourier transformed image (Note: the low pass filter needs to be a circle of 1's within the radius Do)

lowpfilter = zeros(tpa1, tpa2);
Do = 120;

for i = 1:tpa1
    for j = 1:tpa2


       if sqrt((i - tpa1/2)^2 + (j - tpa2/2)^2) <= Do
           lowpfilter(i,j) = 1;
       end

   end
end


filteredmultip = lowpfilter*fftdTestpatt;

3) Take the real components of the inverse 2D FFT and revert the centering.

filteredimagecent = real(ifft2(filteredmultip));


for i = 1:tpa1
    for j = 1:tpa2

        filteredimage(i,j) = (filteredimagecent(i,j))*(-1)^(i + j);
    end
end

Any help or comments would be greatly appreciated.

Answer 1

I'm surprised why this didn't give you an error, but you are performing matrix multiplication, not element-wise multiplication in the frequency domain. Recall that filtering in the frequency domain requires element-wise multiplication of the two transformed signals to perform the equivalent covnvolution operation in spatial domain. You simply need to change the multiplication statement so that it is the element-wise equivalent:

filteredmultip = lowpfilter.*fftdTestpatt;

Also, make sure you convert your image back to the same data type as the original image before display. If your image was uint8 for example, you'll want to use im2uint8 to convert it. This is important primarily for display purposes and for writing images to file. If you left it as double as you've done in your code, showing the image would be visualized as mostly white because displaying images of type double assumes the range of values is from [0,1].

As a side note, I suspect the reason why you aren't getting an error is because your image and filter are square, so the dimensions when it comes to matrix multiplication are valid. If you decide to do this on non-square images, you will definitely get an error.

Warning - Use of an ideal lowpass filter

What you are implementing is filtering with an ideal lowpass filter, so what will happen is that you will see ringing effects when you lowpass filter. The reason why is because this comes straight from signal processing theory. It requires a large (or rather infinite...) number of sinusoids in the spatial domain to realize a hard edge in the frequency domain. Remember that the FFT is a decomposition of your signal in terms of sinusoids. When you use this lowpass filter to filter your image, this is visualized as ringing in the reconstructed image as hard edges need a large summation of sinusoids (hence the fact the ringing) to create them.

To show you these effects, let's demonstrate for an example image. I'm going to use the cameraman image that's part of the image processing toolbox. If I used a radius of Do = 40 and ran your code (corrected), this is the original image and this is what I get after filtering:

As you can see, that's pretty bad. The ringing comes from the hard edges of your filter in the frequency domain. People tend to use a more gradual decrease in magnitude as you move farther away from the centre. A good example is a Gaussian blur. What you would do instead is define a standard deviation and then create a mask that is the same size as your image that represents a 2D Gaussian.

Recall that the 2D Gaussian can be represented as:

^{Source: Wikipedia}

You simply loop over all pixel coordinates and compute the above equation. I didn't multiply by the scaling factor in the front as you really don't need to do it. As such, you can use this code to generate a Gaussian mask, which is equivalent to what you have with your two for loops but doing it more vectorized. I define a grid of 2D coordinates that span the same size as your image, then run through the equation for each location in the image. We of course need to centre the kernel so you must offset by the middle of the image as you have done with your previous lowpass filtering code. I've also decided to use your Do variable as the standard deviation.

Do = 40;
[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = exp(-((X - (tpa2/2)).^2 + (Y - (tpa1/2)).^2) / (2*Do*Do));

So we now get this as the filter response:

... and when I filter, I get:

Compare with the original image:

As you can see, the blur is much better. No ringing. As such, when you filter images, make sure you avoid hard edges in the filter as this will produce ringing artifacts in the spatial domain.

Some more suggestions

I have a few more suggestions to give you to make this code run fast.

Avoid using loops to centre the image

As you already know from signal processing theory, if you pre-multiply the pixel values in the image by (-1)^(i+j) where (i,j) is a spatial location of the pixel you want to filter, this centres your image in the frequency domain. I would actually avoid using loops to do this and take the FFT first then centre the image. You can use a function called fftshift that performs the centering for you. First take the FFT of your image, then invoke this function after:

fftdTestpatt = fft2(dTestpatt);
fftdTestpatt = fftshift(fftdTestpatt); % Centre the FFT

Avoid using loops for generating the filter

I would also avoid using loops to generate your filter. Specifically, for your code with the ideal filter, replace your code with this which follows the same logic as what I had with the Gaussian filter. We can also remove the sqrt operation and make sure that you're comparing with the squared of the radius to avoid any unnecessary computation:

[X,Y] = meshgrid(1:tpa2, 1:tpa1);
lowpfilter = (X - tpa1/2).^2 + (Y - tpa2/2).^2 <= Do^2;
lowpfilter = double(lowpfilter);

Take special care of the element-wise power operation (.^) in the above code. The last statement is important as we need to cast the result back to double as generating the filter first would in fact give you a logical type as the result. We need the double type to ensure the precision of the FFT is respected.

Avoid loops to uncentre the image after filtering

After you're finished filtering, you again multiply by (-1)^(i+j) to uncentre the image. Use the related ifftshift function to uncentre the image after you filter with the FFT:

filteredmultip = ifftshift(filteredmultip); % Uncentre first
filteredimage = real(ifft2(filteredmultip));