OpenCV笔记21 频率域滤波

一、概述

       图像的傅里叶变换及其两个重要的度量：幅度谱和相位谱。了解两个重要的概念：低频和高频。低频指的是图

的傅里叶变换 “ 中心位置 ” 附近的区域。注意，如无特殊说明，后面所提到的图像的傅里叶变换都是中心化后的。高频随着到“ 中心位置 ” 距离的增加而增加，即傅里叶变换中心位置的外围区域，这里的“ 中心位置 ” 指的是傅里叶变换所对应的幅度谱最大值的位置。

频率域滤波器在程序或者数学运算中的呈现可以理解为一个矩阵，该矩阵的宽、高和图像的傅里叶变换的宽、高是相同的，下面所涉及的常用的低通、高通、带通、带阻 等滤波的关键步骤，就是通过一定的准则构造该矩阵的。

二、步骤

步骤如下：

下面通过一个简单的例子来详细解释频率域滤波的整个步骤。

第一步：输入图像矩阵 I 。假设为：

第二步：图像矩阵的每一个像素值乘以 （-1）r+c 得到矩阵 I′ ， I′ =I.* （ -1 ） r+c ，其

中 r 和 c 代表当前像素值在矩阵中的位置索引。

第三步：因为图像矩阵的宽和高均为 7 ，为了利用傅里叶变换的快速算法，对 I′ 补 0 ，使用命令getOptimalDFTSize （ 7 ）得到一个不小于 7 且可以分解为 2p ×3q ×5r 的最小整数，计算结果为8 。所以在矩阵 I′ 的右侧和下侧各补一行 0 ，记为 f ：

 第四步：利用傅里叶变换的快速算法得到复数矩阵F。

 

        OpenCV是将复数矩阵按照双通道存储的，即第一通道存储的是复数矩阵的实部，第二通道存储的是复数矩阵的虚部。

第五步：构建频率域滤波器 Filter 。频率域滤波器本质上是一个和第四步得到的快速傅里叶变换矩阵F 具有相同行数、列数的复数矩阵，一般情况下为实数矩阵，这里假设是 一个全是1的矩阵：

本章提到的频率域滤波，如低通滤波、高通滤波、自定义滤波等，其关键步骤就是 通过一定的标准构造该矩阵以完成图像在频率域上的滤波。

第六步：将第四步得到的快速傅里叶变换矩阵 F 和第五步得到的频率域滤波器 Filter 的对应位置相乘（矩阵的点乘）。当然，如果滤波器是一个实数矩阵，那么在代码实现 中，将傅里叶变换的实部和虚部分别与频率域滤波器进行点乘即可，即Ffilter =F.*Filter，因为这里构造的滤波器是一个全是1 的矩阵，所以 F filter =F 。

第七步：对第六步得到的点乘矩阵 Ffilter 进行傅里叶逆变换，得到复数矩阵 F′ 。

第八步： 取复数矩阵F′的实部 。

第九步：与第二步类似，将第八步得到的矩阵乘以（ -1 ） r+c 。

第十步：因为在快速傅里叶变换的步骤中进行了补 0 操作，所以第九步得到的实部矩 阵的尺寸有可能比原图大，所以要进行裁剪，取该实部矩阵的左上角，尺寸和原图相同。裁剪得到的结果，即为频率域滤波的结果。在该示例中，因为滤波器是一个全是1 的

矩阵，相当于对原图没有做任何处理，即最后滤波的结果和原图是一样的。

频率域滤波算法均是按照上述十个步骤完成的，接下来详细介绍常用滤波器的构建 方法、代码实现及其效果。

三、低通滤波

 针对图像的傅里叶变换，低频信息表示图像中灰度值缓慢变化的区域；而高频信息则正好相反，表示灰度值变化迅速的部分，如边缘。低通滤波，顾名思义，保留傅里叶变换的低频信息；或者削弱傅里叶变换的高频信息；而高通滤波则正好相反，保留傅里叶变换的高频信息，移除或者削弱傅里叶变换的低频信息。

三种常用的低通滤波器：

               H、 W 分别代表图像快速傅里叶变换的高、宽， 傅里叶谱的最大值（中心点）的位置在maxR，maxC , radius 代表截断频率， D （ r ， c ）代表到中心位置的距离

1、理想低通滤波：ilpFilter=[ilpFilter(r,c)]H*w,

 2、巴特沃斯低通滤波器

第二种是巴特沃斯低通滤波器，记 blpFilter=[blpFilter （ r ， c ） ]H×W

 3、高斯低通滤波器

记glpFilter=[glpFilter（r，c）]H×W；

 作用：

滤波器越靠近中心点位置的值越接近于1，越远离中心位置的值就越小于1，与傅里叶变换相乘后，相当于保留了低频信息，消弱或者移除了高频信息

四、低通滤波的代码实现

// 快速傅里叶变换
void  fft2Image(InputArray  image, OutputArray F)
{
	//的到Mat 的类型
	Mat  Image = image.getMat();
	int rows = Image.rows;
	int cols = Image.cols;
	// 获得最佳两列
	int  r = getOptimalDFTSize(rows);
	int  c = getOptimalDFTSize(cols);

	// 在左侧和下面进10
	Mat  f;
	copyMakeBorder(Image, f, 0, r - rows, 0, c - cols, BORDER_CONSTANT,Scalar::all(0));
	// 两个通道，第一个通道用存储实部  第二个通道是存放虚部的
	dft(f, F, DFT_COMPLEX_OUTPUT);

	//傅里叶逆变换 
	//dft(f, F, DFT_INVERSE+DFT_REAL_OUTPUT+DFT_SCALE);//只返回实部
}


// 幅度谱
void amplitudeSpectrum(InputArray  _srcfft, OutputArray  _dstSpectrum) 
{
	if (_srcfft.channels() != 2) return ;
	// 如果是两个通道则开始分离通道  

	vector<Mat> FFT2Channels;
	split(_srcfft, FFT2Channels);

   // 计算傅里叶变换幅度值    sqrt(pow(r,2)+pow(L,2))
	magnitude(FFT2Channels[0], FFT2Channels[1], _dstSpectrum);

}
// 傅里叶谱的灰度级显示
Mat graySpectrum(Mat  s) 
{
	Mat  dst;
	log(s + 1, dst);
	// 归一化
	normalize(dst,dst,0,1,NORM_MINMAX);
	// 为了显示 灰度级  
	dst.convertTo(dst,CV_8UC1,255,0);
	return dst;
}
//  计算傅里叶变换的相位角==phase 算子
Mat  phaseS(Mat sfft)
{
	// 相位谱
	Mat  phase;
	phase.create(sfft.size(), CV_64FC1);
	// 分离通道
	vector<Mat> fft2Channels;
	split(sfft,fft2Channels);
	// 开始计算相位谱
	for (int r = 0; r < phase.rows;r++) {
		for (int c = 0; c < phase.cols; c++) {
			// 实部  虚部
			double  real = fft2Channels[0].at<double>(r,c);
			double  img = fft2Channels[1].at<double>(r,c);
			// atan2 的返回值就是【0,180】，【-180,0】

			phase.at<double>(r, c) = atan2(img,real);
		}
	}
	return phase;
}
// -----------------------------------------------------------

#if  1 //低通滤波

enum typeFilter 
{
	 ILPFILTER,
	 BLPFILTER,
	 GLPFILTER,
};
Mat   createFilter(Size   size, Point  center,float radius, int  type, int n = 2)
{

	Mat  ilpFilter = Mat::zeros(size,CV_32FC1);  // 定义一个size大小的滤波器
	int  rows = size.height;
	int  cols = size.width;

	if (radius <= 0)return ilpFilter;


	// 构建理想的低通滤波器
	if (type== ILPFILTER)
	{
		for (int  r=0;r< rows;r++)
		{
			for (int c = 0; c < cols;c++) 
			{
				// 计算距离
				float  norm2 = pow(abs(float(r - center.y)), 2)+ pow(abs(float(r - center.y)), 2);
				if (norm2 < radius)
				{
					ilpFilter.at<float>(r, c) = 1;
				}
				else 
				{
					ilpFilter.at<float>(r, c) = 0;
				}

			}
		}
	}
	// 
	else if(type == BLPFILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
	/*			float m = sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0));
				float  n = m / radius;
				float  o = pow(n,2*n);
				float  s = 1.0 + o;
				float  i = 1.0 / s;*/
				ilpFilter.at<float>(r, c) = float(1.0 / (1.0 + (pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / radius, 2 * n))));

			}
		}
	}
	// 高斯滤波
	else if (type == GLPFILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
				ilpFilter.at<float>(r, c) = float(exp(-(pow(c-center.x,2.0)+pow(r-center.y,2.0))/(2*pow(radius,2.0))));
			}
		}
	}
	return  ilpFilter;
}

Mat    src;// 输入的图像
Mat    F;// 图像的快速傅里叶变换
Point    maxLocation;//傅里叶谱的最大值的坐标
int   radius = 20; // 截断频率
const   int  max_Radius = 100;// 最大截断频率



Mat  ilpFilter;//低通滤波器
int  filterType=0;//
int  max_FilterType = 2;


Mat   F_ilpFilter;// 低通傅里叶变换
Mat   FlpSpetrum;// 低通傅里叶变换的傅里叶谱的灰度级

Mat  result;// 滤波后的图像

string  lpFilterSpectrum = "低通傅里叶谱";// 




void   callback_lpFilter(int ,void *);


int main()
{
	 src = imread("C:\\Users\\19473\\Desktop\\opencv_images\\611.png",CV_LOAD_IMAGE_GRAYSCALE);
	if (!src.data)
	{
		cout << "原始图读取失败" << endl;
		return -1;
	}
	imshow("原图",src);
	// 数据转换    将数据转换为double 
	Mat  f1; 
	src.convertTo(f1,CV_32FC1,1.0,0.0);
	// 2   每个数  x pow(-1,c+r);
	for (int r = 0; r < f1.rows; r++)
	{
		for (int c = 0; c < f1.cols; c++)
		{
			// 判断奇偶
			if (r+c%2) 
			{
				f1.at<float>(r, c) *= -1;
			}
		}
	}

	// 开始傅里叶变换
	fft2Image(f1,F);// F 就是傅里叶变换之后的
	
	Mat as;
	//得到傅里叶谱
	amplitudeSpectrum(F,as);

	//傅里叶谱的灰度级显示
	Mat  s = graySpectrum(as);
	imshow("原傅里叶谱的灰度级显示",s);

	// 找到傅里叶谱的最大值的坐标
	minMaxLoc(s,NULL,NULL,NULL,&maxLocation);

	//------------------------------------ 低通滤波----------------------------------------------
	namedWindow(lpFilterSpectrum,WINDOW_AUTOSIZE);
	createTrackbar("低通滤波：", lpFilterSpectrum,&filterType, max_FilterType, callback_lpFilter);
	createTrackbar("半径：", lpFilterSpectrum, &radius, max_Radius, callback_lpFilter);
	callback_lpFilter(0,0);
	waitKey(0);

	return 0;
}
// 回调函数
void   callback_lpFilter(int, void*) 
{
	// ----------------- 构建低通滤波器----------------------------------
	ilpFilter = createFilter(F.size(),maxLocation,radius, filterType,2);
	// -------------------- 低通滤波器和图像的傅里叶变换开始点乘=====================================
	F_ilpFilter.create(F.size(),F.type());


	for (int r = 0; r < F_ilpFilter.rows; r++)
	{
		for (int c = 0; c < F_ilpFilter.cols; c++)
		{
			// 分别取出当前位置的快速傅里叶变换和理想低通滤波的值
			Vec2f f_rc = F.at<Vec2f>(r,c);
			float  lpFilter_rc = ilpFilter.at<float>(r,c);
			// 低通滤波器和图像的快速傅里叶变换的对应位置相乘
			F_ilpFilter.at<Vec2f>(r, c) = f_rc * lpFilter_rc;
		}
	}

	// 低通傅里叶变换的傅里叶谱
	amplitudeSpectrum(F_ilpFilter, FlpSpetrum);

	FlpSpetrum = graySpectrum(FlpSpetrum);
	imshow(lpFilterSpectrum, FlpSpetrum);



	// 对傅里叶变换逆变换  并且只要实部
	dft(F_ilpFilter,result,DFT_SCALE+DFT_INVERSE+DFT_REAL_OUTPUT);


	// 乘以  （-1）的（r+c）次方
	for (int r = 0; r < result.rows; r++)
	{
		for (int c = 0; c < result.cols; c++)
		{
			if ((c + r) % 2) result.at<float>(r, c) *= -1;
		}
	}

	// 很重要的一步  将结果转换成 CV_8u 
	result.convertTo(result, CV_8UC1, 1.0, 0);
	result = result(Rect(0,0,src.cols,src.rows)).clone();
	imshow("经过低通滤波后的图像",result);

}
#endif

 五、高通滤波

maxR ， maxC 代表图像傅里叶谱的最大值的位置，D(r,c):

1、 理想高通滤波

2. 巴特沃斯高通滤波器

3. 高斯高通滤波器

 六、带通 带阻滤波器

带通滤波是指只保留某一范围区域的频率带。下面介绍三种常用的带通滤波器。假 设BW 代表带宽， D0 代表带宽的径向中心，其他符号与低通、高通滤波相同。

1. 理想带通滤波器和理想带阻滤波器

 2. 巴特沃斯带通滤波器和巴特沃斯带阻滤波器

 3. 高斯带通滤波器和高斯带阻滤波器

 带阻滤波器和这个上述的恰好相反

七、自定义滤波器

#if  1 //自定义滤波器

Mat image;// 输入的图像
Mat    ImageFFT;// 图像的快速傅里叶变换
Point    maxLocation;//傅里叶谱的最大值的坐标



Mat  ilpFilter;//低通滤波器


Mat   F_ImageSpetrum;// 傅里叶变换的傅里叶谱

Mat  result;// 滤波后的图像

string  windowName = "傅里叶幅度谱的灰度级";// 

bool  drawing_box = false;// 鼠标事件

Point drawPoint;
Rect rectFilter;
bool  gotRectFilter = false;

void  mouseRectHandler(int event,int x,int  y,int,void* ) 
{
	switch (event)
	{
	case CV_EVENT_LBUTTONDOWN:// 鼠标按下
		drawing_box = true;
		// 记录起点
		drawPoint = Point(x, y);
		break;
	case CV_EVENT_MOUSEMOVE:
		if (drawing_box)
		{
			// 将鼠标移动到右下角
			if (x> drawPoint.x&&y>= drawPoint.y)
			{
				rectFilter.x = drawPoint.x;
				rectFilter.y = drawPoint.y;
				rectFilter.width = x - drawPoint.x;
				rectFilter.height= y - drawPoint.y;
			}

			// 将鼠标移动到有右上角
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = drawPoint.x;
				rectFilter.y = y;
				rectFilter.width = x - drawPoint.x;
				rectFilter.height = drawPoint.y- y ;
			}

			// 将鼠标移动到有左上角
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = x;
				rectFilter.y = y;
				rectFilter.width = drawPoint.x - x;
				rectFilter.height = drawPoint.y - y;
			}

			// 将鼠标移动到有左下角
			if (x >= drawPoint.x && y <= drawPoint.y)
			{
				rectFilter.x = x;
				rectFilter.y = drawPoint.y;
				rectFilter.width =  drawPoint.x-x;
				rectFilter.height = y- drawPoint.y ;
			}
		}
		break;
	case CV_EVENT_LBUTTONUP:
		drawing_box = false;
		gotRectFilter = true;
		break;
	default:
		break;
	}
}


int main()
{
	image = imread("C:\\Users\\19473\\Desktop\\opencv_images\\612.png",CV_LOAD_IMAGE_GRAYSCALE);
	if (!image.data)
	{
		cout << "原始图读取失败" << endl;
		return -1;
	}
	imshow("原图", image);
	// 数据转换    将数据转换为double 
	Mat  f1; 
	image.convertTo(f1,CV_32FC1,1.0,0.0);
	// 2   每个数  x pow(-1,c+r);
	for (int r = 0; r < f1.rows; r++)
	{
		for (int c = 0; c < f1.cols; c++)
		{
			// 判断奇偶
			if (r+c%2) 
			{
				f1.at<float>(r, c) *= -1;
			}
		}
	}

	// 开始傅里叶变换
	fft2Image(f1,ImageFFT);// F 就是傅里叶变换之后的
	
	Mat amplSpec;
	//得到傅里叶谱
	amplitudeSpectrum(ImageFFT, amplSpec);

	//傅里叶谱的灰度级显示
	Mat  spectrum = graySpectrum(amplSpec);
	

	// 找到傅里叶谱的最大值的坐标
	minMaxLoc(amplSpec,NULL,NULL,NULL,&maxLocation);

	//------------------------------------ 自定义  滤波----------------------------------------------
	namedWindow(windowName,WINDOW_AUTOSIZE);

	setMouseCallback(windowName, mouseRectHandler,NULL);
	for (;;) 
	{

		spectrum(rectFilter).setTo(0);
		// 自定义 滤波器和傅里叶变换点乘
		ImageFFT(rectFilter).setTo(Scalar::all(0));
		imshow(windowName, spectrum);
		// 按esc 见退出编辑
		if (waitKey(10)==27)
		{
			break;
		}

	}

	Mat  result;

	// 对傅里叶变换逆变换  并且只要实部
	dft(ImageFFT, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);


	// 乘以  （-1）的（r+c）次方
	for (int r = 0; r < result.rows; r++)
	{
		for (int c = 0; c < result.cols; c++)
		{
			if ((c + r) % 2) result.at<float>(r, c) *= -1;
		}
	}

	// 很重要的一步  将结果转换成 CV_8u 
	result.convertTo(result, CV_8UC1, 1.0, 0);
	result = result(Rect(0, 0, image.cols, image.rows)).clone();
	imshow("经过自定义滤波后的图像", result);
	waitKey(0);

	return 0;
}
// 回调函数

#endif