图像的离散傅里叶变换:
傅里叶变换可以将一幅图片分解为正弦和余弦两个分量,换而言之,他可以将一幅图像从其空间域(spatial domain)转换为频域(frequency domain)。这种变换的思想是任何函数可以很精确的接近无穷个sin()函数和cos()函数的和。傅里叶变换提供了这种方法来达到这种效果。
对于数字图像这种离散的信号,频率大小表示信号变化的剧烈程度或者说是信号变化的快慢。频率越大,变化越剧烈,频率越小,信号越平缓,对应到图像中,高频信号往往是图像中的边缘信号和噪声信号,而低频信号包含图像变化频繁的图像轮廓及背景等信号
二维图像的傅里叶变换数学表达式:
式中f(i, j)是图像空间域的值而F是频域的值。傅里叶转换的结果是复数,这也显示出了傅里叶变换是一副实数图像(real image)和虚数图像(complex image)叠加或者是幅度图像(magitude image)和相位图像(phase image)叠加的结果。在实际的图像处理算法中仅有幅度图像(magnitude image)图像能够用到,因为幅度图像包含了我们所需要的所有图像几何结构的信息。但是,如果想通过修改幅度图像或者相位图像来间接修改原空间图像,需要保留幅度图像和相位图像来进行傅里叶逆变换,得到修改后图像。
进行图像傅里叶变换需要用到OpenCV里面的一些函数:
1.dft()
opencv提供的傅里叶变换函数dft(),其定义如下:
C++: 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计算矩阵卷积时非常有效。
2. getOptimalDFTSize()
返回给定向量尺寸经过DFT变换后结果的最优尺寸大小。其函数定义如下:
C++: int getOptimalDFTSize(int vecsize);
参数解释:
- int vecsize: 输入向量尺寸大小(vector size)
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.
这个函数不能直接用于DCT(离散余弦变换)最优尺寸的估计,可以通过getOptimalDFTSize((vecsize+1)/2)*2得到。
3.magnitude()
计算二维矢量的幅值,其定义如下:
C++: void magnitude(InputArray x, InputArray y, OutputArray magnitude);
参数解释:
- InputArray x: 浮点型数组的x坐标矢量,也就是实部
- InputArray y: 浮点型数组的y坐标矢量,必须和x尺寸相同
- OutputArray magnitude: 与x类型和尺寸相同的输出数组
其计算公式如下:
4. copyMakeBorder()
扩充图像边界,其函数定义如下:
C++: void copyMakeBorder(InputArray src, OutputArray dst, int top, int bottom, int left, int right, int borderType, const Scalar& value=Scalar() );
参数解释:
- InputArray src: 输入图像
- OutputArray dst: 输出图像,与src图像有相同的类型,其尺寸应为Size(src.cols+left+right, src.rows+top+bottom)
- int类型的top、bottom、left、right: 在图像的四个方向上扩充像素的值
- int borderType: 边界类型,由borderInterpolate()来定义,常见的取值为BORDER_CONSTANT
- const Scalar& value = Scalar(): 如果边界类型为BORDER_CONSTANT则表示为边界值
5. normalize()
归一化就是把要处理的数据经过某种算法的处理限制在所需要的范围内。首先归一化是为了后面数据处理的方便,其次归一化能够保证程序运行时收敛加快。归一化的具体作用是归纳同意样本的统计分布性,归一化在0-1之间是统计的概率分布,归一化在某个区间上是统计的坐标分布,在机器学习算法的数据预处理阶段,归一化也是非常重要的步骤。其定义如下:
C++: void normalize(InputArray src, OutputArray dst, double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray() )
参数解释:
- InputArray src: 输入图像
- OutputArray dst: 输出图像,尺寸大小和src相同
- double alpha = 1: range normalization模式的最小值
- double beta = 0: range normalization模式的最大值,不用于norm normalization(范数归一化)模式
- int norm_type = NORM_L2: 归一化的类型,主要有
。NORM_INF: 归一化数组的C-范数(绝对值的最大值)
。NORM_L1: 归一化数组的L1-范数(绝对值的和)
。NORM_L2: 归一化数组的L2-范数(欧几里得)
。NORM_MINMAX: 数组的数值被平移或缩放到一个指定的范围,线性归一化,一般较常用。
- int dtype = -1: 当该参数为负数时,输出数组的类型与输入数组的类型相同,否则输出数组与输入数组只是通道数相同,而depth = CV_MAT_DEPTH(dtype)
- InputArray mask = noArray(): 操作掩膜版,用于指示函数是否仅仅对指定的元素进行操作。
示例程序:
//项目:图像离散傅里叶变换(DFT)
#include "stdafx.h"
#include <opencv2/opencv.hpp>
//#include "opencv2/core/core.hpp"
//#include "opencv2/imgproc/imgproc.hpp"
//#include "opencv2/highgui/highgui.hpp"//这三个头文件其实就包含在#include <opencv2/opencv.hpp>下了
#include <iostream>
#include<math.h>
using namespace std;
using namespace cv;
int main(int argc, char** argv)
{
Mat src = imread("C:/Users/59235/Desktop/image/national customs3.jpg",0);//以灰度图像模式读取图像
if (!src.data)
{
printf("could not load image...\n");
return -1;
}
namedWindow("original image", CV_WINDOW_AUTOSIZE);
imshow("original image", src);
//1、将图像扩大到合适的尺寸
int m = getOptimalDFTSize(src.rows);
int n = getOptimalDFTSize(src.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(src, padded, 0, m - src.rows, 0, n - src.cols, BORDER_CONSTANT, Scalar::all(0));
//2、为傅里叶变化的结果(实部和虚部)分配存储空间
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
//将planes数组组合合并成一个多通道的数组complexI
merge(planes, 2, complexI);
//3、进行离散傅里叶变换
dft(complexI, complexI);
//4、将复数转化为幅值,计算公式:log(1+sqrt(Re(DFT(I))^2+Im(DFT(I))^2))
split(complexI, planes);//将多通道数组分离成几个单通道数组,planes[0] = Re(DFT(src), planes[1] = Im(DFT(src))
magnitude(planes[0], planes[1], planes[0]);
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之间的浮点值将矩阵变换为可视的图像格式
normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX);
imshow("频谱幅值", magnitudeImage);
waitKey(0);
return 0;
}
效果图:
(原图)
(傅里叶变换频谱图)
从频谱图上可以看出,当将频谱移频到原点以后,图像中心比较亮。在频谱图中,一个点的亮暗主要与这个频率中点的数目和点的灰度值有关,也就是说在空间域中包含这种频率的点越多。而经过频移后,频率为0的部分,也就是傅里叶变换所得到的常量分量在图像中心,往外扩散,点所代表的频率越来越高。也说明了图像中的“能量”主要集中在低频部分。