图像的傅立叶变换

关于傅立叶变换的原理请看刚萨雷斯和相关博客——傅里叶分析之掐死教程（完整版）更新于2014.06.06。因为博主不是数学专业大佬，只能从代码上面进行讲解。

1. 前言

图像的傅立叶变换是从空间域变换到空间频率域的一个操作。在频率域，我们可以对图像进行滤波、增强等一系列图像处理步骤。相对于空间域的图像处理来说，频率域的图像处理相对复杂，但用途更加的广泛。然而，人们对于图像傅立叶变换的操作依然感到比较陌生，因此本文对图像傅立叶变换进行一个简单的介绍，但不涉及到图像处理的操作，只是简单的介绍和展示一下傅立叶变换和反变换的过程，来给各位提供一点浅薄的参考。

图像进行傅立叶变换后，得到的是一个频谱，该频谱与图像的大小一致，但包含的是空间频率的信息，主要分为高频分量和低频分量。其中，低频分量，也叫直流分量，代表的是原图像的平滑部分，主要分布在频谱图的四角，占据频谱的大部分能量，因此频谱图中低频分量的亮度很高；高频分量，也叫交流分量，代表的是原图像的边缘部分，主要分布在频谱图的中央，占据频谱的少部分能量，因此频谱图中低频分量的亮度很低。

由于直流分量分布在四角，为了便于处理，我们通常将低频分量移动到频谱的中央，而把高频分量移动到频谱的四角。具体的操作如下图所示：也就是将A和D调换，B和C调换，这个过程就是中心化。

下图给出一个直观的未中心化和中心化后的频谱图的幅度谱：可以看到，中心化后的幅度谱的低频分量位于图像中央。

如果对频谱进行了中心化，那么我们在进行傅立叶反变换来恢复到空间域的时候，需要对频谱再进行一次中心化（前后两次中心化相互抵消），然后再傅立叶反变换才能变回原图；如果没有对原图频谱进行中心化，那我们直接对频谱进行傅立叶反变换即可。显而易见，两者的反变换结果是一样的。

2. 代码

test_fft.py

import numpy as np
import cv2
import torch.fft

# 中心化（反中心化）--> 关于图像中心的对称变换。
def FFT_SHIFT(img):
    fimg = img.clone()
    m,n = fimg.shape
    m = int(m / 2)
    n = int(n / 2)

    for i in range(m):
        for j in range(n):
            fimg[i][j] = img[m+i][n+j]
            fimg[m+i][n+j] = img[i][j]
            fimg[m+i][j] = img[i][j+n]
            fimg[i][j+n] = img[m+i][j]

    return fimg

if __name__ == '__main__':
    # 读图
    pattern = cv2.imread("rectangle.png",0)
    pattern = torch.from_numpy(pattern).type(torch.float32)

    # 傅立叶变换
    pattern_fft = torch.fft.fftn(pattern)
    pattern_fft_shift = FFT_SHIFT(pattern_fft) # 中心化
    pattern_amplitude = torch.log(torch.abs(pattern_fft)) # 获取幅度谱
    pattern_amplitude_shift = torch.log(torch.abs(pattern_fft_shift))  # 获取中心化的幅度谱
    pattern_phase = torch.angle(pattern_fft) # 获取相位谱
    pattern_phase = pattern_phase - torch.floor(pattern_phase / (2 * np.pi)) * 2 * np.pi
    pattern_phase_shift = torch.angle(pattern_fft_shift)  # 获取中心化的相位谱
    pattern_phase_shift = pattern_phase_shift - torch.floor(pattern_phase_shift / (2 * np.pi)) * 2 * np.pi

    pattern_amplitude = pattern_amplitude.numpy() # torch -> numpy
    pattern_amplitude_shift = pattern_amplitude_shift.numpy()
    pattern_phase = pattern_phase.numpy()
    pattern_phase_shift = pattern_phase_shift.numpy()

    # 傅立叶逆变换
    # pattern_ifft = torch.fft.ifftn(pattern_fft) # 直接对未中心化的频谱进行傅立叶反变换
    # 或者对中心化的频谱再中心化（称为反中心化），然后再进行傅立叶反变换
    pattern_ifft_shift = FFT_SHIFT(pattern_fft_shift) # 反中心化
    pattern_ifft_shift = torch.fft.ifftn(pattern_ifft_shift) # 傅立叶反变换
    pattern_ifft_shift = torch.abs(pattern_ifft_shift) # 逆变换的结果是复数，要取幅值才能得到原图

    pattern_ifft_shift = pattern_ifft_shift.numpy() # torch -> numpy

    # 保存
    max_p = np.max(pattern_amplitude)
    min_p = np.min(pattern_amplitude)
    pattern_amplitude = 255.0 * (pattern_amplitude - min_p) / (max_p - min_p)
    pattern_amplitude = pattern_amplitude.astype(np.uint8)
    cv2.imwrite("pattern_amplitude.png",pattern_amplitude) # 保存幅度谱
    pattern_amplitude_shift = 255.0 * (pattern_amplitude_shift - min_p) / (max_p - min_p)
    pattern_amplitude_shift = pattern_amplitude_shift.astype(np.uint8)
    cv2.imwrite("pattern_amplitude_shift.png", pattern_amplitude_shift)  # 保存中心化的幅度谱

    max_p = np.max(pattern_phase)
    min_p = np.min(pattern_phase)
    pattern_phase = 255.0 * (pattern_phase - min_p) / (max_p - min_p)
    pattern_phase = pattern_phase.astype(np.uint8)
    cv2.imwrite("pattern_phase.png",pattern_phase)
    pattern_phase_shift = 255.0 * (pattern_phase_shift - min_p) / (max_p - min_p)
    pattern_phase_shift = pattern_phase_shift.astype(np.uint8)
    cv2.imwrite("pattern_phase_shift.png", pattern_phase_shift) # 保存中心化的相位谱

    pattern_ifft_shift = pattern_ifft_shift.astype(np.uint8)
    cv2.imwrite("pattern_ifft.png",pattern_ifft_shift) # 保存傅立叶反变换结果
    

show_image.py

import matplotlib.pyplot as plt
import cv2

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

# 显示结果
if __name__ == '__main__':
    # 读图
    pattern = cv2.imread("1.png",0)
    pattern_amplitude = cv2.imread("pattern_amplitude.png",0)
    pattern_amplitude_shift = cv2.imread("pattern_amplitude_shift.png",0)
    pattern_phase = cv2.imread("pattern_phase.png", 0)
    pattern_phase_shift = cv2.imread("pattern_phase_shift.png", 0)
    pattern_ifft = cv2.imread("pattern_ifft.png", 0)

    plt.subplot(231),plt.imshow(pattern),plt.title("原图")
    plt.subplot(232),plt.imshow(pattern_amplitude),plt.title("幅度谱")
    plt.subplot(233),plt.imshow(pattern_phase),plt.title("相位谱")
    plt.subplot(234),plt.imshow(pattern_ifft),plt.title("傅立叶反变换")
    plt.subplot(235),plt.imshow(pattern_amplitude_shift),plt.title("中心化的幅度谱")
    plt.subplot(236),plt.imshow(pattern_phase_shift),plt.title("中心化的相位谱")
    plt.show()
    

3. 结果展示

下面展示了三张原图的结果，图中的颜色是 plt 默认的 cmap 映射。

4. 附加——numpy.fft实现傅立叶变换

show_image.py

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("rectangle.png", 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)
    print(phase_f)
    # 进行傅里叶逆变换
    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()