离散傅立叶变换与逆变换

一、怎样为一副图像增加一个通道

Mat A=(Mat_<double>(3,3)<<1,2,3,4,5,6,7,8,9);
Mat B=Mat::zeros(A.size(),A.type());
cout<<"the channels numbers of A: "<<A.channels()<<endl;
Mat  planels[]={A,B};
Mat out;
merge(planels,2,out);
cout<<"After changing the channels numbers: "<<out.channels()<<endl;
cout<<out<<endl;
return 0;
结果：
the channels numbers of A: 1
After changing the channels numbers: 2
[1, 0, 2, 0, 3, 0;
4, 0, 5, 0, 6, 0;
7, 0, 8, 0, 9, 0]
请按任意键继续. . .

从上面结果可以看出，我们为原图像增加了一个通道

二、离散傅立叶变换

（1）将图像延扩到最佳尺寸. 离散傅立叶变换的运行速度与图片的尺寸息息相关。当图像的尺寸是2， 3，5的整数倍时，计算速度最快。 因此，为了达到快速计算的目的，经常通过添凑新的边缘像素的方法获取最佳图像尺寸。函数 getOptimalDFTSize() 返回最佳尺寸，而函数 copyMakeBorder() 填充边缘像素:

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

Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);         // 为延扩后的图像增添一个初始化为0的通道

（3）进行离散傅立叶变换. 支持图像原地计算 (输入输出为同一图像):

dft(complexI, complexI);            // 变换结果很好的保存在原始矩阵中

（4）将复数转换为幅度.复数包含实数部分(Re)和复数部分 (imaginary - Im)。 离散傅立叶变换的结果是复数，对应的幅度可以表示为:

M = sqrt{ {Re(DFT(I))}^2 + {Im(DFT(I))}^2}
转化为OpenCV代码:
split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
Mat magI = planes[0];

（5）对数尺度(logarithmic scale)缩放. 傅立叶变换的幅度值范围大到不适合在屏幕上显示。高值在屏幕上显示为白点，而低值为黑点，高低值的变化无法有效分辨。为了在屏幕上凸显出高低变化的连续性，我们可以用对数尺度来替换线性尺度:

M= log{(1 + M)}
转化为OpenCV代码:
magI += Scalar::all(1);                    // 转换到对数尺度
log(magI, magI);

（6）剪切和重分布幅度图象限. 还记得我们在第一步时延扩了图像吗? 那现在是时候将新添加的像素剔除了。为了方便显示，我们也可以重新分布幅度图象限位置(注：将第五步得到的幅度图从中间划开得到四张1/4子图像，将每张子图像看成幅度图的一个象限，重新分布即将四个角点重叠到图片中心)。 这样的话原点(0,0)就位移到图像中心。

magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
int cx = magI.cols/2;
int cy = magI.rows/2;

Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - 为每一个象限创建ROI
Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

Mat tmp;                           // 交换象限 (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);

q1.copyTo(tmp);                    // 交换象限 (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);

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

normalize(magI, magI, 0, 1, CV_MINMAX); // 将float类型的矩阵转换到可显示图像范围(float [0， 1]).

（8）具体实现代码

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

int main(){
Mat I = imread("F:\\opencv_test\\1.jpg" ,0);
if(I.empty())
return -1;
//得到最优的图像大小
int m=getOptimalDFTSize(I.rows);
int n=getOptimalDFTSize(I.cols);
Mat paddImg;
//对图像补零，在图像下方和图像右边不零，上边和左边保持不变
copyMakeBorder(I,paddImg,0,m-I.rows,0,n-I.cols,BORDER_CONSTANT,Scalar::all(0));
//扩充一个通道，方便存储傅立叶变换后的结果
Mat planels[]={Mat_<float>(paddImg),Mat::zeros(paddImg.size(),CV_32F)};//Mat_的用法，为什么paddImg一定要打括号？
Mat complexImg;
merge(planels,2,complexImg);
//进行傅立叶变换
dft(complexImg,complexImg);
//把变换后的实数和复数部分分离开
split(complexImg,planels);
//求幅度
magnitude(planels[0],planels[1],planels[0]);
Mat magI=planels[0];
//压缩动态范围，做对数变换
magI=magI+Scalar::all(0);
log(magI,magI);
//去掉补充的行列
magI = magI(Rect(0, 0, I.cols, I.rows ));//
//平移，使图像的中心变换到原点
int cy=magI.rows/2;
int cx=magI.cols/2;

//把原区域分为四块
Mat q0(magI,Rect(0,0,cx,cy));//Rect的用法，第三个参数到底是列
Mat q1(magI,Rect(cx,0,cx,cy));
Mat q2(magI,Rect(0,cy,cx,cy));//
Mat q3(magI,Rect(cx,cy,cx,cy));

//把左上角和右下角的区域交换
Mat temp;
q0.copyTo(temp);
q3.copyTo(q0);
temp.copyTo(q3);
//把左下角和右上角的区域进行互换
q1.copyTo(temp);
q2.copyTo(q1);
temp.copyTo(q2);

//归一化，方便显示
normalize(magI,magI,1,0,CV_MINMAX);
namedWindow("original image",0);
namedWindow("spectrum magnitude image",0);
imshow("original image",I);
imshow("spectrum magnitude image",magI);

waitKey();
return 0;

}

三、傅立叶逆变换

傅立叶逆变换实现韩式为idft,逆变换后输出的为复数，我们可取其实部显示即可

int main(){
Mat I = imread("F:\\opencv_test\\6.tif",0);
if(I.empty())
return -1;
//傅立叶变换
int m=getOptimalDFTSize(I.rows);
int n=getOptimalDFTSize(I.cols);
Mat paddImg;
copyMakeBorder(I,paddImg,0,m-I.rows,0,n-I.cols,BORDER_CONSTANT,Scalar::all(0));

Mat planels[]={Mat_<float>(paddImg),Mat::zeros(paddImg.size(),CV_32F)};
Mat complexImg;
merge(planels,2,complexImg);
dft(complexImg,complexImg);
//傅立叶逆变换
Mat iImg;
idft(complexImg,iImg);
//实际上，做完逆变换是复数形式，我们取其实部即可
split(iImg,planels);
Mat magI=planels[0];
magI = magI(Rect(0, 0, I.cols, I.rows ));
//归一化，方便显示
normalize(magI,magI,1,0,CV_MINMAX);

imshow("original",I);

namedWindow("spectrum magnitude image",0);

imshow("spectrum magnitude image",magI);

waitKey();
return 0;

}

注：傅立叶变换来自：F:\opencv\OpenCV2.3.2官方文档（中文）\core 模块. 核心功能\离散傅立叶变换 — OpenCV 2.3.2 documentation.htm