图像的二维离散傅里叶变换代码解读以及一些效果展示

二维离散傅里叶变换后的图像的意义：

低频部分反映图像的概貌，高频反映图像的细节，改变换建立了图像在空间域的像素特性与其频率域上的信息强度之间的联系。
——《数字图像处理》贾永红

来自博客 如何理解 图像傅里叶变换的频谱图：

很多人都不了解图像（二维）频谱中的每一点究竟代表了什么，有什么意义?

一句话解释为： 二维频谱中的每一个点都是一个与之一 一对应的二维正弦/余弦波。

代码效果展示

几个图像的频谱图像

代码

我给书上《opencv3编程入门》代码增加了更详细的注释

/*
《opencv3编程入门》p139.离散傅里叶变换
以输入图像为单通道的灰度图像I为例。
    ☆相关概念：
    通道的概念（参考 http://www.xuebuyuan.com/1681143.html） ：
    图像通道在RGB色彩模式下就是指在下就是指那单独的红色R、绿色G、蓝色B部分。
也就是说，一幅完整的图像，是由红色绿色蓝色三个通道组成的。他们共同作用产生了完整的图像。
同样在HSV色系中指的是色调H，饱和度S，亮度V三个通道。
多通道模式是把含有通道的图像分割成单个的通道。
    灰度模式（参考 http://www.xuebuyuan.com/1681143.html） ：
    灰度模式是8位深度的图像模式。也就是28，28=256，在全黑和全白之间插有254个灰度等级的颜色来描绘灰度模式的图像。
所有模式的图像都能换成灰度模式，甚至位图（深度为1，只有黑白两色）也可转换为灰度模式。
*/
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;

int main( )
{

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

	//【2】将输入图像延扩到最佳的尺寸，边界用0补充
	/*
	    傅里叶变换的速度和图像尺寸有关，当尺寸是2、3、5的整数倍，计算速度较快。
	于是进行添凑新的边缘像素。getOptimalDFTSize()用于返回最佳尺寸，copyMakeBorder()用于填充边缘像素
	*/
	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()使几个单通道数组合并成多通道数组
	merge(planes, 2, complexI);         

	//【4】进行就地离散傅里叶变换
	/*
		就地（in-place)含义:
		输入输出为同一图像
	*/
	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)缩放
	/*
		由于幅度范围太大，不适合在屏幕显示。为了在屏幕显示，
	用对数尺度来替换线性尺度
	M1=log(1+M)
	*/
	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】归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
	/*
		幅度值仍然超过可显示范围[0,1],normalize()归一化后可以显示
	*/
	//此句代码的OpenCV2版为：
	normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX); 
	//此句代码的OpenCV3版为:
	//normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX); 

	//【9】显示效果图
	imshow("频谱幅值", magnitudeImage);
	//任意键按下，程序关闭
	waitKey();

	return 0;
}