1.图像的傅里叶频谱的意义
之前的博文其实已经归纳过这方面的内容了。我们常用的图像平滑处理,其实就是一个低通滤波,一定程度上去除高频信号,可以使得图像变得柔和(也就是平滑)。但是,在去除周期性噪声时候,空间域内的滤波(卷积)就不是那么好操作了。所以,这里时候,无论是理解起来方便,还是其他原因,都需要在频域内进行滤波。
详细的叙述还是在下面的博文里面啦!!!!
[数字图像处理]频域滤波(1)–基础与低通滤波器
[数字图像处理]频域滤波(2)–高通滤波器,带阻滤波器与陷波滤波器
2. 傅里叶频谱的计算
这部分的内容,主要就是使用openCV自带的函数
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 )
去求取图像的傅里叶变换。这里,输出结果CvArr* dst由两个通道组成,分别代表了实部与虚部。我们再根据如下算式,就可以得到傅里叶频谱了。
我自己也参考了很多人的代码,然后实现的代码如下。
IplImage* fft2(IplImage* image_input)
{
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
//cout<< " Width" << image_input->width << " " << dftWidth << "\n";
//cout<< "Height" << image_input->height << " " << dftHeight << "\n";
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
//image_Re <--- image_padded
cvConvertScale(image_padded,image_Re);
//image_Im <--- 0
cvZero(image_Im);
//image_Fourier[0] <--- image_Re
//image_Fourier[1] <--- image_Im
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
//image_Fourier[0] ---> image_Re
//image_Fourier[1] ---> image_Im
cvSplit(image_Fourier,image_Re,image_Im,0,0);
//Mag = sqrt(Re^2 + Im^2)
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
// log (1 + Mag)
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
// |-----|-----| |-----|-----|
// | 1 | 3 | | 4 | 2 |
// |-----|-----| ---> |-----|-----|
// | 2 | 4 | | 3 | 1 |
// |-----|-----| |-----|-----|
IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);
int cx = image_Re->width/2;
int cy = image_Re->height/2;
cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy)); // 1
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy)); // 4
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy)); // 4
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy)); // 1
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy)); // 3
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy)); // 2
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy)); // 2
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy)); // 3
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvResetImageROI(image_Re);
cvResetImageROI( Fourier);
cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);
return(Fourier);
}
从这里开始,还是简单的分析一下代码吧。
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
这里参考了文献[2]中的说法,在尺寸数为2,3,5的倍数的场合,计算的速度是最快的。所以使用函数getOptimalDFTSize()
来寻找最匹配的尺寸,然后再同伙cvCopyMakeBorder()
进行多余部分的填充,这里选的配置是将图放在从点(0,0)开始的位置,其余不足的地方,用0进行填充。
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
//image_Re <--- image_padded
cvConvertScale(image_padded,image_Re);
//image_Im <--- 0
cvZero(image_Im);
//image_Fourier[0] <--- image_Re
//image_Fourier[1] <--- image_Im
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
//image_Fourier[0] ---> image_Re
//image_Fourier[1] ---> image_Im
cvSplit(image_Fourier,image_Re,image_Im,0,0);
其实这里的很好理解的,将填充到最适尺寸的图像赋值给image_Re,将image_Im赋值为0。让后将这两层图复制到image_Fourier的两个通道里,然后使用函数cvDFT()
进行傅里叶变换。得到结果还是存在于image_Fourier的两个通道里,分别代表实部与虚部,然后通过cvSplit()
将其抽出到image_Re与image_Im里。
//Mag = sqrt(Re^2 + Im^2)
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
// log (1 + Mag)
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
以上代码,实现了以下计算。
还有就是进行了一个对数变换,这个也没的说,看傅里叶频谱的标配操作。
// |-----|-----| |-----|-----|
// | 1 | 3 | | 4 | 2 |
// |-----|-----| ---> |-----|-----|
// | 2 | 4 | | 3 | 1 |
// |-----|-----| |-----|-----|
IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);
int cx = image_Re->width/2;
int cy = image_Re->height/2;
cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy)); // 1
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy)); // 4
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy)); // 4
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy)); // 1
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy)); // 3
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy)); // 2
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy)); // 2
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy)); // 3
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvResetImageROI(image_Re);
cvResetImageROI( Fourier);
cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);
return(Fourier);
其实重头戏在这里,这里需要一个交换操作。至于为何所求得的傅里叶频谱为什么需要交换的原因是,这个代码求得的结果其实是范围
[数字图像处理]频域滤波(1)–基础与低通滤波器
这里,我使用了ROI操作与cvAddWeighted()
函数进行了实现。其运行的结果如下所示。
恩,基本可以看出来,直流分量也被我移动到了中心,以上代码实现了傅里叶频谱的计算与显示。
3. 不用交换操作的代码
使用MATLAB去求取尺寸为
fft2(f,2*M,2*N)
。使用此函数求得的福利叶变换,其实还是
然后再对函数
[数字图像处理]频域滤波(1)–基础与低通滤波器
为此,实现的代码变为了如下形式。
IplImage* fft2_New(IplImage* image_input)
{
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
cout<< " Width" << image_input->width << " " << dftWidth << "\n";
cout<< "Height" << image_input->height << " " << dftHeight << "\n";
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
//image_Re = image_padded .* (-1)^(x+y);
double pixel;
for(int y=0;y<image_padded->height;y++)
{
for(int x=0;x<image_padded->width;x++)
{
pixel = cvGetReal2D(image_padded,x,y);
pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
cvSetReal2D(image_Re,x,y,pixel);
}
}
//image_Im <--- 0
cvZero(image_Im);
//image_Fourier[0] <--- image_Re
//image_Fourier[1] <--- image_Im
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
//image_Fourier[0] ---> image_Re
//image_Fourier[1] ---> image_Im
cvSplit(image_Fourier,image_Re,image_Im,0,0);
//Mag = sqrt(Re^2 + Im^2)
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
// log (1 + Mag)
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
cvNormalize(image_Re,image_Re,1,0,CV_C,NULL);
return(image_Re);
}
在这里,由于考虑到计算的原因,我将
for(int y=0;y<image_padded->height;y++)
{
for(int x=0;x<image_padded->width;x++)
{
pixel = cvGetReal2D(image_padded,x,y);
pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
cvSetReal2D(image_Re,x,y,pixel);
}
}
其实也就相当于,
原图的傅里叶频谱
使用
从实验结果看来,可以看出以下两点
- 对于原图而言两个函数的结果基本一致,两个函数都得到了正确的结果。
- 使用平滑处理后,频谱的高频成分明显变小,对于空间域的图像而言,图像变得模糊了。
原文发于博客:http://blog.csdn.net/thnh169/
参考文献
[1]opencv 中 傅里叶变换 FFT :http://blog.csdn.net/abcjennifer/article/details/7359952
[2]OpenCV实现基于傅里叶变换的旋转文本校正 : http://johnhany.net/2013/11/dft-based-text-rotation-correction/#imageclose-380
[3]学习OpenCV范例(八)——离散傅立叶变换 : http://blog.csdn.net/chenjiazhou12/article/details/21240647
=============更新日志===================
2015 - 5 - 19 初版
2015 - 5 - 21 修改了某些叙述