學習《Hardware-Efficient Bilateral Filtering for Stereo Matching》一文筆記。，

學習《Hardware-Efficient Bilateral Filtering for Stereo Matching》一文筆記。，

　　本人一直用部落格園的編輯器西部落格，由於CSDN和部落格園模板不一樣，為了您閱讀的方便，建議您通過部落格園的部落格閱讀本文：http://www.cnblogs.com/Imageshop/p/3839879.html

 

      個人收藏了很多香港大學、香港科技大學以及香港中文大學裡專門搞映像研究一些博士的個人網站，一般會不週期性瀏覽他們的作品，最近在看楊慶雄的網點時，發現他又寫了一篇雙邊濾波的文章，並且配有原始碼，於是下載下來研讀了一番，這裡僅僅對一些過程做簡單的記錄，以防時間久了忘記。

　　關於楊慶雄的相關文章可見：Hardware-Efficient Bilateral Filtering for Stereo Matching 以及一篇  Recursive Bilateral Filtering，都配有相關的原始碼。

　　在《Hardware-Efficient Bilateral Filtering for Stereo Matching 》一文中，作者提出了一種新的更加適合於硬體（GPU)實現的高效的雙邊濾波演算法,但作者的參考代碼中只提供了CPU版本的程式，這對於來說已經足夠了，我只玩CPU。

　　和傳統的雙邊濾波不同，這篇文章最佳化的雙邊濾波在空域上使用的是均值濾波器，因此位於待處理像素的中心周邊位置的各像素在空間域的權重都是一樣的，而在值域周邊像素的權重和兩個像素值域的差異有關。用原文的話說就是：This paper investigates into the situation when the spatial distanceis ignored。

　　那麼這樣的雙邊濾波用離散化後的公式表示就是：

                 

　　上式中I(q)表示位於q處像素的強度值，G為衡量兩個像素值域關聯性的一個函數。如果G的參數從另外一個同樣大小的映像T中取樣，則上式就表示在此情況下的聯合雙邊濾波。

　　　　　　

　　實際上，可以把公式（1）即普通的雙邊濾波看成是聯合雙邊濾波（公式2）的一種特殊的情形，即T=I；    

     針對演算法的最佳化，作者提出了金子塔式的演算法，以聯合雙邊濾波為例，如所示。

          

　　以兩層金字塔為例進行說明，第一層稱為fine level，實際上其就等於原始的映像I和導向映像T，第二層稱之為coarse level，在 finle level 中的像素點p， 用

  　　　　         表示在coarse level中P點對應的整數座標

　　　　 　　

                    分別表示在fine level被採樣到的以及沒被採樣到的位置，比如中第一行的54和50

 

      那麼按照演算法那的要求，需要分別按式（6）（7）計算coarse level中T和I的值：

　　　   　

      也就是說coarse level中的映像I採用fine level中 sampled 和 miss-sampled的4個點（二維）的平均值，而導向圖T則直接使用sampled 的值，這主要是為了保留邊緣。這一點很重要，在很多類似的技術中，都需要這樣處理，即有些下採樣的過程必須使用最近鄰取樣（插值），因為其他的比如雙線性都會使得邊緣被羽化。

　　式（7）中，即原始映像。

　　那麼根據式（2），聯合雙邊濾波也可以寫成：

　　             （因為，所以這個式子很容易看懂 ）

     同樣的，在coarse level中，用作為導向圖的雙邊濾波可以寫成：

           （很明顯這裡的累加的數量已經減為上一層的1/4了）

    根據式（4）（5）（7），上式可進一步寫成：

         　　（看不懂？注意式（7）那個平均值的關係，再靜心的想一想就明白了）

     如果在輸入的映像，存在關係並且，對於一維訊號，就是指輸入資料是類似這樣的 10  10 24 24 5 5 7 7 8 8.......那麼下式                

     　　　　（看不懂為什麼非要把那個q!=p提出來）

　　則有  ；

     我們計:

                 　　　　　　（再回頭看看公式（4），你會發現什麼）

　　可以證明：

       　　（這個證明我也懶得驗證了）

　  式（16）可以看出通過在coarse level中利用導向映像計算的結果和fine level中的映像進行插值則可以得到fine level中映像的雙邊濾波的結果。

     但是以上推導全部是基於這個前提條件的，實際映像中很少能還有滿足該條件的。因此，作者提出再對J做一個雙邊濾波，為了覆蓋miss-sampled的那個點，雙邊濾波自小的半徑必須是1，同時為了考慮效果和效率，一般半徑取2比較合適了。

     這樣做，其實仔細想起來似乎是個很有趣的迴圈，對J再做一個雙邊濾波，又可以用上面的方式，對第二層coarse level在做一次coarse level，如此迴圈下去。要知道每進行一次coarse level處理的資料量只是原始的1/4了。當然，我們不能無限下去，適當的次數是必要的。當達到底層的coarse level後，在反向的進行每一層fine level的操作，直到最頂層。

     如果每層coarse level的雙邊濾波的半徑都為r，如果進行了N層的處理，最相當於對最頂層的映像進行半徑為 r * 2N的雙邊濾波。

     從理論上分析，當N取無限大, r =2 時， 該演算法的執行時間為 = 4/3 *原始頂層映像5×5 地區的雙邊濾波的直接實現。因此演算法的執行時間於range kernel的參數無關。 

     更多的細節可能需要讀者參考文章附帶的代碼進行體會。

     作者提出該演算法非常適合於GPU實現，並且給出了GPU和CPU版本程式的速度比較，如下表：

    

　　其中的HEBF(Hardware-Efficient Bilateral Filtering)即為本文的演算法，加速比達到了282，我對此有些疑問的。

     第一：作者的計時標準是什麼，是從使用者提供了輸入資料I，導向資料T，以及調節參數後，包括記憶體配置等等處理整個過程的時間嗎？

     第二：如果是第一條所說的，那麼 那些從作者提供的CPU代碼上看有些部分的代碼必須在CPU上完成，個人認為就那些代碼的執行時間也不止2.4ms的，那作者如何做到的。

     第三：作者提供的CPU代碼可以看出，整個程式的並行性並不是特彆強，前後的執行也有嚴重的依賴關係。只有一些內部的迴圈有很強的並行性，但那些可並行性的代碼的計算量特別小，如果是在多核上用多線程來並發，可能線程切換的時間比本身的計算時間還長，GPU雖然是輕量級的核，是否會好點呢，本人不瞭解GPU的特性，有待高人解答。

　　CPU版本的程式如果想利用多線程實現，有一個建議就是把每通道的資料單獨來處理，雖然會多了前後的拆分和合并的過程，以及一些迴圈變數的多次計量外，在四核的電腦上估計速度會有一倍的提升。（R/G/B三通道三個線程必行執行）。    

      總的來說，文章雖然採用的不是傳統意義上標準的雙邊濾波，但是效果和速度還是相對來說比較不錯的，在很多應用場合是完全可以替代標準雙邊濾波使用的。

      關於彩色映像的值域相似度的測量，原文使用的是如下的方式：

/// <summary>/// 計算兩個像素的距離，此處可以使用多種表達方式，內嵌函式。/// </summary>/// <param name="PixelA">像素A的記憶體位址。</param>/// <param name="PixelB">像素B的記憶體位址。</param>/// <returns>返回兩個像素的距離測度。</returns>inline unsigned char MaxRGBEuclideanDistance (unsigned char *PixelA ,unsigned char *PixelB){    int DiffR, DiffG, DiffB;    DiffB = abs(PixelA[0] - PixelB[0]);     DiffG = abs(PixelA[1] - PixelB[1]);     DiffR = abs(PixelA[2] - PixelB[2]);     return max(max(DiffR, DiffG), DiffB);}

　　個人認為也可以從其他方面考慮，比如歐式距離，棋盤距離，甚至每個通道單獨想減求絕對值都可以。

     在coarse level層面的雙邊濾波上，作者使用的是brute-force 方式，即最原始的暴力迴圈：

/// <summary>/// 對輸入資料按照導向資料值進行雙邊濾波處理。/// </summary>/// <param name="WegightedImage">記錄加權乘積後的映像資料。</param>/// <param name="ImageWeight">記錄加權值。</param>/// <param name="Table">權值的尋找表。</param>/// <param name="DownInput">輸入資料。</param>/// <param name="Guide">導向資料。</param>/// <param name="Width">輸入資料的一維尺寸。</param>/// <param name="Height">輸入資料的二位尺寸。</param>/// <param name="Radius">半徑值。</param>void DoBilateralFilter(double *WegightedImage, double *ImageWeight, double *Table, unsigned char *DownInput, unsigned char *Guide, int Width, int Height, int Radius){    int X, Y, XX, YY, MinX, MinY, MaxX, MaxY, Index, IndexXY, IndexXXYY;    double SumR, SumG, SumB, SumWeight, Weight;    for (Y = 0; Y < Height; Y++)    {        MinY = max(Y - Radius, 0);                //    防止訪問超出映像範圍的像素        MaxY = min(Y + Radius, Height - 1);        Index = Y * Width;        IndexXY = Y * Width * 3;        for (X = 0; X < Width; X++)        {            MinX = max(X - Radius, 0);            MaxX = min(X + Radius, Width - 1);            SumR = 0; SumG = 0; SumB = 0; SumWeight = 0;            for (YY = MinY; YY < MaxY; YY++)            {                IndexXXYY = YY * Width * 3 + MinX * 3;                for (XX = MinX; XX < MaxX; XX++)                {                    Weight = Table[MaxRGBEuclideanDistance(Guide + IndexXXYY, Guide + IndexXY)];    　//    Guide[XX,YY]  VS  Guide[X,Y]                    SumB += DownInput[IndexXXYY] * Weight;                                            //    DownInput[XX,YY,0]                    SumG += DownInput[IndexXXYY + 1] * Weight;                                        //    DownInput[XX,YY,1]                    SumR += DownInput[IndexXXYY + 2] * Weight;                                        //    DownInput[XX,YY,2]                    SumWeight += Weight;                    IndexXXYY += 3;                }            }            ImageWeight[Index] = SumWeight;                                                            //    ImageWeight[XX,YY]            WegightedImage[IndexXY] = SumB;                                                            //    WegightedImage[XX,YY,0]            WegightedImage[IndexXY + 1] = SumG;                                                        //    WegightedImage[XX,YY,1]            WegightedImage[IndexXY + 2] = SumR;                                                        //    WegightedImage[XX,YY,2]            Index++;            IndexXY += 3;        }    }}

　　當半徑比較小時，這種方式的實現似乎也是沒有辦法的事情了，有興趣的朋友可以搜尋下這篇文章：Reshuffling: A Fast Algorithm for Filtering with Arbitrary Kernels看看有沒有什麼課最佳化的地方了。

　　原始的論文已經提供了原始碼，我在其基礎上改成了我習慣的方式，並做了一個UI供朋友們測試，需要代碼的朋友請直接下載論文代碼並自己動手改寫，不要找我要代碼。

     測試程式介面即：http://files.cnblogs.com/Imageshop/HEBF.rar

   

    

    原文作者提供了這一組圖：

     

     有誰知道是那一篇論文有提到用雙邊濾波實現灰階映像的上色演算法，麻煩告之一下了，我覺得這個很有意思。

 

****************************基本上我不提供原始碼，但是我會盡量用文字把對應的演算法描述清楚或提供參考文檔************************

*************************************因為靠自己的努力和實踐寫出來的效果才真正是自己的東西，人一定要靠自己*******************

****************************作者： laviewpbt   時間： 2014.7.12    聯絡QQ:  33184777 轉載請保留本行資訊**********************