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

第5章 Core组件进阶

5.5 离散傅里叶变换&＃xff08;Discrete Fourier Transform&＃xff0c;DFT&＃xff09;

5.5.1 离散傅里叶变换原理

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

5.5.2 dft()函数

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

void dft(InputArray src,OutputArray dst,int flags&＃61;0,int nonzeroRows&＃61;0)

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

4.示例&＃xff1a;用dft函数计算两个二维实矩阵卷积

void consolveDFT(InputArray A, InputArray B, OutputArray C)
{//【1】初始化输出矩阵C.create(abs(A.rows - B.rows) &＃43; 1, abs(A.cols - B.cols) &＃43; 1, A.type());Size dftSize;//【2】计算DFT变换尺寸dftSize.width &＃61; getOptimalDFTSize(A.cols &＃43; B.cols - 1);dftSize.height &＃61; getOptimalDFTSize(A.rows &＃43; 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】就地操作&＃xff08;in_place&＃xff09;&＃xff0c;进行快速傅里叶变换&＃xff0c;并将nonzeroRows参数置为非零&＃xff0c;以进行更快速的处理dft(tempA, tempA, 0, A.rows);dft(tempB, tempB, 0, B.rows);//【6】将得到的频谱相乘&＃xff0c;结果存放于tempA中mulSpectrums(tempA, tempB, tempA);//计算两个傅里叶频谱的每个元素的乘法//【7】将结果变换为频域&＃xff0c;且尽管结果行&＃xff08;result rows&＃xff09;都为非零&＃xff0c;我们只需要其中的C.rows的第一行&＃xff0c;所以采用nonzeroRows&＃61;&＃61;C.rowsdft(tempA, tempA, DFT_INVERSE &＃43; DFT_SCALE, C.rows);//【8】将结果复制到C中tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);//所有的临时缓冲区将被自动释放&＃xff0c;所以无须收尾操作
}


5.5.3 返回DFT最优尺寸大小&＃xff1a;getOptimalDFTSize()函数

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

5.5.4 扩充图像边界&＃xff1a;copyMakeBorder()函数

1.原型&＃xff1a;

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

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

5.5.5 计算二维矢量的幅值&＃xff1a;magnitude()函数

1.原型&＃xff1a;

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

2.参数说明&＃xff1a;
&＃xff08;1&＃xff09;表示矢量的浮点型X坐标值&＃xff08;实部&＃xff09;
&＃xff08;2&＃xff09;表示矢量的浮点型Y坐标值&＃xff08;虚部&＃xff09;
&＃xff08;3&＃xff09;输出的幅值&＃xff0c;它和第一个参数x有着同样的尺寸和类型
3.原理&＃xff1a;
                

5.5.6 计算自然对数&＃xff1a;log()函数

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

void log(InputArray src,OutputArray dst)

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

1.原型&＃xff1a;

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

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

5.5.8 示例程序&＃xff1a;离散傅里叶变换

#include
#include
#include
#include
#include
using namespace cv;
using namespace std;int main()
{//[1]以灰度模式读取原始图像并显示Mat srcImage &＃61; imread("love.jpg", 0);if (!srcImage.data) { printf("读取图片错误&＃xff0c;请确定目录下是否有imread函数指定图片存在~&＃xff01;\n"); return false; }imshow("原始图像", srcImage);//[2]将输入图像延扩到最佳尺寸&＃xff0c;边界用0扩充(扩大图像尺寸可以提高计算速度)int m &＃61; getOptimalDFTSize(srcImage.rows);int n &＃61; getOptimalDFTSize(srcImage.cols);//将添加的像素初始化为0Mat padded;copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));//[3]为傅里叶变换的结果&＃xff08;实部和虚部&＃xff09;分配存储空间//傅里叶变换结果是复数&＃xff0c;即每个原图像值结果会有两个图像值,且频域值范围至少存储在float格式中&＃xff0c;所以将输入图像转换成浮点类型&＃xff0c;并多加一个额外通道存储复数部分//将planes数组组合合并成一个多通道的数组complexIMat planes[] &＃61; { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };Mat complexI;merge(planes, 2, complexI);//[4]进行就地离散傅里叶变换dft(complexI, complexI);//[5]将复数转换为幅值&＃xff0c;即&＃61;>log(1&＃43;sqrt(Re(DFT(I)^2&＃43;IM(DFT(I))^2))split(complexI, planes); //将多通道数组complexI分离成几个单通道数组,planes[0]&＃61;Re(DFT(I),planes[1]&＃61;Im(DFT(I))magnitude(planes[0], planes[1], planes[0]);//planes[0]&＃61;magnitudeMat magnitudeImage &＃61; planes[0];//[6]进行对数尺寸&＃xff08;logarithmic scale&＃xff09;缩放//傅里叶变换的幅度值范围大到不合适再屏幕显示&＃xff0c;高值在屏幕上显示为白点&＃xff0c;低值为黑点&＃xff0c;高低值的变化无法有效分辨&＃xff0c;为了在屏幕上凸显高低变化的连续性&＃xff0c;可以用对数尺度来替换线性尺度&＃xff1a;M&＃39;&＃61;log(1&＃43;M&＃xff09;magnitudeImage &＃43;&＃61; Scalar::all(1);log(magnitudeImage, magnitudeImage);//求自然对数//[7]剪切和重分布幅度图象限//因为[2]中进行了图像延扩&＃xff0c;现在需要剔除&＃xff0c;为了方便显示&＃xff0c;可以重新分布幅度图象限位置//若有奇数行或奇数列&＃xff0c;进行频谱裁剪magnitudeImage &＃61; magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));//重新排列傅里叶图像中的象限&＃xff0c;使得原点位于图像中心int cx &＃61; magnitudeImage.cols / 2;int cy &＃61; 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区域的右下//交换象限&＃xff08;左上和右下进行交换&＃xff09;Mat tmp;q0.copyTo(tmp);q3.copyTo(q0);tmp.copyTo(q3);//交换象限&＃xff08;右上与左下进行交换&＃xff09;q1.copyTo(tmp);q2.copyTo(q1);tmp.copyTo(q2);//[8]归一化&＃xff0c;用0到1之间的浮点值将矩阵变换为可视的图像格式//有了重分布后的幅度图&＃xff0c;但是幅度值仍然超过可显示范围[0&＃xff0c;1]&＃xff0c;使用normalize函数将幅度归一化到可显示范围normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//[9]显示效果图imshow("频谱幅值", magnitudeImage);waitKey(0);return 0;
}