Cross-Scale Cost Aggregation for Stereo Matching立體匹配演算法介紹
阿新 • • 發佈：2019-01-31
最近，研究了下CVPR2014上的一篇基於多尺度代價聚合的立體匹配演算法，這個作者提供了原始碼，運行了下，發現效果真心不錯，不開後端處理的話，時間在0.4s左右。這個演算法比較牛逼的有兩點：
1：結合多尺度思想，對原始影象進行下采樣，然後在每層影象上計算匹配代價，進行代價聚合，然後多尺度得到的視差進行結合，作為最終的代價聚合值。
2：提供了一個框架，裡面包含立體匹配很多常用的經典的方法，可以在每一步使用不同的方式完成。
參考一篇部落格：http://blog.csdn.net/wsj998689aa/article/details/44411215
本文結合原始碼，對這個演算法的流程進行詳細分析。
首先進行匹配代價計算：
原始碼裡有三種方法，CEN（就是census），CG（Census + TAD + 梯度）,GRD（梯度+TAD） ，這三種方式幾乎是當下進行第一步匹配代價計算的常用方式。
1：CEN：選定一個視窗，對視窗中每一個畫素與中心畫素進行比較，大於中心畫素即為0，否則為1。從而得到一個二進位制系列，分別對左圖和右圖進行計算每個畫素的匹配代價。並得到初步的代價聚合，計算每個畫素在視差範圍內每個可能視差的代價聚合值。左右匹配，即兩個二進位制系列相似值，每個對應位是否相等。
void CenCC::buildCV( const Mat& lImg, const Mat& rImg, const int maxDis, Mat* costVol )
{
// for TAD + Grd input image must be CV_64FC3//輸入的影象必須是64FC3
CV_Assert( lImg.type() == CV_64FC3 && rImg.type() == CV_64FC3 );
int hei = lImg.rows;//影象寬高
int wid = lImg.cols;
Mat lGray, rGray;
Mat tmp;
lImg.convertTo( tmp, CV_32F );//64FC3影象資料轉化成32F 下面是從右向左進行匹配的，就是在尋找視差方向上時不同而已
// build cost volume
bitset<CENCUS_BIT> lB;
bitset<CENCUS_BIT> rB;
pRCode = rCode;
for( int y = 0; y < hei; y ++ ) {
int index = y * wid;
for( int x = 0; x < wid; x ++ ) {
rB = *pRCode;
for( int d = 0; d < maxDis; d ++ ) {
double* cost = ( double* ) rCostVol[ d ].ptr<double>( y );
cost[ x ] = CENCUS_BIT;
if( x + d < wid ) {
lB = lCode[ index + x + d ];
cost[ x ] = ( rB ^ lB ).count();
}
}
pRCode ++;
}
}這裡就完成了利用census原理計算匹配代價，然後進行初始代價聚合。
2：TAD + 梯度
這個演算法運用影象三通道顏色值的差的絕對值AD，結合一個閾值，做成TAD樣式，求三通道均值; 然後利用引數為1的sobel運算元求解x方向的影象梯度，同樣加一個閾值，做成TGRD的樣式。這兩個特徵進行加權相加。這裡還對影象邊緣進行判斷，是邊緣的話，就把左右TAD和TGRD的右換成邊緣閾值。
// X Gradient
// sobel size must be 1
Sobel( lGray, lGrdX, CV_64F, 1, 0, 1 );
Sobel( rGray, rGrdX, CV_64F, 1, 0, 1 );
lGrdX += 0.5;
rGrdX += 0.5;inline double myCostGrd( double* lC, double* rC,
double* lG, double* rG )
{
double clrDiff = 0;
// three color
for( int c = 0; c < 3; c ++ ) {
double temp = fabs( lC[ c ] - rC[ c ] );
clrDiff += temp;
}
clrDiff *= 0.3333333333;
// gradient diff
double grdDiff = fabs( lG[ 0 ] - rG[ 0 ] );
clrDiff = clrDiff > TAU_1 ? TAU_1 : clrDiff;
grdDiff = grdDiff > TAU_2 ? TAU_2 : grdDiff;
return ALPHA * clrDiff + ( 1 - ALPHA ) * grdDiff;
}3:census + TAD + 梯度
就是把上面兩個加權疊加，下式的前面是TAD+TGRD的匹配代價，後面是census的匹配代價，得到最終的匹配代價。
cost[ x ] = 2 - exp( - tmpCost / 35 ) - exp( - 255 * cost[ x ] / 15 );//這就是TAD+TGD+Census的匹配代價計算接下來就是進行代價聚合了。
代價聚合：原始碼提供了多種方式，雙邊濾波器bilateral filter，引導濾波器guided image filter，box filter，還有非區域性的代價聚合方法（《A Non-Local Cost Aggregation Method for Stereo Matching》），基於分割樹的代價聚合方法《Segment-Tree based Cost Aggregation for Stereo Matching》。基於濾波器的方法是設定一個視窗，在視窗內進行代價聚合。雙邊濾波器就是對視窗內畫素進行距離加權和亮度加權。
後兩種都拋棄了固定視窗模式,基於NL的代價聚合是使用最小生成樹代替視窗。基於ST的代價聚合是使用分割塊代替視窗。
這裡主要分析下基於分割樹的代價聚合方法。
這個方式是CVPR2013的文章，是NL方法的升級版。原理類似。
這個方式中，可以不採用第一步計算的匹配代價，不過它自身計算匹配代價和TAD+TGRD類似。這裡就使用第一步的匹配代價。
1：初始化一個圖G（V，E），每個畫素就是頂點V，邊是每兩個畫素之間進行三通道分別求差的絕對值，取最大的那個差值作為邊的權值。
float CColorWeight::GetWeight(int x0, int y0, int x1, int y1) const {//得到三通道兩個畫素值差值最大的那個通道的差值的絕對值
return (float)std::max(
std::max( abs(imgPtr(y0,x0)[0] - imgPtr(y1,x1)[0]), abs(imgPtr(y0,x0)[1] - imgPtr(y1,x1)[1]) ),
abs(imgPtr(y0,x0)[2] - imgPtr(y1,x1)[2])
);
}2：利用邊的權值對原影象進行分割，構建分割樹。
依次對每個邊連線的畫素，判斷邊權值是否大於閾值，大於則不在同一個分割塊，否則就在同一個分割塊。這裡先對邊按權值大小進行升序排列，因為權值小的更容易是一個分割塊。這樣就構成了一個個分割塊，每個分割塊就是一個小的樹。
然後再次判斷每個邊連線的兩個畫素，意在連線每個分割塊，使每個小樹作為影象整體樹的子樹連線在一起。
這樣就構成了整幅影象的樹。
int cnt = 0;
universe *segment_graph(int num_vertices, int num_edges, edge *edges,
float c, unsigned char *edges_mask) {
// sort edges by weight 按邊的權重大小升序排列邊
std::sort(edges, edges + num_edges);
// make a disjoint-set forest 產生一個並查集森林，每棵樹中的成員都可由根結點所代表，意思是將圖中的點分割成子樹，每個子樹的根節點代表子樹。然後再把各子樹連線起來
universe *u = new universe(num_vertices);
// init thresholds
float *threshold = new float[num_vertices];
for (int i = 0; i < num_vertices; i++) //對每個畫素點，初始化閾值為c/1
threshold[i] = THRESHOLD(1,c);
// for each edge, in non-decreasing weight order...對每個邊
//這裡就是從第0個畫素開始，若邊權重小於閾值，就劃為同一個分割塊，若不小於，則不是同一個分割塊，這裡就得到整幅影象的若干個分割塊。
//每個分割塊的任意畫素的elts[].p都是最開始不是上一個分割塊的那個畫素的座標一維值。size記錄當前分割塊所包含的畫素的數量。
//當一個邊的一個端點首次加入比較時，此時的邊有個mask值為255,分割塊之間的邊的mask都為0，重複比較的兩點之間的邊的mask也為0。這裡是避免樹上節點的重複連線。
for (int i = 0; i < num_edges; i++) {
edge *pedge = &edges[i];
// components connected by this edge找到這個邊連線的兩個畫素
int a = u->find(pedge->a);
int b = u->find(pedge->b);
if (a != b) //如果這兩個畫素不在同一個分割塊內
{
if (pedge->w <= threshold[a] && pedge->w <= threshold[b])//如果這兩個畫素連的邊的權值小於等於這兩個畫素的閾值
{//就把這兩個畫素合併，a作為b的根節點，b作為子節點，代表a，以後u->find(a)得到的就是b的屬性p這樣就可以每個節點只有兩個子節點
edges_mask[i]=255;
u->join(a, b);
a = u->find(a);
threshold[a] = pedge->w + THRESHOLD(u->size(a), c);//THRESHOLD(u->size(a)得到此時這個分塊的大小，其實是b的大小。。。
}
}
}
//cv::Mat sg_img = cv::Mat(256,320,CV_8UC3);
char sg_img_1[256*320*3];
cv::Mat sg_img = cv::Mat(256,320,CV_8UC3,sg_img_1);
if((cnt++)%3 == 0)
{
for(int i=0;i<num_edges;i++)
{
int a = u->find(edges[i].a);
int b = u->find(edges[i].b);
sg_img_1[edges[i].a * 3] = a%(256*256);
sg_img_1[edges[i].a * 3+1 ] =256 - a%256;
sg_img_1[edges[i].a * 3+2 ] = a;
sg_img_1[edges[i].b * 3] = b%(256*256);
sg_img_1[edges[i].b * 3+1 ] =256 - b%256;
sg_img_1[edges[i].b * 3+2 ] = b;
}
cv::imshow("sg_img",sg_img);
cv::waitKey(1);
}
//added by X. Sun: re-organizing the structures to be a single tree
for (int i = 0; i < num_edges; i++)
{
int a = u->find(edges[i].a);
int b = u->find(edges[i].b);
if (a != b)
{
int size_min = MIN(u->size(a), u->size(b));//求分塊a和b的最少大小
u->join(a, b);
//record
edges_mask[i]=255;
if(size_min > MIN_SIZE_SEG)
edges[i].w += PENALTY_CROSS_SEG;
}
}
// free up
delete []threshold;
return u;
}3：對整顆樹的所有節點計算其父節點和子節點，並進行排序。組成一個整體樹，從根節點到葉節點進行儲存，便於下面基於樹結構進行迴圈代價聚合。
//step 2: build node based graph 這裡求出每個點與周圍四鄰域的點的父子關係，並給出這兩個點之間的距離。
//因為mask並不全為255,所以只有當這個點的鄰域點是根據當前點判斷進入分割塊，或者分割塊之間，才算作子節點。
TreeNode *AdjTable = new TreeNode[pixelsNum];
for(int i = 0;i < pixelsNum;i++) AdjTable[i].id = i;//初始化每個畫素的id為畫素的位置
for(int i = 0;i < edgeNum;i++) {
if(!edges_mask[i]) continue;//如果該邊沒有被計算是否分割，則繼續
int pa = edges[i].a;
int pb = edges[i].b;
int dis = std::min(int(edges[i].w * weightProvider.GetScale() + 0.5f), 255);
int x0, y0, x1, y1;
x0 = pa % m_imgSize.width; y0 = pa / m_imgSize.width;//求出x0在影象中的確切位置
x1 = pb % m_imgSize.width; y1 = pb / m_imgSize.width;
TreeNode &nodeA = AdjTable[pa];
TreeNode &nodeB = AdjTable[pb];
nodeA.children[nodeA.childrenNum].id = pb;
nodeA.children[nodeA.childrenNum].dist = (uchar)dis;
nodeA.childrenNum++;
nodeB.children[nodeB.childrenNum].id = pa;
nodeB.children[nodeB.childrenNum].dist = (uchar)dis;
nodeB.childrenNum++;
}
//step 3: build ordered tree
if(!m_tree.empty()) m_tree.clear();
m_tree.resize(pixelsNum);
bool *isVisited = new bool[pixelsNum];
memset(isVisited, 0, sizeof(bool) * pixelsNum);
m_tree[0] = AdjTable[0];
isVisited[0] = true;
int start = 0, end = 1;
while(start < end) {
TreeNode &p = m_tree[start++];
for(int i = 0;i < p.childrenNum;i++) {
if(isVisited[p.children[i].id]) continue;
isVisited[p.children[i].id] = true;
TreeNode c;
c.id = p.children[i].id;
c.father.id = p.id;
c.father.dist = p.children[i].dist;
TreeNode &t = AdjTable[c.id];
for(int j = 0;j < t.childrenNum;j++) {
if(t.children[j].id != p.id) {
c.children[c.childrenNum++] = t.children[j];
}
}
m_tree[end++] = c;
}
}4：從葉節點到根節點，然後再從根節點到葉節點進行代價聚合。論文中的聚合圖很好的說明了這個問題。
從葉節點到根節點就是每個畫素點的代價聚合值為當前畫素的匹配代價 + 該畫素所有葉節點的匹配代價聚合值的加權和，權值就是邊權值。從下到上，逐層計算。
從根節點到葉節點就是對每個畫素點在上一步計算的代價聚合值 + 權值 ×（其父節點的代價聚合值 - 其本身畫素對父節點代價聚合值的貢獻）。
void CSegmentTree::Filter(cv::Mat costVol, int maxLevel) {
cv::Mat costBuffer = costVol.clone();
KIdx_<float, 3> costPtr((float *)costVol.data, m_imgSize.height, m_imgSize.width, maxLevel);
KIdx_<float, 3> bufferPtr((float *)costBuffer.data, m_imgSize.height, m_imgSize.width, maxLevel);
int pixelsNum = m_imgSize.area();
//first pass: from leaf to root
for(int i = pixelsNum - 1;i >= 0;i--) {
TreeNode &node = m_tree[i];
float *cost = &bufferPtr(node.id * maxLevel);
for(int z = 0;z < node.childrenNum;z++) {
float *child_cost = &bufferPtr(node.children[z].id * maxLevel);
float weight = m_table[node.children[z].dist];
for(int k = 0;k < maxLevel;k++) {
cost[k] += child_cost[k] * weight;
}
}
}
//second pass: from root to leaf
memcpy(&costPtr(0), &bufferPtr(0), sizeof(float) * maxLevel);//這個就是最開始的父節點。
for(int i = 1;i < pixelsNum;i++) {//注意上面的迴圈順序和此處的不同
TreeNode &node = m_tree[i];
float *final_cost = &costPtr(node.id * maxLevel);
float *cur_cost = &bufferPtr(node.id * maxLevel);
float *father_cost = &costPtr(node.father.id * maxLevel);
float weight = m_table[node.father.dist];
for(int k = 0;k < maxLevel;k++) {
final_cost[k] = weight * (father_cost[k] - weight * cur_cost[k]) + cur_cost[k];
}
}
}這樣就完成了每個畫素的代價聚合。同理對右圖匹配左圖也一樣。
由於這篇文章是基於多尺度的代價聚合，那麼就對原始影象進行下采樣，實驗設定為5層金字塔。在每一層上進行匹配代價計算，代價聚合。這裡要注意的是每次下采樣，影象縮小一倍，那麼視差範圍變為maxDis_c = maxDis_c / 2 + 1; 求出的視差值也要縮小1倍，disSc_c *= 2;//視差值乘以2.因為影象變小了一倍，然後視差d肯定比原視差差一倍。
求出每層影象的代價聚合值之後，合併到原始影象層。
這裡先計算每層影象求得的代價聚合值相對於原始層影象代價聚合值的加權因子。
Mat regMat = Mat::zeros( PY_LVL, PY_LVL, CV_64FC1 );
for( int s = 0; s < PY_LVL; s ++ ) { //對金字塔的每層之間，及每層與上下層之間的相關因子進行計算
if( s == 0 ) {
regMat.at<double>( s, s ) = 1 + REG_LAMBDA;
regMat.at<double>( s, s + 1 ) = - REG_LAMBDA;
} else if( s == PY_LVL - 1 ) {
regMat.at<double>( s, s ) = 1 + REG_LAMBDA;
regMat.at<double>( s, s - 1 ) = - REG_LAMBDA;
} else {
regMat.at<double>( s, s ) = 1 + 2 * REG_LAMBDA;
regMat.at<double>( s, s - 1 ) = - REG_LAMBDA;
regMat.at<double>( s, s + 1 ) = - REG_LAMBDA;
}
}
Mat regInv = regMat.inv( );//對因子矩陣求反
double* invWgt = new double[ PY_LVL ];
for( int s = 0; s < PY_LVL; s ++ ) {
invWgt[ s ] = regInv.at<double>( 0, s );
}然後從原始影象的每個畫素開始，依次對每個畫素在每一層的代價聚合進行加權和。最後得到的和就是原始層影象的代價聚合值。
//
// Left Cost Volume
//
for( int d = 1; d < smPyr[ 0 ]->maxDis; d ++ ) {
// printf( ".s.v." );
for( int y = 0; y < hei; y ++ ) {
for( int x = 0; x < wid; x ++ ) {
int curY = y;
int curX = x;
int curD = d;
double sum = 0;
for( int s = 0; s < PY_LVL; s ++ ) {
double curCost = smPyr[ s ]->costVol[ curD ].at<double>( curY, curX );
#ifdef _DEBUG
if( y == 160 && x == 160 ) {
printf( "\ns=%d(%d,%d)\td=%d\tcost=%.4lf", s, curY, curX, curD, curCost );
}
#endif
sum += invWgt[ s ] * curCost;//這裡就體現了多尺度的特點，就是多尺度的畫素代價聚合值加權求和
curY = curY / 2;//這是為了求下一層的代價聚合值
curX = curX / 2;
curD = ( curD + 1 ) / 2;
}
smPyr[ 0 ]->costVol[ d ].at<double>( y, x ) = sum;//再把所有層的代價聚合值加權和得到的總的代價聚合值返給原始影象層，作為最終的代價聚合值
}
}求右圖匹配左圖也類似。
其實多尺度方法也算是比較簡單了，運用多尺度進行代價聚合，作者說符合生物學習慣，先遠後近，遠處看個大概，近處就能看清細節了。思路很好，但含金量並不算太高，原來發CVPR包裝也是很必要的。開玩笑了，實測這種方法效果還真心不錯。
代價聚合完之後，那就WTA求最優視差值吧。WTA也能叫做演算法……
void SSCA::Match( void )//WTA演算法求左右視差
{
//printf( "\n\tMatch" );
for( int y = 0; y < hei; y ++ ) {
uchar* lDisData = ( uchar* ) lDis.ptr<uchar>( y );
for( int x = 0; x < wid; x ++ ) {
double minCost = DOUBLE_MAX;
int minDis = 0;
for( int d = 1; d < maxDis; d ++ ) {
double* costData = ( double* )costVol[ d ].ptr<double>( y );
if( costData[ x ] < minCost ) {
minCost = costData[ x ];
minDis = d;
}
}
lDisData[ x ] = minDis * disSc;
}
}
#ifdef COMPUTE_RIGHT
for( int y = 0; y < hei; y ++ ) {
uchar* rDisData = ( uchar* ) rDis.ptr<uchar>( y );
for( int x = 0; x < wid; x ++ ) {
double minCost = DOUBLE_MAX;
int minDis = 0;
for( int d = 1; d < maxDis; d ++ ) {
double* costData = ( double* )rCostVol[ d ].ptr<double>( y );
if( costData[ x ] < minCost ) {
minCost = costData[ x ];
minDis = d;
}
}
rDisData[ x ] = minDis * disSc;
}
}
#endif
}求出初始視差後，接下來就是比較繁瑣的視差求精了。
原始碼中給出了兩種方式，
1：加權中值濾波，進行左右一致性檢測; 進行不可靠點進行找左右最近的有效點選視差值較小的那個; 然後進行對不可靠點進行加權中值濾波;
2：分割後，進行左右一致性檢測; 進行不可靠點進行找左右最近的有效點選視差值較小的那個; 然後進行對不可靠點進行加權中值濾波; 然後求分割塊的平面引數，對分割塊內的每個畫素進行視差求精。
後處理一般都是這樣，很多方法都是前面隨便求個初始值，重要的就是後處理。這塊程式碼比較多，有時間再說了。
本文只是簡單的從程式碼上進行分析這篇論文的原理，沒有放實驗圖片，我也是太懶了….。