码迷，mamicode.com
首页 > 其他好文 > 详细

[openCV]图像的傅里叶频谱

时间：2015-05-19 19:07:08      阅读：295      评论：0      收藏：0      [点我收藏+]

标签：opencv   图像合成   c   

1.图像的傅里叶频谱的意义

之前的博文其实已经归纳过这方面的内容了。我们常用的图像平滑处理，其实就是一个低通滤波，一定程度上去除高频信号，可以使得图像变得柔和（也就是平滑）。但是，在去除周期性噪声时候，空间域内的滤波（卷积）就不是那么好操作了。所以，这里时候，无论是理解起来方便，还是其他原因，都需要在频域内进行滤波。
详细的叙述还是在下面的博文里面啦！！！！
[数字图像处理]频域滤波(1)–基础与低通滤波器
[数字图像处理]频域滤波(2)–高通滤波器，带阻滤波器与陷波滤波器

2. 傅里叶频谱的计算

这部分的内容，主要就是使用openCV自带的函数
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 )
去求取图像的傅里叶变换。这里，输出结果CvArr* dst由两个通道组成，分别代表了实部与虚部。我们再根据如下算式，就可以得到傅里叶频谱了。
|F(u,v)|=R2(u,v)+F2(u,v)???????????????√2
我自己也参考了很多人的代码，然后实现的代码如下。

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);

以上代码，实现了以下计算。
|F(u,v)|=R2(u,v)+F2(u,v)???????????????√2
还有就是进行了一个对数变换，这个也没的说，看傅里叶频谱的标配操作。

//  |-----|-----|           |-----|-----|   
//  |  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);

其实重头戏在这里，这里需要一个交换操作。至于为何所求得的傅里叶频谱为什么需要交换的原因是，这个代码求得的结果其实是范围[0,2π]内的傅里叶变换。而我们所需要的是[?π,π]内的结果。详细的原因，还是在我以前的博客里有说明。
[数字图像处理]频域滤波(1)–基础与低通滤波器

这里，我使用了ROI操作与cvAddWeighted()函数进行了实现。其运行的结果如下所示。
技术分享
恩，基本可以看出来，直流分量也被我移动到了中心，以上代码实现了傅里叶频谱的计算与显示。

3. 不用交换操作的代码

其实使用Matlab的fft2()函数所求的，就是范围[0,2π]内的傅里叶变换。当时，解决的办法是，现将图像进行以下变换。
f(x,y)=f(x,y)×(?1)x+y
然后再对函数f(x,y)进行福利叶变换，所得到结果，就是所需要的是[?π,π]内的结果。为此，实现的代码变为了如下形式。

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);
}

在这里，由于考虑到计算的原因，我将f(x,y)=f(x,y)×(?1)x+y这个计算，用下面代码去进行了实现。

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);
    }
}

其实也就相当于，x+y是偶数的话f(x,y)不变，x+y是奇数的话f(x,y)变为负。用这样一个简单的判断去实现。其运行的结果，如下所示。

原图的傅里叶频谱
技术分享
使用11×11高斯滤波器后的傅里叶频谱
技术分享

从实验结果看来，可以看出以下两点

  1. 对于原图而言两个函数的结果基本一致，两个函数都得到了正确的结果。
  2. 使用平滑处理后，频谱的高频成分明显变小，对于空间域的图像而言，图像变得模糊了。

原文发于博客：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

[openCV]图像的傅里叶频谱

标签：opencv   图像合成   c   

原文地址：http://blog.csdn.net/zhoufan900428/article/details/45845957
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！