(九)OpenCV-Python學習—影象傅立葉變換

阿新 • • 發佈：2020-10-10

　　對於二維圖片，可以對其進行傅立葉變換，獲取圖片的頻譜資訊。頻譜有很多應用，包括顯著性檢測，卷積定理，頻率域濾波等，下面是圖片傅立葉變換的一些基本概念：

1. 影象傅立葉變換

　　對於M行N列的影象矩陣f(x,y)，f(x, y)表示第x行y列的畫素值，則存在複數矩陣F，有以下公式：

　　F(u,v)稱為f(x, y)的傅立葉變換，f(x,y)稱為F(u,v)的傅立葉逆變換

　　opencv提供函式dft()可以對影象進行傅立葉變換和傅立葉逆變換，函式引數如下：

dst =cv.dft(src, flags=0, nonzeroRows=0)
    src: 輸入影象矩陣，只支援CV_32F或者CV_64F的單通道或雙通道矩陣(單通道的代表實數矩陣，雙通道代表複數矩陣)
    flags: 轉換的識別符號，取值包括DFT_COMPLEX_OUTPUT,DFT_REAL_OUTPUT,DFT_INVERSE,DFT_SCALE, DFT_ROWS，通常組合使用
        DFT_COMPLEX_OUTPUT: 輸出複數形式
        DFT_REAL_OUTPUT: 只輸出複數的實部
        DFT_INVERSE:進行傅立葉逆變換
        DFT_SCALE:是否除以M 
*N （M行N列的圖片，共有有M*N個畫素點）
        DFT_ROWS:輸入矩陣的每一行進行傅立葉變換或者逆變換
        (傅立葉正變換一般取flags=DFT_COMPLEX_OUTPUT，
         傅立葉逆變換一般取flags= DFT_REAL_OUTPUT+DFT_INVERSE+DFT_SCALE)     
    nonzerosRows: 當設定為非0時，對於傅立葉正變換，表示輸入矩陣只有前nonzerosRows行包含非零元素；對於傅立葉逆變換，表示輸出矩陣只有前nonzerosRows行包含非零元素

返回值：
    dst: 單通道的實數矩陣或雙通道的複數矩陣

2. 快速傅立葉變化

　　對於M行N列的影象矩陣，傅立葉變換理論上需要 $次運算，非常耗時。而當和時，傅立葉變換可以通過O(MNlog(M*N)) 次運算就能完成運算，即快速傅立葉變換。當圖片矩陣的行數和列數都可以分解成時，opencv中的dft()會進行傅立葉變換快速演算法，所以在計算時一般先對二維矩陣進行擴充補0，已滿足規則，opencv提供函式getOptimalDFTSize()函式來計算需要補多少行多少列的0，其引數如下：$

retval=cv.getOptimalDFTSize(vecsize)
    vecsize: 整數，圖片矩陣的行數或者列數，函式返回一個大於或等於vecsize的最小整數N，且N滿足N 
=2^p*3^q*5^r

　　快速傅立葉變換dft()使用示例程式碼如下：

#coding:utf-8

import cv2
import numpy as np

img_gray = cv2.imread(r"D:\data\receipt_rotate.jpg", 0)
rows, cols = img_gray.shape[:2]

#快速傅立葉變換，輸出複數
rPadded = cv2.getOptimalDFTSize(rows)
cPadded = cv2.getOptimalDFTSize(cols)
imgPadded = np.zeros((rows, cols), np.float32)  # 填充
imgPadded[:rows, :cols] = img_gray
fft_img = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)
print(fft_img)

#快速傅立葉逆變換，只輸出實數部分
ifft_img = cv2.dft(fft_img, flags=cv2.DFT_REAL_OUTPUT+cv2.DFT_INVERSE+cv2.DFT_SCALE)
ori_img = np.copy(ifft_img[:rows, :cols])  # 裁剪
print(img_gray)
print(ori_img)
print(np.max(ori_img-img_gray))   #9.1552734e-05，接近於0

new_gray_img = ori_img.astype(np.uint8)

cv2.imshow("img_gray", img_gray)
cv2.imshow("new_gray_img", new_gray_img)
cv2.waitKey(0)
cv2.destroyAllWindows()

快速傅立葉變換demo

3. 傅立葉幅度譜和相位譜

　　影象矩陣進行快速傅立葉變換後得到複數矩陣F，記Real為矩陣F的實部，Imaginary為矩陣的虛部，則F可以表示為：

　　幅度譜(Amplitude Spectrum)，又稱相位譜，則可以表示為：

　　相位譜(phase spectrum)可以表示為：

　　因此，F也可以表示為：

4. 幅度譜和相位譜的灰度化

　　進行快速傅立葉變換後，一般會計算圖片矩陣的幅度譜和相位譜，為了觀察，一般會將幅度譜和相位譜進行灰度化，即將幅度譜和相位譜轉變為灰度圖，從而進行圖片視覺化。

　　因為幅度譜的最大值在左上角(0, 0)處，通常將其移動到幅度譜的中心位置(w/2, h/2)，因此需要在進行影象傅立葉變換前，將影象矩陣乘以 $, 即傅立葉譜的中心化(-1)^(r+c) 。綜上所述，標準的傅立葉變換處理步驟如下：$

　　計算幅度譜和相位譜的程式碼和結果如下：

#coding:utf-8

import cv2
import numpy as np

def fftImage(gray_img, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), np.float32)
    imgPadded[:rows, :cols] = gray_img
    fft_img = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)  #輸出為複數，雙通道
    return fft_img

#計算幅度譜
def amplitudeSpectrum(fft_img):
    real = np.power(fft_img[:, :, 0], 2.0)
    imaginary = np.power(fft_img[:, :, 1], 2.0)
    amplitude = np.sqrt(real+imaginary)
    return amplitude

#幅度譜的灰度化
def graySpectrum(amplitude):
    amplitude = np.log(amplitude+1)  #增加對比度
    spectrum = cv2.normalize(amplitude,  0, 1, cv2.NORM_MINMAX, dtype=cv2.CV_32F)
    spectrum *= 255

    #歸一化，幅度譜的灰度化
    # spectrum = amplitude*255/(np.max(amplitude))
    # spectrum = spectrum.astype(np.uint8)

    return spectrum

#計算相位譜並灰度化
def phaseSpectrum(fft_img):
    phase = np.arctan2(fft_img[:,:,1], fft_img[:, :, 0])
    spectrum = phase*180/np.pi  #轉換為角度，在[-180,180]之間
    # spectrum = spectrum.astype(np.uint8)
    return spectrum

def stdFftImage(img_gray, rows, cols):
    fimg = np.copy(img_gray)
    fimg = fimg.astype(np.float32)
    # 1.影象矩陣乘以（-1）^(r+c), 中心化
    for r in range(rows):
        for c in range(cols):
            if(r+c)%2:
                fimg[r][c] = -1*img_gray[r][c]
    fft_img = fftImage(fimg, rows, cols)
    amplitude = amplitudeSpectrum(fft_img)
    ampSpectrum = graySpectrum(amplitude)
    return ampSpectrum


if __name__ == "__main__":
    img_gray = cv2.imread(r"D:\data\receipt_rotate.jpg", 0)
    # img_gray = cv2.imread(r"D:\data\carbrand.png", 0)
    rows, cols = img_gray.shape[:2]
    fft_img = fftImage(img_gray, rows, cols)
    amplitude = amplitudeSpectrum(fft_img)
    ampSpectrum = graySpectrum(amplitude)   #幅度譜灰度化
    phaSpectrum = phaseSpectrum(fft_img)    #相位譜灰度化

    ampSpectrum2 = stdFftImage(img_gray, rows, cols) # 幅度譜灰度化並中心化
    cv2.imshow("img_gray", img_gray)
    cv2.imshow("ampSpectrum", ampSpectrum)
    cv2.imshow("ampSpectrum2", ampSpectrum2)
    cv2.imshow("phaSpectrum", phaSpectrum)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

幅度譜灰度級

5. 卷積定理

　　卷積定理指：空間域上的卷積等於頻率域上的乘積，空間域的乘積等於頻率域的卷積。即空間域中兩個函式的卷積的傅立葉變換等於兩個函式的傅立葉變換在頻率域中的乘積，反之兩個函式的傅立葉變換的卷積等於兩個函式的乘積，表示式如下：

　　對應影象上，若I是M行N列的影象矩陣，k是m行n列的卷積核，I與k的全卷積I* $k具有M+m-1行N+n-1列，這是因為卷積計算過程中採用0擴充邊界，假設I_padded 是對 I 填充邊界後的矩陣，k_padded是對k填充邊界後卷積核，FTIp和FTkp分別是I_padded和k_padded的傅立葉變換，則：$

　　卷積定理的一個運用就是能通過快速傅立葉變換來計算影象處理中的卷積，其大致流程如下：

將圖片I和卷積核k，進行邊界填充為I_padded和k_padded
計算I_padded和k_padded的快速傅立葉變換，得到FTIp和FTkp
將FTIp和FTkp的乘積進行傅立葉逆變換，並只取實部，得到影象矩陣
將上一步得到的影象矩陣進行裁剪，即為最終的卷積影象矩陣

6. 普殘差顯著性檢測

　　視覺注意性機制是一種選擇性注意，總是選擇圖片中最顯著的，與其他內容差異較大的成分(大多為圖片前景的某部分)。視覺顯著性檢測有很多方法，基於幅度譜殘差，是一種簡單有效的方法，其演算法步驟如下：

1.計算影象的快速傅立葉變換矩陣F，並計算幅度譜的灰度級graySpectrum
2.計算相位譜phaseSpectrum，然後根據相位譜計算對應的正弦譜和餘弦譜
3.對前面計算的幅度譜灰度級graySpectrum進行均值平滑，記為
4.計算普殘差(spectralResidual)，普殘差的計算公式如下：

5.對普殘差進行冪指數運算，即exp(spectrualResidual)
6.將第5步得到的冪指數作為信的幅度譜，仍然使用原來的相位譜，通過傅立葉逆變換，得到複數矩陣F2
7.對複數矩陣F2，計算其實部和虛部的平方，再進行高斯平滑，灰度級轉換，即得到顯著性

　　程式碼使用如下：

#coding:utf-8
import cv2
import numpy as np

def fftImage(gray_img, rows, cols):
    rPadded = cv2.getOptimalDFTSize(rows)
    cPadded = cv2.getOptimalDFTSize(cols)
    imgPadded = np.zeros((rPadded, cPadded), np.float32)
    imgPadded[:rows, :cols] = gray_img
    fft_img = cv2.dft(imgPadded, flags=cv2.DFT_COMPLEX_OUTPUT)  #輸出為複數，雙通道
    return fft_img


if __name__ == "__main__":
    img_gray = cv2.imread(r"C:\Users\silence_cho\Desktop\Messi.jpg", 0)
    rows, cols = img_gray.shape[:2]
    fft_image = fftImage(img_gray, rows, cols)

    # 計算幅度譜及其log值
    amplitude = np.sqrt(np.power(fft_image[:, :, 0], 2.0) + np.power(fft_image[:, :, 1], 2.0)) #計算幅度譜
    ampSpectrum = np.log(amplitude+1)   # 幅度譜的log值

    #計算相位譜, 餘弦譜，正弦譜
    phaSpectrum = np.arctan2(fft_image[:,:,1], fft_image[:, :, 0])  # 計算相位譜，結果為弧度
    cosSpectrum = np.cos(phaSpectrum)
    sinSpectrum = np.sin(phaSpectrum)

    #幅度譜灰度級均值平滑
    meanAmpSpectrum= cv2.boxFilter(ampSpectrum, cv2.CV_32FC1, (3, 3))

    #殘差
    spectralResidual = ampSpectrum - meanAmpSpectrum
    expSR = np.exp(spectralResidual)

    #實部，虛部，複數矩陣
    real = expSR*cosSpectrum
    imaginary = expSR*sinSpectrum
    new_matrix = np.zeros((real.shape[0], real.shape[1], 2), np.float32)
    new_matrix[:, :, 0] = real
    new_matrix[:, :, 1] = imaginary
    ifft_matrix = cv2.dft(new_matrix, flags=cv2.DFT_COMPLEX_OUTPUT+cv2.DFT_INVERSE)

    #顯著性
    # saliencymap = np.sqrt(np.power(ifft_matrix[:, :, 0], 2) + np.power(ifft_matrix[:, :, 1], 2))
    saliencymap = np.power(ifft_matrix[:, :, 0], 2) + np.power(ifft_matrix[:, :, 1], 2)
    saliencymap = cv2.GaussianBlur(saliencymap, (11, 11), 2.5)
    saliencymap = saliencymap/np.max(saliencymap)
    #伽馬變換，增加對比度
    saliencymap = np.power(saliencymap, 0.5)
    saliencymap = np.round(saliencymap*255)
    saliencymap = saliencymap.astype(np.uint8)

    cv2.imshow("img_gray", img_gray)
    cv2.imshow("saliencymap", saliencymap)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

普殘差顯著性檢測

參考： https://zhuanlan.zhihu.com/p/115002897?utm_source=ZHShareTargetIDMore

　　 https://blog.csdn.net/kena_m/article/details/49406687

(九)OpenCV-Python學習—影象傅立葉變換

(九)OpenCV-Python學習—影象傅立葉變換

python用opencv 影象傅立葉變換

python 影象的離散傅立葉變換例項

Python 實現影象快速傅立葉變換和離散餘弦變換

OpenCV學習(5)--離散傅立葉變換、濾波、侵蝕、擴張

[Python影象處理] 三十二.傅立葉變換（影象去噪）與霍夫變換（特徵識別）萬字詳細總結

python影象的頻域(傅立葉變換、低通濾波)

FFT快速傅立葉變換的python實現過程解析

使用python實現離散時間傅立葉變換的方法

學習總結-快速傅立葉變換（FastFourierTransform）

數字影象處理（二）傅立葉變換

opencv-dft傅立葉變換(把空域變成頻域)

OpenCV 離散傅立葉變換

【筆記】傅立葉變換學習筆記

【學習筆記】快速離散傅立葉變換（FFT）（遞迴版）

【影象處理筆記】傅立葉變換【影象處理筆記】總目錄圖傅立葉變換

快速傅立葉變換FFT

快速傅立葉變換（FFT）

【數論】快速傅立葉變換 - 多項式相乘

多項式乘法與離散傅立葉變換

(九)OpenCV-Python學習—影象傅立葉變換

相關推薦