DFT原理:(单变量离散傅里叶变换)
数学基础:
任何一个函数都可以转换成无数个正弦和余弦函数的和的形式。
通常观察傅里叶变换后的频域函数可以获得两个重要的信息:幅频曲线和相频曲线。
在数字图像处理中的作用:
在数字图像处理中,对一张图片进行傅里叶变换后我们获得的是:实数图像(幅度图像)+虚数图像(相位图像)
傅里叶变换在数字图像处理中将空间域信息转为频域信息。
如果需要得到图像中的几何结构信息,就要用到离散傅里叶变换
在频域中高频和低频意义:
高频代表了图像的细节、纹理信息,(噪声)
低频代表图像的轮廓信息。
所以低通滤波器可以得到轮廓。
傅里叶变换可以应用的场景:
图像的增强与图像去噪
图像分割的边缘检测
图像特征提取
图像压缩
Opencv中离散傅里叶变换函数:
dft(input,output,flags,nonzeroRows);
//flags:标识符
//nonzeroRows:默认值为0,非零时,表示你想处理的那一行C.rows,计算时更加高效
dft函数应用实例
例1:用dft函数计算两个二维实矩阵卷积
例子包含的小知识点
1、Size类型
CvSize结构表示矩形尺寸的结构,结构体中分别定义了矩形的宽度和高度,具体定义如下:
typedef struct CvSize {
int width; /* 矩形宽度,单位为象素 /
int height; / 矩形高度,单位为象素 */
}CvSize;
2、getOptimalDFTSize()函数
返回DFT最优尺寸大小:getOpimalDFTSize()
函数返回给定向量尺寸的傅里叶最优尺寸大小、
input:向量尺寸,即图像的rows\cols
3、mulSpectrums()函数
void cvMulSpectrums( const CvArr* src1, const CvArr* src2, CvArr* dst, int flags );
src1
第一输入数组
src2
第二输入数组
dst
输出数组,和输入数组有相同的类型和大小。
flags
下面列举的值的组合:
DFT_COMPLEX_OUTPUT- 把数组的每一行视为一个单独的频谱 (参见 cvDFT 的参数讨论).
DFT_REAL_OUTPUT - 在做乘法之前取第二个输入数组的共轭.
函数 cvMulSpectrums 执行两个 CCS-packed 或者实数或复数傅立叶变换的结果复数矩阵的每个元素的乘法。
4、dft只能处理浮点数,所以需要将输入图像转为float类型
全部代码:
void convolveDFT(Mat& A,Mat& B, Mat& C)
{//【1】初始化输出矩阵C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type());Size dftSize; //???什么意思//【2】计算DFT变换的尺寸dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);dftSize.height = getOptimalDFTSize(A.rows + 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】就地操作,进行快速傅里叶变换,并将nonzeroRows参数置为非零,以进行更加快速的处理???为什么dft(tempA, tempA, 0, A.rows);dft(tempB, tempB, 0, B.rows);//【6】将得到的频谱相乘,结果存放于tempA中mulSpectrums(tempA,tempB,tempA, DFT_COMPLEX_OUTPUT);//DFT_REAL_OUTPUT//【7】将结果变换为频域且尽管行结果都为非零,我们只需要其中C.rows的第一行,所以采用nonzeroRows==C.rowsdft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows);//【8】将结果复制到C中tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);//所有的临时缓冲区将被自动释放,所以无须收尾操作
}
int main()
{SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), FOREGROUND_INTENSITY | FOREGROUND_GREEN); //字体为绿色//载入原图Mat srcImage = imread("D:\\opencv_picture_test\\形态学操作\\coin_inv.png", 0); //读取灰度图Mat kernel = (Mat_<float>(3, 3) << 1, 1, 1, 1, 1, 1, 1, 1, 1);cout << kernel;Mat floatI = Mat_<float>(srcImage);// change image type into floatMat filteredI;convolveDFT(floatI, kernel, filteredI);normalize(filteredI, filteredI, 0, 1,NORM_MINMAX); // Transform the matrix with float values into a// viewable image form (float between values 0 and 1).//imshow("image", srcImage);imshow("filtered", filteredI); //这里显示报错,但是可以用ImageWatch查看,暂时不知道原因waitKey(0);return 0;
}
参考链接:
https://blog.csdn.net/lichengyu/article/details/18848281
利用傅里叶变换卷积和利用核游走整个图像进行卷积运算的区别
一般求法中,利用核游走整个图像进行卷积运算,实际上进行的是相关运算,真正意义上的卷积,应该首先把核翻转180度,再在整个图像上进行游走。OpenCV中的filter2D实际上做的也只是相关,而非卷积。
例2:显示一幅图像傅里叶变换后的频域图像
需要用的函数:
1、扩充图像边界:copyMakeBorder()
C++: void copyMakeBorder (InputArray src, OutputArray dst, int top, int bottom, int left, int right,int borderType, const Scalar&
value=Scalar () );
top
bottom
left
right
分别表示在四个方向上扩充多少个像素。
第七个参数:边界类型,常见取值BORDER_CONSTANT;
第八个参数:默认值为0,当边界类型取值为BORDER_CONSTANT时,这个参数表示边界值。
2、计算二维矢量的幅值:magnitude()
C++: void magnitude (InputArray X,InputArray Y, OutputArray magnitude)
X:矢量实部
Y:矢量虚部
magnitude:输出幅值,与x有相同的尺寸和类型
3、log()函数
C++: void 1og (InputArray srC, OutputArray dst)
计算每个数组元素绝对值的自然对数。
傅里叶变换的幅值范围大到不适合在屏幕上显示。
为了凸显出高低变化的连续性,可以用对数尺度来替换线性尺度。
3、矩阵归一化:normalize()函数
void cv::normalize(InputArry src,InputOutputArray dst,double alpha=1,double beta=0,int norm_type=NORM_L2,int dtype=-1,InputArray mark=noArry())
alpha 用来规范值或者规范范围,并且是下限;
beta 只用来规范范围并且是上限,因此只在NORM_MINMAX中起作用;
1.NORM_L1、NORM_INF、NORM_L2模式下归一化结果与beta无关,只与alpha有关,详见第4部分的公式说明;
2.NORM_MINMAX中alpha、beta都起作用,同时需要注意的是alpha和beta的取值顺序与归一化结果无关。即alpha=255,beta=0和alpha=0,beta=255最后的归一化结果是相同的。
归一化公式:
而其中的dtype为负数时,输出数组的type与输入数组的type相同;
否则,输出数组与输入数组只是通道数相同,而tpye=CV_MAT_DEPTH(dtype).
参考链接
实现流程:
1、载入原图
2、将图像扩大到合适的尺寸(当图像的尺寸是2.3.5的整数倍时,运行速度最快)
//将输入图像延扩到最佳尺寸,边界用0补充
3、为傅里叶变换的结果(实部和虚部)分配存储空间
4、进行离散傅里叶变化
5、将复数转化为幅值
6、进行对数尺度缩放(傅里叶变换的幅值范围大到不适合在屏幕上显示。为了凸显出高低变化的连续性,可以用对数尺度来替换线性尺度。)
公式:M1 = log(1+M);
7、剪切和重分布幅度图像限
因为在第二部中延扩了图像,现在需要将添加的像素剔除
为了方便显示,也可以重新分布幅度图像象限位置(将第五步得到的幅度图从中间划开得到4张1/4子图像,将每张子图都看成幅度图的一个象限,重新分布即将4个交点重叠到图像中心)这样的话原点(0,0)就为一道图像中心了。
//剪切和重分布幅度图像限
//如有奇数行或奇数列,进行频谱裁剪
//重新排列傅里叶图像中的象限,使得原点位于图像中心。
//交换象限(左上与右下进行交换)
//交换象限(左下与右上进行交换)
8、归一化
这一步仍然是为了显示。现在有了重分布后的幅度图,但是幅度值仍然超过了可显示范围【0,1】。这里使用归一化函数。
9、显示效果图
全部代码:
#include <opencv2/opencv.hpp>
#include <iostream>
#include "windows.h"
#include <stdio.h>
//#include "My_ImageProssing_base.h"using namespace cv;
using namespace std;void convolveDFT(Mat& A,Mat& B, Mat& C)
{//【1】初始化输出矩阵C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type());Size dftSize; //???什么意思//【2】计算DFT变换的尺寸dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);dftSize.height = getOptimalDFTSize(A.rows + 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】就地操作,进行快速傅里叶变换,并将nonzeroRows参数置为非零,以进行更加快速的处理???为什么dft(tempA, tempA, 0, A.rows);dft(tempB, tempB, 0, B.rows);//【6】将得到的频谱相乘,结果存放于tempA中mulSpectrums(tempA,tempB,tempA, DFT_COMPLEX_OUTPUT);//DFT_REAL_OUTPUT//【7】将结果变换为频域且尽管行结果都为非零,我们只需要其中C.rows的第一行,所以采用nonzeroRows==C.rowsdft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows);//【8】将结果复制到C中tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);//所有的临时缓冲区将被自动释放,所以无须收尾操作
}
int main()
{SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), FOREGROUND_INTENSITY | FOREGROUND_GREEN); //字体为绿色//1、载入原图Mat srcImage = imread("D:\\opencv_picture_test\\形态学操作\\coin_inv.png",0); //读取灰度图//2、将图像扩大到合适的尺寸(当图像的尺寸是2.3.5的整数倍时,运行速度最快)//【2】将输入图像延扩到最佳尺寸,边界用0补充int m = getOptimalDFTSize(srcImage.rows);int n = getOptimalDFTSize(srcImage.cols);//将添加的像素初始化为0Mat padded;copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, n - srcImage.cols, BORDER_CONSTANT,0);//3、为傅里叶变换的结果(实部和虚部)分配存储空间Mat planes[] = { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };Mat complexI;merge(planes, 2, complexI);//4、进行离散傅里叶变化dft(complexI, complexI);//5、将复数转化为幅值split(complexI, planes);//将多通道数组complexI分离成几个单通道数//planes[0] = Re(DFT(I));//planes[1] = Im(DFT(I));//计算矢量幅值magnitude(planes[0], planes[1], planes[0]);//将幅值存入planes[0] Mat magnitudeImage = planes[0];//6、进行对数尺度缩放magnitudeImage += Scalar::all(1);log(magnitudeImage, magnitudeImage);//就地操作,求自然对数//7、剪切和重分布幅度图像限magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));//这个&-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); //将q0与q3互换q3.copyTo(q0);tmp.copyTo(q3);//交换象限(左下与右上进行交换)q1.copyTo(tmp); //将q1与q2互换q2.copyTo(q1);tmp.copyTo(q2);//8、归一化//这一步仍然是为了显示。现在有了重分布后的幅度图,但是幅度值仍然超过了可显示范围【0, 1】。这里使用归一化函数。normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//9、显示效果图imshow("原图", srcImage);imshow("频谱幅值", magnitudeImage);waitKey(0);return 0;
}
srcImage:
padded:
magnitudeImage:
每个图像的具体性质: