OpenCV实现频率域滤波——以高斯低通滤波去噪为例

  最近由于作业原因，试着用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.csdn.net/dang_boy/article/details/76150067