图像的傅里叶变换
图像的傅里叶变换
备份自：http://blog.rainy.im/2015/11/03/fourier-transform-in-image-processing/
傅里叶变换（Fourier Transform）是非常重要的数学分析工具，同时也是一种非常重要的信号处理方法。我记得本科课程电路原理中有提到过，但由于计算过于复杂，好像是超出考试范围了，所以并没有深入学习。最近实验中需要对图像进行滤波处理，文献中提到的方法通常是经过傅里叶变换之后对频域进行过滤，将图像中的低频信息与高频信息区分开来。
理解傅里叶变换对非数学专业的人来说比较难以理解的原因主要有两方面。首先，由于涉及到比较复杂繁琐的数学操作，看到下面的这两个公式，一般人可能当时就蒙了：
$$
\widehat{f}(t) = \int_{-\infty}^{\infty} f(x) e^{-2\pi ixt}dx;
$$
$$
f(x) = \int_{-\infty}^{\infty}\widehat f (t) e^{2\pi itx}dt;
$$
另一方面的原因则是由于变换过程比较抽象，很难从直觉上去把握在傅里叶变换过程中到底发生了什么。关于傅里叶变换的科普文章，有一篇是知乎上传阅较广的《傅里叶分析之掐死教程》，作者用了尽可能少的数学公式和图形分解来解决这两方面的问题。在此之前，国外有个专门科普数学概念的网站（Better Explained）也写了一篇类似的科普文章，但更彻底的是，全文都不涉及到任何数学公式和推断，完全用英语来向读者解释傅里叶变换过程。第一篇文章以音乐和乐谱为例，第二篇作者用“奶昔”作为例子，我想了半天终于找到一个更通俗的例子：
我们来想象一下，假设这块面饼的厚度是$N$层面条，每一层都是由一根弯曲成正弦曲线形状的面条排列而成，有些面条波浪较大，也就是排列较为稀疏，而密有些排列较密集。我们把这$N$层面条挤压到一起，就得到上图这一块杂乱无序、世间独一无二的面饼。傅里叶变换所做的事就是把上面的过程反过来，我们可以从一块完整的面饼得到最初的$N$层面条。如果有些人的口味比较特殊，喜欢波浪大又稀疏的面，于是我们就将排列太紧凑（高频面）剔除之后再重新压制成一块新的面饼，这就是我们最终想要的滤波（Filter）的过程。
带着对面饼的想象，我们来看一种更为抽象、优雅的描述（from Wikipedia）：
一般的波形或者说信号（Signal）都是基于时间尺度上的采样结果，因此也称为时域（Time Domain），而上面泡面的例子和我们将要处理的图像信号则是基于空间尺度上的采样，但好像并没有“空域（Space Domain）”这一说，毕竟我们对空间的感知仍然依赖于时间。不过在空间尺度上我们可以更直观地认为信号是静止，例如下面这张图像（灰度图），其实是由250x250个像素点组成，每个像素点的灰度值（$[0, 255]$）就是基于像素坐标的空间采样的结果：
右边的3D Fourier就是一块长相奇怪的面饼。
实践
对傅里叶变换有了大概的了解之后可以先动手尝试一下，来更加直观地感受一下（实际上完全可以在不理解的情况下，直接上手）。这里用到的是OpenCV + Numpy，实际上OpenCV和Numpy都提供了快速傅里叶变换（FFT）算法：
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('Joseph_Fourier_250.jpg', 0)
f = np.fft.fft2(img) # 快速傅里叶变换算法得到频率分布
fshift = np.fft.fftshift(f) # 默认结果中心点位置是在左上角，转移到中间位置
fimg = np.log(np.abs(fshift)) # fft 结果是复数，求绝对值结果才是振幅
# 展示结果
plt.subplot(121), plt.imshow(img, 'gray'), plt.title('Original Fourier')
plt.subplot(122), plt.imshow(fimg, 'gray'), plt.title('Fourier Fourier')
plt.show()
右边图就是频率分布图谱，其中越靠近中心的位置频率越低，越亮（灰度值越高）的位置代表该频率的信号振幅越大。fft的结果是复数形式，保留了图像的全部信息，但去绝对值得到的频谱图只表现了振幅而没有体现相位。
回想一下高中时候学过的三角函数：
$$
f(x) = A sin(\omega x+\varphi) = A sin(2\pi fx + \varphi)
$$
一个正弦波是由下面三个参数决定的：
- 角速度（频率）$\omega = 2\pi f$ ；
- 振幅$A$；
- 相位$\varphi$。
除了上面这个公式之外，还可以用另外一种形式来（唯一地）表示一个正弦波（from BetterExplained）：
即：
$$
cos(x) + i sin(x) \Leftrightarrow a + i b
$$
所以说，fft的复数结果保留了正弦波成分的所有信息，但频谱图只展现了频率和振幅的分布。因此可以根据fft的结果还原原始图像，但是我们做傅里叶变换的目的并不是为了观察图像的频率分布（至少不是最终目的），更多情况下是为了对频率进行过滤。过滤的方法一般有三种：低通（Low-pass）、高通（High-pass）、带通（Band-pass）。所谓低通就是保留图像中的低频成分，过滤高频成分，可以把过滤器想象成一张渔网，根据上文对频谱图的解读，想要低通过滤器，就是将高频区域的信号全部拉黑，而低频区域全部保留：
img = cv.imread('Joseph_Fourier_250.jpg', 0)
f = np.fft.fft2(img) # 快速傅里叶变换算法得到频率分布
fshift = np.fft.fftshift(f) # 默认结果中心点位置是在左上角，转移到中间位置
lpButterMask = butterWorthLPF(img.shape[:2], 24, 2)
hpButterMask = butterWorthHPF(img.shape[:2], 36, 2)
lpFshift = fshift * lpButterMask
maskedInvf = np.fft.ifft2(np.fft.ifftshift(lpFshift))
lpfImg = np.abs(maskedInvf)
hpFshift = fshift * hpButterMask
maskedInvf = np.fft.ifft2(np.fft.ifftshift(hpFshift))
hpfImg = np.abs(maskedInvf)
plt.subplot(221), plt.imshow(lpButterMask, 'gray'), plt.title('Butterworth LPF')
plt.subplot(222), plt.imshow(lpfImg, 'gray'), plt.title('LPF Image')
plt.subplot(223), plt.imshow(hpButterMask, 'gray'), plt.title('HPF Spectrum')
plt.subplot(224), plt.imshow(hpfImg, 'gray'), plt.title('HPF Image')
plt.show()
很显然，滤波器的选择也是很重要的，这里用到的是 Butterworth 滤波器，有兴趣的可以自己实现一下:P
总结
傅里叶变换真是又伟大、又深奥、又方便！至少对于一个微积分已经忘记差不多、三角函数公式都要搜索半天才能回忆起来的“文科生”来说，从头学习一遍简直是又有虐脑又有惊喜！这篇文章可能更加偏重于记录我自己的消化过程，如果想要更加细致、深入深刻以及深入浅出地介绍，请参考下方参考链接。
这两天被虐脑的感觉真是酸爽，打算接续这一篇下去写一个专题，就叫#BlowYourMind3000#，哈哈哈，下篇预告《贝叶斯》。
参考
-
Schyns, P. G., & Oliva, A. (1994). From blobs to boundary edges: evidence for time- and spatial-scale-dpendent scene recognition.Psychological Science*, 5*(5), 195-200.
