Python image processing: frequency domain filtering, noise reduction and image enhancement

Click “Xiaobai Xue Vision” above and choose to add “star” or “pin“

Heavy stuff, delivered as soon as possible

This article is about 4300 words, it is recommended to read for 8 minutes
This article will discuss the various stages involved in the frequency transformation of an image from FFT to inverse FFT, combined with the use of FFT shift and inverse FFT shift.

Image processing has become an integral part of our daily lives, involving various fields such as social media and medical imaging. Images obtained from digital cameras or other sources such as satellite photos and medical scans may require preprocessing to remove or enhance noise. Frequency domain filtering is a possible solution that can remove noise while enhancing image sharpening.

Fast Fourier Transform (FFT) is a mathematical technique that transforms images from the spatial domain to the frequency domain. It is a key tool for frequency transformation in image processing. By leveraging a frequency domain representation of an image, we can effectively analyze the image based on its frequency content, thereby simplifying the application of filtering procedures to remove noise.

This article uses three Python libraries, namely openCV, Numpy and Matplotlib.

import cv2
 import numpy as np
 from matplotlib import pyplot as plt
 img = cv2.imread('sample.png',0) # Using 0 to read image in grayscale mode
 plt.imshow(img, cmap='gray') #cmap is used to specify imshow that the image is in greyscale
 plt.xticks([]), plt.yticks([]) # remove tick marks
 plt.show()

1. Fast Fourier Transform (FFT)

Fast Fourier transform (FFT) is a widely used mathematical technique that allows images to be transformed from the spatial domain to the frequency domain and is a basic component of frequency transformation. Using FFT analysis, the periodicity of the image can be obtained and divided into different frequency components to generate an image spectrum that displays the amplitude and phase of each image’s respective frequency component.

f = np.fft.fft2(img) #the image 'img' is passed to np.fft.fft2() to compute its 2D Discrete Fourier transform f
 mag = 20*np.log(np.abs(f))
 plt.imshow(mag, cmap = 'gray') #cmap='gray' parameter to indicate that the image should be displayed in grayscale.
 plt.title('Magnitude Spectrum')
 plt.xticks([]), plt.yticks([])
 plt.show()

The code above calculates the magnitude of the Fourier transform f using np.abs(), converts to a logarithmic scale with np.log(), and then multiplies by 20 to get the magnitude in decibels. This is done to make the amplitude spectrum easier to visualize and interpret.

2. FFT displacement

In order for the filtering algorithm to be applied to the image, the zero frequency component of the image is moved to the center of the spectrum using FFT shifting.

fshift = np.fft.fftshift(f)
 mag = 20*np.log(np.abs(fshift))
 plt.imshow(mag, cmap = 'gray')
 plt.title('Centered Spectrum'), plt.xticks([]), plt.yticks([])
 plt.show()

3. Filtering

One purpose of frequency transformation is to use various filtering algorithms to reduce noise and improve image quality. The two most commonly used image sharpening filters are the Ideal high-pass filter and the Gaussian high-pass filter. These filters all use the frequency domain representation of the image obtained through the Fast Fourier Transform (FFT) method.

An ideal high-pass filter is an infinitely long filter with infinite frequency bandwidth and ideal passband and stopband responses. Signals of all frequencies within its passband are completely transmitted, while signals of all frequencies within its stopband are completely suppressed.

In the frequency domain, the amplitude-frequency response of an ideal filter is:

Within the passband, the amplitude-frequency response is 1
Within the stopband, the amplitude-frequency response is 0

In the time domain, the impulse response of an ideal filter is:

Within the passband, the impulse response is an infinitely long sequence of unit impulse functions.
Within the stopband, the impulse response is zero

Since an ideal filter has infinite bandwidth in the frequency domain, it cannot be implemented in practical applications. The digital filters used in practice are usually based on the approximation of the ideal filter, so they are just an Ideal.

Gaussian high-pass filter is a filter commonly used in digital image processing. Its function is to retain high-frequency detail information in the image and suppress low-frequency signals. This filter is based on a Gaussian function and has a smooth frequency response that can adapt to a variety of image details.

The frequency response of a Gaussian high-pass filter can be expressed as:

H(u,v) = 1 – L(u,v)

Among them, L(u,v) is a low-pass filter, which can be represented by a Gaussian function. The response of a high-pass filter is obtained by subtracting the response of the low-pass filter from 1. In practice, different parameter settings are usually used to adjust the Gaussian function to achieve different filtering effects.

Disk-shaped images are used to define the frequency components to be retained or suppressed when Fourier transforming the image. In this context, an ideal filter usually refers to an ideal low-pass or high-pass filter that can be selected in the frequency domain to retain or suppress signals within a specific frequency range. After applying this ideal filter to the Fourier transform of the image, and then performing the inverse transformation, the image processed by the filter can be obtained.

We won’t introduce the specific details, let’s look directly at the code:

The following functions are circular masks that generate ideal high-pass and low-pass filters.

import math
 def distance(point1,point2):
     return math.sqrt((point1[0]-point2[0])**2 + (point1[1]-point2[1])**2)


 def idealFilterLP(D0,imgShape):
     base = np.zeros(imgShape[:2])
     rows, cols = imgShape[:2]
     center = (rows/2,cols/2)
     for x in range(cols):
         for y in range(rows):
             if distance((y,x),center) < D0:
                 base[y,x] = 1
     return base


 def idealFilterHP(D0,imgShape):
     base = np.ones(imgShape[:2])
     rows, cols = imgShape[:2]
     center = (rows/2,cols/2)
     for x in range(cols):
         for y in range(rows):
             if distance((y,x),center) < D0:
                 base[y,x] = 0
     return base

The following function generates a Gaussian high and low pass filter with the required circular mask

import math
 def distance(point1,point2):
     return math.sqrt((point1[0]-point2[0])**2 + (point1[1]-point2[1])**2)


 def gaussianLP(D0,imgShape):
     base = np.zeros(imgShape[:2])
     rows, cols = imgShape[:2]
     center = (rows/2,cols/2)
     for x in range(cols):
         for y in range(rows):
             base[y,x] = math.exp(((-distance((y,x),center)**2)/(2*(D0**2))))
     return base


 def gaussianHP(D0,imgShape):
     base = np.zeros(imgShape[:2])
     rows, cols = imgShape[:2]
     center = (rows/2,cols/2)
     for x in range(cols):
         for y in range(rows):
             base[y,x] = 1 - math.exp(((-distance((y,x),center)**2)/(2*(D0**2))))
     return base

Filter example

fig, ax = plt.subplots(2, 2) # create a 2x2 grid of subplots
 fig.suptitle('Filters') # set the title for the entire figure


 # plot the first image in the top-left subplot
 im1 = ax[0, 0].imshow(np.abs(idealFilterLP(50, img.shape)), cmap='gray')
 ax[0, 0].set_title('Low Pass Filter of Diameter 50 px')
 ax[0, 0].set_xticks([])
 ax[0, 0].set_yticks([])


 # plot the second image in the top-right subplot
 im2 = ax[0, 1].imshow(np.abs(idealFilterHP(50, img.shape)), cmap='gray')
 ax[0, 1].set_title('High Pass Filter of Diameter 50 px')
 ax[0, 1].set_xticks([])
 ax[0, 1].set_yticks([])


 # plot the third image in the bottom-left subplot
 im3 = ax[1, 0].imshow(np.abs(gaussianLP(50,img.shape)), cmap='gray')
 ax[1, 0].set_title('Gaussian Filter of Diameter 50 px')
 ax[1, 0].set_xticks([])
 ax[1, 0].set_yticks([])


 # plot the fourth image in the bottom-right subplot
 im4 = ax[1, 1].imshow(np.abs(gaussianHP(50,img.shape)), cmap='gray')
 ax[1, 1].set_title('Gaussian Filter of Diameter 50 px')
 ax[1, 1].set_xticks([])
 ax[1, 1].set_yticks([])


 # adjust the spacing between subplots
 fig.subplots_adjust(wspace=0.5, hspace=0.5)


 # save the figure to a file
 fig.savefig('filters.png', bbox_inches='tight')

Multiply the filter and shifted image to get the filtered image

In order to obtain the final filtered image with the desired frequency response, the key is to perform a point-wise multiplication of the shifted image with the filter in the frequency domain.

This process multiplies the corresponding pixels of two image elements. For example, when applying a low-pass filter, we multiply the shifted Fourier transform image point-by-point with the low-pass filter.

This operation suppresses high frequencies and preserves low frequencies, and vice versa for a high-pass filter. This multiplication process is essential to remove unwanted frequencies and enhance desired frequencies, resulting in a sharper and sharper image.

It allows us to obtain the desired frequency response and obtain the final filtered image in the frequency domain.

4. Multiplying Filter and Shifted Image

The multiplicative filter is a filter with pixel values as weights. It obtains the filtered pixel values by multiplying the weight of the filter with the pixel values of the image. Specifically, assuming that the weight of the multiplication filter is h(i,j) and the pixel value of the image is f(m,n), then the filtered pixel value g(x,y) can be expressed as:

g(x,y) = ∑∑ f(m,n)h(x-m,y-n)

Among them, ∑∑ means summing all (m,n).

The translated image refers to the result of a translation operation on the image. The translation operation usually refers to translating the pixels of the image along the x-axis and y-axis directions. The translated images have the same size and resolution as the original image, but their pixel positions have changed. In the filtering operation, the translated image can be used for convolution with the filter to implement the filtering operation.

It allows us to obtain the desired frequency response and obtain the final filtered image in the frequency domain.

fig, ax = plt.subplots()
 im = ax.imshow(np.log(1 + np.abs(fftshifted_image * idealFilterLP(50,img.shape))), cmap='gray')
 ax.set_title('Filtered Image in Frequency Domain')
 ax.set_xticks([])
 ax.set_yticks([])


 fig.savefig('filtered image in freq domain.png', bbox_inches='tight')

When visualizing the Fourier spectrum, the choice between using np.log(1 + np.abs(x)) and 20*np.log(np.abs(x)) is a matter of personal preference and can depend on the specific s application.

Typically Np.log (1 + np.abs(x)) is used because it helps visualize the spectrum more clearly by compressing the dynamic range of the data. This is accomplished by taking the logarithm of the absolute value of the data and adding 1 to avoid taking the logarithm of zero.

And 20*np.log(np.abs(x)) scales the data by 20 times and takes the logarithm of the absolute value of the data, which makes it easier to see the smaller amplitude differences between different frequencies. But it doesn’t compress the dynamic range of the data like np.log(1 + np.abs(x)) does.

Both methods have their own advantages and disadvantages, which ultimately depend on the specific application and personal preference.

5. Inverse FFT displacement

After frequency domain filtering, we need to move the image back to its original position and then apply the inverse FFT. To achieve this, an inverse FFT shift is used, which reverses the shift process performed previously.

fig, ax = plt.subplots()
 im = ax.imshow(np.log(1 + np.abs(np.fft.ifftshift(fftshifted_image * idealFilterLP(50,img.shape)))), cmap='gray')
 ax.set_title('Filtered Image inverse fft shifted')
 ax.set_xticks([])
 ax.set_yticks([])


 fig.savefig('filtered image inverse fft shifted.png', bbox_inches='tight')

6. Inverse Fast Fourier Transform (IFFT)

Inverse Fast Fourier Transform (IFFT) is the final step in image frequency transformation. It is used to transfer images from the frequency domain back to the spatial domain. The result of this step is an image with reduced noise and improved clarity compared to the original image in the spatial domain.

fig, ax = plt.subplots()
 im = ax.imshow(np.log(1 + np.abs(np.fft.ifft2(np.fft.ifftshift(fftshifted_image * idealFilterLP(50,img.shape))))), cmap='gray')
 ax.set_title('Final Filtered Image In Spatial Domain')
 ax.set_xticks([])
 ax.set_yticks([])


 fig.savefig('final filtered image.png', bbox_inches='tight')

Summary

We then string together all the operations and display them, and the function draws all the images.

def Freq_Trans(image, filter_used):
     img_in_freq_domain = np.fft.fft2(image)
 
     # Shift the zero-frequency component to the center of the frequency spectrum
     centered = np.fft.fftshift(img_in_freq_domain)
 
     # Multiply the filter with the centered spectrum
     filtered_image_in_freq_domain = centered * filter_used
 
     # Shift the zero-frequency component back to the top-left corner of the frequency spectrum
     inverse_fftshift_on_filtered_image = np.fft.ifftshift(filtered_image_in_freq_domain)
 
     # Apply the inverse Fourier transform to obtain the final filtered image
     final_filtered_image = np.fft.ifft2(inverse_fftshift_on_filtered_image)
 
     return img_in_freq_domain,centered,filter_used,filtered_image_in_freq_domain,inverse_fftshift_on_filtered_image,final_filtered_image

Use high-pass, low-pass ideal filters and Gaussian filters with diameters of 50, 100 and 150 pixels respectively to demonstrate their impact on enhancing image clarity.

fig, axs = plt.subplots(12, 7, figsize=(30, 60))
 
 filters = [(f, d) for f in [idealFilterLP, idealFilterHP, gaussianLP, gaussianHP] for d in [50, 100, 150]]
 
 for row, (filter_name, filter_diameter) in enumerate(filters):
     # Plot each filter output on a separate subplot
     result = Freq_Trans(img, filter_name(filter_diameter, img.shape))
 
     for col, title, img_array in zip(range(7),
    ["Original Image", "Spectrum", "Centered Spectrum", f"{filter_name.__name__} of Diameter {filter_diameter} px",
     f"Centered Spectrum multiplied by {filter_name.__name__}", "Decentralize", "Processed Image"],
    [img, np.log(1 + np.abs(result[0])), np.log(1 + np.abs(result[1])), np.abs(result[2]), np.log (1 + np.abs(result[3])),
      np.log(1 + np.abs(result[4])), np.abs(result[5])]):
         
         axs[row, col].imshow(img_array, cmap='gray')
         axs[row, col].set_title(title)
 
 plt.tight_layout()
 plt.savefig('all_processess.png', bbox_inches='tight')
 plt.show()

It can be seen that when changing the diameter of the circular mask, the impact on image sharpness is different. As the diameter increases, more frequencies are suppressed, resulting in smoother images with less detail. Reducing the diameter allows more frequencies to pass through, resulting in sharper images and more detail. In order to achieve the desired effect, it is important to choose the right diameter, as using a diameter that is too small will result in a filter that is not effective enough, while using a diameter that is too large will result in too much detail being lost.

Generally speaking, Gaussian filters are more commonly used in image processing tasks due to their smoothness and robustness. In some applications, where a sharper cutoff is required, an ideal filter may be more suitable.

Using FFT to modify the image frequency is an effective method to reduce noise and improve image sharpness. This involves using an FFT to convert the image to the frequency domain, using appropriate techniques to filter the noise, and using an inverse FFT to convert the modified image back to the spatial domain. By understanding and implementing these techniques, we can improve image quality for a variety of applications.

Editor: Yu Tengkai

Proofreading: Lin Yilin

Download 1: OpenCV-Contrib extension module Chinese version tutorial

Reply in the background of the "Xiaobai Xue Vision" public account: Chinese tutorial on extension module, you can download the first Chinese version of OpenCV extension module tutorial on the entire network, covering extension module installation, SFM algorithm, stereo vision, target tracking, biological vision, ultrasonic vision Resolution processing and more than twenty chapters.


Download 2: Python visual practical project 52 lectures
Reply in the background of the "Xiaobai Xue Vision" public account: Python visual practical projects, you can download them, including image segmentation, mask detection, lane line detection, vehicle counting, adding eyeliner, license plate recognition, character recognition, emotion detection, text content extraction, 31 practical vision projects, including facial recognition, support rapid school computer vision.


Download 3: 20 lectures on OpenCV practical projects
Reply in the background of the "Xiaobai Xue Vision" public account: OpenCV practical projects 20 lectures, you can download 20 practical projects based on OpenCV to achieve advanced OpenCV learning.


Communication group

Welcome to join the public account reader group to communicate with peers. There are currently WeChat groups for SLAM, 3D vision, sensors, autonomous driving, computational photography, detection, segmentation, recognition, medical imaging, GAN, algorithm competitions, etc. (will be gradually subdivided in the future). Please scan the WeChat ID below to join the group, and note: "nickname + school/company + research direction", for example: "Zhang San + Shanghai Jiao Tong University + visual SLAM". Please note according to the format, otherwise it will not be approved. After successful addition, you will be invited to join the relevant WeChat group according to the research direction. Please do not send advertisements in the group, otherwise you will be asked to leave the group, thank you for your understanding~

The knowledge points of the article match the official knowledge archives, and you can further learn relevant knowledge. Python entry skill treeArtificial intelligenceMachine learning toolkit Scikit-learn384379 people are learning the system