個人收藏了很多香港大學、香港科技大學以及香港中文大學裡專門搞影象研究一些博士的個人網站，一般會不定期的瀏覽他們的作品，最近在看楊慶雄的網點時，發現他又寫了一篇雙邊濾波的文章，並且配有原始碼，於是下載下來研讀了一番，這裡僅僅對一些過程做簡單的記錄，以防時間久了忘記。
關於楊慶雄的相關文章可見：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 轉載請保留本行資訊**********************