Python圖像處理：頻域?yàn)V波降噪和圖像增強(qiáng)

作者：Muhammad Ahmed 2023-03-09 15:25:49

圖像處理已經(jīng)成為我們?nèi)粘Ｉ钪胁豢苫蛉钡囊徊糠?，涉及到社交媒體和醫(yī)學(xué)成像等各個領(lǐng)域。通過數(shù)碼相機(jī)或衛(wèi)星照片和醫(yī)學(xué)掃描等其他來源獲得的圖像可能需要預(yù)處理以消除或增強(qiáng)噪聲。頻域?yàn)V波是一種可行的解決方案，它可以在增強(qiáng)圖像銳化的同時(shí)消除噪聲。

快速傅里葉變換(FFT)是一種將圖像從空間域變換到頻率域的數(shù)學(xué)技術(shù)，是圖像處理中進(jìn)行頻率變換的關(guān)鍵工具。通過利用圖像的頻域表示，我們可以根據(jù)圖像的頻率內(nèi)容有效地分析圖像，從而簡化濾波程序的應(yīng)用以消除噪聲。本文將討論圖像從FFT到逆FFT的頻率變換所涉及的各個階段，并結(jié)合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)是一種廣泛應(yīng)用的數(shù)學(xué)技術(shù)，它允許圖像從空間域轉(zhuǎn)換到頻率域，是頻率變換的基本組成部分。利用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()計(jì)算傅里葉變換f的幅度，np.log()轉(zhuǎn)換為對數(shù)刻度，然后乘以20得到以分貝為單位的幅度。這樣做是為了使幅度譜更容易可視化和解釋。

2、FFT位移

為了使濾波算法應(yīng)用于圖像，利用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、過濾

頻率變換的的一個目的是使用各種濾波算法來降低噪聲和提高圖像質(zhì)量。兩種最常用的圖像銳化濾波器是Ideal high-pass filter 和Gaussian high-pass filter。這些濾波器都是使用的通過快速傅里葉變換(FFT)方法獲得的圖像的頻域表示。

Ideal high-pass filter（理想濾波器）是一種無限長的、具有無限頻帶寬和理想通帶和阻帶響應(yīng)的濾波器。其通帶內(nèi)所有頻率的信號都被完全傳遞，而阻帶內(nèi)所有頻率的信號則完全被抑制。

在頻域上，理想濾波器的幅頻響應(yīng)為：

在通帶內(nèi)，幅頻響應(yīng)為 1
在阻帶內(nèi)，幅頻響應(yīng)為 0

在時(shí)域上，理想濾波器的沖激響應(yīng)為：

在通帶內(nèi)，沖激響應(yīng)為一個無限長的單位沖激函數(shù)序列
在阻帶內(nèi)，沖激響應(yīng)為零

由于理想濾波器在頻域上具有無限帶寬，因此它無法在實(shí)際應(yīng)用中實(shí)現(xiàn)。實(shí)際中使用的數(shù)字濾波器通常是基于理想濾波器的逼近，所以才被成為只是一個Ideal。

高斯高通濾波器（Gaussian high-pass filter）是一種在數(shù)字圖像處理中常用的濾波器。它的作用是在圖像中保留高頻細(xì)節(jié)信息，并抑制低頻信號。該濾波器基于高斯函數(shù)，具有光滑的頻率響應(yīng)，可以適應(yīng)各種圖像細(xì)節(jié)。

高斯高通濾波器的頻率響應(yīng)可以表示為：

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

其中，L(u,v)是一個低通濾波器，它可以用高斯函數(shù)表示。通過將低通濾波器的響應(yīng)從1中減去，可以得到一個高通濾波器的響應(yīng)。在實(shí)際中，通常使用不同的參數(shù)設(shè)置來調(diào)整高斯函數(shù)，以達(dá)到不同的濾波效果。

圓形掩膜（disk-shaped images）是用于定義在圖像中進(jìn)行傅里葉變換時(shí)要保留或抑制的頻率分量。在這種情況下，理想濾波器通常是指理想的低通或高通濾波器，可以在頻域上選擇保留或抑制特定頻率范圍內(nèi)的信號。將這個理想濾波器應(yīng)用于圖像的傅里葉變換后，再進(jìn)行逆變換，可以得到經(jīng)過濾波器處理后的圖像。

具體的細(xì)節(jié)我們就不介紹了，直接來看代碼：

下面函數(shù)是生成理想高通和低通濾波器的圓形掩膜

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

下面函數(shù)生成一個高斯高、低通濾波器與所需的圓形掩膜

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')

相乘過濾器和移位的圖像得到過濾圖像

為了獲得具有所需頻率響應(yīng)的最終濾波圖像，關(guān)鍵是在頻域中對移位后的圖像與濾波器進(jìn)行逐點(diǎn)乘法。

這個過程將兩個圖像元素的對應(yīng)像素相乘。例如，當(dāng)應(yīng)用低通濾波器時(shí)，我們將對移位的傅里葉變換圖像與低通濾波器逐點(diǎn)相乘。

此操作抑制高頻并保留低頻，對于高通濾波器反之亦然。這個乘法過程對于去除不需要的頻率和增強(qiáng)所需的頻率是必不可少的，從而產(chǎn)生更清晰和更清晰的圖像。

它使我們能夠獲得期望的頻率響應(yīng)，并在頻域獲得最終濾波圖像。

4、乘法濾波器（Multiplying Filter）和平移后的圖像（Shifted Image）

乘法濾波器是一種以像素值為權(quán)重的濾波器，它通過將濾波器的權(quán)重與圖像的像素值相乘，來獲得濾波后的像素值。具體地，假設(shè)乘法濾波器的權(quán)重為h(i,j)，圖像的像素值為f(m,n)，那么濾波后的像素值g(x,y)可以表示為：

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

其中，∑∑表示對所有的(m,n)進(jìn)行求和。

平移后的圖像是指將圖像進(jìn)行平移操作后的結(jié)果。平移操作通常是指將圖像的像素沿著x軸和y軸方向進(jìn)行平移。平移后的圖像與原始圖像具有相同的大小和分辨率，但它們的像素位置發(fā)生了變化。在濾波操作中，平移后的圖像可以用于與濾波器進(jìn)行卷積運(yùn)算，以實(shí)現(xiàn)濾波操作。

它使我們能夠獲得期望的頻率響應(yīng)，并在頻域獲得最終濾波圖像。

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')

在可視化傅里葉頻譜時(shí)，使用np.log(1+np.abs(x))和20*np.log(np.abs(x))之間的選擇是個人喜好的問題，可以取決于具體的應(yīng)用程序。

一般情況下會使用Np.log (1+np.abs(x))，因?yàn)樗ㄟ^壓縮數(shù)據(jù)的動態(tài)范圍來幫助更清晰地可視化頻譜。這是通過取數(shù)據(jù)絕對值的對數(shù)來實(shí)現(xiàn)的，并加上1以避免取零的對數(shù)。

而20*np.log(np.abs(x))將數(shù)據(jù)按20倍縮放，并對數(shù)據(jù)的絕對值取對數(shù)，這可以更容易地看到不同頻率之間較小的幅度差異。但是它不會像np.log(1+np.abs(x))那樣壓縮數(shù)據(jù)的動態(tài)范圍。

這兩種方法都有各自的優(yōu)點(diǎn)和缺點(diǎn)，最終取決于具體的應(yīng)用程序和個人偏好。

5、逆FFT位移

在頻域?yàn)V波后，我們需要將圖像移回原始位置，然后應(yīng)用逆FFT。為了實(shí)現(xiàn)這一點(diǎn)，需要使用逆FFT移位，它反轉(zhuǎn)了前面執(zhí)行的移位過程。

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)是圖像頻率變換的最后一步。它用于將圖像從頻域傳輸回空間域。這一步的結(jié)果是在空間域中與原始圖像相比，圖像減少了噪聲并提高了清晰度。

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')

總結(jié)

我們再把所有的操作串在一起顯示，

函數(shù)繪制所有圖像

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像素，展示它們對增強(qiáng)圖像清晰度的影響。

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()

可以看到，當(dāng)改變圓形掩膜的直徑時(shí)，對圖像清晰度的影響會有所不同。隨著直徑的增加，更多的頻率被抑制，從而產(chǎn)生更平滑的圖像和更少的細(xì)節(jié)。減小直徑允許更多的頻率通過，從而產(chǎn)生更清晰的圖像和更多的細(xì)節(jié)。為了達(dá)到理想的效果，選擇合適的直徑是很重要的，因?yàn)槭褂锰〉闹睆綍?dǎo)致過濾器不夠有效，而使用太大的直徑會導(dǎo)致丟失太多的細(xì)節(jié)。

一般來說，高斯濾波器由于其平滑性和魯棒性，更常用于圖像處理任務(wù)。在某些應(yīng)用中，需要更尖銳的截止，理想濾波器可能更適合。

利用FFT修改圖像頻率是一種有效的降低噪聲和提高圖像銳度的方法。這包括使用FFT將圖像轉(zhuǎn)換到頻域，使用適當(dāng)?shù)募夹g(shù)過濾噪聲，并使用反FFT將修改后的圖像轉(zhuǎn)換回空間域。通過理解和實(shí)現(xiàn)這些技術(shù)，我們可以提高各種應(yīng)用程序的圖像質(zhì)量。

責(zé)任編輯：華軒來源： DeepHub IMBA