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