（图像处理之滤波）OpenCV实现频率域的低通高斯滤波（C++）

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

本文链接： https://blog.youkuaiyun.com/weixin_40647819/article/details/80600918

  最近由于作业原因，试着用OpenCV实现频率域滤波，但是OpenCV中并没有像MATLAB中fftshift这样的中心化操作，所以我写了一个频率域滤波的函数，以后用频率域滤波的时候在主函数中调用即可。当然，水平有限，编写的代码并不优美，有问题请大家留言批评指正。

  在这里我不介绍傅里叶变换，频率域滤波和高斯低通滤波器的原理，想必大家已经有了大概了解，本文关注于OpenCV中的代码实现。

 废话少说，先上一张例图：

这幅图像的噪声为周期性噪声，带用阻滤波器或者陷波滤波器去噪较好，这里用高斯低通去噪，效果一般，周期性没有完全去除。

下面是频率域滤波的流程

1、计算原图像的DFT，得到F（U，V）;

2、将频谱F（U，V）的零频点移到中心位置;

3、计算滤波器函数H（U，V）与上步结果F（U，V）的乘积G（U，V）;

4、将上一步结果G（U，V）的零频点移回;（在OpenCV中没有这步也没关系，原因待解）

5、计算上步的IDFT，得到g（X，Y）;

6、取g（X，Y）的模作为结果。

下面根据上述步骤编写代码如下：

    
    

      
      

       
       

      
      
      
      

       
       

        
        //**************************************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        //频率域滤波——以高斯低通为例
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        //**************************************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        #include<opencv2/opencv.hpp>
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        #include<iostream>
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        using namespace std;
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        using namespace cv;
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        Mat gaussianlbrf(Mat scr,
        
        float sigma);
        
        //高斯低通滤波器函数
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        Mat freqfilt(Mat scr,Mat blur);
        
        //频率域滤波函数
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        int main()
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        {
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat input=imread(
        
        "imageHill.jpg",CV_LOAD_IMAGE_GRAYSCALE);
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	imshow(
        
        "原图",input);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        int w=getOptimalDFTSize(input.cols);  
        
        //获取进行dtf的最优尺寸
       
       
      
      

      
      

       
       

      
      
      
      

       
       
    
        
        int h=getOptimalDFTSize(input.rows); 
        
        //获取进行dtf的最优尺寸
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat padded;  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            copyMakeBorder(input,padded,
        
        0,h-input.rows,
        
        0,w-input.cols,  BORDER_CONSTANT  , Scalar::all(
        
        0));  
        
        //边界填充
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            padded.convertTo(padded,CV_32FC1); 
        
        //将图像转换为flaot型
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat gaussianKernel=gaussianlbrf(padded,
        
        45);
        
        //高斯低通滤波器
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat 
        
        out=freqfilt(padded,gaussianKernel);
        
        //频率域滤波
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	imshow(
        
        "结果图 sigma=40",
        
        out);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	waitKey(
        
        0);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        return 
        
        0;
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        }
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        //*****************高斯低通滤波器***********************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        Mat gaussianlbrf(Mat scr,
        
        float sigma)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        {
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
          Mat gaussianBlur(scr.size(),CV_32FC1); 
        
        //，CV_32FC1
       
       
      
      

      
      

       
       

      
      
      
      

       
       
  
        
        float d0=
        
        2*sigma*sigma;
        
        //高斯函数参数，越小，频率高斯滤波器越窄，滤除高频成分越多，图像就越平滑
       
       
      
      

      
      

       
       

      
      
      
      

       
       
  
        
        for(
        
        int i=
        
        0;i<scr.rows ; i++ )
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
          {
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	  
        
        for(
        
        int j=
        
        0; j<scr.cols ; j++ )
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	  {
       
       
      
      

      
      

       
       

      
      
      
      

       
       
		  
        
        float d=pow(
        
        float(i-scr.rows/
        
        2),
        
        2)+pow(
        
        float(j-scr.cols/
        
        2),
        
        2);
        
        //分子,计算pow必须为float型
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	      gaussianBlur.at<
        
        float>(i,j)=expf(-d/d0);
        
        //expf为以e为底求幂（必须为float型）
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	  }
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
          }
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
           imshow(
        
        "高斯低通滤波器",gaussianBlur);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
        
        return gaussianBlur;
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        }
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        //*****************频率域滤波*******************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        Mat freqfilt(Mat scr,Mat blur)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        {
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        //***********************DFT*******************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat plane[]={scr, Mat::zeros(scr.size() , CV_32FC1)}; 
        
        //创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数）
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat complexIm;  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            merge(plane,
        
        2,complexIm);
        
        //合并通道 （把两个矩阵合并为一个2通道的Mat类容器）
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	dft(complexIm,complexIm);
        
        //进行傅立叶变换，结果保存在自身  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        //***************中心化********************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            split(complexIm,plane);
        
        //分离通道（数组分离）
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       
    
        
        int cx=plane[
        
        0].cols/
        
        2;
        
        int cy=plane[
        
        0].rows/
        
        2;
        
        //以下的操作是移动图像  (零频移到中心)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part1_r(plane[
        
        0],Rect(
        
        0,
        
        0,cx,cy));  
        
        //元素坐标表示为(cx,cy)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part2_r(plane[
        
        0],Rect(cx,
        
        0,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part3_r(plane[
        
        0],Rect(
        
        0,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part4_r(plane[
        
        0],Rect(cx,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat temp; 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part1_r.copyTo(temp);  
        
        //左上与右下交换位置(实部)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part4_r.copyTo(part1_r);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part4_r);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
      
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part2_r.copyTo(temp);  
        
        //右上与左下交换位置(实部)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part3_r.copyTo(part2_r);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part3_r);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat part1_i(plane[
        
        1],Rect(
        
        0,
        
        0,cx,cy));  
        
        //元素坐标(cx,cy)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part2_i(plane[
        
        1],Rect(cx,
        
        0,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part3_i(plane[
        
        1],Rect(
        
        0,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part4_i(plane[
        
        1],Rect(cx,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	 part1_i.copyTo(temp);  
        
        //左上与右下交换位置(虚部)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part4_i.copyTo(part1_i);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part4_i);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
      
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part2_i.copyTo(temp);  
        
        //右上与左下交换位置(虚部)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part3_i.copyTo(part2_i);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part3_i);
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        //*****************滤波器函数与DFT结果的乘积****************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat blur_r,blur_i,BLUR;  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	multiply(plane[
        
        0], blur, blur_r); 
        
        //滤波（实部与滤波器模板对应元素相乘）
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	multiply(plane[
        
        1], blur,blur_i);
        
        //滤波（虚部与滤波器模板对应元素相乘）
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	Mat plane1[]={blur_r, blur_i};
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            merge(plane1,
        
        2,BLUR);
        
        //实部与虚部合并
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	  
        
        //*********************得到原图频谱图***********************************
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	magnitude(plane[
        
        0],plane[
        
        1],plane[
        
        0]);
        
        //获取幅度图像，0通道为实部通道，1为虚部，因为二维傅立叶变换结果是复数  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	plane[
        
        0]+=Scalar::all(
        
        1);  
        
        //傅立叶变o换后的图片不好分析，进行对数处理，结果比较好看 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            log(plane[
        
        0],plane[
        
        0]);    
        
        // float型的灰度空间为[0，1])
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	normalize(plane[
        
        0],plane[
        
        0],
        
        1,
        
        0,CV_MINMAX);  
        
        //归一化便于显示
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            imshow(
        
        "原图像频谱图",plane[
        
        0]);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        //******************IDFT*******************************
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        /*
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part111(BLUR,Rect(0,0,cx,cy));  //元素坐标(cx,cy)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part222(BLUR,Rect(cx,0,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part333(BLUR,Rect(0,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            Mat part444(BLUR,Rect(cx,cy,cx,cy));  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	 part111.copyTo(temp);  //左上与右下交换位置(虚部)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part444.copyTo(part111);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part444);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
              
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part222.copyTo(temp);  //右上与左下交换位置
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            part333.copyTo(part222);  
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            temp.copyTo(part333);
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	*/
       
       
      
      

      
      

       
       

      
      
      
      

       
       
 
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	idft( BLUR, BLUR);	
        
        //idft结果也为复数
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        	split(BLUR,plane);
        
        //分离通道，主要获取通道
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            magnitude(plane[
        
        0],plane[
        
        1],plane[
        
        0]);  
        
        //求幅值(模)
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
            normalize(plane[
        
        0],plane[
        
        0],
        
        1,
        
        0,CV_MINMAX);  
        
        //归一化便于显示
       
       
      
      

      
      

       
       

      
      
      
      

       
       
	
        
        return plane[
        
        0];
        
        //返回参数
       
       
      
      

      
      

       
       

      
      
      
      

       
       

        
        }
       
       
      
      
    
    

效果图：

注：1，结果的黑边为填充效果。

      2，务必注意矩阵数据类型，做DFT，必须为浮点型，计算滤波器函数H（U，V）与DFT的乘积G（U，V）时，数据类型务必一致，包括通道数，单通道类型不能与多通道类型相乘，关于数据类型，可以参看这里。

参考：https：//blog.youkuaiyun.com/dang_boy/article/details/76150067