图像处理之傅里叶变换和小波变换

最近在看物体识别论文摘要，好多论文中涉及到使用离散余弦傅里叶变换DFT（Discrete Fourier Transform）对图像进行处理，因此特地看了这部分的内容，傅里叶变换和小波变换。

一、DFT的原理：

以二维图像为例，归一化的二维离散傅里叶变换可以写成如下形式：

其中f(x,y)表示图像的空间域的值，而F表示频域的值，傅里叶转换的结果为复数，这也表明，傅里叶变换其实是一副实数图像和虚数图像叠加或幅度图像和相位图像叠加的结果，在实际的图像处理算法中，仅有幅度图像能够用得到，因为其包含了图像的所有几何结构信息。但是如果想通过修改幅度图和相位图来修改原空间图像，需保留幅度图和相位图来进行傅里叶变换，从而得到修改后的图像。

在频域里面，高频部分代表了图像边缘、线条以及纹理等细节信息，低频部分代表了图像的轮廓信息。在这里首先介绍下空间域和频率域：

空间域：

一般情况下，空间域的图像为f(x,y),形象一点就是一个二维矩阵，每个坐标对应一个颜色值。

频率域：

频率：对于图像来说，可以指图像颜色值的梯度，即灰度级的变化速度。

幅度：频率的权，即该频率所占的比例。

二、代码实现与效果

#include<iostream>
#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
using namespace std;
using namespace cv;
int main()
{
	//(1)读取原图像
	Mat src = imread("101200.jpg",0);
	if (!src.data)
	{
		cout << "Reading image error!" << endl;
		return false;
	}
	imshow("src",src);
	//（2）将输入图像扩展到最佳尺寸，将添加的像素扩展为0
	int m = getOptimalDFTSize(src.rows);
	int n = getOptimalDFTSize(src.cols);
	Mat padded;
	copyMakeBorder(src,padded,0,m-src.rows,0,n-src.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)将复数转换为幅值，即=> 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];
	//(5)进行对数尺度(logarithmic scale)缩放
	//由于幅度范围太大，不适合在屏幕显示。为了在屏幕显示，用对数尺度来替换线性尺度M1=log(1+M)
	magnitudeImage += Scalar::all(1);//转换到对数尺度
	log(magnitudeImage, magnitudeImage);//求自然对数
	//(6)剪切和重分布幅度图象限
	/*
	剔除第二步添加的像素。重分布是把四个象限的四张图像拼接到一起
	*/
	//若有奇数行或奇数列，进行频谱裁剪      
	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);
	//(7)归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
	/*
	幅度值仍然超过可显示范围[0,1],normalize()归一化后可以显示
	*/
	normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX); 
	//【9】显示效果图
	imshow("dft", magnitudeImage);
	//逆变换
	Mat _complexim;
	complexI.copyTo(_complexim);//把变换结果复制一份，进行逆变换，也就是恢复原图
	Mat iDft[] = { Mat::zeros(planes[0].size(), CV_32F), Mat::zeros(planes[0].size(), CV_32F) };//创建两个通道，类型为float，大小为填充后的尺寸
	idft(_complexim, _complexim);//傅立叶逆变换
	split(_complexim, iDft);//结果貌似也是复数
	magnitude(iDft[0], iDft[1], iDft[0]);//分离通道，主要获取0通道
	normalize(iDft[0], iDft[0], 1, 0, CV_MINMAX);//归一化处理，float类型的显示范围为0-1,大于1为白色，小于0为黑色
	imshow("idft",iDft[0]);
	waitKey(0);
	return 0;
}

实现结果：

二、小波变换

# include<opencv2/opencv.hpp>
# include<iostream>

using namespace std;
using namespace cv;

int main()
{
	Mat img = imread("101200.jpg", 0);
	int Height = img.cols;
	int Width = img.rows;
	int depth = 1;    //定义分解深度,也就是几级分解
	int depthcount = 1;
	Mat tmp = Mat::ones(Width, Height, CV_32FC1);//CV_32FC1表示32位float，这在32位编译器上是32位float，也就是单精度。CV_64FC1在32位编译器上是64位float，也就是双精度。
	Mat wavelet = Mat::ones(Width, Height, CV_32FC1);
	Mat imgtmp = img.clone();
	imgtmp.convertTo(imgtmp, CV_32FC1);
	while (depthcount <= depth){
		Width = img.rows / depthcount;
		Height = img.cols / depthcount;
		for (int i = 0; i < Width; i++)
		{
			for (int j = 0; j < Height / 2; j++)
			{
				tmp.at<float>(i, j) = (imgtmp.at<float>(i, 2 * j) + imgtmp.at<float>(i, 2 * j + 1)) / 2;//整体信息
				tmp.at<float>(i, j + Height / 2) = (imgtmp.at<float>(i, 2 * j) - imgtmp.at<float>(i, 2 * j + 1)) / 2;//细节信息
			}
		}
		for (int i = 0; i < Width / 2; i++)
		{
			for (int j = 0; j < Height; j++)
			{
				wavelet.at<float>(i, j) = (tmp.at<float>(2 * i, j) + tmp.at<float>(2 * i + 1, j)) / 2;//整体信息
				wavelet.at<float>(i + Width / 2, j) = (tmp.at<float>(2 * i, j) - tmp.at<float>(2 * i + 1, j)) / 2;//细节信息
			}
		}
		imgtmp = wavelet;
		depthcount++;
	}
	wavelet.convertTo(wavelet, CV_8UC1);
	wavelet += 100;            //图像暗度过低，所以这里我加了50
	imshow("小波变换", wavelet);
	waitKey(0);
	return 0;
}

实现效果：

.