《OpenCV3编程入门》学习笔记5 Core组件进阶（五）离散傅里叶变换（DFT）

第5章 Core组件进阶

5.5 离散傅里叶变换（Discrete Fourier Transform，DFT）

5.5.1 离散傅里叶变换原理

1.对一张图像使用傅里叶变换就是把它分解成正弦和余弦，将图像从空间域（spatial domain）转换到频域（frequency domain）
2.理论基础：任一函数都可以表示成无数个正弦和余弦函数的和的形式
3.二维图像的傅里叶变换数学公式表达：
              ,  
  f是空间域值（spatial domain），F是频域值（frequency domain）,转换之后的频域值是复数，因此，显示傅里叶变换之后的结果需要使用实数图像（real image）加虚数图像（complex image），或幅度图像（magitude image）加相位图像（phase image）。
4.频域里面，高频部分代表图像的细节、纹理信息；低频部分代表图像的轮廓结构，如果对一幅精细的图像使用低通滤波器，那么滤波后的结果会只剩下轮廓，如果图像受到的噪声正好位于某个特定的“频率”范围内，则可以通过滤波器来恢复原来的图像。
5.傅里叶变换在图像处理中可以做到图像增强与图像去噪、图像分割之边缘检测、图像特征提取、图像压缩等。

5.5.2 dft()函数

1.作用：对一维或二维浮点数组进行正向或反向离散傅里叶变换
2.原型：

void dft(InputArray src,OutputArray dst,int flags=0,int nonzeroRows=0)

3.参数说明：
（1）输入矩阵
（2）运算结果
（3）转换标识符（默认0）
            
（4）参数设为非零时，函数会假设只有输入矩阵的第一个非零行包含非零元素（没有设置DFT_INVERSE标识符）或只有输出矩阵的第一个非零行包含非零元素（设置了DFT_INVERSE标识符），如此函数可以对其他行进行更高效的处理，以节省时间开销，尤其在采用DFT计算矩阵卷积时非常有效。

4.示例：用dft函数计算两个二维实矩阵卷积

void consolveDFT(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 - 1);

	//【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);

	//所有的临时缓冲区将被自动释放，所以无须收尾操作
}

5.5.3 返回DFT最优尺寸大小：getOptimalDFTSize()函数

  返回给定向量尺寸的傅里叶最优尺寸大小，为了提高离散傅里叶变换的运行速度，需要扩充图像，具体扩充多少，由此函数计算得到：
                 int getOptimalDFTSize(int vecsize)

5.5.4 扩充图像边界：copyMakeBorder()函数

1.原型：

void copyMakeBorder(InputArray src,OutputArray dst,int top,int bottom,int left,int right,int borderType,const Scalar&value=Scalar())

2.参数说明：
（1）源图像
（2）运算结果，需和源图片有一样的尺寸和类型，Size(src.cols+left+right,src.rows+top+bottom)
（3）top，bottom，left，right，分别表示在源图像的四个方向上扩充多少像素
（4）边界类型，常见取值为BORDER_CONSTANT
（5）默认值Scalar(),可理解为默认值为0，borderType取值为BORDER_CONSTANT，这个参数表示边界值

5.5.5 计算二维矢量的幅值：magnitude()函数

1.原型：

void magnitude(InputArray x,InputArray y,OutputArray magnitude)

2.参数说明：
（1）表示矢量的浮点型X坐标值（实部）
（2）表示矢量的浮点型Y坐标值（虚部）
（3）输出的幅值，它和第一个参数x有着同样的尺寸和类型
3.原理：
                

5.5.6 计算自然对数：log()函数

1.作用：计算每个数组元素绝对值的自然对数
2.原型：

void log(InputArray src,OutputArray dst)

3.参数说明：输入图像，得到的对数值
4.原理：
                
5.5.7 矩阵归一化：normalize()函数

1.原型：

void normalize(InputArray src,OutputArray dst,double alpha=1,double beta=0,int norm_type=NORM_L2,int dtype=-1,InputArray mask=noArray() )

2.参数说明：
（1）输入图像
（2）运算结果，和源图片有一样的尺寸和类型
（3）归一化后的最大值，默认值1
（4）归一化后的最大值，默认值0
（5）归一化类型，有NORM_INF、NORM_L1、NORM_L2（默认）和NORM_MINMAX等参数可选
（6）默认值-1，当此参数取负值时，输出矩阵和src有同样的类型，否则它和src有同样的通道数，且此时图像深度为CV_MAT_DEPTH（dtype）
（7）可选的操作掩膜，默认值noArray()

5.5.8 示例程序：离散傅里叶变换

#include<opencv2/opencv.hpp>
#include<opencv2/core/core.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<iostream>
using namespace cv;
using namespace std;

int main()
{
	//[1]以灰度模式读取原始图像并显示
	Mat srcImage = imread("love.jpg", 0);
	if (!srcImage.data) { printf("读取图片错误，请确定目录下是否有imread函数指定图片存在~！\n"); return  false; }
	imshow("原始图像", srcImage);

	//[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));

	//[3]为傅里叶变换的结果（实部和虚部）分配存储空间
	//傅里叶变换结果是复数，即每个原图像值结果会有两个图像值,且频域值范围至少存储在float格式中，所以将输入图像转换成浮点类型，并多加一个额外通道存储复数部分
	//将planes数组组合合并成一个多通道的数组complexI
	Mat planes[] = { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);

	//[4]进行就地离散傅里叶变换
	dft(complexI, complexI);

	//[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];

	//[6]进行对数尺寸（logarithmic scale）缩放
	//傅里叶变换的幅度值范围大到不合适再屏幕显示，高值在屏幕上显示为白点，低值为黑点，高低值的变化无法有效分辨，为了在屏幕上凸显高低变化的连续性，可以用对数尺度来替换线性尺度：M'=log(1+M）
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);//求自然对数

	//[7]剪切和重分布幅度图象限
	//因为[2]中进行了图像延扩，现在需要剔除，为了方便显示，可以重新分布幅度图象限位置
	//若有奇数行或奇数列，进行频谱裁剪
	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);

	//[8]归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
	//有了重分布后的幅度图，但是幅度值仍然超过可显示范围[0，1]，使用normalize函数将幅度归一化到可显示范围
	normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);

	//[9]显示效果图
	imshow("频谱幅值", magnitudeImage);

	waitKey(0);
	return 0;
}