opencv实现频域平滑图像

本文是参考opencv低通滤波器这篇文章改进得到的，地址：http://blog.csdn.net/huangli19870217/article/details/1285661。

关于高通滤波，其实只要在createfilter（）函数的基础上改写滤波器的函数表示就行，这样就可以实现了，有兴趣的可以自己试试。

#include<opencv2\core\core.hpp>
#include <opencv2\opencv.hpp>
#include <opencv2\highgui\highgui.hpp>
using namespace cv;
using namespace std;
//创建低通滤波器
void createfilter(Mat& lowTemp,float fc)
{
	//参数准备
	int i, j;
	int state = -1;
	double tempD;
	float width, height;
	width = lowTemp.cols;
	height =lowTemp.rows;
	float x, y;
	x = width / 2;
	y = height / 2;

	//D0截止频率
	double D0 = fc * min(lowTemp.rows, lowTemp.cols) / 2.0;

	for (i = 0; i < height; i++)
	{
		float* ptr = lowTemp.ptr<float>(i);
		for (j = 0; j < width; j++)
		{
			if (i > y && j > x)
			{
				state = 3;
			}
			else if (i > y)
			{
				state = 1;
			}
			else if (j > x)
			{
				state = 2;
			}
			else
			{
				state = 0;
			}

			switch (state)
			{
			case 0:
				tempD = (double)(i * i + j * j); tempD = sqrt(tempD); break;
			case 1:
				tempD = (double)((height - i) * (height - i) + j * j); tempD = sqrt(tempD); break;
			case 2:
				tempD = (double)(i * i + (width - j) * (width - j)); tempD = sqrt(tempD); break;
			case 3:
				tempD = (double)((height - i) * (height - i) + (width - j) * (width - j)); tempD = sqrt(tempD); break;
			default:
				break;
			}

			//二维高斯低通滤波器传递函数
			//tempD = exp(-0.5 * pow(tempD / D0, 2));
			//ptr[j*2]= tempD;
			//ptr[2*j+1]= tempD;



			//衰减系数为2的二维指数低通滤波器传递函数
			//tempD = exp(-pow(tempD / D0, 2));
			//ptr[j*2]= tempD;
			//ptr[2*j+1]= tempD;

			//2阶巴特沃思低通滤波器传递函数
			//tempD = 1 / (1 + pow(tempD / D0, 2 * 2));
			//ptr[j*2]= tempD;
			//ptr[2*j+1]= tempD;


			//二维理想低通滤波器传递函数
				if(tempD <= D0)
				{
					ptr[j*2] = 1.0;
			    	ptr[j*2+1] = 1.0;
				}
				else
				{
					ptr[j*2] = 0.0;
			    	ptr[j*2+1] = 0.0;
				}
		}
	}
}

int main(int argc, char ** argv)
{
	//【1】以灰度模式读取原始图像并显示
	Mat srcImg = imread("test1.jpg", 0);
	imshow("原始图像",srcImg);
	//【2】将输入图像延扩到最佳的尺寸，边界用0补充
	int m = getOptimalDFTSize(srcImg.rows);
	int n = getOptimalDFTSize(srcImg.cols);
	//将添加的像素初始化为0.
	Mat padded;
	copyMakeBorder(srcImg, padded, 0, m - srcImg.rows, 0, n - srcImg.cols, BORDER_CONSTANT, Scalar::all(0));
	//【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
	//将planes数组组合合并成一个多通道的数组complexI
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);
	//【4】进行就地离散傅里叶变换
	dft(complexI, complexI);
	//【5】将复数转换为幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
	split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组
	magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude  
	Mat magnitudeImage = planes[0];
	//【6】进行对数尺度(logarithmic scale)缩放
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);//求自然对数
	//【7】剪切和重分布幅度图象限
	//若有奇数行或奇数列，进行频谱裁剪      
	magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
	//重新排列傅立叶图像中的象限，使得原点位于图像中心  
	int cx = magnitudeImage.cols / 2;
	int cy = magnitudeImage.rows / 2;
	Mat q0(magnitudeImage, Rect(0, 0, cx, cy));   // ROI区域的左上
	Mat q1(magnitudeImage, Rect(cx, 0, cx, cy));  // ROI区域的右上
	Mat q2(magnitudeImage, Rect(0, cy, cx, cy));  // ROI区域的左下
	Mat q3(magnitudeImage, Rect(cx, cy, cx, cy)); // ROI区域的右下
	//交换象限（左上与右下进行交换）
	Mat tmp;
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	//交换象限（右上与左下进行交换）
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
	//【8】归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
	normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX);

	//【10】创建低通滤波器
	Mat lowTemp = Mat(complexI.size(),complexI.type());
	createfilter(lowTemp, 0.3);
	//将复数转换为幅值
	split(lowTemp, planes); // 将多通道数组complexI分离成几个单通道数组
	magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude  
	Mat magFilter = planes[0];
	//进行对数尺度(logarithmic scale)缩放
	magFilter += Scalar::all(1);
	log(magFilter, magFilter);//求自然对数
	//剪切和重分布幅度图象限
	//若有奇数行或奇数列，进行频谱裁剪      
	magFilter = magFilter(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
	//重新排列傅立叶图像中的象限，使得原点位于图像中心  
	int cx1 = magFilter.cols / 2;
	int cy1 = magFilter.rows / 2;
	Mat q01(magFilter, Rect(0, 0, cx1, cy1));   // ROI区域的左上
	Mat q11(magFilter, Rect(cx, 0, cx1, cy1));  // ROI区域的右上
	Mat q21(magFilter, Rect(0, cy1, cx1, cy1));  // ROI区域的左下
	Mat q31(magFilter, Rect(cx1, cy1, cx1, cy1)); // ROI区域的右下
	//交换象限（左上与右下进行交换）
	Mat tmp1;
	q01.copyTo(tmp1);
	q31.copyTo(q01);
	tmp1.copyTo(q31);
	//交换象限（右上与左下进行交换）
	q11.copyTo(tmp1);
	q21.copyTo(q11);
	tmp1.copyTo(q21);
	//归一化，用0到1之间的浮点值将矩阵变换为可视的图像格式
	normalize(magFilter, magFilter, 0, 1, CV_MINMAX);
	imshow("滤波器频谱图", magFilter);

	//【11】频域滤波
	mulSpectrums(complexI, lowTemp, complexI,0);
	//【12】傅里叶反变换
	idft(complexI, complexI, cv::DFT_SCALE);
	//【13】还原原始图像
	split(complexI, planes);
	Mat dft_filter_img;
	dft_filter_img.create(complexI.size(), CV_8UC1);
	planes[0].convertTo(dft_filter_img, CV_8UC1);
	dft_filter_img = dft_filter_img(cv::Rect(0, 0, srcImg.cols, srcImg.rows));//剪切图像

	imshow("低通滤波结果", dft_filter_img);
	waitKey(0);
	return 0;
}

原始图像：

滤波器频谱图：

滤波结果：