码上实战立体匹配系列经典PatchMatch: 迭代传播

Posted 2022-08-03 李迎松~

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

算法效果图镇楼：

上一篇博客代价计算中，我们讲解了PatchMatchStereo（PMS）的代价计算器，有了代价计算器，我们将在迭代传播步骤中频繁调用它计算某个像素的聚合代价。本篇，博主要介绍的是PMS最关键的核心：迭代传播。迭代传播有多重要呢？相比读过前面几篇的同学一定提前认识了未经迭代传播，也就是只做随机初始化的结果：

左影像

右影像

左视差图

右视差图

迭代传播可以让这个完全看不懂的麻点图变成无限接近于文章开头展示的结果。

来来来，代码拿来！

【码上实战】【立体匹配系列】经典PatchMatch: （5）迭代传播

• • DoPropagation
  • • SpatialPropagation
    • PlanePropagation
    • ViewPropagation
  • 3次迭代实验

我将迭代传播写到两个独立的文件里：pms_propagation.h和pms_propagation.cpp，这样结构会比较清晰，两个文件只实现了一个类：PMSPropagation，它将完成迭代传播的所有步骤，从前文可以得知步骤有三：（1）空间传播（2）视图传播（3）平面优化。

DoPropagation

我们先看看类PMSPropagation的唯一公有函数：

public:
	/** \\brief 执行传播一次 */
	void DoPropagation();

它是PMS主类调用的唯一接口，调用它一次就可以完成一次迭代，调用多次就完成多次迭代，PMS的原文中是推荐迭代3次。

DoPropagation方法里面的实现，大家或许能猜到是依次执行上面3个步骤，这3个步骤我放到私有函数里：

/**
* \\brief 空间传播
 * \\param x 像素x坐标
 * \\param y 像素y坐标
 * \\param direction 传播方向
 */
void SpatialPropagation(const sint32& x, const sint32& y, const sint32& direction) const;

/**
 * \\brief 视图传播
 * \\param x 像素x坐标
 * \\param y 像素y坐标
 */
void ViewPropagation(const sint32& x, const sint32& y) const;

/**
 * \\brief 平面优化
 * \\param x 像素x坐标
 * \\param y 像素y坐标
 */
void PlaneRefine(const sint32& x, const sint32& y) const;

看注释大家应该没有太多疑问，唯一需要多说一句的是传播是对于每个像素做的，也就是每个像素都要依次做完3个步骤才轮到下一个像素，这里博主也走了一点弯路，我开始是全图做完空间传播，再全图做视图传播，跑出来的结果打了我face。

我们来看具体实现。

首先，我们看看开放的唯一接口DoPropagation都干了些什么呢？

void PMSPropagation::DoPropagation()

	if(!cost_cpt_left_|| !cost_cpt_right_ || !img_left_||!img_right_||!grad_left_||!grad_right_ ||!cost_left_||!plane_left_||!plane_right_||!disparity_map_||
		!rand_disp_||!rand_norm_) 
		return;
	

	// 偶数次迭代从左上到右下传播
	// 奇数次迭代从右下到左上传播
	const sint32 dir = (num_iter_%2==0) ? 1 : -1;
	sint32 y = (dir == 1) ? 0 : height_ - 1;
	for (sint32 i = 0; i < height_; i++) 
		sint32 x = (dir == 1) ? 0 : width_ - 1;
		for (sint32 j = 0; j < width_; j++) 

			// 空间传播
			SpatialPropagation(x, y, dir);

			// 平面优化
			if (!option_.is_fource_fpw) 
				PlaneRefine(x, y);
			

			// 视图传播
			ViewPropagation(x, y);

			x += dir;
		
		y += dir;
	
	++num_iter_;

开始的指针检查咱们就不说了，从后面开始看，遍历每个像素进行迭代传播没什么多说的，这里我要重点讲解的有两点：

1. 传播方向。原文中说的比较清楚，偶数次迭代就从左上角像素传播到右下角像素，奇数次迭代就调个头从右下角像素传播到左上角像素。代码里num_iter_表示迭代序号，初始化为0，也就是偶数次迭代，每执行完一次迭代num_iter_就加1，下一次自然就变成了奇数。
2. 传播顺序。在原文中传播顺序应该是空间传播-视图传播-平面优化，这里我做了一个小调整，次序是空间传播-平面优化-视图传播。原因是第一次空间传播往往得不到很好的效果，按照以前的顺序，左视图做完第一次空间传播后再通过视图传播到右视图的结果并不会很好，效果会打折扣，而调整后的顺序，做完空间传播和平面优化，视差结果会变得更好，从而使接下来的视图传播更有质量。这样做还有另一个原因，我下面单独介绍视图传播的时候会讲解。

进入每个子步骤的实现吧。

SpatialPropagation

首先，来看空间传播SpatialPropagation。看看代码：

void PMSPropagation::SpatialPropagation(const sint32& x, const sint32& y, const sint32& direction) const

	// ---
	// 空间传播

	// 偶数次迭代从左上到右下传播
	// 奇数次迭代从右下到左上传播
	const sint32 dir = direction;

	// 获取p当前的视差平面并计算代价
	auto& plane_p = plane_left_[y * width_ + x];
	auto& cost_p = cost_left_[y * width_ + x];
	auto* cost_cpt = dynamic_cast<CostComputerPMS*>(cost_cpt_left_);

	// 获取p左(右)侧像素的视差平面，计算将平面分配给p时的代价，取较小值
	const sint32 xd = x - dir;
	if (xd >= 0 && xd < width_) 
		auto& plane = plane_left_[y * width_ + xd];
		if (plane != plane_p) 
			const auto cost = cost_cpt->ComputeA(x, y, plane);
			if (cost < cost_p) 
				plane_p = plane;
				cost_p = cost;
			
		
	

	// 获取p上(下)侧像素的视差平面，计算将平面分配给p时的代价，取较小值
	const sint32 yd = y - dir;
	if (yd >= 0 && yd < height_) 
		auto& plane = plane_left_[yd * width_ + x];
		if (plane != plane_p) 
			const auto cost = cost_cpt->ComputeA(x, y, plane);
			if (cost < cost_p) 
				plane_p = plane;
				cost_p = cost;
			
		
	

总的来说，函数体里有两步：

第一步，将左边（如果是反向，则为右边）像素的视差平面赋给当前像素，计算新的代价，判断新的代价是否比当前代价小，如果更小，则接受该视差平面为新的视差平面，代价也会更新。

第二步，将上边（如果是反向，则为下边）像素的视差平面赋给当前像素，计算新的代价，判断新的代价是否比当前代价小，如果更小，则接受该视差平面为新的视差平面，代价也会更新。

逻辑是比较清晰的，总共也没多少行代码，理解起来不难。

我们来看看只做一次空间传播，结果如何：（大家可以屏蔽掉其他两个步骤的代码来实验）

左影像

右影像

左视差图

右视差图

结果并不太好。如果我选择前端平行窗口模型呢？（把option的is_fource_fpw设置为true，即为Frontal-Parallel Window）

左影像

右影像

左视差图

右视差图

看起来好多了，这是为什么呢？同学们知道原因吗？我们可以就用Frontal-Parallel Window吗？大家可以在留言区讨论下。我后面会给答案大家。

PlanePropagation

平面优化显得要复杂一些，但是理解起来也不难。先看代码：

void PMSPropagation::PlaneRefine(const sint32& x, const sint32& y) const

	// --
	// 平面优化
	const auto max_disp = static_cast<float32>(option_.max_disparity);
	const auto min_disp = static_cast<float32>(option_.min_disparity);

	// 随机数生成器
	std::random_device rd;
	std::mt19937 gen(rd());
	const auto& rand_d = *rand_disp_;
	const auto& rand_n= *rand_norm_;

	// 像素p的平面、代价、视差、法线
	auto& plane_p = plane_left_[y * width_ + x];
	auto& cost_p = cost_left_[y * width_ + x];
	auto* cost_cpt = dynamic_cast<CostComputerPMS*>(cost_cpt_left_);

	float32 d_p = plane_p.to_disparity(x, y);
	PVector3f norm_p = plane_p.to_normal();

	float32 disp_update = (max_disp - min_disp) / 2.0f;
	float32 norm_update = 1.0f;
	const float32 stop_thres = 0.1f;

	// 迭代优化
	while (disp_update > stop_thres) 

		// 在 -disp_update ~ disp_update 范围内随机一个视差增量
		float32 disp_rd = rand_d(gen) * disp_update;
		if (option_.is_integer_disp) 
			disp_rd = static_cast<float32>(round(disp_rd));
		

		// 计算像素p新的视差
		const float32 d_p_new = d_p + disp_rd;
		if (d_p_new < min_disp || d_p_new > max_disp) 
			disp_update /= 2;
			norm_update /= 2;
			continue;
		

		// 在 -norm_update ~ norm_update 范围内随机三个值作为法线增量的三个分量
		PVector3f norm_rd;
		if (!option_.is_fource_fpw) 
			norm_rd.x = rand_n(gen) * norm_update;
			norm_rd.y = rand_n(gen) * norm_update;
			float32 z = rand_n(gen) * norm_update;
			while (z == 0.0f) 
				z = rand_n(gen) * norm_update;
			
			norm_rd.z = z;
		
		else 
			norm_rd.x = 0.0f; norm_rd.y = 0.0f;	norm_rd.z = 0.0f;
		

		// 计算像素p新的法线
		auto norm_p_new = norm_p + norm_rd;
		norm_p_new.normalize();

		// 计算新的视差平面
		auto plane_new = DisparityPlane(x, y, norm_p_new, d_p_new);

		// 比较Cost
		if (plane_new != plane_p) 
			const float32 cost = cost_cpt->ComputeA(x, y, plane_new);

			if (cost < cost_p) 
				plane_p = plane_new;
				cost_p = cost;
				d_p = d_p_new;
				norm_p = norm_p_new;
			
		

		disp_update /= 2.0f;
		norm_update /= 2.0f;
	

一开始，我获取了两个随机数生成器，一个用来生成视差的随机值，一个用来生成法线的随机值（随机数生成器的初始化是放在PMSPropagation类的构造函数里的）。并获取了像素$p$的平面、代价、视差和法线以做后用。

平面传播的原理，有必要再介绍下：

PMS将$f_p$转换为点加法向量的表达方式，并设置两个参数：${\Delta }_{z_0}^{max}$和${\Delta }_n^{max}$。${\Delta }_{z_0}^{max}$为点$P(x_0,y_0,z_0)$的z-坐标的可变化范围，${\Delta }_n^{max}$为法向量$\vec{n}$各分量的可变化范围。在$[-{\Delta }_{z_0}^{max},{\Delta }_{z_0}^{max}]$范围内随机一个值${\Delta }_{z_0}$加到$z_0$上得$z_0^{\prime }=z_0+{\Delta }_{z_0}$，由此得到新的点$P^{\prime }(x_0,y_0,z以上是关于码上实战立体匹配系列经典PatchMatch:\; 迭代传播的主要内容，如果未能解决你的问题，请参考以下文章$

码上实战立体匹配系列经典PatchMatch: 迭代传播

码上实战立体匹配系列经典PatchMatch: 后处理

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

码上实战立体匹配系列经典PatchMatch: 随机初始化

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

码上实战立体匹配系列经典AD-Census: 多步骤视差优化