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

     在此示例中，将展示如何计算以及显示傅里叶变换后的幅度图像。由于数字图像的离散性，像素值的取值范围也是有限的。比如在一张灰度图像中，像素灰度值一般在0到255之间。因此,我们这里讨论的也仅仅是离散傅里叶变换(DFT)。如果需要得到图像中的几何结构信息，那么就要用到离散傅里叶变换了。下面的步骤将以输入图像为单通道的灰度图像Ⅰ为例，进行分步说明。
1.【第一步】载入原始图像
    我们在这一步以灰度模式读取原始图像，进行是否读取成功的检测，并显示出读取到的图像。代码如下。

//[1]以灰度模式读取原始图像并显示
	Mat srcImage = imread("E:/pictures/2.jpg", 0);
	imshow("原始图像", srcImage);

2.【第二步】将图像扩大到合适的尺寸
     离散傅里叶变换的运行速度与图片的尺寸有很大关系。当图像的尺寸是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));

3.【第三步】为傅里叶变换的结果（实部和虚部）分配存储空间
    傅里叶变换的结果是复数,这就是说对于每个原图像值,结果会有两个图像值。此外，频域值范围远远超过空间值范围，因此至少要将频域储存在 float格式中。所以我们将输入图像转换成浮点类型，并多加一个额外通道来储存复数部分。

	//[3]为傅里叶变换的结果（实部和虚部）分配存储空间
	//将planes数组组合合并成一个多通道的数组complexI
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);

4.【第四步】进行离散傅里叶变换
这里的离散傅里叶变换为图像就地计算模式（ in-place，输入输出为同一图像)。

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

5.【第五步】将复数转换为幅值
复数包含实数部分（Re）和虚数部分(imaginary - Im)。离散傅里叶变换的结果是复数，对应的幅度可以表示为:
M= √Re(DET(I))3+ Im(DFT(T))'
那么，转化为OpenCV代码，就是下面这样的。

	//[5]将复数转换为幅值
	split(complexI, planes);//将多通道数组complexI分离成几个单通道数组
	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage = planes[0];

6.【第六步】进行对数尺度（( logarithmic scale）缩放.
傅里叶变换的幅度值范围大到不适合在屏幕上显示。高值在屏幕上显示为白点，而低值为黑点，高低值的变化无法有效分辨。为了在屏幕上凸显出高低变化的连续性，我们可以用对数尺度来替换线性尺度，公式如下。
M=log(1+M)
而写成OpenCV代码，就是下面这样的。

	//[6]进行对数尺度缩放
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);

7.【第七步】剪切和重分布幅度图象限
因为在第二步中延扩了图像，那现在是时候将新添加的像素剔除了。为了方便显示，也可以重新分布幅度图象限位置（注:将第五步得到的幅度图从中间划开，得到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);

8.【第八步】归一化
这一步仍然是为了显示。现在有了重分布后的幅度图，但是幅度值仍然超过可显示范围[0,1]。我们使用normalize()函数将幅度归一化到可显示范围。

//[8]归一化
	normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);

9.【第九步】显示效果图
处理完成，最后进行显示即可。

源代码如下：
 

#include "opencv2\core\core.hpp"
#include "opencv2\imgproc\imgproc.hpp"
#include "opencv2\highgui\highgui.hpp"
#include <opencv2/core/utils/logger.hpp>
#include <iostream>

using namespace std;
using namespace cv;

int main()
{
	cv::utils::logging::setLogLevel(utils::logging::LOG_LEVEL_SILENT);//控制台不在输出日志文件
	//[1]以灰度模式读取原始图像并显示
	Mat srcImage = imread("E:/pictures/2.jpg", 0);
	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]为傅里叶变换的结果（实部和虚部）分配存储空间
	//将planes数组组合合并成一个多通道的数组complexI
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);
	//[4]进行就地离散傅里叶变换
	dft(complexI, complexI);
	//[5]将复数转换为幅值
	split(complexI, planes);//将多通道数组complexI分离成几个单通道数组
	magnitude(planes[0], planes[1], planes[0]);
	Mat magnitudeImage = planes[0];
	//[6]进行对数尺度缩放
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);
	//[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);
	//[8]归一化
	normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);
	//[9]显示效果图
	imshow("频谱幅值", magnitudeImage);
	waitKey();
	return 0;
}

运行结果如图：