图像傅里叶变换基础在本节中，我们将介绍离散傅里叶变换( Discrete Fourier Transform, DFT

开启掘金成长之旅！这是我参与「掘金日新计划 · 12 月更文挑战」的第18天，点击查看活动详情

前言

在本节中，我们将介绍离散傅里叶变换( Discrete Fourier Transform, DFT )相关基础概念，并使用 Numpy 的 fft 模块在图像中应用 DFT，或者我们也可以使用 Scipy 的 fftpack 模块。

图像傅里叶变换

我们首先介绍 2D 傅里叶变换及 2D 逆离散傅里叶变换。(灰度)图像可以被定义为 2D 函数 $f(x,y)$ ，其中 $(x,y)∈\{0,...,M-1\}×\{0,...,N-1\}$ 。DFT 可以将图像从其空间表示 $f(x,y)$ 改变为频域表示 $f(u,v)$ ，其中 $(u,v)∈\{0,...,M-1\}×\{0,...,N-1\}$ 表示频率分量/傅里叶基向量：

绘图7.png

(1) 首先，导入所有必需的库：

import numpy as np
import numpy.fft as fp
from skimage.io import imread
from skimage.color import rgb2gray
from skimage.metrics import peak_signal_noise_ratio
import matplotlib.pyplot as plt

(2) 定义函数 plot_image()，使用 matplotlib.pyplot 的 imshow() 函数绘制图像：

def plot_image(im, title):
    plt.imshow(im, cmap='gray')
    plt.axis('off')
    plt.title(title, size=10)

(3) 接下来，定义函数 plot_freq_spectrum() 以绘制图像的频(功率)谱，该函数可以通过可选参数控制是否显示颜色映射图和轴刻度。我们需要使用 numpy.fft.fftshift() 函数将零功率谱系数 (0,0) 移动到功率谱中心，然后可视化功率谱：

numpy.fft.fftshift(x, axes=None)

傅里叶变换会返回一复数数组；复数部分对应于相位，实数部分对应于我们感兴趣的幅度：

def plot_freq_spectrum(F, title, cmap=plt.cm.gray, show_axis=True, colorbar=False):
    plt.imshow((20*np.log10(0.1 + fp.fftshift(F))).real.astype(int), cmap=cmap)
    if not show_axis:
        plt.axis('off')
    if colorbar:
        plt.colorbar()
    plt.title(title, size=10)

(4) 接下来，创建一些简单图像并获取 DFT 的功率谱。首先分别生成一对大小为 100x100 且具有(等间距)水平和垂直条纹的周期图像：

h, w = 100, 100
images = list()
im = np.zeros((h,w))
for x in range(h):
    im[x,:] = np.sin(x)
images.append(im)
im = np.zeros((h,w))
for y in range(w):
    im[:,y] = np.sin(y)
images.append(im)

(5) 接下来，我们生成对角线周期图像，然后使用半径为 10 的圆形掩码(以图像的中心为圆心)，从根据此图像中创建另一图像：

im = np.zeros((h,w))
for x in range(h):
    for y in range(w):
        im[x,y] = np.sin(x + y)
images.append(im)
im = np.zeros((h,w))
for x in range(h):
    for y in range(w):
        if (x-h/2)**2 + (y-w/2)**2 < 100:
            im[x,y] = np.sin(x + y)
images.append(im)

(7) 最后，生成另一对图像，第一个图像带有一个实心圆，第二个带有以图像中心为中心的实心正方形：

im = np.zeros((h,w))
for x in range(h):
    for y in range(w):
        if (x-h/2)**2 + (y-w/2)**2 < 25:
            im[x,y] = 1
images.append(im)
im = np.zeros((h,w))
im[h//2 -5:h//2 + 5, w//2 -5:w//2 + 5] = 1
images.append(im)

(8) 我们将所有图像都存储在列表中，使用 numpy.fft 的 fft2() 函数使用快速傅里叶变换算法来计算 DFT，并绘制生成的每个输入图像的输出功率谱：

numpy.fft.fft2(a, s=None, axes=(-2, -1), norm=None)

下图显示了相应图像及其 DFT：

plt.figure(figsize=(25,10))
i = 1
for im in images:
    plt.subplot(2,6,i), plot_image(im, 'image {}'.format(i))
    plt.subplot(2,6,i+6), plot_freq_spectrum(fp.fft2(im), 'DFT {}'.format(i), show_axis=False)
    i += 1
plt.tight_layout()
plt.show()

理想情况下，对于图像 1 和 2，水平和垂直周期模式的 DFT 应分别为垂直和水平对齐的点(对应于频率)。但是，如上图所示，由于边缘效应，结果为一条垂直和水平线(线上的亮点对应于相应的频率)。

同样，对于图像 3，周期对角线模式应在垂直对角线方向上产生两个对角线点，但由于边缘效应，得到了其他许多线。如果使用二值圆形掩码处理图像 3 中心的圆，则会降低边缘效应，我们可以看到沿主对角线方向的亮点，对应于频率分量。