码上实战立体匹配系列经典PatchMatch: 代价计算

Posted 2022-08-03 李迎松~

下载完整源码，点击进入: https://github.com/ethan-li-coding/PatchMatchStereo
欢迎同学们在Github项目里讨论，如果觉得博主代码质量不错，右上角star一下！感谢！

算法效果图镇楼：

上一篇博客随机初始化中，我们讲解了算法的第一个步骤：随机初始化 的实现代码，按照步骤，后面一步应该是迭代传播优化，但是我们不得不先讲另外一部分，它非常关键，贯穿于迭代优化步骤全程，它就是：PatchMatchStereo（PMS）的代价计算。

【码上实战】【立体匹配系列】经典PatchMatch: （4）代价计算

PMS代价计算的理论知识，我们早在理论恒叨的第一篇中有所介绍：

【理论恒叨】【立体匹配系列】经典PatchMatch: （1）Slanted support windows倾斜支持窗模型

细节再不多言。

PMS计算两个像素的匹配代价公式为：（$q$为左视图像素，$q^{\prime }$为右视图像素）

式中，$||I_p-I_q||$为像素$q$和$q^{\prime }$在RGB空间的L1-distance（就是$|r-r|+|g-g|+|b-b|$），
$||▽I_p-▽I_{q^{\prime }}||$为$q$和$q^{\prime }$的梯度之差的绝对值，该式的含义就是两个像素的颜色相差越大、梯度相差越大，则不相似的程度越大，匹配代价就越大，$\alpha$参数为自定义参数，在颜色和梯度之间做一个平衡。注意到根据平面算出来的视差值是小数值，所以像素$q^{\prime }$的位置不是整数而是小数，它的颜色值和梯度值通过线性内插得到。

如何计算左视图像素 $p$ 在视差平面为$f(a_f,b_f,c_f)$时的代价呢？

以$p$为中心，作一个适当大小的窗口，通过窗口内每个像素在视差平面为$f(a_f,b_f,c_f)$时的代价加权聚合得到：

式中，$W_p$是一个以$p$为中心的方形窗口，$w(p,q)$是自适应的权值，计算公式如下：

总结一下就是，像素$p$在视差平面为$f(a_f,b_f,c_f)$时的聚合代价，是把同样的视差平面赋给窗口内所有其他像素后得到的代价加权得到的。换句话说：最适合我的，也是适合你们的！

我们来看代码。

我把代价计算相关的代码都放到一个hpp文件中：cost_computer.hpp，文件名一看便知其意。我写了一个代价计算基类：

/**
 * \\brief 代价计算器基类
 */
class CostComputer 
public:
	/** \\brief 代价计算器默认构造 */
	CostComputer(): img_left_(nullptr), img_right_(nullptr), width_(0), height_(0), patch_size_(0), min_disp_(0),
	                max_disp_(0) 

以派生子类的方式实现具体的代价计算类，这样便于扩展其他不同的代价计算方法，而不用过多的修改频繁使用代价计算类的迭代传播函数体。

而代码里实现了PMS论文中的方法，也就是上面提到的代价计算方法，放到类CostComputerPMS 中：

/**
 * \\brief 代价计算器：PatchMatchStero原文代价计算器
 */
class CostComputerPMS : public CostComputer 
public:

	/** \\brief PMS代价计算器默认构造 */
	CostComputerPMS(): grad_left_(nullptr), grad_right_(nullptr), gamma_(0), alpha_(0), tau_col_(0), tau_grad_(0)  ;

我们重点介绍它。

一些简单的不必要重点介绍的代码本文就略过了，比如成员函数、构造函数之类的，重点介绍下它的成员函数。

首先是Compute函数，它就是我上面最先说到的计算两个像素之间的匹配代价。我们看代码：

/**
 * \\brief 计算左影像p点视差为d时的代价值，未做边界判定
 * \\param x		p点x坐标
 * \\param y		p点y坐标
 * \\param d		视差值
 * \\return 代价值
 */
inline float32 Compute(const sint32& x, const sint32& y, const float32& d) override

	const float32 xr = x - d;
	if (xr < 0.0f || xr >= static_cast<float32>(width_)) 
		return (1 - alpha_) * tau_col_ + alpha_ * tau_grad_;
	
	// 颜色空间距离
	const auto col_p = GetColor(img_left_, x, y);
	const auto col_q = GetColor(img_right_, xr, y);
	const auto dc = std::min(abs(col_p.b - col_q.x) + abs(col_p.g - col_q.y) + abs(col_p.r - col_q.z), tau_col_);

	// 梯度空间距离
	const auto grad_p = GetGradient(grad_left_, x, y);
	const auto grad_q = GetGradient(grad_right_, xr, y);
	const auto dg = std::min(abs(grad_p.x - grad_q.x)+ abs(grad_p.y - grad_q.y), tau_grad_);

	// 代价值
	return (1 - alpha_) * dc + alpha_ * dg;

函数输入看注释就明白，将像素的坐标和视差值作为输入。函数体代码前三行是右视图的$x$坐标计算（因为是核线像对所以$y$坐标和输入一致）以及边界处理。首先通过两个像素的位置获取像素的颜色值和梯度值，并计算两者的L1-distance，然后代价值等于两空间距离的加权和，权值alpha_是PMS的参数，一开始就被指定了。

以上只是计算了两个像素之间的匹配代价，再来看聚合代价的计算代码：

/**
 * \\brief 计算左影像p点视差平面为p时的聚合代价值
 * \\param x		p点x坐标
 * \\param y 	p点y坐标
 * \\param p		平面参数
 * \\return 聚合代价值
 */
inline float32 ComputeA(const sint32& x, const sint32& y, const DisparityPlane& p) const

	const auto pat = patch_size_ / 2;
	const auto& col_p = GetColor(img_left_, x, y);
	float32 cost = 0.0f;
	for (sint32 r = -pat; r <= pat; r++) 
		const sint32 yr = y + r;
		for (sint32 c = -pat; c <= pat; c++) 
			const sint32 xc = x + c;
			if (yr < 0 || yr > height_ - 1 || xc < 0 || xc > width_ - 1) 
				continue;
			
			// 计算视差值
			const float32 d = p.to_disparity(xc,yr);
			if (d < min_disp_ || d > max_disp_) 
				cost += COST_PUNISH;
				continue;
			

			// 计算权值
			const auto& col_q = GetColor(img_left_, xc, yr);
			const auto dc = abs(col_p.r - col_q.r) + abs(col_p.g - col_q.g) + abs(col_p.b - col_q.b);
#ifdef USE_FAST_EXP
			const auto w = fast_exp(double(-dc / gamma_));
#else
			const auto w = exp(-dc / gamma_);
#endif

			// 聚合代价
			const auto grad_q = GetGradient(grad_left_, xc, yr);
			cost += w * Compute(col_q, grad_q, xc, yr, d);
		
	
	return cost;

如前所述，聚合代价是以$(x,y)$为中心的窗口内所有像素在同一个视差平面$p$下的代价加权和，所以代码里对窗口内的所有像素进行遍历，计算（1）权值（2）每个像素在视差平面p下的视差d，及视差为d时的匹配代价。匹配代价计算是调用上面的代价计算函数（代码里用的是另外一个等价函数，为了避免重复获取在函数外先获取了左视图像素的颜色和梯度）。最后将每个像素计算的代价进行加权求和，就得到像素$(x,y)$的聚合代价。

代码里的USE_FAST_EXP是一种快速计算exp的算法，是为了给算法加速，和算法主题无关，不作多说。

另外，还有几个关于获取像素颜色和梯度数据的函数，他们是基于线性内插法，原理十分简单，这里也不做细说，大家看看代码即可。

/**
* \\brief 获取像素点的颜色值
* \\param img_data	颜色数组
* \\param x			像素x坐标，实数，线性内插得到颜色值
* \\param y			像素y坐标
* \\return 像素(x,y)的颜色值
*/
inline PVector3f GetColor(const uint8* img_data, const float32& x,const sint32& y) const
/**
* \\brief 获取像素点的梯度值
* \\param grad_data	梯度数组
* \\param x			像素x坐标，实数，线性内插得到梯度值
* \\param y			像素y坐标
* \\return 像素(x,y)的梯度值
*/
inline PVector2f GetGradient(const PGradient* grad_data, const float32& x, const sint32& y) const

好了，今天就到这里啦，我知道大家可能更想学习迭代传播，但还是下一篇再介绍吧，代价计算也是PMS的核心，要先掌握才能往后。放心，不会等太久！

下载完整源码，点击进入: https://github.com/ethan-li-coding/PatchMatchStereo.git
欢迎同学们在Github项目里讨论，如果觉得博主代码质量不错，给颗小星星，Follow 我！感激不尽！

码上实战系列

【码上实战】【立体匹配系列】经典PatchMatch: （1）框架
【码上实战】【立体匹配系列】经典PatchMatch: （2）主类
【码上实战】【立体匹配系列】经典PatchMatch: （3）随机初始化
【码上实战】【立体匹配系列】经典PatchMatch: （4）代价计算
【码上实战】【立体匹配系列】经典PatchMatch: （5）迭代传播
【码上实战】【立体匹配系列】经典PatchMatch: （6）后处理

博主简介：
Ethan Li 李迎松（知乎：李迎松）
武汉大学 摄影测量与遥感专业博士

主方向立体匹配、三维重建

2019年获测绘科技进步一等奖（省部级）

爱三维，爱分享，爱开源
GitHub： https://github.com/ethan-li-coding
邮箱：[e-mail]

个人微信：

欢迎交流！

关注博主不迷路，感谢！
博客主页：https://ethanli.blog.csdn.net