离散傅里叶变换
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;
void ShowHelpText()
{
//输出欢迎信息和OpenCV版本
printf("\n\n\t\t\t非常感谢购买《OpenCV3编程入门》一书!\n");
printf("\n\n\t\t\t此为本书OpenCV2版的第28个配套示例程序\n");
printf("\n\n\t\t\t 当前使用的OpenCV版本为:" CV_VERSION );
printf("\n\n ----------------------------------------------------------------------------\n");
}
int main( )
{
//【1】以灰度模式读取原始图像并显示
Mat srcImage = imread("1.jpg", 0);
if(!srcImage.data ) { printf("读取图片错误,请确定目录下是否有imread函数指定图片存在~! \n"); return false; }
imshow("原始图像" , srcImage);
ShowHelpText();
//【2】将输入图像延扩到最佳的尺寸,边界用0补充
int m = getOptimalDFTSize( srcImage.rows );
int n = getOptimalDFTSize( srcImage.cols );
//该函数返回给定向量尺寸的傅里叶最优尺寸大小。为了提高离散傅里叶变换的运行速度,需要扩充图像,而具体扩充多少,就由这个函数来计算得到。
//DFT变换在一个向量尺寸上不是一个单调函数,当计算两个数组卷积或对一个数组进行光学分析,它常常会用0扩充一些数组来得到稍微大点的数组以达到比原来数组计算更快的目的。一个尺寸是2阶指数(2,4,8,16,32…)的数组计算速度最快,一个数组尺寸是2、3、5的倍数(例如:300 = 5*5*3*2*2)同样有很高的处理效率。
//getOptimalDFTSize()函数返回大于或等于vecsize的最小数值N,这样尺寸为N的向量进行DFT变换能得到更高的处理效率。在当前N通过p,q,r等一些整数得出N = 2^p*3^q*5^r.
//将添加的像素初始化为0.
Mat padded;
copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));
//该函数作用是扩充图像边界,第一个参数,输入图像;第二个参数,输出图像(与源图像有一样的类型,且size为Size(src.cols+left+right,src.rows+top+bottom);接下来
//四个参数分别为int类型的top,bottom,left,right分别表示在源图像的四个方向扩充多少像素;第七个参数,borderType类型的,边界类型,常见取值为BORDER_CONSTANT Various border types, image boundaries are denoted with '|'
/* BORDER_REPLICATE: aaaaaa|abcdefgh|hhhhhhh
* BORDER_REFLECT: fedcba|abcdefgh|hgfedcb
* BORDER_REFLECT_101: gfedcb|abcdefgh|gfedcba
* BORDER_WRAP: cdefgh|abcdefgh|abcdefg
* BORDER_CONSTANT: iiiiii|abcdefgh|iiiiiii with some specified 'i';具体可见网页收藏*/
//第八个参数const Scalar&类型的value,由默认值Scalar(),可以理解为默认值为0.当borderType取值为BORDER_CONSTANT时,这个参数表示边界值。
//【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
//将planes数组组合合并成一个多通道的数组complexI
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); //void merge(const Mat* mv, size_t count, OutputArray dst);融合有两种原型
//傅里叶变换结果是复数,这就是说对于每个原图像值,结果会有两个图像值。此外,频域值范围远远超过空间值范围,因此至少要将频域存储在float格式中,所以我们要将输入图像转换成
//浮点类型,并多加一个额外通道来储存复数部分。
//【4】进行就地离散傅里叶变换
dft(complexI, complexI);
/*.void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0);
InputArray src: 输入图像,可以是实数或虚数
OutputArray dst: 输出图像,其大小和类型取决于第三个参数flags
int flags = 0: 转换的标识符,有默认值0.其可取的值如下所示:
。DFT_INVERSE: 用一维或二维逆变换取代默认的正向变换
。DFT_SCALE: 缩放比例标识符,根据数据元素个数平均求出其缩放结果,如有N个元素,则输出结果以1/N缩放输出,常与DFT_INVERSE搭配使用。
。DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换;此标识符可在处理多种适量的的时候用于减小资源的开销,这些处理常常是三维或高维变换等复杂操作。
。DFT_COMPLEX_OUTPUT: 对一维或二维的实数数组进行正向变换,这样的结果虽然是复数阵列,但拥有复数的共轭对称性(CCS),可以以一个和原数组尺寸大小相同的实数数组进行填充,这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组(像简单光谱分析等等),通过设置标志位可以使函数生成一个全尺寸的复数输出数组。
。DFT_REAL_OUTPUT: 对一维二维复数数组进行逆向变换,这样的结果通常是一个尺寸相同的复数矩阵,但是如果输入矩阵有复数的共轭对称性(比如是一个带有DFT_COMPLEX_OUTPUT标识符的正变换结果),便会输出实数矩阵。
. int nonzeroRows = 0: 当这个参数不为0,函数会假设只有输入数组(没有设置DFT_INVERSE)的第一行或第一个输出数组(设置了DFT_INVERSE)包含非零值。这样的话函数就可以对其他的行进行更高效的处理节省一些时间,这项技术尤其是在采用DFT计算矩阵卷积时非常有效。*/
//【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
//magnitude()函数用于计算二维矢量的幅值。第一个参数,InputArray类型的x,表示矢量的浮点型X坐标值,也就是实部,第二个参数,InputArray类型的y,
//表示矢量的浮点型Y坐标值,也就是虚部,第三个参数,OutputArray类型的magnitude,输出的幅值,它和第一个参数x有着同样的尺寸和类型。dst(I)=根号(x(I)^2+y(I)^2)
Mat magnitudeImage = planes[0];
//【6】进行对数尺度(logarithmic scale)缩放;对数变换,原因是幅值太大,必须缩小到合适显示的区间,加一防止出现负数
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage);//求自然对数e为底dst(I)=log|src(I)|ifsrc(I)!=0;=C otherwise
//【7】剪切和重分布幅度图象限
//若有奇数行或奇数列,进行频谱裁剪
magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
//magnitudeImage.cols & -2 是要进行位的和运算。magnitudeImage.cols和-2需要分别转化为无符整型,再转化为二进制数进行和运算。由于-2是负数,再转换的过程中存在溢出等问题,
//以8位二进制为例子,-2(补码)实际转化为二进制的数为11111110。当一个二进制数与11111110进行位的和运算时,该二进制数的最小数位就为0,转化为无符整型时,就是一个偶数。
//重新排列傅立叶图像中的象限,使得原点位于图像中心
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(magnitudeImage, magnitudeImage,0,1, CV_MINMAX);
//矩阵归一化:第一个参数输入对象;第二个参数输出对象;第三个参数double类型的alpha,归一化后的最大值,由默认值1;第四个参数double类型的beta,归一化后的最大值有默认值0;
//底五个参数int类型的norm_type。归一化类型,有NORM_INF、NORM_L1和NORM_L2和NORM_MINMAX等参数可选,有默认值NORM_L2;第六个参数,int类型的dtype,有默认值-1.当此参数取负值时,
//输出矩阵和src有同样的类型,否则,它和src有同样的通道数,且此时图像深度为CV_MAT_DEPTH(dtype);第七个参数,InputArray类型的mask,可选的操作掩膜
//之所以归一化到0到1因为该图形是32位浮点型
//【9】显示效果图
imshow("频谱幅值", magnitudeImage); //在imshow中在把【0,1】映射到【0,255】
waitKey();
return 0;
}