学习《Hardware-Efficient Bilateral Filtering for Stereo Matching》一文笔记。

　　个人收藏了很多中国香港大学、中国香港科技大学以及中国香港中文大学里专门搞图像研究一些博士的个人网站，一般会不定期的浏览他们的作品，最近在看杨庆雄的网点时，发现他又写了一篇双边滤波的文章，并且配有源代码，于是下载下来研读了一番，这里仅仅对一些过程做简单的记录，以防时间久了忘记。

　　关于杨庆雄的相关文章可见：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 转载请保留本行信息**********************