OpenCV实现图像基础频率域滤波

写在前面：刚接触数字图像处理的频域滤波时，很陌生。感觉这个技术的范围很窄，还不如直接处理空域那么真实。最近发现有些情况不得不在频域操作。感官图像的恢复和重建是最大的应用点。特此实现一些基本的频域滤波操作，为后期学习打好基础…

1. 频率域滤波步骤

前处理：包括对图像边界填充，使之达到OpenCV傅里叶变换最佳尺寸，然后就是将f(x,y)乘以(-1)^(x+y)，使傅里叶变换位于填充后图像大小的频率矩形的中心。

后处理：包括从反变换回来的图的左上象限提取一个大小为MXN(原图尺寸)的区域，得到与输入图像大小相同的滤波后的结果g(x, y)。

具体操作步骤详见《数字图像处理》第4章。

2.总结

1. H(u, v)是实函数，即虚部为0。
2. G(u, v) = H(u, v)*F(u, v) = H(u, v) *(R(u, v) + jI(u, v)) =H(u, v)*R(u, v) + H(u, v)*jI(u, v)。
3. 反变换获得图像g(x, y)取实部即可。
4. 在进行高通滤波时，滤波后的图像中会含有负值，这时直接将像素值归一化到[0, 1]，得到的结果同时标定了正值和负值，图像呈现深灰色，虽然也对，但是为了呈现另一种效果（底色为黑）我将负值进行了手动截断为0，然后在进行归一化。

3.程序实现

//***************************************tool***************************************//


//function:	使傅里叶频谱中心化
void fftshift(const cv::Mat& src, cv::Mat& dst)
{
	dst = src.clone();
	if (dst.type() != CV_32F) dst.convertTo(dst, CV_32F);
	
	for (int r = 0; r < dst.rows; ++r)
	{
		for (int c = 0; c < dst.cols; ++c)
		{
			if ((r + c) % 2 != 0)	//r + c 不为偶数，当前点变为负数
				dst.at<float>(r, c) = -dst.at<float>(r, c);
		}
	}
}


//function:	创建最佳傅里叶变换尺寸
cv::Size Make_BetterSize(const cv::Mat src)
{
	int w = cv::getOptimalDFTSize(src.cols);
	int h = cv::getOptimalDFTSize(src.rows);

	//width, height
	return Size(w, h);
}


//function:	将图像中负数截断为0
//opencv中应该有现成的工具函数，需要查一下文档
void MinusToZero(const cv::Mat& src, cv::Mat& dst)
{
	dst = src.clone();
	
	for (int r = 0; r < dst.rows; ++r)
	{
		for (int c = 0; c < dst.cols; ++c)
		{
			if (dst.at<float>(r, c) < 0)
				dst.at<float>(r, c) = 0.0;
		}
	}
}


//*****************************calculate filtering kernel*****************************//


//function:	创建理想低通滤波器
void calcILPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	Mat H(filterSize, CV_32FC1, Scalar::all(0));

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//当距离小于截止频率时，设置为1
			if (d <= (double)D0)
				H_data[x] = 1.f;
		}
	}
	//改为双通道图像便于与F(u,v)相乘
	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	巴特沃斯低通滤波器
void calcButterWorth(const cv::Mat& src, int D0, Mat& dst, float n)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据巴特沃斯低通滤波器计算H(u, v)
			H_data[x] = 1.0 / (1 + pow(d / D0, 2 * n));
		}
	}
	
	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	创建高斯低通滤波器
void calcGLPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据gauss低通滤波器计算H(u, v)
			H_data[x] = expf((-d * d) / (2 * D0 * D0));
		}
	}

	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	创建高斯高通滤波器
//这里未做简化，只需要用1减去高斯低通核即可创建高通核
void calcGHPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::ones(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据gauss低通滤波器计算H(u, v)
			H_data[x] = 1.f - expf((-d * d) / (2 * D0 * D0));
		}
	}

	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}

//********************************执行G(u,v) = F(u,v) * H(u,v)*******************************//


//function:		将输入图像在频率域进行处理，然后输出
//InputImage:	输入单通道图像
//H:			频率域滤波核
void frequency_filter(const cv::Mat& InputImage, const cv::Mat& H, cv::Mat& OutputImage)
{
	cv::Mat padded;
	int m = getOptimalDFTSize(InputImage.rows);
	int n = getOptimalDFTSize(InputImage.cols);
	//on the border add zero values
	cv::copyMakeBorder(InputImage, padded, 0, m - InputImage.rows, 0, n - InputImage.cols, BORDER_CONSTANT, Scalar::all(0));
	cv::Mat planes[] = { Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
	//使傅里叶变换中心化
	fftshift(planes[0], planes[0]);
	//创建一个复数图像F(u, v)
	cv::Mat complexF;
	cv::merge(planes, 2, complexF);
	//执行傅里叶变换
	cv::dft(complexF, complexF);

	cv::Mat complexFH;
	cv::Mat ifft;
	//采用对应像素相乘G(u, v) = H(u, v) * F(u, v)
	cv::multiply(complexF, H, complexFH);
	cv::idft(complexFH, ifft, DFT_REAL_OUTPUT);
	
	//最后不要忘了再完成一次移动
	fftshift(ifft, ifft);
	//为了显示效果将负数截断为0
	MinusToZero(ifft, ifft);
	normalize(ifft, ifft, 0, 1, NORM_MINMAX);
	ifft.convertTo(ifft, CV_8UC1, 255);
	//裁剪回原尺寸输出
	ifft(cv::Rect(0, 0, InputImage.cols, InputImage.rows)).copyTo(OutputImage);
}


//*********************************main()*******************************//

int main()
{
	string path = "blurring_test.tif";

	Mat SrcImage = imread(path, IMREAD_GRAYSCALE);
	if (!SrcImage.data) {
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}

	//高通滤波结果和高通核
	Mat H_outimg, H_kernel_H;
	//低通滤波结果和低通核
	Mat L_outimg, L_kernel_H;

	calcGLPF(SrcImage, 60, L_kernel_H);
	frequency_filter(SrcImage, L_kernel_H, L_outimg);

	calcGHPF(SrcImage, 80, H_kernel_H);
	frequency_filter(SrcImage, H_kernel_H, H_outimg);
	
	imwrite("low_blurring_test.png", L_outimg);
	imwrite("high_blurring_test.png", H_outimg);
	
	cv::waitKey(0);
	return 0;
}

4. 测试结果

原图：

高斯低通滤波器：

高斯高通滤波器：