(九)OpenCV-Python學習—圖像傅里葉變換

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

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列的圖像矩陣，傅里葉變換理論上需要 (m*n)²次運算，非常耗時。而當 M=2^m 和N=2ⁿ 時，傅里葉變換可以通過O(MNlog(M*N)) 次運算就能完成運算，即快速傅里葉變換。當圖片矩陣的行數和列數都可以分解成 2^p*3^q*5^r時，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()

 

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