(OpenCV — 10)离散傅里叶变换
离散傅里叶变换 ( Discrete Fourier Transfonn, 缩写为 DFT ) , 是指傅里叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅里叶变换 ( DTFT ) 频域的采样 。 在形式上 ， 变换两端（ 时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号 的主值序列 。 即使对有限长的离散信号做 DFT, 也应当对其经过周期延拓成为周期信号再进行变换 。 在实际应用中，通常采用快速傅里叶变换来高效计算 DFT 。
离散傅里叶变换的原理
简单来说 ， 对一张图像使用傅里叶变换就是将它分解成正弦和余弦两部分 ，也就是将图像从空间域 ( spatial domain ) 转换到频域 ( frequency domain ) 。
这一转换的理论基础为 ： 任一 函数都可以表示成无数个正弦和余弦函数的和的形式。傅里叶变换就是一个用来将函数分解的工具。
二维图像的傅里叶变换可以用以下数学公式表达。
式中 是空间域 ( spatial domain ) 值， F 是频域 ( frequency domain ) 值 。转换之后 的 频域值是复数，因 此 ，显示傅里叶变换之后的结果需要使用实数图像 ( real image ) 加虚数图像 ( complex image ) , 或者幅度图像 ( magitude image )加相位图像 ( phase image ) 的形式 。 在实际 的 图像处理过程中，仅仅使用了幅度图像，因为幅度图像包含了原图像的几乎所有我们需要的几何信息 。 然而，如果想通过修改幅度图像或者相位图像的方法来间接修改原空间图像 ， 需要使用逆傅里叶变换得到修改后的空间图像，这样就必须同时保留幅度图像和相位图像了。
(时域)实数图像 + 虚数图像 = (频域)幅度图像 + 相位图像
在此示例中，我们将展示如何计算以及显示傅里叶变换后的幅度图像 。 由于数字图像 的离散性 ， 像素值的取值范围也是有限的。比如在一张灰度图像中 ， 像素灰度值一般在 0 到 255 之间 。 因此，我们这里讨论的也仅仅是离散傅里叶变换( DFT ) 。如果需要得到图像中的几何结构信息 ， 那你就要用到它了 。
在频域里面 ， 对于一幅图像，高频部分代表了图像的细节 、 纹理信息；低频部分代表了图像的轮廓信息 。 如果对一幅精细的图像使用低通滤波器，那么滤波后的结果就只剩下轮廓了 。 这与信号处理的基本思想是相通的 。 如果图像受到的噪声恰好位千某个特定的＂频率 ”范围内，则可以通过滤波器来恢复原来的图像 。傅里叶变换在图像处理中可以做到图像增强与图像去噪、图像分割之边缘检测、图像特征提取、图像压缩等。
dft()函数详解
dft 函数的作用是对一 维或二维浮点数数组进行正向或反向离散傅里叶变换。
void dft(InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0);- 第一个参数， lnputArray 类型的 src 。 输入矩阵，可以为实数或者虚数 。
- 第二 个参数，OutputArray 类型的 dst。函数调用后的运算结果存在这里，其尺寸和类型取决于标识符 ，也就是第三个参数 flags 。
- 第三个参数 ， int 类型的 flags。转换 的标识符，有默认值 0 , 取值可以为表中标识符的结合。
dft 标识符取值列表
第四个参数， int 类型的 nonzeroRows, 有默认值 0。当此参数设为非零时（最好是取值为想要处理的那一行的值，比如 C.rows ), 函数会假设只有输入矩阵的第一个非零行包含非零元素（没有设置 DFT_INVERSE 标识符） ，或只有输出矩阵的第一个非零行包含非零元素（设置了 DFT_INVERSE 标识符） 。 这样的话 ， 函数就可对其他行进行更高效的处理，以节省时间开销 。这项技术尤其是在采用 DFT 计算矩阵卷积时非常有效 。
讲解完函数的参数含义，下面我们看一个用 dft 函数计算两个二维实矩阵卷积的示例核心片段 。
void convolveDFT (InputArray A, InputArray B, OutputArray C)
{
//[1]初始化输出矩阵
C.create(abs(A.rows - B.rows)+1 , abs(A.cols - B.cols )+1 , A.type());
Size dftSize ;
//[2]计算DFT变换的尺寸
dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
dftSize.height = getOptimalDFTSize(A.rows + B.rows - l) ;
//[3]分配临时纹冲区并初始化置零
Mat tempA(dftSize, A.type() , Scalar::all(0));
Mat tempB(dftSize, B.type() , Scalar::all(0));
//[4]分别复制 A 和 B 到 tempA 和 tempB 的左上角
Mat roiA(tempA, Rect(0 , 0, A.cols, A.rows));
A.copyTo(roiA);
Mat roiB(tempB, Rect(0 , 0, B.cols, B.rows));
B.copyTo(roiB);
//[5]就地操作 (in-place), 进行快速傅里叶变换 ， 并将 nonzeroRows 参数置为非零，以进行更加快速地处理
dft(tempA, tempA, 0 , A.rows);
dft(tempB, tempB, 0 , B.rows);
//[6] 将得到的频谱相乘，结果存放于 tempA 中
mulSpectrums (tempA, tempB, tempA);
//[7]将结果文换为频域 ， 且尽管结果行 (result rows) 都为非零 ， 我们只需要其中的C.rows 的 第 一行，所以采 用 nonzeroRows == C . rows
dft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows) ;
//[8]将结果复制到 C 中
tempA(Rect(0 , 0 , C.cols, C.rows)).copyTo(C);
//所有的临时缓冲区将被自动释放，所以无须收尾操作
} 此示例中的注释 已经十分详尽。其中出现了新的 函数 MulSpectrums , 它 的作用是计算两个傅里叶频谱的每个元素的乘法 ， 前两个参数为输入的参加乘法运算的两个矩阵，第三个参数为得到的乘法结果矩阵。
返回 DFT 最优尺寸大小： getOptimalDFTSize()函数
getOptimalDFTSize 函数返回给定向量尺寸的傅里叶最优尺寸大小 。 为了提高离散傅里叶变换的运行速度 ， 需要扩充图像，而具体扩充多少，就由这个函数来计算得到 。
int getOptimalDFTSize(int vecsize);此函数的唯一一个参数为 int 类型的 vecsize, 向量尺寸，即图片的 rows 、 cols 。
扩充图像边界： copyMakeBorder()函数
void copyMakeBorder(InputArray src,
OutputArray dst,
int top,
int bottom,
int left,
int right,
int borderType,
const Scalar& value = Scalar() );- 第一个参数 ， lnputAnay 类型的 src , 输入图像 ，即源图像，填 Mat 类的对象即可。
- 第二个参数 ， OutputArray 类型的 dst, 函数调用后的运算结果存在这里 ，即这个参数用千存放函数调用后 的输出 结果，需和源 图片有一样的尺寸和类型 ， 且 s ize 应该为 Size ( src.cols+left+right, src.rows+top+bottom ) 。
- 接下来的 4 个参数分别为 int 类型的 top 、 bottom 、 left 、 right , 分别表示在源图像的四个方 向 上扩充多少像素 ，例如 top=2, bottorn=2, left=2, right=2就意味着在源图像的上下左右各扩充两个像素宽度的边界 。
- 第七个参数 ， borderType 类型的 ， 边界类型，常见取值为BORDER_CONSTANT, 可参考 borderlnter-polate()得到更多的细节。
- 第八个参数 ， const Scalar&类型 的 value , 有默认值 Scalar(), 可以理解为默认值为 0。当 borderType 取值为 BORDER_CONSTANT 时， 这个参数表示边界值 。
计算二维矢量的幅值： magnitude()函数
void magnitude(InputArray x, InputArray y, OutputArray magnitude);第一个参数 ， InputArray 类型 的 X, 表示矢量 的浮点型 X 坐标值 ，也就是实部 。
第一个参数 ， LnputArray 类型的 y, 表示矢量的浮点型 Y 坐标值，也就是虚部 。
第三 次参数 ， OutputArray 类型的 magnitude , 输出的幅值，它和第一个参数 x 有着同样的尺寸和类型。
下式可以表示 magnitude()函数的原理 ：
计算自 然对数 ： log()函数
void log (InputArray src, OutputArray dst）第一个参数为输入图像，第二个参数为得到的对数值 。其原理如下 。
矩阵归一化： normalize()函数
void normalize(InputArray src,
OutputArray dst,
double alpha=1 ,
double beta=0 ,
int norm_type=NORM_L2 ,
int dtype=-1,
InputArray mask=noArray());- 第一个参数 ， lnputArray 类型的 src 。 输入图像 ，即 源图像，填 Mat 类的对象即可 。
- 第二个参数， OutputArray 类型的 dst。 函数调用后的运算结果存在这里，和源图片有一样的尺寸和类型。
- 第三个参数 ， doubl e 类型的 alpha 。 归一化后的最大值，有默认值 1 。
- 第四个参数， double 类型的 beta 。 归一化后的最大值，有默认值 0 。
- 第五个参数， int 类型 的 nonn_type。归一化类型 ， 有 NORM_INF 、NORM_L l 、NORM_L2 和 NORM_MfNMAX 等参数可选 ，有默认值 NORM_L2 。
- 第六个参数 ， int 类型 的 dtype, 有默认值4 。当此参数取负值时，输出矩阵和 src 有同样的类型 ， 否则 ，它和 src 有同样的通道数，且此时图像深度为CV_MAT_DEPTH ( dtype ) 。
- 第七个参数 ， .lnputArray 类型的 mask , 可选的操作掩膜，有默认值 noArray() 。
示例程序 ： 离散傅里叶变换
下面我们将学习一个 以 dft()函数为核心，对图像求傅里叶变换的 ， 有详细注释的示例程序 。
在此示例中，将展示如何计算以及显示傅里叶变换后的幅度图像 。 由于数字图像的离散性，像素值的取值范围也是有限的。比如在一张灰度图像中 ， 像素灰度值一般在 0 到 255 之间 。因此 ， 我们这里讨论的也仅仅是离散傅里叶变换( DFT) 。如果需要得到图像中的几何结构信息，那么就要用到离散傅里叶变换了 。 下面的步骤将以输入图像为单通道 的灰度图像 I 为例，进行分步说明 。
【 第一步 】 载入原始图像
我们在这一步 以灰度模式读取原始图像，进行是否读取成功的检测，并显示出读取到的图像 。 代码如下 。
//【1】以灰度模式读取原始图像并显示
Mat srcImage = imread("1.jpg", 0);
if(!srcImage.data ) { printf("读取图片错误，请确定目录下是否有imread函数指定图片存在~！ \n"); return false; }
imshow("原始图像" , srcImage); 【第二步】将图像扩大到合适的尺寸
离散傅里叶变换的运行速度与图片的尺寸有很大关系。当图像的尺寸是 2 、 3 、5 的 整数倍时，计算速度最快 。 因此，为了达到快速计算的目的， 经常通过添凑新的边缘像素的方法获取最佳图像尺寸。函数 getOptimalDFTSize()用于返回最佳尺寸，而函数 copyMakeBorder()用于填充边缘像素 ，这一步代码如下 。
//【2】将输入图像延扩到最佳的尺寸，边界用0补充
int m = getOptimalDFTSize( srcImage.rows );
int n = getOptimalDFTSize( srcImage.cols );
//将添加的像素初始化为0.
Mat padded;
copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));【第三步】为傅里叶变换的结果 （ 实部和虚部 ） 分配存储空间
傅里叶变换的结果是复数，这就是说对于每个原图像值，结果会有两个图像值。此外，频域值范围远远超过空间值范围， 因此至少要将频域储存在 float 格式中。所以我们将输入图像转换成浮点类型，并多加一个额外通道来储存复数部分。
//【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
//将planes数组组合合并成一个多通道的数组complexI
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); 【第四步】进行离散傅里叶变换
这里的离散傅里叶变换为图像就地计算模式（ in-place , 输入输出为同一 图像） 。
//【4】进行就地离散傅里叶变换
dft(complexI, complexI); 【第五步】将复数转换为幅值
复数包含实数部分 ( Re ) 和虚数部分 ( imaginary - Im ) 。离散傅里叶变换的结果是复数 ， 对应的幅度可以表示为：
那么 ， 转化为 OpenCV 代码，就是下面这样的。
//【5】将复数转换为幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组，planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
Mat magnitudeImage = planes[0];【第六步】进行对数尺度 ( logarithmic scale ) 缩放 ．
傅里叶变换的幅度值范围大到不适合在屏幕上显示。高值在屏幕上显示为白点 ，而低值为黑点，高低值的变化无法有效分辨。为了在屏幕上凸显出高低变化的连续性，我们可 以用对数尺度来替换线性尺度，公式如下。
而写成 OpenCV 代码，就是下面这样的 。
//【6】进行对数尺度(logarithmic scale)缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage);//求自然对数【第七步】 剪切和重分布幅度图象限
因为在第二步中延扩了图像，那现在是时候将新添加的像素剔除了。为了方便显示，也可以重新分布幅度图象限位置（注 ： 将第五步得到 的幅度图从中间划开 ， 得到 4 张 1 /4 子图像 ， 将每张子图像看成幅度图的一个象限，重新分布，即将 4 个角点重叠到图片中心） 。 这样的话原点 ( 0,0 ) 就位移到图像中心 了。 OpenCV代码如下 。
//【7】剪切和重分布幅度图象限
//若有奇数行或奇数列，进行频谱裁剪
magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
//重新排列傅立叶图像中的象限，使得原点位于图像中心
int cx = magnitudeImage.cols/2;
int cy = magnitudeImage.rows/2;
Mat q0(magnitudeImage, Rect(0, 0, cx, cy)); // ROI区域的左上
Mat q1(magnitudeImage, Rect(cx, 0, cx, cy)); // ROI区域的右上
Mat q2(magnitudeImage, Rect(0, cy, cx, cy)); // ROI区域的左下
Mat q3(magnitudeImage, Rect(cx, cy, cx, cy)); // ROI区域的右下
//交换象限（左上与右下进行交换）
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
//交换象限（右上与左下进行交换）
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);【第八步】归一化
这一步仍然是为了显示 。 现在有了重分布后的幅度图，但是幅度值仍然超过可显示范围 [0, 1] 。 我们使用 normalizeO函数将幅度归一化到可显示范围。
//【8】归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
//此句代码的OpenCV2版为：
//normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX);
//此句代码的OpenCV3版为:
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX); 【第九步】显示效果图
处理完成 ， 最后进行显示即 可。
//【9】显示效果图
imshow("频谱幅值", magnitudeImage);
waitKey();将上述代码组合到一起 ， 便得到了程序的完整源代码。
原始图：
离散傅里叶变换效果图