【碼上實戰】【立體匹配系列】經典AD-Census: （4）十字交叉域代價聚合

Ethan Li 李迎鬆發表於2020-09-30

抱歉讓大家久等了！

下載完整原始碼，點選進入: https://github.com/ethan-li-coding/AD-Census
歡迎同學們在Github專案裡討論！

上篇代價計算中，帶大家深入瞭解了ADCensusStereo代價計算步驟的實現程式碼，並做了實驗展示了演算法的結果。這裡先貼一下代價計算的實驗結果圖：

Census

AD-Census

正如上篇結尾所述，AD-Census的代價計算方法結合了AD和Census匹配策略的優點，是一個優秀的策略。

而無論如何，代價計算的結果都無法作為最終的視差結果，後續的優化步驟對ADCensus來說都是不可少的步驟。本篇即介紹第一個優化步驟：基於十字交叉域的代價聚合 的程式碼實現。

文章目錄

演算法

演算法原理請看博文：

經典AD-Census: （2）十字交叉域代價聚合（Cross-based Cost Aggregation）

不再贅述，以免喧賓奪主。

程式碼實現

類設計

成員函式

我將十字交叉域代價聚合寫成一個類CrossAggregator，放在檔案cross_aggregator.h/cross_aggregator.cpp中。

/**
 * \brief 十字交叉域代價聚合器
 */
class CrossAggregator {
public:
	CrossAggregator();
	~CrossAggregator();
}

我們需要為CrossAggregator類設計一些公有介面，從而可以呼叫它們來達到代價聚合的目的。

我們沿用代價計算器設計的思路。

第一個介面是 初始化函式Initialize ，為代價聚合做一些準備工作，例如很重要的一點我需要開闢一塊記憶體用來存放每個畫素的十字交叉臂資料，以便後面進行聚合。

第二類介面是必不可少的 設定資料SetData 以及 設定引數SetParam ，脫離資料和引數，演算法啥也不是。

第三個就是關鍵介面 Aggregate 了，最終還是要靠呼叫它來完成代價聚合。

最後兩個額外的介面是 get_arms_ptr 和 get_cost_ptr，前者獲取十字交叉臂陣列，這純粹是因為後處理步驟需要；後者是獲取聚合後的代價資料指標，它很重要，因為後面的掃描線優化步驟需要它。

我們來看看最終的類介面設計：

/**
 * \brief 初始化代價聚合器
 * \param width		影像寬
 * \param height	影像高
 * \return true:初始化成功
 */
bool Initialize(const sint32& width, const sint32& height, const sint32& min_disparity, const sint32& max_disparity);

/**
 * \brief 設定代價聚合器的資料
 * \param img_left		// 左影像資料，三通道
 * \param img_right		// 右影像資料，三通道
 * \param cost_init		// 初始代價陣列
 */
void SetData(const uint8* img_left, const uint8* img_right, const float32* cost_init);

/**
 * \brief 設定代價聚合器的引數
 * \param cross_L1		// L1
 * \param cross_L2		// L2
 * \param cross_t1		// t1
 * \param cross_t2		// t2
 */
void SetParams(const sint32& cross_L1, const sint32& cross_L2, const sint32& cross_t1, const sint32& cross_t2);

/** \brief 聚合 */
void Aggregate(const sint32& num_iters);

/** \brief 獲取所有畫素的十字交叉臂資料指標 */
CrossArm* get_arms_ptr();

/** \brief 獲取聚合代價陣列指標 */
float32* get_cost_ptr();

如上介面為類的全部公有函式，呼叫者可按邏輯次序呼叫。具體介面的引數含義，註釋寫的很清楚，大家對著看就行了。

需要提到的一點是，get_arms_ptr函式範圍的型別CrossArm，它是我自定義的十字交叉臂結構體，儲存了畫素的上下左右臂的長度。如下所示：

/**
* \brief 交叉十字臂結構
* 為了限制記憶體佔用，臂長型別設定為uint8，這意味著臂長最長不能超過255
*/
struct CrossArm {
	uint8 left, right, top, bottom;
	CrossArm() : left(0), right(0), top(0), bottom(0) { }
};

為了完成代價聚合，我們需要一些子步驟，它們包括：

構建十字交叉臂
計算畫素的支援區畫素數量
聚合所有畫素在某固定視差下的代價

它們被設計為類的私有成員函式：

/** \brief 構建十字交叉臂 */
void BuildArms();
/** \brief 搜尋水平臂 */
void FindHorizontalArm(const sint32& x, const sint32& y, uint8& left, uint8& right) const;
/** \brief 搜尋豎直臂 */
void FindVerticalArm(const sint32& x, const sint32& y, uint8& top, uint8& bottom) const;
/** \brief 計算畫素的支援區畫素數量 */
void ComputeSupPixelCount();
/** \brief 聚合某個視差 */
void AggregateInArms(const sint32& disparity, const bool& horizontal_first);

在聚合函式 Aggregate 中，這些私有函式被依次呼叫從而完成代價聚合。

成員變數

最後是成員函式類不可或缺的成員變數，比如影像尺寸、交叉臂資料（左影像）、影像資料、代價資料、演算法引數等，正是它們在類的作用域內始終保持著演算法計算所需要的值，才能達到最後的計算目的。得用私有型別把它們保護起來，僅限類的內部使用，我們來看看吧！

/** \brief 影像尺寸 */
sint32	width_;
sint32	height_;

/** \brief 交叉臂 */
vector<CrossArm> vec_cross_arms_;

/** \brief 影像資料 */
const uint8* img_left_;
const uint8* img_right_;

/** \brief 初始代價陣列指標 */
const float32* cost_init_;
/** \brief 聚合代價陣列 */
vector<float32> cost_aggr_;

/** \brief 臨時代價資料 */
vector<float32> vec_cost_tmp_[2];
/** \brief 支援區畫素數量陣列 0：水平臂優先 1：豎直臂優先 */
vector<uint16> vec_sup_count_[2];
vector<uint16> vec_sup_count_tmp_;

sint32	cross_L1_;			// 十字交叉視窗的空間域引數：L1
sint32  cross_L2_;			// 十字交叉視窗的空間域引數：L2
sint32	cross_t1_;			// 十字交叉視窗的顏色域引數：t1
sint32  cross_t2_;			// 十字交叉視窗的顏色域引數：t2
sint32  min_disparity_;			// 最小視差
sint32	max_disparity_;			// 最大視差

/** \brief 是否成功初始化標誌	*/
bool is_initialized_;

說明一下，在成員變數中，我們的初始代價資料cost_init_是指標型，因為是直接使用的代價計算器返回的資料，而聚合代價資料cost_aggr_是vector容器，它將在本類的初始化函式中開闢記憶體，並存在於代價聚合器的整個生命週期，將通過介面get_cost_ptr輸出其地址，給後續步驟使用。

類實現

我們先看看代價聚合器類CrossAggregator的初始化函式實現：

bool CrossAggregator::Initialize(const sint32& width, const sint32& height, const sint32& min_disparity, const sint32& max_disparity)
{
	width_ = width;
	height_ = height;
	min_disparity_ = min_disparity;
	max_disparity_ = max_disparity;
	
	const sint32 img_size = width_ * height_;
	const sint32 disp_range = max_disparity_ - min_disparity_;
	if (img_size <= 0 || disp_range <= 0) {
		is_initialized_ = false;
		return is_initialized_;
	}

	// 為交叉十字臂陣列分配記憶體
	vec_cross_arms_.clear();
	vec_cross_arms_.resize(img_size);

	// 為臨時代價陣列分配記憶體
	vec_cost_tmp_[0].clear();
	vec_cost_tmp_[0].resize(img_size);
	vec_cost_tmp_[1].clear();
	vec_cost_tmp_[1].resize(img_size);

	// 為儲存每個畫素支援區畫素數量的陣列分配記憶體
	vec_sup_count_[0].clear();
	vec_sup_count_[0].resize(img_size);
	vec_sup_count_[1].clear();
	vec_sup_count_[1].resize(img_size);
	vec_sup_count_tmp_.clear();
	vec_sup_count_tmp_.resize(img_size);

	// 為聚合代價陣列分配記憶體
	cost_aggr_.resize(img_size * disp_range);

	is_initialized_ = !vec_cross_arms_.empty() && !vec_cost_tmp_[0].empty() && !vec_cost_tmp_[1].empty() 
					&& !vec_sup_count_[0].empty() && !vec_sup_count_[1].empty() 
					&& !vec_sup_count_tmp_.empty() && !cost_aggr_.empty();
	return is_initialized_;
}

前面變數賦值不必多說，初始化主要的任務是為交叉十字臂陣列分配記憶體，為臨時代價陣列分配記憶體，以及為儲存每個畫素支援區畫素數量的陣列分配記憶體。（如果讀懂了演算法理論，那麼這些操作應該不難理解，都是為了最後計算聚合後的代價）。此外，聚合代價陣列的記憶體也在此被開闢，並由本類終生維護。（自己計算出的結果自己來維護，合情合理！）

初始化之後，SetData和SetParam函式都比較好理解，我也就不佔用篇幅了。

當以上準備工作做好後，我們就開始真正實現代價聚合的功能。

從原理可知，要完成代價聚合，我們必須經過兩個步驟：

構建聚合區域，也就是構建十字交叉臂
執行聚合

我們知道對於每個畫素，都有兩個方向的臂，一個水平方向臂，一個豎直方向臂，而構建臂的約束是一模一樣的，分為距離約束和顏色約束，簡單來說是距離不能超過設定閾值，顏色差也不能超過設定閾值。這麼做是為了讓畫素聚合區域內的畫素都和其視差相近。AD-Census雖然不是簡單的在臂延伸的過程中計算畫素與中心畫素的距離和色差，但是也不復雜，只是多計算一個延伸方向上與上一個畫素的距離和色差，讓構造更加穩健一些。具體實現原理，請看博文：

經典AD-Census: （2）十字交叉域代價聚合

我寫了兩個函式分別構造水平臂和豎直臂，其實除了方向不一樣，構造方法及引數都是一模一樣的。這裡只貼出畫素 $(x, y)$ 水平臂的構造程式碼：

void CrossAggregator::FindHorizontalArm(const sint32& x, const sint32& y, uint8& left, uint8& right) const
{
	// 畫素資料地址
	const auto img0 = img_left_ + y * width_ * 3 + 3 * x;
	// 畫素顏色值
	const ADColor color0(img0[0], img0[1], img0[2]);
	
	left = right = 0;
	//計算左右臂,先左臂後右臂
	sint32 dir = -1;
	for (sint32 k = 0; k < 2; k++) {
		// 延伸臂直到條件不滿足
		// 臂長不得超過cross_L1
		auto img = img0 + dir * 3;
		auto color_last = color0;
		sint32 xn = x + dir;
		for (sint32 n = 0; n < std::min(cross_L1_, MAX_ARM_LENGTH); n++) {

			// 邊界處理
			if (k == 0) {
				if (xn < 0) {
					break;
				}
			}
			else {
				if (xn == width_) {
					break;
				}
			}

			// 獲取顏色值
			const ADColor color(img[0], img[1], img[2]);

			// 顏色距離1（臂上畫素和計算畫素的顏色距離）
			const sint32 color_dist1 = ColorDist(color, color0);
			if (color_dist1 >= cross_t1_) {
				break;
			}

			// 顏色距離2（臂上畫素和前一個畫素的顏色距離）
			if (n > 0) {
				const sint32 color_dist2 = ColorDist(color, color_last);
				if (color_dist2 >= cross_t1_) {
					break;
				}
			}

			// 臂長大於L2後，顏色距離閾值減小為t2
			if (n + 1 > cross_L2_) {
				if (color_dist1 >= cross_t2_) {
					break;
				}
			}

			if (k == 0) {
				left++;
			}
			else {
				right++;
			}
			color_last = color;
			xn += dir;
			img += dir * 3;
		}
		dir = -dir;
	}
}

搞清楚是哪兩個畫素參與距離計算就差不多弄懂了。其一是在臂延伸過程中新的畫素與中心畫素的距離；其二是新的畫素與其前一個畫素的距離。

豎直臂的構造方法和水平臂一樣，只是方向不同而已。

接下來，我們需要計算出每個畫素的支援區（十字交叉臂覆蓋的區域）畫素的數量。為什麼要有這一步呢？原因是我們在原理部分了解到，代價聚合是將支援區的代價值累加，再除以支援區的畫素數量，也就是計算支援區代價的均值，賦給中心畫素的代價值。

同時，我們必須注意到的是，不同的聚合方向，支援區並不相同，即先水平再豎直的聚合方向和先豎直再水平的聚合方向，兩者的支援區是不同的。所以需要分別用兩種方向的陣列來儲存支援區的畫素數。我們來看程式碼：

void CrossAggregator::ComputeSupPixelCount()
{
	// 計算每個畫素的支援區畫素數量
	// 注意：兩種不同的聚合方向，畫素的支援區畫素是不同的，需要分開計算
	bool horizontal_first = true;
	for (sint32 n = 0; n < 2; n++) {
		// n=0 : horizontal_first; n=1 : vertical_first
		const sint32 id = horizontal_first ? 0 : 1;
		for (sint32 k = 0; k < 2; k++) {
			// k=0 : pass1; k=1 : pass2
			for (sint32 y = 0; y < height_; y++) {
				for (sint32 x = 0; x < width_; x++) {
					// 獲取arm數值
					auto& arm = vec_cross_arms_[y*width_ + x];
					sint32 count = 0;
					if (horizontal_first) {
						if (k == 0) {
							// horizontal
							for (sint32 t = -arm.left; t <= arm.right; t++) {
								count++;
							}
						}
						else {
							// vertical
							for (sint32 t = -arm.top; t <= arm.bottom; t++) {
								count += vec_sup_count_tmp_[(y + t)*width_ + x];
							}
						}
					}
					else {
						if (k == 0) {
							// vertical
							for (sint32 t = -arm.top; t <= arm.bottom; t++) {
								count++;
							}
						}
						else {
							// horizontal
							for (sint32 t = -arm.left; t <= arm.right; t++) {
								count += vec_sup_count_tmp_[y*width_ + x + t];
							}
						}
					}
					if (k == 0) {
						vec_sup_count_tmp_[y*width_ + x] = count;
					}
					else {
						vec_sup_count_[id][y*width_ + x] = count;
					}
				}
			}
		}
		horizontal_first = !horizontal_first;
	}
}

在計算的過程中，我們通過模擬構造支援區的過程來計算支援區的畫素數，即像原文中所述，先統計畫素在一個方向的支援區畫素數量，儲存於臨時陣列裡，再沿另一個方向的臂累加臨時陣列裡的儲存數值。由於兩個構造方向的支援區不同，所以vec_sup_count_是一個容量為2的陣列，分別儲存著兩個構造方向每個畫素支援區數量。

最後，我們可以進行迭代聚合了，迭代次數是作為引數傳入的，我們用一個for迴圈來完成迭代，每執行一次迭代，我們會轉換下支援區構造順序（偶數次迭代是先水平再豎直構造，奇數次是先豎直再水平構造）。之後演算法將對候選視差下的每一個視差 $d$ 進行獨立的全圖聚合（視差 $d_1$ 與視差 $d_2$ 之間的代價聚合是完全獨立的，是可以高度並行的）。

我們先來看單視差 $d$ 下的聚合程式碼實現：

void CrossAggregator::AggregateInArms(const sint32& disparity, const bool& horizontal_first)
{
	// 此函式聚合所有畫素當視差為disparity時的代價

	if (disparity < min_disparity_ || disparity >= max_disparity_) {
		return;
	}
	const auto disp = disparity - min_disparity_;
	const sint32 disp_range = max_disparity_ - min_disparity_;
	if (disp_range <= 0) {
		return;
	}

	// 將disp層的代價存入臨時陣列vec_cost_tmp_[0]
	// 這樣可以避免過多的訪問更大的cost_aggr_,提高訪問效率
	for (sint32 y = 0; y < height_; y++) {
		for (sint32 x = 0; x < width_; x++) {
			vec_cost_tmp_[0][y * width_ + x] = cost_aggr_[y * width_ * disp_range + x * disp_range + disp];
		}
	}

	// 逐畫素聚合
	const sint32 ct_id = horizontal_first ? 0 : 1;
	for (sint32 k = 0; k < 2; k++) {
		// k==0: pass1
		// k==1: pass2
		for (sint32 y = 0; y < height_; y++) {
			for (sint32 x = 0; x < width_; x++) {
				// 獲取arm數值
				auto& arm = vec_cross_arms_[y*width_ + x];
				// 聚合
				float32 cost = 0.0f;
				if (horizontal_first) {
					if (k == 0) {
						// horizontal
						for (sint32 t = -arm.left; t <= arm.right; t++) {
							cost += vec_cost_tmp_[0][y * width_ + x + t];
						}
					} else {
						// vertical
						for (sint32 t = -arm.top; t <= arm.bottom; t++) {
							cost += vec_cost_tmp_[1][(y + t)*width_ + x];
						}
					}
				}
				else {
					if (k == 0) {
						// vertical
						for (sint32 t = -arm.top; t <= arm.bottom; t++) {
							cost += vec_cost_tmp_[0][(y + t) * width_ + x];
						}
					} else {
						// horizontal
						for (sint32 t = -arm.left; t <= arm.right; t++) {
							cost += vec_cost_tmp_[1][y*width_ + x + t];
						}
					}
				}
				if (k == 0) {
					vec_cost_tmp_[1][y*width_ + x] = cost;
				}
				else {
					cost_aggr_[y*width_*disp_range + x*disp_range + disp] = cost / vec_sup_count_[ct_id][y*width_ + x];
				}
			}
		}
	}
}

程式碼按照原文的思路，先對每個畫素，將某一方向的代價值累加到臨時儲存陣列；再對每個畫素沿另一方向的臂展對臨時儲存值累加，最後累加值除以支援區數量，得到最終的代價值。

為了加快速度，在聚合之前，我們將視差為 $d$ 時的初始代價，轉存至一個臨時的二維陣列，這樣聚合過程中無需頻繁訪問三維的代價陣列，減少訪問耗時。

以上我們就將聚合的子步驟都實現完畢。剩下的工作就簡單了，在類的公有聚合介面 Aggregate 中，我們只需要一次呼叫以上子步驟完成聚合即可。程式碼如下：

void CrossAggregator::Aggregate(const sint32& num_iters)
{
	if (!is_initialized_) {
		return;
	}

	const sint32 disp_range = max_disparity_ - min_disparity_;

	// 構建畫素的十字交叉臂
	BuildArms();

	// 代價聚合
	// horizontal_first 代表先水平方向聚合
	bool horizontal_first = true;

	// 計算兩種聚合方向的各畫素支援區畫素數量
	ComputeSupPixelCount();

	// 先將聚合代價初始化為初始代價
	memcpy(&cost_aggr_[0], cost_init_, width_*height_*disp_range*sizeof(float32));

	// 多迭代聚合
	for (sint32 k = 0; k < num_iters; k++) {
		for (sint32 d = min_disparity_; d < max_disparity_; d++) {
			AggregateInArms(d, horizontal_first);
		}
		// 下一次迭代，調換順序
		horizontal_first = !horizontal_first;
	}
}

函式的引數是聚合的迭代次數，原文推薦為4。

最後在主類ADCensusStereo的代價聚合函式CostAggregation中，我們呼叫代價聚合器CrossAggregator的公有成員函式完成代價聚合：

void ADCensusStereo::CostAggregation()
{
	// 設定聚合器資料
	aggregator_.SetData(img_left_, img_right_, cost_computer_.get_cost_ptr());
	// 設定聚合器引數
	aggregator_.SetParams(option_.cross_L1, option_.cross_L2, option_.cross_t1, option_.cross_t2);
	// 代價聚合
	aggregator_.Aggregate(4);
}

從繁至簡，也是C++的特性之一。

實驗

按照慣例，我們會提供一些實驗結果。

和上篇代價計算一樣，我們依然採用Cone資料。

左檢視

右檢視

代價計算的結果，我在開頭已經貼出，代價聚合的實驗結果如下：

代價計算

代價聚合

一個字，妙！通過代價聚合，視差圖質量得到了顯著的提升！

好了今天我們就到此結束，下一篇為大家帶來的是 掃描線優化。拜拜！

下載AD-Census完整原始碼，點選進入: https://github.com/ethan-li-coding/AD-Census
歡迎同學們在Github專案裡討論，如果覺得博主程式碼質量不錯，右上角給顆星！感謝！

博主簡介：
Ethan Li 李迎鬆（知乎：李迎鬆）
武漢大學攝影測量與遙感專業博士
主方向立體匹配、三維重建
2019年獲測繪科技進步一等獎（省部級）

愛三維，愛分享，愛開源
GitHub： https://github.com/ethan-li-coding （歡迎follow和star）

個人微信：

歡迎交流！

關注博主不迷路，感謝！
部落格主頁：https://ethanli.blog.csdn.net

AAAI 2020 開源論文 | 用於深度立體匹配的自適應單峰匹配代價體濾波
2020-02-20
AI
4個Python經典專案實戰，練手必備哦
2019-09-10
Python
DL4J實戰之四：經典卷積例項(GPU版本)
2021-10-15
卷積GPU
軟體測試經典測試題（4）
2018-11-02
經典演算法研究系列：九、影象特徵提取與匹配之SIFT演算法
2013-09-22
演算法特徵
DL4J實戰之三：經典卷積例項(LeNet-5)
2021-10-14
卷積
實戰分享丨MySQL 與Django版本匹配相關經驗
2020-09-08
MySqlDjango
經典老遊戲，在獨立開發者手中復甦
2019-07-10
遊戲
30個關於Shell指令碼的經典案例（上）
2019-09-06
指令碼
[Jira] Jira4系列、Jira6系列破解實戰
2016-11-03
Java經典例項：在文字中匹配換行符
2016-11-10
Java
UVA 673 括號的匹配——經典棧的應用
2014-07-25
python怎麼讀_4個Python經典專案實戰，練手必備哦！
2020-12-07
Python
現代化SharePoint經典網站
2018-04-23
網站
面試70問經典回答 - 上
2010-12-31
面試
視訊教程-Python機器學習經典案例實戰-Python
2020-05-28
Python機器學習
cocos2dx實現經典飛機大戰
2016-02-25
Mac經典戰略策略遊戲:
2022-05-05
Mac遊戲
Java經典例項：正規表示式，找到匹配的文字
2016-11-09
Java
交叉熵代價函式定義及其求導推導（讀書筆記）
2019-03-04
熵函式求導筆記
人工智慧領域經典資料集
2018-06-12
人工智慧
經典n皇后問題java程式碼實現
2011-10-19
Java
安卓入門---安卓開發實戰經典1-3章
2014-11-24
安卓
微軟經典面試100題系列（部分）
2014-07-10
微軟面試
Scala 深入淺出實戰經典第45講： scala中context bounds程式碼例項
2015-07-30
Context
Java經典例項：正規表示式，替換匹配的文字
2016-11-09
Java
家庭區域網組建經典例項分享
2016-07-29
rust實戰系列-base64編碼
2022-07-12
Rust
清華尹成帶你實戰GO案例（42)Go 經典hello world
2018-05-22
Go
戲說領域驅動設計（十七）——實體實戰
2022-03-23
經典T-SQL程式碼
2016-07-15
SQL
【Java】經典示例程式碼
2014-02-08
Java
【80/90年代共同記憶！】用容器技術重溫經典，復活坦克大戰！
2018-08-10
網路管理實戰寶典
2017-11-11
Eclipse經典實用操作
2017-09-12
Eclipse
RC4演算法：流密碼演算法的經典之作
2024-03-11
演算法密碼
Java基礎系列22：有關連結串列的經典面試題目解析與程式碼實現
2017-06-06
Java面試題
移動端開發——關於區域性區域滾動總結 | 實戰系列
2019-01-02

【碼上實戰】【立體匹配系列】經典AD-Census: （4）十字交叉域代價聚合

文章目錄

演算法

程式碼實現

類設計

成員函式

成員變數

類實現

實驗

相關文章