密集比對之動态規劃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);
}
           

本人水準有限，如有錯誤，還望不吝指正，代碼有一定删減，沒有重新編譯，如有錯誤，請自行調試，有問題請郵箱聯系[email protected]。