图像的傅里叶变换（opencv）

对于非数学专业的人来说，理解图像的傅里叶变化真的是一件很困难的事情，刚开始满头雾水的我疯狂的在网上找文章，功夫不负有心人，我找到了这个：

https://www.cnblogs.com/h2zZhou/p/8405717.html

真的是非数学专业的救命草啊！！！！

理解了傅里叶变换的原理之后就需要自己动手实现了

图像经过基础的傅里叶变换得到的矩阵元素是复数类型，从复数矩阵中得到每个元素的模组成的矩阵即为该图像的幅度谱，再从复数矩阵中得到每个元素的相位组成的矩阵即为该图像的相位谱，对复数矩阵做傅里叶逆变换即可得到原来的图像。

刚开始我用cv2.dft实现，但是它的返回的是一个两层的三维矩阵，第一层是复数的实数部分，第二层是虚数部分，比较麻烦。只实现了它的幅度谱和相位谱，代码如下：

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


def fft2_image(image):
    """
    傅里叶变换理论上需要0((MN)^2)次运算，这是非常耗时的，并极大地降低了傅里叶在图像处理中的应用。
    幸运的是，当M = 2^m和N = 2^n时，或者对于任意的M和N，傅里叶变换通过O(MN log(MN))次运算就可以完成，
    称为“快速傅里叶变换”，快速傅里叶变换时针对于行数和列数均满足可以分解为2^p×3^q×5^r的情况的
    利用函数getOptimalDFTSize实现快速傅里叶变换：
    该函数可在矩阵的右侧和左侧补0，使得该矩阵的行数和列数均满足可以分解为为2^p×3^q×5^r
    :param image:输入图像
    :return:傅里叶变换的结果
    """
    # 得到行和列
    r, c = image.shape[:2]
    # 得到快速傅里叶变换的最优扩充
    rPadded = cv2.getOptimalDFTSize(r)
    cPadded = cv2.getOptimalDFTSize(c)
    # 边缘扩充，下边缘和右边缘的扩充值为0
    fft2 = np.zeros((rPadded, cPadded, 2))
    fft2[:r, :c, 0] = image
    # 快速傅里叶变换
    cv2.dft(fft2, fft2, cv2.DFT_COMPLEX_OUTPUT)
    return fft2


def amplitude_spectrum(fft2):
    """
    幅度谱函数，即就是求每个(u, v)的模
    :param fft2:傅里叶变换的复数形式的输出
    :return:幅度的灰度值数组，但是其中很大一部分值大于255，所以还需要一个幅度谱显示函数
    """
    real2 = np.power(fft2[:, :, 0], 2.0)
    imag2 = np.power (fft2[:, :, 1], 2.0)
    amplitude = np.sqrt(real2 + imag2)
    return amplitude


def gray_spectrum(amplitude):
    """
    幅度谱显示函数，一般采用对数函数对幅度谱进行数值压缩，
    再进行归一化，这样的到的幅度谱的灰度级显示的对比度会比较高。
    :param amplitude: 取傅里叶变换幅度之后的矩阵
    :return: 归一化后的矩阵
    """
    # 对比度拉伸
    amplitude = np.log(amplitude + 1.0)
    # 归一化，傅里叶谱的灰度级数显示（归一化的对象灰度值必须是numpy.float类型，否则图像像素值全为0）
    spectrum = np.zeros(amplitude.shape)
    cv2.normalize(amplitude, spectrum, 0, 1, cv2.NORM_MINMAX)
    return spectrum


def phase_spectrum(fft2):
    """
    相位谱函数，利用Numpy中的函数arctan2，
    该函数的第一个参数是输入的虚部矩阵，第二个参数是输入参数的实部矩阵，
    返回值为对应位置的相位角，数值范围为[-π,π]，
    可以将返回值的矩阵除以π再乘以180，规格化到[-180,180]
    :param fft2:输入矩阵，fft2[:, :, 1],为虚部矩阵，fft2[:, :, 1]为实部矩阵
    :return:
    """
    # 得到行数列数
    rows, cols = fft.shape[:2]
    # 计算相位角
    phase = np.arctan2(fft2[:, :, 0], fft2[:, :, 1])
    # 将相位角转换为[-180, 180]
    spectrum = phase/math.pi*180
    return spectrum


def show_images(images):
    """
    显示函数
    :param images: 原始图像列表
    :return: none
    """
    for i in range(len(images)):
        img = images[i]
        #行，列，索引
        x = 2
        plt.subplot(x, math.ceil(len(images)/x), i+1)
        plt.imshow(img, cmap="gray")
        title = "("+str(i+1)+")"
        plt.title(title,fontsize=10)
        plt.xticks([])
        plt.yticks([])
    plt.show()


if __name__ == "__main__":
    image = cv2.imread("..\\images\\rectangle.jpg", 0)
    # 计算图像矩阵的傅里叶变换（这里由于尺寸问题没用快速傅里叶变换）
    fft = fft2_image(image)
    # fft = np.zeros((image.shape[0], image.shape[1], 2))
    # fft[:, :, 0] = image
    # cv2.dft(fft, fft, cv2.DFT_COMPLEX_OUTPUT)
    # 求幅度谱
    amplitude = amplitude_spectrum(fft)
    # 幅度谱的灰度级数显示
    amp_spectrum = gray_spectrum(amplitude)
    cv2.imshow("0", amp_spectrum)

    # 相位谱的灰度级数显示
    phase_spe = phase_spectrum(fft)
    cv2.imshow("1", phase_spe)
    # 以下是傅里叶变换的幅度谱的中心化
    # 第一步：图像矩阵乘以(-1)^(r+c)
    rows, cols = image.shape
    fimg = np.copy(image)
    for r in range(rows):
        for c in range(cols):
            if (r+c)%2:
                fimg[r][c] = -1*image[r][c]
    # 第二步，快速傅里叶变换
    imgfft2 = fft2_image(fimg)
    # 第三步，求傅里叶变换的幅度谱
    am_spe = amplitude_spectrum(imgfft2)
    # 幅度谱的灰度级显示
    gray_spe1 = gray_spectrum(am_spe)
    cv2.imshow("2", gray_spe1)
    gray_spe2 = gray_spe1*255
    gray_spe2 = gray_spe2.astype(np.uint8)
    cv2.imshow("3", gray_spe2)
    images = [image, amp_spectrum, phase_spe, gray_spe1, gray_spe2]
    show_images(images)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
    """
    images = [phase_spe]
    show_images(images)
    
    cv2.imshow("1", phase_spe)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
    
    file = open(r'..\\array.txt', 'w')
    np.savetxt(file, amp_spectrum)
    file.close()
    
    # 傅里叶逆变换
    ifft2 = np.zeros(fft2.shape[:2])
    cv2.dft(fft2, ifft2, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)
    # 裁剪
    img = np.copy(ifft2[:image.shape[0], :image.shape[1]])
    # 裁剪后的结果image等于img
    # 通过判断原矩阵减去逆变换裁剪后的矩阵是否为零矩阵，验证两者是否相同
    # print("{:.20f}".format(np.max(img - image)))
    """

但是，用numpy包中带的傅里叶变换的函数就非常简单了。

import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

if __name__ == "__main__":
    image = cv2.imread("..\\images\\rectangle.jpg", 0)
    #对image进行快速傅里叶变换，输出为复数形式
    fft = np.fft.fft2(image)
    # np.fft.fftshift将零频率分量移动到频谱中心
    fft_shift = np.fft.fftshift(fft)
    # 取绝对值，将复数形式变为实数
    # 取绝对值目的为了将数据变化到较小的范围
    s1 = np.log(np.abs(fft))
    s2 = np.log(np.abs(fft_shift))
    # np.angle能够直接根据复数的虚部和实部求出角度，默认的角度是弧度
    phase_f = np.angle(fft)
    phase_f_shift = np.angle(fft_shift)
    # 进行傅里叶逆变换
    ifft_shift = np.fft.ifftshift(fft_shift)
    ifft = np.fft.ifft2(ifft_shift)
    #出来的是复数，无法显示
    image_back = np.abs(ifft)
    # 显示
    plt.subplot(231), plt.imshow(image, 'gray'), plt.title("原图"), plt.xticks([]), plt.yticks([])
    plt.subplot(232), plt.imshow(s1, 'gray'), plt.title("幅度谱中心化以前"), plt.xticks([]), plt.yticks([])
    plt.subplot(233), plt.imshow(s2, 'gray'), plt.title("幅度谱中心化以后"), plt.xticks([]), plt.yticks([])
    plt.subplot(234), plt.imshow(phase_f, 'gray'), plt.title("相位谱中心化以前"), plt.xticks([]), plt.yticks([])
    plt.subplot(235), plt.imshow(phase_f_shift, 'gray'), plt.title("相位谱中心化以后"), plt.xticks([]), plt.yticks([])
    plt.subplot(236), plt.imshow(image_back, 'gray'), plt.title("傅里叶逆变换结果"), plt.xticks([]), plt.yticks([])
    plt.show()

运行结果：

问题：为什么中心化以前幅度谱的四周较亮？

因为傅里叶频谱是对称的，以一维傅里叶为例（用numpy函数绘制了一个非周期信号对其进行傅里叶变换）：

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
plt.rcParams['axes.unicode_minus']=False       #显示负号

if __name__ == "__main__":
    x = np.linspace(0, 10)
    y = np.zeros(len(x))
    for i in range(len(x)):
        if x[i] < 1:
            y[i] = 4
        elif 4 < x[i] < 8:
            y[i] = 3
        else:
            y[i] = 0
    plt.subplot(211), plt.plot(y)
    transformed = np.fft.fft(y)
    plt.subplot(212), plt.plot(transformed)
    plt.show()

结果：

我也不知道为什么变换结果上50附近没有绘制出来，若是有人看到能不能指教一下？

但是，忽略这一点点小瑕疵，非周期函数的傅里叶变换频谱是关于频率中心对称的。

根据二维傅里叶变换的可分离性（一个二维傅里叶变换可分解为两个一维傅里叶变换），可得二维傅里叶变换是关于中点对称的，四个角落为低频分量。又因为图像大多由低频分量组成，所以会呈现出四个角落比较亮的情况。

在别人的博客上看到的一张图：