Digital Image Processing

Author: Sandipan Dey

Convolution and Frequency Domain Filtering

In this chapter we will continue with 2D convolution and understand how convolution can be done faster in the frequency domain (with basic concepts of the convolution theorem). We will see the basic differences between correlation and convolution with an example on an image. We will also describe an example from SciPy that will show how to find the location of specific patterns in an image with a template image using cross-correlationtion using kernels, such as box-kernel or Gaussian kernel) in the frequency domain, such as high-pass, low-pass, band-pass, and band-stop filters, and how to implement them with Python libraries by using examples. We shall demonstrate with examples how some filters can be used for image denoising (for example, the band-reject or notch filters to remove periodic noise from an image, or the inverse or Wiener filters to deblur an image that's blurred with a Gaussian / motion-blur kernel).

The topics to be covered in this chapter are as follows:

  1. Convolution theorem and frequency domain Gaussian blur
  2. Filtering in the frequency domain (with the SciPy ndimage module and scikit-image)

Convolution theorem and frequency domain Gaussian blur Filtering in the frequency domain (with the SciPy ndimage module and scikit-image)

Convolution theorem and frequency domain Gaussian blur

In this section, we will see more applications of convolution on images using Python modules such as scipy signal and ndimage. Let's start with convolution theorem and see how the convolution operation becomes easier in the frequency domain.

Application of the convolution theorem

The convolution theorem says that convolution in an image domain is equivalent to a simple multiplication in the frequency domain:

png

 

Following diagram shows the application of fourier transforms:

png

The next diagram shows the basic steps in frequency domain filtering. We have the original image, F, and a kernel (a mask or a degradation/enhancement function) as input. First, both input items need to be converted into the frequency domain with DFT, and then the convolution needs to be applied, which by convolution theorem is just an (element-wise) multiplication. This outputs the convolved image in the frequency domain, on which we need to apply IDFT to obtain the reconstructed image (with some degradation or enhancement on the original image):

png

 

Import libraries

some note

Frequency domain Gaussian blur filter with numpy fft

 

png

Gaussian kernel in the frequency domain

In this section, we will see how the Gaussian kernel looks like in the frequency domain in 2D and 3D plot.

The Gaussian LPF kernel spectrum in 2D

The next code block shows how to plot the spectrum of a Gaussian kernel in 2D with the log transform:

 

Since the Gaussian kernel is a low-pass filter, its frequency spectrum has higher values for the frequencies in the center (it allows more low-frequency values), and gradually decreases as one moves away from the center to the higher frequency values

png

Frequency domain Gaussian blur filter with scipy signal.fftconvolve()

The following code block shows how the SciPy signal module's fftconvolve() function can be used to run the convolution in the frequency domain (internally by just a multiplication and the convolution theorem):

 

png

The following code block shows how to plot the original and the blurred image spectrum after convolution:

 

png

The same code block as the previous one is used to compute a blurred version of the Victoria Memorial Hall image, but this time the fftconvolve() function is invoked with a larger (50 x 50) Gaussian kernel with σ=10. The following screenshot shows the output of the same code with the new image and the modified kernel

 

png

 

 

png

Comparing the runtimes of SciPy convolve() and fftconvolve() with the Gaussian blur kernel

We can use the Python timeit module to compare the runtimes of the image domain and the frequency domain convolution functions. Since the frequency domain convolution involves a single matrix multiplication instead of a series of sliding window arithmetic computations, it is expected to be much faster. The following code compares the runtimes:

 

 

png

The following code visualizes the difference in between runtimes. Each of the functions has been run on the same input image with the same Gaussian kernel 100 times and then the boxplot of the times taken for each function is plotted:

png

Filtering in the frequency domain (HPF, LPF, BPF, and notch filters)

If we remember from the image processing pipeline described in Chapter 1, Getting Started with Image Processing, the immediate next step after image acquisition is image pre-processing. Images are often **corrupted by random variations in intensity and illumination or have poor contrast**, and therefore can't be used directly and need to be enhanced. This is where the filters are used.

What is a filter?

High-Pass Filter (HPF)

This filter allows only high frequencies from the frequency domain representation of the image (obtained with DFT) and blocks all low frequencies beyond a cut-off value. The image is reconstructed with inverse DFT, and since the *high-frequency components correspond to edges, details, noise, and so on, HPFs tend to extract or enhance them. The next few sections demonstrate how to use different functions from the numpy, scipy, and scikit-image libraries to implement an HPF and the impact of an HPF on an image.

The Python code to implement an HPF is shown in the following code. As can be seen, the high frequency components correspond more to the edges of an image, and the average (flat) image information is lost as we get rid of more and more low frequency components: Mode in PIL, image.convert('RGB')

png

plot the original rhino image spectrum in the logarithm domain, as well, show the original rhino image spectrum on following screen

png

The following screenshot shows the spectrum of the rhino image after applying the HPF,the output of the previous code:

png

The following code block shows how to obtain the image back from the previous spectrum by applying ifft2():

png

The screenshot above shows the output of the preceding code: the rhino image after the HPF is applied. As can be seen, the edges in the image become more prominent—HPF finds the edges (which correspond to higher frequencies) in the image:

The following code block shows the application of the HPF on the cameraman grayscale image, with different frequency cut-off values, F:

png

How SNR changes with frequency cut-off

The following code block shows how to plot the change in the signal to noise ratio (SNR) with the cutoff frequency (F) for the HPF:

png

Low-pass filter (LPF)

This filter allows only the low frequencies from the frequency domain representation of the image (obtained with DFT), and blocks all high frequencies beyond a cut-off value. The image is reconstructed with inverse DFT, and since the high-frequency components correspond to edges, details, noise, and so on, LPF tends to remove these. The next few sections demonstrate how to use different functions from the numpy, scipy, and scikit-image libraries to implement an LPF and the impact of an LPF on an image.

The next code block demonstrates how to use the LPF (weighted average filter) to blur the lena grayscale image:

LPF with fourier_gaussian

This function from the scipy ndimage module implements a multi-dimensional Gaussian Fourier filter. The frequency array is multiplied with the Fourier transform of a Gaussian kernel of a given size.

 

png

The preceding code block displays the frequency spectrum of the image after the box kernel is applied:

png

LPF with scipy fftpack

We can implement a LPF on an image with the following steps:

  1. Perform a 2D FFT with scipy.fftpack fft2 and obtain the frequency domain representation of the image
  2. Keep only the low frequency components (get rid of the high frequency components)
  3. Perform an inverse FFT to reconstruct the image

The Python code to implement an LPF is shown in the following code. As can be seen from the next screenshots, the high frequency components correspond more to the average (flat) image information and the detailed information (for example, the edges) of an image is lost as we get rid of more and more high-frequency components.

For example, if we keep only the first-frequency component and discard all others, in the resulting image obtained after inverse FFT we can hardly see the rhinos, but as we keep more and more higher frequencies, they become prominent in the final image:

png

The following code block shows how to plot the spectrum of the image in the logarithm domain after blocking the high frequencies; in other words, only allowing the low frequencies:

png

The following code block shows the application of LPF on the cameraman grayscale image, with different frequency cutoff values, F:

png

How SNR changes with frequency cutoff

The following code block shows how to plot the change in signal to noise ratio (SNR) with cutoff frequency (F) for LPF:

 

png

The picture above shows how SNR of the output image decreases with the increase in the cut-off frequency of LPF. The red horizontal line indicates the original image's SNR, plotted for comparison

Band-pass filter (BPF) with DoG

The Difference of Gaussian (DoG) kernel can be used as a BPF to allow the frequencies in a certain band and discard all other frequencies. The following code block shows how the DoG kernel can be used with fftconvolve() to implement a BPF:

 

png

png

Band-stop (notch) filter

This filter blocks/rejects a few chosen frequencies from the frequency domain representation of the image (obtained with DFT), and hence the name. It is useful for removing periodic noise from images, as discussed in the next section.

Using a notch filter to remove periodic noise from images

In this example, we will first add some periodic (sinusoidal) noise to the parrot image to create a noisy parrot image (this can happen because of interference with some electrical signal) and then observe the effect of the noise in the frequency domain of the image using the following code block:

 

png

Now let's design a band-stop/band-reject (notch) filter to eliminate the frequencies that are responsible for noise by setting the corresponding frequency components to zero in the next code block:

png

The screenshot shows the output of the preceding code block, the restored image by applying the notch filter. As can be seen, the original image looks sharper than the restored one, since some true frequencies from the original image are also rejected by the band-reject filter along with the noise.

Image restoration

In image restoration, the degradation is modeled. This enables the effects of the degradation to be (largely) removed. The challenge is the loss of information and noise. The following diagram shows the basic image degradation model:

png

 

In the next few sections, we will describe a couple of degradation models (namely _***inverse***_ and Wiener filters).

Deconvolution and inverse filtering with FFT

Given a blurred image with a known (assumed) blur kernel, a typical image processing task is to get back (at least an approximation of) the original image. This particular task is known as deconvolution. One of the naive filters that can be applied in the frequency domain to achieve this is the inverse filter that we are going to discuss in this section.

Let's first start by blurring the grayscale lena image with Gaussian blur using the following code:

png

Now we can use the inverse filter on the blurred image (using the same H) to get back the original image. The following code block demonstrates how to do it:

png

Please refer textbook for an advanced discussion, exercise and question 4 in the Question section.

Skip next example for motion deblur

Similarly, we can deblur an image blurred with a known motion blur kernel using the inverse filter. The code remains the same; only the kernel changes, as shown in the following code. Note that we need to create a zero-padded kernel of a size equal to the size of the original image before we can apply the convolution in the frequency domain (using np.pad(); the details are left as an exercise for you):

 

png

 

Image deconvolution with the Wiener filter

We already saw how to to obtain the (approximate) original image from the blurred image (with a known blur kernel) using the inverse filter in the last section. Another important task in image processing is the removal of noise from a corrupted signal. This is also known as image restoration. The following code block shows how the scikit-image restoration module's unsupervised Wiener filter can be used for image denoising with deconvolution:

png

The screenshot shows the output of the previous code—the original, blurred-noisy, and restored image using the unsupervised Wiener filter

Image denoising with FFT

The next example is taken from plot_fft_image_denoise. This example demonstrates how to denoise an image first by blocking Fourier elements with high frequencies using a LPF with FFT. Let's first display the noisy grayscale image with the following code block:

png

 

The screenshot shows the output of the preceding code block, the original noisy image

The following code block displays the frequency spectrum of the noisy image:

png

The screenshot shows the output of the preceding code, the Fourier spectrum of the original noisy image.

Filter in FFT

The following code block shows how to reject a bunch of high frequencies and implement an LPF to attenuate noise (which corresponds to high-frequency components) from the image:

 

png

The above screenshot shows the output of the preceding code, the filtered spectrum with the LPF.

Reconstructing the final image

The following code block shows how to reconstruct the image with IFFT from the filtered Fourier coefficients:

 

png

The above screenshot shows the output of the preceding code—a much cleaner output image obtained from the original noisy image with filtering in the frequency domain

Summary

In this chapter, we discussed a few important concepts primarily related to 2D convolution and its related applications in image processing, such as filtering in the spatial domain. We also discussed a few different frequency domain filtering techniques and illustrated them with quite a few examples with the scikit-image numpy fft, scipy, fftpack, signal, and ndimage modules. We started with the convolution theorem and its application in frequency domain filtering, various frequency domain filters such as LPF, HPF, and notch filters, and finally deconvolution and its application in designing filters for image restoration, such as inverse and Wiener filters.

On completion of this chapter, the reader should be able to write Python code to do 2D convolution/filtering, and should also be able to write Python code to implement time/frequency domain filters with/without convolution.

In the next chapter, we will start on different image enhancement techniques based on the concepts introduced in the last two chapters.

Further reading

  1. https://www.cs.cornell.edu/courses/cs1114/2013sp/sections/S06_convolution.pdf
  2. http://www.aip.de/groups/soe/local/numres/bookcpdf/c13-3.pdf
  3. http://www.cse.usf.edu/~r1k/MachineVisionBook/MachineVision.files/MachineVision_Chapter4.pdf
  4. https://web.stanford.edu/class/ee367/slides/lecture6.pdf
  5. https://pdfs.semanticscholar.org/presentation/50e8/fb095faf6ed51e03c85a2fcb7eb1ae1b1009.pdf
  6. http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf