Fourier Transform in Medical Imaging with Python Implementation

Medical image processing and signal processing are huge tasks in radiology. You receive a lot of signals with great variability. These signals are used to test the health of the patient and uncover any ailment. All this is done by applying the Fourier transform.

This article will cover the Fourier transform concept and observe its implementation in the Python programming language. We will also observe how image processing is done using the Fourier transform.

Fourier transform is a mathematical tool widely used in medical imaging and signal processing to decompose signals into their constituent frequencies. It enables the analysis and manipulation of complex signals, facilitating tasks such as denoising, filtering, and feature extraction. Python provides powerful libraries like NumPy and SciPy that make it easy to implement Fourier transform algorithms and apply them to real-world problems in radiology and beyond.

What is Fourier Transform?

The Fourier transform is used to understand and decompose signals into their components and frequencies. After applying the Fourier transform, we receive a sinusoidal curve. Let us look at the formula for the Fourier transform.

The Fourier transform formula may look intimidating at first glance, but it essentially represents the relationship between a signal in the time domain and its representation in the frequency domain. The formula integrates the product of the signal and a complex exponential term over all time, yielding the Fourier transform of the signal. This transformation allows us to analyze the signal’s frequency components and perform various operations in the frequency domain.

Now, let us try to understand this concept using Python. We have decomposed a sample frequency of 1000Hz into Frequency Domain signal and magnitude.

import numpy as np
import matplotlib.pyplot as plt

# Generate a sample signal
fs = 1000  # Sampling frequency (Hz)
t = np.arange(0, 1, 1/fs)  # Time vector from 0 to 1 second
f1 = 10  # Frequency of the first sine wave (Hz)
f2 = 50  # Frequency of the second sine wave (Hz)
signal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)  # Sum of two sine waves

# Plot the sample signal
plt.figure(figsize=(10, 4))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('Time Domain Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

# Compute the Fourier transform
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), 1/fs)

# Plot the Fourier transform (magnitude spectrum)
plt.subplot(2, 1, 2)
plt.plot(frequencies, np.abs(fft_result))
plt.title('Frequency Domain Signal (Magnitude Spectrum)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, fs/2)  # Display only positive frequencies
plt.tight_layout()
plt.show()

Let us look at the output of the code above.

From the output above, we have decomposed it into the time domain and frequency domain signal. Now, let us move on to the next section, which explains how to apply the Fourier transform to image processing.

Fourier Transform for Image processing

We have applied the Fourier transform for image processing using the code below. Essentially, we have an image with some disturbance or noise in it. By applying the Fourier transform, we have removed the noise and thus received a cleaner image. For the code below, we can get images from Google.

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, color

# Step 1: Load the image
image = io.imread('lena.png')  # Load the image (you can replace 'lena.png' with your image file)

# Step 2: Convert the image to grayscale
gray_image = color.rgb2gray(image)

# Step 3: Compute the Fourier transform of the image
fft_image = np.fft.fft2(gray_image)

# Step 4: Shift the zero frequency component to the center
fft_image_shifted = np.fft.fftshift(fft_image)

# Step 5: Visualize the magnitude spectrum (optional)
magnitude_spectrum = np.abs(fft_image_shifted)
plt.imshow(np.log(1 + magnitude_spectrum), cmap='gray')
plt.title('Magnitude Spectrum')
plt.colorbar()
plt.show()

# Step 6: Create a mask to filter out high-frequency noise
rows, cols = gray_image.shape
center_row, center_col = rows // 2, cols // 2
mask_size = 50
mask = np.zeros((rows, cols))
mask[center_row - mask_size:center_row + mask_size, center_col - mask_size:center_col + mask_size] = 1

# Step 7: Apply the mask to the shifted Fourier transform
fft_image_filtered = fft_image_shifted * mask

# Step 8: Inverse Fourier transform to obtain the denoised image
denoised_image = np.abs(np.fft.ifft2(np.fft.ifftshift(fft_image_filtered)))

# Step 9: Plot the original and denoised images
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(gray_image, cmap='gray')
plt.title('Original Image')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(denoised_image, cmap='gray')
plt.title('Denoised Image')
plt.axis('off')

plt.show()

The above code is similar to the code we witnessed in the last section, with a slight change. Rather than generating random data, we take an image as an input.

Conclusion

Who knew that Fourier Transform could be so fascinating? We’ve seen how this mathematical tool can be useful for medical imaging and signal processing.

With Python by our side, we can easily put Fourier Transform into action. We can take a signal, break it into frequency components, and then manipulate it as required. As we continue to explore medical imaging and signal processing, Fourier transform can be a great companion.

Hope you enjoyed reading it!!