图像的傅里叶频谱

た入场券 2022-10-01 02:50 196阅读 0赞

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

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

# 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);
    }
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34 
       35 
       36 
       37 
       38 
       39 
       40 
       41 
       42 
       43 
       44 
       45 
       46 
       47 
       48 
       49 
       50 
       51 
       52 
       53 
       54 
       55 
       56 
       57 
       58 
       59 
       60 
       61 
       62 
       63 
       64 
       65 
       66 
       67 
       68 
       69 
       70 
       71 
       72 
       73 
       74 
       75 
       76 
       77 
       78 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34 
       35 
       36 
       37 
       38 
       39 
       40 
       41 
       42 
       43 
       44 
       45 
       46 
       47 
       48 
       49 
       50 
       51 
       52 
       53 
       54 
       55 
       56 
       57 
       58 
       59 
       60 
       61 
       62 
       63 
       64 
       65 
       66 
       67 
       68 
       69 
       70 
       71 
       72 
       73 
       74 
       75 
       76 
       77 
       78

从这里开始，还是简单的分析一下代码吧。

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)); 
     
      
       1 
       2 
       3 
       4 
       5 
       6 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6

这里参考了文献\[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);
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19

其实这里的很好理解的，将填充到最适尺寸的图像赋值给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 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9

以上代码，实现了以下计算。  
|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);
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34

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

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

# 3. 不用交换操作的代码 #

使用MATLAB去求取尺寸为M×N图像f的傅里叶频谱时候，通常会用`fft2(f,2*M,2*N)`。使用此函数求得的福利叶变换，其实还是\[0,2π\]范围内的傅里叶变换。要想使得傅里叶频谱的DC分量位于中心的话，其实还是需要一些别的操作的。冈萨雷斯的《数字图像处理》一书中，对于这个问题，是利用了傅里叶变换的评议特性，即  
f(x,y)=f(x,y)×(−1)x\+y  
然后再对函数f(x,y)进行福利叶变换，所得到结果，就是所需要的是\[−π,π\]内傅里叶变换的结果。具体的原理与推导，还是参看冈萨雷斯的《数字图像处理》英文原本的p258页左右，中文译本的p148左右。当然，嫌麻烦的话，还可以看我的博文，博文中也简明推导了平移特性。  
[\[数字图像处理\]频域滤波(1)–基础与低通滤波器][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);
    }
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34 
       35 
       36 
       37 
       38 
       39 
       40 
       41 
       42 
       43 
       44 
       45 
       46 
       47 
       48 
       49 
       50 
       51 
       52 
       53 
       54 
       55 
       56 
       57 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
       10 
       11 
       12 
       13 
       14 
       15 
       16 
       17 
       18 
       19 
       20 
       21 
       22 
       23 
       24 
       25 
       26 
       27 
       28 
       29 
       30 
       31 
       32 
       33 
       34 
       35 
       36 
       37 
       38 
       39 
       40 
       41 
       42 
       43 
       44 
       45 
       46 
       47 
       48 
       49 
       50 
       51 
       52 
       53 
       54 
       55 
       56 
       57

在这里，由于考虑到计算的原因，我将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);
        }
    }
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9 
      
     
      
       1 
       2 
       3 
       4 
       5 
       6 
       7 
       8 
       9

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

原图的傅里叶频谱  
![原图的结果][20150519171301897]  
使用11×11高斯滤波器后的傅里叶频谱  
![这里写图片描述][20150519171446167]

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

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

[OpenCV]: http://lib.csdn.net/base/opencv
[1]: http://blog.csdn.net/zhoufan900428/article/details/17194289
[20150519165535662]: /images/20220705/42a0df1899024973bb34c2c341d3420c.png
[20150519171301897]: /images/20220705/f4f4e2b13f2b43028656c41600c3e812.png
[20150519171446167]: /images/20220705/7d32ccf535db4b6eb3935971a93ce973.png