密集匹配之动态规划DP

       动态规划通常应用于最优化问题的求解中，Baker、Ohta 等将动态规划引入立体匹配中来获取视差图。动态规划匹配的过程可以分为两个阶段，建立视差空间和动态规划优化，将立体匹配问题转化为视差空间内寻找最优路径的问题。

       密集匹配通常会充分利用影像间的核线约束条件，对立体像对进行核线纠正，这样同名像点肯定位于对应的同名核线上，降低了匹配的难度。视差空间影像DSI（Disparity Space Image）是一种反映核线立体像对间同一扫描线视差关系的数据结构，将视差分析转化为一种类似的谱分析。其数学表达为：

其中w 为图像的宽度， NaN 表示在该视差值为d 时不可能有对应的同名点，L((* ))  是对应点之间的相似性测度函数。在匹配核线影像时，逐行计算像素的DSI，将每一行的DSI 依次叠加起来构成一个三维空间，称之为视差空间，视差空间包含了每个像素所有可能的视差取值。如下图所示，图中的x、y 代表图像的横纵坐标，d 表示视差搜索范围。

在建立好视差空间后，立体匹配问题就会转化为在视差空间内寻找最佳路径的问题，该路径上的视差累积的能量最小，相似性最大。一般采用从图像左边到右边的顺序，逐行进行。

实现过程如下：

void DynamicProgramMatch(const cv::Mat& img1, const cv::Mat& img2,
	int min_disp, int max_disp, int P1, int P2,
	int m_size, cv::Mat& disp)
{
	if (img1.empty() || img1.type() != CV_8UC1
		|| img2.empty() || img2.type() != CV_8UC1
		|| img1.size() != img2.size() || m_size < 1) return;

	int disp_range = max_disp - min_disp + 1;///视差范围
	disp = cv::Mat(img1.size(), CV_8UC1, cv::Scalar::all(0));///左视差图
	for (int r = 0; r < img1.rows; ++r)
	{
		printf("动态规划...%d%%\r", (int)(r * 100.0 / (img1.rows))); // 方便查看进度
		
		cv::Mat diff(cv::Size(disp_range, img1.cols), CV_64FC1, cv::Scalar::all(0)); //视差空间
		cv::Mat p_mat(cv::Size(disp_range, img1.cols), CV_8UC1, cv::Scalar::all(0)); // 路径走向
		for (int c1 = 0; c1 < img1.cols; ++c1)
		{
			for (int d = 0; d <= disp_range; ++d)
			{
				int c2 = c1 - d - min_disp;///当前视差下左像素对应对的右像素坐标
				if (c2 >= 0)
				{
					diff.ptr<double>(c1)[d] =
						SAD(img1, c1, r, img2, c2, r, m_size);//SAD匹配代价
				}
				else
				{
					diff.ptr<double>(c1)[d] =
						c1 < min_disp ? 0 : diff.ptr<double>(c1)[d - 1]; // 防止边界超限处代价对整体路径的影响
				}
				
			}
		}

		cv::Mat e_mat_cur(cv::Size(1, disp_range), CV_64FC1, cv::Scalar::all(0));
		for (int c = 1; c < img1.cols; ++c)
		{
			cv::Mat e_mat_pre = e_mat_cur.clone(); // 能量空间
			for (int cur = 0; cur < disp_range; ++cur)
			{
				double cost_min = FLT_MAX;//当前状态下路径代价
				int p_min = 0; //路径走向
				double e_cur =  diff.ptr<double>(c)[cur]; //数据项
				for (int pre = 0; pre < disp_range; pre++)
				{
					int deta_d = abs(cur - pre);
					double e_smooth = deta_d > 1 ? P2 : (deta_d == 0 ? 0 : P1); //视差等于1与大于1时候的平滑惩罚不同
					double e_data = e_cur + e_smooth + e_mat_pre.ptr<double>(pre)[0]; // 路径能量值
					if (e_data < cost_min)
					{
						cost_min = e_data;
						p_min = pre;
					}
				}
				e_mat_cur.ptr<double>(cur)[0] = cost_min;
				p_mat.ptr<uchar>(c)[cur] = p_min;
			}
		}

		int p_min = 0;
		double e_min = e_mat_cur.ptr<double>(0)[0];
		for (int i = 0; i < disp_range; ++i)
		{
			if (e_mat_cur.ptr<double>(i)[0] < e_min)
			{
				p_min = i;
				e_min = e_mat_cur.ptr<double>(i)[0];
			}
		}
		//取得视差
		for (int c = img1.cols - 1; c >= 0; --c)
		{
			int d = p_mat.ptr<uchar>(c)[p_min];
			p_min = d;
			disp.ptr<uchar>(r)[c] = d + min_disp; 
		}
	}
}

实现样例如下：

void DynamicProgrammingTest()
{
	cv::Mat img1 = cv::imread("cones\\im2.png");
	cv::Mat img2 = cv::imread("cones\\im6.png");
	cv::imshow("img1", img1);
	cv::imshow("img2", img2);
	cv::Mat gray1, gray2;
	cv::cvtColor(img1, gray1, cv::COLOR_BGR2GRAY);
	cv::cvtColor(img2, gray2, cv::COLOR_BGR2GRAY);
	cv::Mat disp;
	dense_matching::DynamicProgramMatch(gray1, gray2, 0, 100, 1, 4, 3, disp);
	cv::Mat disp_show;
	cv::normalize(disp, disp_show, 0, 255, cv::NORM_MINMAX);
	cv::imshow("disp", disp_show);
	cv::waitKey(0);
}

本人水平有限，如有错误，还望不吝指正，代码有一定删减，没有重新编译，如有错误，请自行调试，有问题请邮箱联系1299771369@qq.com。