Stereo Matching文獻筆記之（八）：《On Building an Accurate Stereo Matching System on Graphics Hardware》讀後感~

周末喽~！又有時間寫寫部落格了，繼續之前的論文讀後感，我發現寫部落格有好幾個好處：1.分享。2.備忘。3.鍛煉表達能力。不知道有沒有童鞋和我有一樣的結論，今天寫寫今年四月份精讀過的一篇文章《On Building an Accurate Stereo Matching System on Graphics Hardware》，文章名咋看起來有點像硬體相關文獻，為什麼叫做一個系統，我想可能是由于作者來自于企業的研究院才這麼起名的，但說出它的别名大家就都知道了，就是AD-Census，這是2011年提出來的算法，作者與SegmentTree是同一人，引用率頗高哈！言歸正傳，本文說一說我對AD-Census的了解，了解不正确的地方，還請各位童鞋批評指正！！

(轉載請注明：http://blog.csdn.net/wsj998689aa/article/details/49403353, 作者：迷霧forest)

算法詳解

這篇文章的亮點我認為有三處：1. 适合并行；2. 基于特征融合的代價計算；3. 基于自适應區域的代價聚合；下面就這幾個方面詳細說一說。

1. 适合并行

這點是相當吸引人的，也是作者的出發點，衆所周知，全局算法不适合并行，為啥？因為建立了複雜而且漂亮的能量函數，需要用同樣複雜而且漂亮的疊代優化算法求解，可惜的是這樣的優化算法如果想并行處理難度太高，并且也快不到哪裡去，是以導緻現有很多全局算法無法得到應用，我們隻能眼睜睜的看着middlebery上的頂級算法，一邊淌着哈喇子，一邊望洋興歎。不過本文關于并行這塊的解釋就不多說了，不親自動手實作是無法體會出來并行帶來的快感。是以主要說說對下面兩點的了解，要知道AD-Census就算拼效果也是得過頭把交椅的。。。

2. 特征融合

我看過的文獻中，很少有在代價計算這一步中融合多種特征的，一般隻采用一種特征而已并且對這塊内容的研究也偏少，作者另辟蹊徑的融合了現有特征AD和Census，AD就是最普通的顔色差的絕對值，立體比對算法中經常使用，其公式如下所示：

其中，i代表不同的通道，這個公式相信大家肯定都一目了然，根據左右兩圖的顔色差的大小來定義cost，這裡值得一提的是幾個名詞：reference image，base image，guidance image和match image，我之前就被繞暈過，因為不同的論文總喜歡用不同的詞彙來代表左圖和右圖，reference image一般翻譯成為參考圖像，base image和guidance image就是reference image，一般被翻譯成為基準圖像，match image一般被翻譯成為比對圖像，那麼參考圖像和比對圖像的關系是什麼呢？我真的好想說，不就是一個代表左圖一個代表右圖嘛！！！但是咱們是搞研究的，必須要嚴謹一些才行，我們回憶一下立體比對的流程，左圖和右圖先做代價計算，怎麼做的呢，周遊左圖中的每一個像素，然後根據視差範圍中的每一個視內插補點，來找到對應右圖的像素，然後根據公式計算代價，然後再針對左圖，周遊每個像素進行代價聚合計算，這就是重點，如果你在左圖上計算代價聚合，那麼左圖就叫做參考圖像，右圖就是比對圖像，反之，就反過來叫。。。

話說回來，Census指的是一種代價計算方法，其屬于非參數代價計算方法中的一種（另外一個代表是rank transform），準确的說它是一種距離度量，它的計算過程的前半部很像經典的紋理特征LBP，就是在給定的視窗内，比較中心像素與周圍鄰居像素之間的大小關系，大了就為1，小了就為0，然後每個像素都對應一個二值編碼序列，然後通過海明距離來表示兩個像素的相似程度，Census代碼如下所示：

int hei = lImg.rows;
    int wid = lImg.cols;
    Mat lGray, rGray;
    Mat tmp;
    lImg.convertTo( tmp, CV_32F );
    cvtColor( tmp, lGray, CV_RGB2GRAY );
    lGray.convertTo( lGray, CV_8U, 255 );
    rImg.convertTo( tmp, CV_32F );
    cvtColor( tmp, rGray, CV_RGB2GRAY );
    rGray.convertTo( rGray, CV_8U, 255 );

    // prepare binary code 
    int H_WD = CENCUS_WND / 2;
    bitset<CENCUS_BIT>* lCode = new bitset<CENCUS_BIT>[ wid * hei ];
    bitset<CENCUS_BIT>* rCode = new bitset<CENCUS_BIT>[ wid * hei ];
    bitset<CENCUS_BIT>* pLCode = NULL;
    bitset<CENCUS_BIT>* pRCode = NULL;

    // 代價計算
    // 計算左圖
    for(int i = 0; i < reflect_pts_num; i++)
    {
        int repeated = reflect[i].repeated;
        if(repeated == 1)
        {
            continue;
        }

        int x     = reflect[i].x;
        int y     = reflect[i].y;
        int index = reflect[i].index;

        uchar pLData = lGray.at<uchar>( y, x );
        pLCode       = &(lCode[index]);

        int bitCnt = 0;
        for( int wy = - H_WD; wy <= H_WD; wy ++ )
        {
            int qy = ( y + wy + hei ) % hei;
            for( int wx = - H_WD; wx <= H_WD; wx ++ )
            {
                if( wy != 0 || wx != 0 ) 
                {
                    int qx = ( x + wx + wid ) % wid;
                    uchar qLData = lGray.at<uchar>( qy, qx );
                    ( *pLCode )[bitCnt] = ( pLData > qLData );
                    bitCnt ++;
                }
            }
        }     
    }

    // 計算右圖
    pRCode = rCode;
    for( int y = 0; y < hei; y ++ ) {
        uchar* pRData = ( uchar* ) ( rGray.ptr<uchar>( y ) );
        for( int x = 0; x < wid; x ++ ) {
            int bitCnt = 0;
            for( int wy = - H_WD; wy <= H_WD; wy ++ ) {
                int qy = ( y + wy + hei ) % hei;
                uchar* qRData = ( uchar* ) ( rGray.ptr<uchar>( qy ) );
                for( int wx = - H_WD; wx <= H_WD; wx ++ ) {
                    if( wy != 0 || wx != 0 ) {
                        int qx = ( x + wx + wid ) % wid;
                        ( *pRCode )[ bitCnt ] = ( pRData[ x ] > qRData[ qx ] );
                        bitCnt ++;
                    }
                }
            }
            pRCode ++;
        }
    }

    // 代價體計算
    bitset<CENCUS_BIT> lB;
    bitset<CENCUS_BIT> rB;
    for(int i = 0; i < reflect_pts_num; i++)
    {
        int repeated = reflect[i].repeated;
        if(repeated == 1)
        {
            continue;
        }

        int x = reflect[i].x;
        int y = reflect[i].y;
        int index = reflect[i].index;

        lB = lCode[index];
        for(int d = 0; d < maxDis; d ++)
        {
            if(x - d >= 0)
            {
                rB = rCode[ index - d ];
                costVol[d].at<double>(y, x) = ( lB ^ rB ).count();                
            }
            else
            {
                costVol[d].at<double>(y, x) = CENCUS_BIT;
            }
        }
    }
           

那為什麼融合起來效果就會好呢？

這個是重點，Census具有灰階不變的特性，所謂灰階不變指的就是像素灰階值的具體大小和編碼之間的相關性不是很強，它隻關心像素之間的大小關系，即使你從5變成了10，但隻要中心像素是15，就一點事情都木有，這樣的性質我們肯定會想到它一定對噪聲和魯棒，的确是這樣。但是它的缺點也很明顯，按照文章的說法， 對于結構重複的區域這個特征就不行了，那基于顔色的特征AD呢？它是對顔色值很敏感的，一旦區域内顔色相近（低紋理）或者有噪聲那麼挂的妥妥的，但是對重複結構卻不會這樣，基于這種互補的可能性，作者嘗試将二者進行融合，這是一種很簡單的線性融合但是卻取得了很好的效果：

其中，下面公式的目的就是歸一化，我們注意看，兩種計算方法得到的取值範圍可能會相差很大，比如一個取值1000，另一個可能隻取值10，這樣便有必要将兩者都歸一化到[0，1]區間，這樣計算出來的C更加有說服力，這是一種常用的手段。論文中給出了AD，Census，AD-Census對一些細節處理的效果圖，可以看得出來各自的優缺點，第一行是重複紋理區域，第二行是低紋理區域，白色與黑色都說明計算的結果很差。

下面分别給出AD、Census、AD-Census的效果圖，正如上一大段的分析，我們會發現AD往往比Census在物體邊緣上的處理更好一些，邊緣明顯清晰，但是AD得到的噪聲太多，并且在低紋理區域，比如中間那個燈罩，AD出現了很大的空洞，在這一點上Census做的相對較好，AD-Census在物體邊緣上的效果是二者的折中，但噪聲更少，整體效果更加理想。

3. 自适應區域

這塊内容是該文獻的重點，再說之前我們先回顧一下一般的代價聚合思想，局部算法采用一個固定或者自适應視窗來代價聚合，全局算法采用整幅圖像得到的抽象結構來代價聚合，例如MST，馬爾科夫圖模型等等，這些在之前的部落格都描述過，我們總結一下二者的共同點，不就是事先确定一個有意義的區域嗎？确定好了之後便可直接在這個區域内進行代價聚合，OK！統一了思想，就可以說說這篇文章是怎麼做的了。

ADCensus建立了一個灰常有意思的，工程化的區域結構，分割算法耗時不舍得用，幹脆我直接用相對暴力的方案對圖像進行分割好了，看圖說話，下圖就是作者采用的分割方法：

采用方法的思想很簡單，目前像素假設是p，我對p先進行垂直方向的周遊，如果像素q滿足兩個限制，那麼就算同一分割區域，否則周遊停止，然後再在得到的N個q和p的水準方向根據同樣的規則進行周遊，于是就得到了對應的分割區域。然而，作者從來都沒有說對圖像進行了分割，隻是說确定p的cross，其實就是分割的意思。限制如下所示，這個限制是經過改造過的（cross region不是作者提出來的，也是引用他人的方法，作者對限制從兩個擴充為三個），這麼改的原因是作者考慮到了之前的限制方案會導緻将邊緣點也包括進去，這樣會對邊緣的視差計算十分不利，于是提出了一個更加嚴格的限制形式，考慮到了相鄰像素色彩上的差異。

下面就可以在像素所屬的區域内進行代價聚合了，就是簡單地将區域内各個像素的ADCensus值相加，但是這裡有個很大的問題作者說的很模糊，作者在文章中強調，為了保證代價聚合的穩定性，需要進行“先水準後垂直”，“先垂直後水準”兩種代價聚合方案各兩次。我當時很迷糊，因為我認為一旦每個像素所屬的區域确定了，這兩種方案得到的代價聚合值肯定是一樣的，除非作者在建構區域的時候，也是采用“先水準後垂直”，“先垂直後水準”的方案，根據我這個猜想還寫了代碼做了實驗，結果發現效果一般，并且反複的看文獻，發現作者說的很明确，就是代價聚合的時候采用兩種方案，而不是區域建構的時候，這到底是怎麼回事？

忽然有一天秋葉旁落，我終于明白了作者的意圖，别忘記了，我們還有一幅圖像呢，就是右圖！！！這麼重要的線索我竟然忽略了，答案是要對左圖和右圖分别進行區域建構，然後代價聚合的時候，如果采用“先水準後垂直”的方案，那麼就先取左右兩個對應區域的交集，然後在将交集中的代價值都加起來，進一步計算垂直方向的代價值的和。另一種方案就是先垂直方向的區域相交，再水準求和，這樣就能得到不同的代價聚合結果。兩種方案各自執行兩次，每一次都用之前新得到的代價聚合值，注意這裡作者隻是簡單的将區域内各像素對應的代價值相加，沒有考慮到權值，可能是為了速度吧，當然加上權值效果肯定會更加好一些。

上圖是代價聚合的過程，是先水準後垂直的方案，代價聚合之後，對于一般的局部算法而言，基本上就到此為止了，但是ADCensus還有一個大招呼之欲來，那就是大名鼎鼎的“掃描線優化”，這個掃描線優化是動态規劃的一種方法，在史上最經典立體比對文獻SGM中首次被使用，具體的思想本文就不詳細說了，由于代價聚合的結果不大靠譜，可以考慮将其視作資料項，建立全局能量函數（公式如下所示），這樣便直接過渡到了全局算法。

其中，第一項C就是代價聚合項，後面兩項分别考慮到了視差的微變（低紋理區域）和劇烈變化（物體邊緣），優化這個能量函數做采用的方法就是“掃描線優化法”，其公式如下所示，這裡和SGM一模一樣，不做過多的解釋，因為以後還會和大家聊聊我對SGM的了解。

然而，ADCensus在P1和P2的設定上不同于SGM，進行了調整，具體的公式就不再粘貼了。做了實驗，如下三幅圖所示，左圖是直接利用代價聚合得到的視差圖，中間一副是進一步通過掃描線優化之後得到的視差圖，第三幅圖是二者的差異，藍色代表差異微弱，紅色代表差異較大，黑色代表沒有差異，可以明顯看到，經過掃描線優化處理之後，視差圖在細節上明顯處理的更好，邊緣更加平滑，但是出現了拖尾現象（燈杆子那裡）。

後續

和其他文獻相同，作者也采用left-right-check的方法将像素點分類成為三種：遮擋點，非穩定點，穩定點。對于遮擋點和非穩定點，隻能基于穩定點的視差傳播來解決，本文在視差傳播這一塊采用了兩種方法，一個是疊代區域投票法，另一個是16方向極線插值法，後者具體來說：沿着點p的16個方向開始搜尋，如果p是遮擋點，那麼遇到的第一個穩定點的視差就是p的視差，如果p是非穩定點，那麼遇到第一個與p點顔色相近的穩定點的視差作為p的視差。針對于視差圖的邊緣，作者直接提取兩側的像素的代價，如果有一個代價比較小，那麼就采用這個點的視內插補點作為邊緣點的視內插補點，至于邊緣點是如何檢測出來的，很簡單，對視差圖随便應用于一個邊緣檢測算法即可。做完這些之後，别忘記亞像素求精，這和WTA一樣，是必不可少的。流程圖如下，忘記說了，再來一個中值濾波吧，因為大家都這麼玩，屬于後處理潛規則。。

這裡重點說說疊代區域投票法，這是我自己的翻譯，英文稱呼是“Iterative Region Voting”，它的目的是對outlier進行填充處理，一般來說outlier遮擋點居多，之前的部落格也介紹過，填充最常用的方法就是用附近的穩定點就行了，省時省力，就是不利于并行處理，作者要設計的是一個完全适合GPU程式設計的算法，是以采用了疊代區域投票法，具體做法是對之前區域建構所得到的每個區域求取視差直方圖（不要歸一化），例如，得到的直方圖共有15個bin，最大的bin值是8，那麼outlier的視差就由這個8來決定，但是穩定點的個數必須得比較多，比較多才有統計穩定性，數學形式化就是如下公式：

其中，Sp就是穩定點的個數，Hp就是最大bin值，一般Ts，Th兩個參數都經驗設定，說白了就是得好好調試一番。此外，這個方法是疊代的，這隻是針對穩定點個數具有統計意義的區域，有些outlier由于區域内穩定點個數不滿足公式，這樣的區域用此方法是處理不來的，隻能進一步通過16方向極線插值來進一步填充，二者配合起來能夠取得不錯的效果，自己做了實驗，這兩種方法的順序也必須一先一後，否則效果也不行，說明一個是大迂回戰略，目的是消滅有生力量，一個是殲滅戰，打的是漏網之魚，二者珠聯璧合，可喜可賀。

總結

ADCensus是個好算法！簡單易于實作，完全有利于并行處理，具有實用化價值，大家可以動手編碼試一試，文章裡面也提供了并行計算指導。這就是我的了解，如有錯誤請不吝賜教哈~

不過我對這篇論文有個疑問，就是在自适應區域那裡，為什麼作者偏偏要加上一個自适應區域内代價聚合呢？總所周知，SGM中沒有這一個步驟，并且自适應區域其實就是分割的目的，在區域内代價聚合并沒有設定每個像素比對代價之前的權值，如果我們将自适應區域這一步去掉，其實ADCensus和SGM唯一的差別就隻有比對代價的計算方式了，前者不用多說，本文解釋了是一種融合代價，後者采用的或是BT或是互資訊，莫非融合特征就必須伴随自适應區域？？未來有空閑我便會重新編碼去究其原因，最近琅琊榜實在是看多了，抱着懷疑一切的态度也開始懷疑上了作者的動機，根據我的猜測，原因無外乎有三種。

1. 自适應區域确實有用，沒有了它視差圖效果就會大打折扣，所謂正途。

2. 自适應區域是為了發表論文才加上去的，如果沒有這一塊，論文的表面含金量就縮水了，所謂雞肋。

3. 自适應區域的真實目的是為了後續的疊代區域求精，是以索性也用于代價聚合，所謂醉翁之意不在酒。

不知是否有童鞋也有同樣的疑問，還請不吝賜教。