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