OpenCV实现同态滤波

同态滤波是属于图像增强的一个小算法，其原理和代码实现在众多博客中均有提及，再此，只对学习中一些自认为有用的知识点进行总结。

实现和学习过程中的一些总结：

• 同态滤波类似于灰度变换，都是对灰度级进行动态范围压缩和对比度增强来改善图像的外观。
• 图像$f(x,y)$可以表示为其照射分量$i(x,y)$和反射分量$r(x,y)$的乘积，即：$$f(x,y)=i(x,y)r(x,y)$$
• 照射分量对应图像中的低频部分，在光源中来说对应漫反射光，反射分量对应图像中高频部分，在光源中对于镜面反射光。
• 在算法开始，先要对图像进行对数处理，若$f(x,y)$中有任何零值，则必须将图像加1，以避免$In(0)$的出现，然后从最终结果中减去1。

一些未出现的工具函数可以参考我之前的博客：图像基础频率域滤波

---

代码实现如下：

//func:		创建滤波传递函数
//size:		扩展后图像的大小
//gamma_L:	控制低频下限
//gamma_H:	控制高频上限
//c,D0:		控制滤波函数剖面的偏斜度
cv::Mat calc_Homomorphic(cv::Size size, double gamma_L, double gamma_H, double c, double D0)
{
	cv::Mat result(size, CV_64FC1);

	int cx = size.width / 2;
	int cy = size.height / 2;

	for (int i = 0; i < result.rows; ++i)
	{
		double* p = result.ptr<double>(i);
		for (int j = 0; j < result.cols; ++j)
		{
			double d_2 = std::pow(i - cy, 2) + std::pow(j - cx, 2);
			//同态滤波传递函数
			p[j] = (gamma_H - gamma_L) * (1 - std::exp(-c * d_2 / (D0 * D0))) + gamma_L;
		}
	}

	//创建双通道图像，便于与复数图像相乘
	cv::Mat planes[] = { result.clone(), result.clone() };
	cv::Mat Homomorphic;
	merge(planes, 2, Homomorphic);

	return Homomorphic;
}


//************************main入口*************************//

int main()
{
	std::string path = "PET_image.tif";

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

	//源图像转化为对数形式
	cv::Mat srclogMat;
	src.convertTo(srclogMat, CV_64FC1);
	srclogMat = srclogMat + cv::Scalar::all(1);
	cv::log(srclogMat, srclogMat);

	//扩展图像，扩展部分使用零填充
	cv::Mat paddedMat;
	int m = cv::getOptimalDFTSize(src.rows);
	int n = cv::getOptimalDFTSize(src.cols);
	//on the border add zero values
	cv::copyMakeBorder(srclogMat, paddedMat, 0, m - src.rows, 0, n - src.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
	cv::Mat paddedMatdouble;
	paddedMat.convertTo(paddedMatdouble, CV_64FC1);
	//使傅里叶变换中心化
	Dftshift(paddedMatdouble, paddedMatdouble);

	//创建一个复数图像F(u, v)
	cv::Mat complexF;
	cv::dft(paddedMatdouble, complexF, cv::DFT_COMPLEX_OUTPUT);

	//创建同态滤波传递函数
	double gamma_L = 0.25, gamma_H = 2.0;
	double c = 1.0, D0 = 160;

	cv::Mat Homomorphic = calc_Homomorphic(paddedMat.size(), gamma_L, gamma_H, c, D0);

	cv::Mat complexFH;
	cv::Mat iDft;
	//采用对应像素相乘G(u, v) = H(u, v) * F(u, v)
	cv::multiply(complexF, Homomorphic, complexFH);
	cv::idft(complexFH, iDft, cv::DFT_REAL_OUTPUT | cv::DFT_SCALE);
	//最后不要忘了再完成一次移动
	Dftshift(iDft, iDft);

	//之前进行了对数变换，现在要变换回来
	cv::exp(iDft, iDft);
	iDft = iDft - cv::Scalar::all(1);
	cv::Mat result = iDft(cv::Rect(0, 0, src.cols, src.rows)).clone();
	//截断负值，归一化并转换为uchar
	MinusToZero(result, result);
	normalize(result, result, 0, 1, cv::NORM_MINMAX);
	result.convertTo(result, CV_8UC1, 255);

	cv::imwrite("ascend_ET.png", result);
	cv::waitKey(0);
	return 0;
}

测试结果：