离散傅里叶变换
- 傅里叶变换
- 代码实现
- 扩展图片
- 创建储存实部和虚部值的矩阵
- 进行离散傅里叶变换
- 将复数转换成振幅
- 对数转换
- 裁剪和重排
- 归一化
- 离散傅里叶变换在图像处理中的应用
- 参考
傅里叶变换
在图片处理中，傅里叶变化会将对图片的时域分析转变为频域分析。
傅里叶的基本思路就是，任何函数都可以近似地变成无限个
sin
\sin
sin和
cos
\cos
cos函数的和。对于具有
N
−
1
N-1
N−1个采样点的离散信号，可以将其中每个采样点的幅值与频率为
k
k
k的
sin
\sin
sin或
cos
\cos
cos函数中对应点的幅值相乘，就得到了该采样点在频率为
k
k
k的
sin
\sin
sin或
cos
\cos
cos函数上的变换结果。将频率
k
k
k进行变化，又可以得到新的结果，将所有这些结果相加，就得到了在该采样点上的离散傅里叶变换结果。用数学公式来表示：
F
(
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
e
−
i
∗
k
n
(
2
π
N
)
F(k) = \sum\limits_{n=0}^{N-1}f(n)e^{-i*kn(\frac{2\pi}{N})}
F(k)=n=0∑N−1f(n)e−i∗kn(N2π)
根据欧拉公式(
e
i
x
=
cos
x
+
i
sin
x
e^{ix}=\cos x +i\sin x
eix=cosx+isinx)展开得：
F
(
k
)
=
∑
n
=
0
N
−
1
f
(
n
)
[
cos
k
n
(
2
π
N
)
−
i
sin
k
n
(
2
π
N
)
]
F(k) = \sum\limits_{n=0}^{N-1}f(n)[\cos{kn(\frac{2\pi}{N})} -i\sin{kn(\frac{2\pi}{N})}]
F(k)=n=0∑N−1f(n)[coskn(N2π)−isinkn(N2π)]
- k k k为 sin \sin sin或 cos \cos cos函数的角频率
- F ( k ) F(k) F(k)是频率为 k k k时傅里叶变换的结果
- n n n为第n个采样点
- f ( n ) f(n) f(n)为离散信号在第 n n n个采样点上的幅值
- i i i为复数中的i，即 − 1 \sqrt{-1} −1
- N N N为采样点的总数+1
对于图片数据来说，像素就是采样点，像素上的值就相当于采样点的幅值。所以在图片上的傅里叶变换，有以下公式：
F
(
k
,
l
)
=
∑
m
=
0
N
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
e
−
i
2
π
(
k
m
N
+
l
n
N
)
F(k,l) = \sum\limits_{m=0}^{N-1} \sum\limits_{n=0}^{N-1} f(m,n)e^{-i2\pi(\frac{km}{N}+\frac{ln}{N})}
F(k,l)=m=0∑N−1n=0∑N−1f(m,n)e−i2π(Nkm+Nln)
根据欧拉公式展开可得：
F
(
k
,
l
)
=
∑
m
=
0
N
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
cos
2
π
(
k
m
N
+
l
n
N
)
−
i
sin
2
π
(
k
m
N
+
l
n
N
)
F(k,l) = \sum\limits_{m=0}^{N-1} \sum\limits_{n=0}^{N-1} f(m,n)\cos{2\pi(\frac{km}{N}+\frac{ln}{N})} -i \sin {2\pi(\frac{km}{N}+\frac{ln}{N})}
F(k,l)=m=0∑N−1n=0∑N−1f(m,n)cos2π(Nkm+Nln)−isin2π(Nkm+Nln)
- 由于图片数据为二维数组， f ( m , n ) f(m,n) f(m,n)代表m行n列的像素值
- k k k和 l l l分别代表应用在 k k k行和 l l l列上的频率
- F(k, l)表示频率分别为 k k k和 l l l时求得的傅里叶变换结果
- i i i仍然为 − 1 \sqrt{-1} −1
也可以说 f f f是时域上的像素值，而 F F F是频域上的像素值。傅里叶变换的结果是一个复数。在图片处理算法中，为了方便查看变换结果，一般需要将复数转换成振幅图片(magnitude image)。它虽然只能展示每个像素值的信息，不能展示频率或相位的信息，但这些也不是我们想在图片中看到的。所以，这里还是用振幅图片来展示傅里叶变换。
电子图片是离散信号，其中的像素值都是有特定值域的。比如说一张灰度图片中所有的像素值都在0到255之间。所以对图片进行的傅里叶转换用到的是离散傅里叶变换（DFT）。
将一张灰度图片进行灰度转换需要以下步骤：
- 扩展图片
- 创建储存实部和虚部的矩阵
- 进行离散傅里叶变换
- 将复数转换为振幅
- 对数转换
- 裁剪和重排
- 归一化
代码实现
扩展图片
图片的尺寸会影响DFT的运算。当图片尺寸是2、3或5的倍数的时候，DFT的运算速度最快。所以对原始图片进行适当的扩展，将会提高运算速度。getOptimalDFTSize()函数能根据输入的图片尺寸计算最佳尺寸，copyMakeBorder()函数则可以扩展图片（新增的像素全设为0）。
Mat padded;
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
copyMakeBorder(I, padded, 0, m-I.rows, o, n-I.cols, BORDER_CONSTANT, Scalaar::all(0));
copyMakeBorder()函数的API如下：
void cv::copyMakeBorder(InputArray src,
OutputArray dst,
int top,
int bottom,
int left,
int right,
int borderType,
const Scalar& value = Scalar())
src原始图片dst输出图片，其尺寸为(src.cols+left+right, src.rows+top+bottom) .top顶部扩展的像素个数bottom底部扩展的像素个数left左边……right右边……borderType边框类型；扩展的像素像边框一样包围着原始图片value边框类型为BORDER_CONSTANT时，边框内的像素的值
创建储存实部和虚部值的矩阵
傅里叶变换的结果是复数，所以其中每个像素值都有两部分——实部和虚部。而且，频域的范围要比时域大很多，所以，其数据类型至少得是浮点型。以下代码，将单通道的浮点数矩阵扩展成双通道的矩阵，用来同时储存结果的实部和虚部。
Mat planes[]{ Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
merge函数将数组中的矩阵合并为一个多通道的矩阵，其API如下：
void cv::merge( const Mat * mv,
size_t count,
OutputArray dst)
mv矩阵数组；数组中的所有矩阵必须有相同的尺寸。count数组中矩阵的数量，必须大于0。dst输出矩阵与矩阵数组中的第一个矩阵的尺寸相同；通道数量与count相同。
傅里叶变换的复数结果矩阵示意图如下：
进行离散傅里叶变换
OpenCV中已经有现成的DFT函数：
dft(complexI, complexI);
第1个参数为输入矩阵，第2个参数为输出矩阵。这里用同一个矩阵对象来储存计算结果。
将复数转换成振幅
对于复数
z
=
x
+
i
y
z=x+iy
z=x+iy，它的模，即
∣
z
∣
=
x
2
+
y
2
|z|=\sqrt{x^2+y^2}
∣z∣=x2+y2。因此在图片数据中，将DFT算出的复数结果转换成振幅需要以下计算：
M
=
R
e
(
D
F
T
(
I
)
)
2
+
I
M
(
D
F
T
(
I
)
)
2
M=\sqrt{Re(DFT(I))^2+IM(DFT(I))^2}
M=Re(DFT(I))2+IM(DFT(I))2
- R E RE RE为复数的实部
- I M IM IM为复数的虚部
- D F T ( I ) DFT(I) DFT(I)为矩阵I的DFT结果
- M M M为振幅结果
在OpenCV中可用以下代码实现：
split(complexI, planes); //planes[0] = Re(DFT(I)), planes[1] = Img(DFT(I))
magnitude(planes[0], planes[1], planes[0]);
Mat magI = planes[0];
其中，split函数将多通道矩阵变成几个单通道矩阵。如上面的代码中，complexI为需要进行分裂的矩阵，planes为接收分裂结果的矩阵数组中，有2个矩阵，一个储存了DFT结果的实部，一个储存了DFT结果的虚部。
magnitude函数计算振幅结果，planes[0]为实部矩阵，planes[1]为虚部矩阵，后面一个planes[0]是用来储存计算结果的矩阵。
对数转换
由于傅里叶系数的动态值域太宽，无法在屏幕上显示，太大的值会变成白点，太小的值会变成黑点。所以要对计算结果进行对数转换，以缩小其值域：
M
1
=
ln
(
1
+
M
)
M_1=\ln(1+M)
M1=ln(1+M)
magI += Scalar::all(1);
log(magI, magI);
log函数对第一个参数进行取自然对数的运算，然后将结果储存在第二个参数中。
裁剪和重排
因为一开始我们将原始图片进行了扩展，所以现在要进行相应的裁剪。
而且为了更好地呈现结果，还需要将结果分成4个大小相等的矩形区域，并进行重排，好让原来的原点能够在中心。
//通过按位与运算来确定裁剪后的列数和行数
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
//将傅里叶图片的4个1/4进行重新排列，以使原点处于图片中心
int cx = magI.cols/2;
int cy = magI.rows/2;
Mat q0(magI, Rect(0, 0, cx, cy)); // 左上角1/4
Mat q1(magI, Rect(cx, 0, cx, cy)); // 右上角1/4
Mat q2(magI, Rect(0, cy, cx, cy)); // 左下角1/4
Mat q3(magI, Rect(cx, cy, cx, cy)); // 右下角1/4
Mat tmp; // 左上角和右下角交换
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 右上角和左下角交换
q2.copyTo(q1);
tmp.copyTo(q2);
归一化
傅里叶图片中的数值仍然是超出显示范围的浮点数，所以最后还需要对其进行归一化，使其变成0-1之间的浮点数
normalize(magI, magI, 0, 1, NORM_MINMAX);
这里使用normalize()函数进行归一化操作，其中第一个参数为输入图片，第二个参数为输出图片，第三、四个参数分别为归一化的下限和上限，最后一个参数为归一化类型。
离散傅里叶变换在图像处理中的应用
离散傅里叶变换可以用来呈现图片中的几何方向。例如，检测图片中的文字或其他对象是否是水平的。
- 当文字是水平的时候，傅里叶变换的结果如下：
- 当文字有一定的倾斜时，傅里叶变换的结果如下：
频域中的主要成分（亮点部分）与图片中的文本对象的倾斜方向是一致的。这样就可以计算倾斜的角度，从而进行相应的对齐操作。
再比如，下面的原图中，山和岸都有点向右倾斜，傅里叶变换后得到的频域图中的主要成分也像右倾斜：
- 原图：
- 傅里叶变换频域图：
参考
- Discrete Fourier Transform, The Core Functionality (core module), OpenCV Tutorials
- 《深度实践OCR：基于深度学习的文字识别》3.1.1.2 傅里叶特征算子
