頻率域濾波基本操作

九钟發表於2024-03-30

1.基礎知識

  在頻率域中，影像的低頻分量對應著影像中較為緩慢變化的部分，而高頻分量則對應著影像中較為快速變化的部分。舉例來說，對於一幅房間影像，牆壁和地板通常具有較為平滑的灰度變化，這些平滑的變化可以被認為是低頻分量。而在影像中的物體邊緣、紋理或其他細節部分，則表現出較快速的灰度變化，這些快速變化可以被認為是高頻分量。因此，透過傅立葉變換，我們可以將影像分解為不同頻率的分量，從而更好地理解影像中的灰度變化模式，並針對不同頻率的分量進行處理，以實現各種影像處理任務，如平滑、銳化、邊緣檢測等。

  頻率域中的濾波過程如下：

image

  中心化操作是將影像的低頻分量（影像的平均亮度或顏色漸變）位於頻譜的中心，而高頻分量（影像的邊緣和細節）分佈在邊緣，這使得頻率分佈更加直觀和易於分析。

  濾波過程用公式表達為：

\[g(x,y) = Real\left \{ {\Im^{-1}[H(u,v)F(u,v)]} \right \} \]

  其中\(\Im^{-1}\)是IDFT(傅立葉逆變換), F(u,v)是輸出影像f(x,y)的DFT, H(u,v)是濾波器傳遞函式，g(x,y)是濾波後的影像。函式F,H和g是大小為PxQ的陣列，即經過填充的輸入影像的大小。

  變換中的低頻與影像中緩慢變化的灰度分量有關，如房間的牆壁和室外少雲的天空等。另一方面，高頻由灰度的急劇過渡造成，如邊緣和噪聲。因此，我們可以認為衰減高頻而透過低頻的函式H(u,v)(稱為低通濾波器)將影像模糊；具有相反性質的濾波器（稱為高通濾波器）將使影像的細節更加清晰，但會降低影像的對比度。

2.低通濾波

2.1 理想低通濾波器（ILPF）

  以原點為中心的一個圓無衰減的透過所有頻率，而圓外頻率被截止

\[H(u,v) = \begin{cases}1, D(u,v) \le D_0 \\0,D(u,v) > D_0 \end{cases}\]

式中，\(D_0\)是一個正常數（即圓的半徑），\(D(u,v)\)是頻率中點\(（u,v）\)到\(P\times Q\)頻率域中心的距離，即：

\[D(u,v) = [(u-P/2)^2+(v-Q/2)^2]^{1/2} \]

注意\(P\)與\(Q\)都是填充零後的大小。

示意圖：

image

2.2 高斯低通濾波器（GLPF）

\[H(u,v) = e^{-D(u,v)^2/{2D_0}^2} \]

示意圖：

image

2.3 巴特沃斯低通濾波器（BLPF）

\[H(u,v) = \frac{1}{1+[D(u,v)/D_0]^{2n}} \]

示意圖:

image

3.高通濾波

高通濾波可以由低通濾波簡單變換得到，如下表所示：

理想高通濾波 高斯高通濾波 巴特沃斯高通濾波
\(H(u,v) = \begin{cases}0,D(u,v)\le D_0 \\1,D(u,v)>D_0\end{cases}\) \(H(u,v) = 1 - e^{-D(u,v)^2/{2D_0}^2}\) $H(u,v) = \frac{1}{1+[D_0/D(u,v)]^{2n}} $

示意圖：

image

4.例項

下面我們對一張飛機照片進行頻率域的變換來體會頻率域濾波的特點

image

效果圖：

image

  可以看到，上面三幅圖對圖片進行了低通濾波，截止頻率\(D_0\)設定為20,下面三幅圖對圖片進行了高通濾波，截止頻率設定為2。我們可以發現高斯濾波與巴特沃斯濾波的效果幾乎相同，這是由於這裡我們對巴特沃斯濾波的階數\(n\)設定為1所造成的，較高的n值可以使BLPF函式逼近ILPF的特性，而使用較低的n值可以使BLPF逼近GLPF函式的特性。
  其次還需要注意的是使用理想低通濾波器會產生振鈴效應（由於理想濾波器在頻率域中的截止邊界是突變的，即從完全透過到完全阻止的過渡是突然的），振鈴效應會在濾波後的影像中產生明顯的振鈴或震盪現象，表現為影像中的明暗交替區域出現週期性的光暗環或邊緣，這是由於頻率域中的截止邊界導致的頻率分量的洩漏。這種振鈴效應會使得影像的細節部分產生不自然的震盪，降低影像的質量。而觀察GLPF的示意圖可以其截止邊界較為平滑，因此不會產生振鈴效應。
  ILPF振鈴效應示意圖如下：

image

 使用BLPF可以用截止頻率來控制低頻與高頻的過渡情況（比如醫學影像），其在低階時只有較小的振鈴效應和較小的負值，但遠沒有ILPF那樣的明顯，但是高階振鈴效應非常明顯。

 BLPF振鈴效應示意圖如下：

image

5.總結

 在實際應用中，理想低通濾波器常用於頻域影像平滑，高斯低通濾波器用於頻域影像平滑並緩解振鈴效應，巴特沃斯低通濾波器可調節截止頻率及斜率以更精確地控制頻域影像平滑效果。我們需要根據具體的問題來選擇方法。

6.程式碼

 對三種濾波器進行基本操作的程式碼如下：

import numpy as np
import matplotlib.pyplot as plt
import cv2

def make_transform_fun(image,d,n,T):
    """
    :param image: 原始影像
    :param d:  截止頻率
    :param n:  階數
    :param T:  高通、低通或高斯
    :return:   傳遞函式
    """
    # 獲取影像尺寸
    rows, cols = image.shape

    # 生成中心在原點的 x 和 y 座標網格
    x = np.arange(cols) - cols // 2
    y = np.arange(rows) - rows // 2
    x, y = np.meshgrid(x, y)

    # 計算每個點到中心點的距離 D(u,v)
    # 建立兩個網格 x 和 y 來實現，它們分別包含了每個點相對於中心點的橫縱座標值
    D = np.sqrt(x ** 2 + y ** 2)

    # 根據濾波器型別應用不同的轉換,利用廣播機制
    if T == 'L':  # 低通濾波器
        transform_fun = 1 / (1 + (D / d) ** (2 * n))
    elif T == 'H':  # 高通濾波器
        transform_fun = 1 / (1 + (d / (D + 1e-9)) ** (2 * n))  # 新增小值避免除以零
    elif T == 'GL':  # 高斯低通
        transform_fun = np.exp(-(D ** 2) / (2 * d ** 2))
    elif T == 'GH':  # 高斯高通
        transform_fun = 1 - np.exp(-(D ** 2) / (2 * d ** 2))

    return transform_fun

def IdealLowPassFilter(image):
    """
    理想低通濾波器
    :param img: 頻域圖
    :return: 幅度影像
    """
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    rows,cols = image.shape
    crow,ccol = int(rows/2),int(cols/2)     # 中心位置
    mask = np.zeros((rows,cols),np.uint8)     # 2個通道，分別儲存實部和虛部
    mask[crow - 20:crow + 20,ccol - 20:ccol + 20] = 1

    # 掩膜影像和頻譜影像乘積
    f = fshift * mask

    # 傅立葉逆變換
    ishift = np.fft.ifftshift(f)
    iimg = np.fft.ifft2(ishift)
    iimg = np.abs(iimg)     # 返回幅度值(將逆變換後的複數影像轉換為幅度影像)

    return iimg

def IdealHighPassFilter(image):
    """
    理想高通濾波器
    :param image: 原圖
    :return: 幅度影像
    """
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    rows,cols = image.shape
    crow,ccol = int(rows / 2),int(cols / 2)
    fshift[crow - 2:crow + 2,ccol - 2:ccol + 2] = 0

    # 傅立葉逆變換
    ishift = np.fft.ifftshift(fshift)
    iimg = np.fft.ifft2(ishift)
    iimg = np.abs(iimg)  # 返回幅度值(將逆變換後的複數影像轉換為幅度影像)

    return iimg

def GaussianPassfilter(image,d,T):
    """
    高斯濾波器
    :param image: 原影像
    :param d:  截止頻率
    """
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)
    T = 'G' + T
    d_fun = make_transform_fun(image,d,1,T)
    new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * d_fun)))
    return new_img

def ButterworthPassfilter(image,d,n,T):
    """
    巴特沃爾濾波器
    :param image: 原影像
    :param d: 截止頻率
    :param n: 階數
    :param T: 高通或低通
    :return: 幅度圖
    """
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    d_fun = make_transform_fun(image,d,n,T)
    new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * d_fun)))
    return new_img

if __name__ == "__main__":
    image = cv2.imread("Airport.jpg",cv2.IMREAD_GRAYSCALE)

    # 這裡設定低通半徑為20，高通為2
    result1 = IdealLowPassFilter(image)
    result2 = IdealHighPassFilter(image)
    result3 = GaussianPassfilter(image,20,'L')     # D0方差為20，越小波形越尖
    result4 = GaussianPassfilter(image,2,'H')
    result5 = ButterworthPassfilter(image,20,1,'L')   # n為1類似於高斯,D0方差為20
    result6 = ButterworthPassfilter(image,2,1,'H')

    plt.figure(figsize=(10,5))
    plt.subplot(231), plt.axis('off'), plt.title('ILPF Image D20'), plt.imshow(result1,cmap='gray')
    plt.subplot(234), plt.axis('off'), plt.title('IHPF Image D2'), plt.imshow(result2, cmap='gray')
    plt.subplot(232), plt.axis('off'), plt.title('GLPF Image D20'), plt.imshow(result3, cmap='gray')
    plt.subplot(235), plt.axis('off'), plt.title('GHPF Image D2'), plt.imshow(result4, cmap='gray')
    plt.subplot(233), plt.axis('off'), plt.title('BLPF Image D20'), plt.imshow(result5, cmap='gray')
    plt.subplot(236), plt.axis('off'), plt.title('BLPF Image D2'), plt.imshow(result6, cmap='gray')
    plt.savefig("result.jpg")
    plt.show()

7.參考

1.https://blog.csdn.net/qq_38463737/article/details/118682500
2.數字影像處理第四版（岡薩雷斯）
3.https://zh.wikipedia.org/wiki/%E4%BD%8E%E9%80%9A%E6%BB%A4%E6%B3%A2%E5%99%A8

相關文章