圖像處理已經成為我們日常生活中不可或缺的一部分，涉及到社交媒體和醫學成像等各個領域。通過數位相機或衛星照片和醫學掃描等其他來源獲得的圖像可能需要預處理以消除或增強噪聲。頻域濾波是一種可行的解決方案，它可以在增強圖像銳化的同時消除噪聲。
快速傅立葉變換(FFT)是一種將圖像從空間域變換到頻率域的數學技術，是圖像處理中進行頻率變換的關鍵工具。通過利用圖像的頻域表示，我們可以根據圖像的頻率內容有效地分析圖像，從而簡化濾波程序的應用以消除噪聲。本文將討論圖像從FFT到逆FFT的頻率變換所涉及的各個階段，並結合FFT位移和逆FFT位移的使用。
本文使用了三個Python庫，即openCV、numpy和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、快速傅立葉變換(FFT)
快速傅立葉變換(FFT)是一種廣泛應用的數學技術，它允許圖像從空間域轉換到頻率域，是頻率變換的基本組成部分。利用FFT分析，可以得到圖像的周期性，並將其劃分為不同的頻率分量，生成圖像頻譜，顯示每個圖像各自頻率成分的振幅和相位。
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()上面代碼使用np.abs()計算傅立葉變換f的幅度，np.log()轉換為對數刻度，然後乘以20得到以分貝為單位的幅度。這樣做是為了使幅度譜更容易可視化和解釋。
2、FFT位移
為了使濾波算法應用於圖像，利用FFT移位將圖像的零頻率分量被移動到頻譜的中心
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、過濾
頻率變換的的一個目的是使用各種濾波算法來降低噪聲和提高圖像質量。兩種最常用的圖像銳化濾波器是Ideal high-pass filter 和Gaussian high-pass filter。這些濾波器都是使用的通過快速傅立葉變換(FFT)方法獲得的圖像的頻域表示。
Ideal high-pass filter（理想濾波器） 是一種無限長的、具有無限頻帶寬和理想通帶和阻帶響應的濾波器。其通帶內所有頻率的信號都被完全傳遞，而阻帶內所有頻率的信號則完全被抑制。
在頻域上，理想濾波器的幅頻響應為：
- 在通帶內，幅頻響應為 1
- 在阻帶內，幅頻響應為 0
在時域上，理想濾波器的衝激響應為：
- 在通帶內，衝激響應為一個無限長的單位衝激函數序列
- 在阻帶內，衝激響應為零
由於理想濾波器在頻域上具有無限帶寬，因此它無法在實際應用中實現。實際中使用的數字濾波器通常是基於理想濾波器的逼近，所以才被成為只是一個Ideal。
高斯高通濾波器（Gaussian high-pass filter）是一種在數字圖像處理中常用的濾波器。它的作用是在圖像中保留高頻細節信息，並抑制低頻信號。該濾波器基於高斯函數，具有光滑的頻率響應，可以適應各種圖像細節。
高斯高通濾波器的頻率響應可以表示為：
H(u,v) = 1 - L(u,v)
其中，L(u,v)是一個低通濾波器，它可以用高斯函數表示。通過將低通濾波器的響應從1中減去，可以得到一個高通濾波器的響應。在實際中，通常使用不同的參數設置來調整高斯函數，以達到不同的濾波效果。
圓形掩膜（disk-shaped images）是用於定義在圖像中進行傅立葉變換時要保留或抑制的頻率分量。在這種情況下，理想濾波器通常是指理想的低通或高通濾波器，可以在頻域上選擇保留或抑制特定頻率範圍內的信號。將這個理想濾波器應用於圖像的傅立葉變換後，再進行逆變換，可以得到經過濾波器處理後的圖像。
具體的細節我們就不介紹了，直接來看代碼：
下面函數是生成理想高通和低通濾波器的圓形掩膜
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下面函數生成一個高斯高、低通濾波器與所需的圓形掩膜
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過濾示例
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')相乘過濾器和移位的圖像得到過濾圖像
為了獲得具有所需頻率響應的最終濾波圖像，關鍵是在頻域中對移位後的圖像與濾波器進行逐點乘法。
這個過程將兩個圖像元素的對應像素相乘。例如，當應用低通濾波器時，我們將對移位的傅立葉變換圖像與低通濾波器逐點相乘。
此操作抑制高頻並保留低頻，對於高通濾波器反之亦然。這個乘法過程對於去除不需要的頻率和增強所需的頻率是必不可少的，從而產生更清晰和更清晰的圖像。
它使我們能夠獲得期望的頻率響應，並在頻域獲得最終濾波圖像。
4、乘法濾波器（Multiplying Filter）和平移後的圖像（Shifted Image）
乘法濾波器是一種以像素值為權重的濾波器，它通過將濾波器的權重與圖像的像素值相乘，來獲得濾波後的像素值。具體地，假設乘法濾波器的權重為h(i,j)，圖像的像素值為f(m,n)，那麼濾波後的像素值g(x,y)可以表示為：
g(x,y) = ∑∑ f(m,n)h(x-m,y-n)
其中，∑∑表示對所有的(m,n)進行求和。
平移後的圖像是指將圖像進行平移操作後的結果。平移操作通常是指將圖像的像素沿著x軸和y軸方向進行平移。平移後的圖像與原始圖像具有相同的大小和解析度，但它們的像素位置發生了變化。在濾波操作中，平移後的圖像可以用於與濾波器進行卷積運算，以實現濾波操作。
此操作抑制高頻並保留低頻，對於高通濾波器反之亦然。這個乘法過程對於去除不需要的頻率和增強所需的頻率是必不可少的，從而產生更清晰和更清晰的圖像。
它使我們能夠獲得期望的頻率響應，並在頻域獲得最終濾波圖像。
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')在可視化傅立葉頻譜時，使用np.log(1+np.abs(x))和20*np.log(np.abs(x))之間的選擇是個人喜好的問題，可以取決於具體的應用程式。
一般情況下會使用Np.log (1+np.abs(x))，因為它通過壓縮數據的動態範圍來幫助更清晰地可視化頻譜。這是通過取數據絕對值的對數來實現的，並加上1以避免取零的對數。
而20*np.log(np.abs(x))將數據按20倍縮放，並對數據的絕對值取對數，這可以更容易地看到不同頻率之間較小的幅度差異。但是它不會像np.log(1+np.abs(x))那樣壓縮數據的動態範圍。
這兩種方法都有各自的優點和缺點，最終取決於具體的應用程式和個人偏好。
5、逆FFT位移
在頻域濾波後，我們需要將圖像移回原始位置，然後應用逆FFT。為了實現這一點，需要使用逆FFT移位，它反轉了前面執行的移位過程。
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、快速傅立葉逆變換(IFFT)
快速傅立葉逆變換(IFFT)是圖像頻率變換的最後一步。它用於將圖像從頻域傳輸回空間域。這一步的結果是在空間域中與原始圖像相比，圖像減少了噪聲並提高了清晰度。
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')總結
我們再把所有的操作串在一起顯示，
函數繪製所有圖像
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使用高通、低通理想濾波器和高斯濾波器的直徑分別為50、100和150像素，展示它們對增強圖像清晰度的影響。
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()可以看到，當改變圓形掩膜的直徑時，對圖像清晰度的影響會有所不同。隨著直徑的增加，更多的頻率被抑制，從而產生更平滑的圖像和更少的細節。減小直徑允許更多的頻率通過，從而產生更清晰的圖像和更多的細節。為了達到理想的效果，選擇合適的直徑是很重要的，因為使用太小的直徑會導致過濾器不夠有效，而使用太大的直徑會導致丟失太多的細節。
一般來說，高斯濾波器由於其平滑性和魯棒性，更常用於圖像處理任務。在某些應用中，需要更尖銳的截止，理想濾波器可能更適合。
利用FFT修改圖像頻率是一種有效的降低噪聲和提高圖像銳度的方法。這包括使用FFT將圖像轉換到頻域，使用適當的技術過濾噪聲，並使用反FFT將修改後的圖像轉換回空間域。通過理解和實現這些技術，我們可以提高各種應用程式的圖像質量。
作者：Muhammad Ahmed