双边滤波

1. 简介

图像平滑是一个重要的操作，而且有多种成熟的算法。这里主要简单介绍一下Bilateral方法（双边滤波）。Bilateral blur相对于传统的高斯blur来说很重要的一个特性即可可以保持边缘（Edge Perseving），这个特点对于一些图像模糊来说很有用。一般的高斯模糊在进行采样时主要考虑了像素间的空间距离关系，但是却并没有考虑像素值之间的相似程度，因此这样我们得到的模糊结果通常是整张图片一团模糊。Bilateral blur的改进就在于在采样时不仅考虑像素在空间距离上的关系，同时加入了像素间的相似程度考虑，因而可以保持原始图像的大体分块进而保持边缘。在于游戏引擎的post blur算法中，bilateral blur常常被用到，比如对SSAO的降噪。

 

双边滤波（Bilateral filter）是一种可以保边去噪的滤波器。之所以可以达到此去噪效果，是因为滤波器是由两个函数构成。一个函数是由几何空间距离决定滤波器系数。另一个由像素差值决定滤波器系数。可以与其相比较的两个filter：高斯低通滤波器(http://en.wikipedia.org/wiki/Gaussian_filter)和α-截尾均值滤波器（去掉百分率为α的最小值和最大之后剩下像素的均值作为滤波器），后文中将结合公式做详细介绍。

 

2.原理

滤波算法中，目标点上的像素值通常是由其所在位置上的周围的一个小局部邻居像素的值所决定。在2D高斯滤波中的具体实现就是对周围的一定范围内的像素值分别赋以不同的高斯权重值，并在加权平均后得到当前点的最终结果。而这里的高斯权重因子是利用两个像素之间的空间距离（在图像中为2D）关系来生成。通过高斯分布的曲线可以发现，离目标像素越近的点对最终结果的贡献越大，反之则越小。其公式化的描述一般如下所述：

其中的c即为基于空间距离的高斯权重，而用来对结果进行单位化。

高斯滤波在低通滤波算法中有不错的表现，但是其却有另外一个问题，那就是只考虑了像素间的空间位置上的关系，因此滤波的结果会丢失边缘的信息。这里的边缘主要是指图像中主要的不同颜色区域（比如蓝色的天空，黑色的头发等），而Bilateral就是在Gaussian blur中加入了另外的一个权重分部来解决这一问题。

 

像素化的图像中当然无法进行无线积分，而且也没必要如此做，因而在使用前需要对其进行离散化。而且也不需要对于每个局部像素从整张图像上进行加权操作，距离超过一定程度的像素实际上对当前的目标像素影响很小，可以忽略的。限定局部子区域后的离散化公就可以简化，即双边滤波器中，输出像素的值依赖于邻域像素的值的加权组合g(i,j)为：

其中权重系数w(i,j,k,l)取决于定义域核和值域核的乘积。

 

定义域核d(i,j,k,l)为：

值域核r(i,j,k,l)为：

那么权重系数w(i,j,k,l)为定义域核和值域核的乘积，即：

 

双边滤波同时考虑了空间域与值域的差别，而Gaussian Filter和α均值滤波分别只考虑了空间域和值域差别。

 

 

 

上述理论公式就构成了Bilateral滤波实现的基础。为了直观地了解高斯滤波与双边滤波的区别，我们可以从下列图示中看出依据。假设目标源图像为下述左右区域分明的带有噪声的图像（由程序自动生成），蓝色框的中心即为目标像素所在的位置，那么当前像素处所对应的高斯权重与双边权重因子3D可视化后的形状如后边两图所示：

 

  左图为原始的噪声图像；中间为高斯采样的权重；右图为Bilateral采样的权重。从图中可以看出Bilateral加入了相似程度分部以后可以将源图像左侧那些跟当前像素差值过大的点给滤去，这样就很好地保持了边缘。

 

3. 代码实现

在opencv中双边滤波是平滑函数cvSmooth中的一个子方法。

void cvSmooth( const void* srcarr, void* dstarr, int smooth_type,   int param1, int param2, double param3, double param4 )

{

      cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;

 

    CV_Assert( dst.size() == src.size() && (smooth_type == CV_BLUR_NO_SCALE || dst.type() == src.type()) );

    if( param2 <= 0 )

        param2 = param1;

    if( smooth_type == CV_BLUR || smooth_type == CV_BLUR_NO_SCALE )

        cv::boxFilter( src, dst, dst.depth(), cv::Size(param1, param2), cv::Point(-1,-1),smooth_type == CV_BLUR, cv::BORDER_REPLICATE);

    else if( smooth_type == CV_GAUSSIAN )

        cv::GaussianBlur( src, dst, cv::Size(param1, param2), param3, param4, cv::BORDER_REPLICATE );

    else if( smooth_type == CV_MEDIAN )

        cv::medianBlur( src, dst, param1 );

    else

        cv::bilateralFilter( src, dst, param1, param3, param4, cv::BORDER_REPLICATE ); 

}

 

其中bilateralFilter即是双边滤波的函数。

void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,  double sigmaColor, double sigmaSpace, int borderType )

{

    Mat src = _src.getMat();

    _dst.create( src.size(), src.type() );

    Mat dst = _dst.getMat();

 

    if( src.depth() == CV_8U )

        bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace,borderType );

    else if( src.depth() == CV_32F )

        bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace,borderType );

    else

        CV_Error(CV_StsUnsupportedFormat,"only for 8u and 32f images");

}

 

 

下面只列出bilateralFilter_8u函数实现。

static void bilateralFilter_8u(const Mat& src, Mat& dst, int d,  double sigma_color, double sigma_space,int borderType)

{

    int cn = src.channels();

    int i, j, maxk, radius;

    Size size = src.size();

    CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && src.type() == dst.type() && src.size() == dst.size() &&  src.data != dst.data );

 

    if( sigma_color <= 0 )        sigma_color = 1;

    if( sigma_space <= 0 )        sigma_space = 1;

    double gauss_color_coeff = -0.5/(sigma_color*sigma_color);

    double gauss_space_coeff = -0.5/(sigma_space*sigma_space);

 

    if( d <= 0 )   

            radius = cvRound(sigma_space*1.5);

    else           

           radius = d/2;

    radius = MAX(radius, 1);

    d = radius*2 + 1;

 

    Mat temp;

    copyMakeBorder(src,temp,radius,radius,radius,radius,borderType);

    vector<float> _color_weight(cn*256);

    vector<float> _space_weight(d*d);

    vector<int> _space_ofs(d*d);

    float* color_weight = &_color_weight[0];

    float* space_weight = &_space_weight[0];

    int* space_ofs = &_space_ofs[0];

 

    // initialize color-related bilateral filter coefficients

    for( i = 0; i < 256*cn; i++ )

        color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);

 

    // initialize space-related bilateral filter coefficients

    for( i = -radius, maxk = 0; i <= radius; i++ )

    {

        j = -radius;

        for( ;j <= radius; j++ )

        {

            double r = std::sqrt((double)i*i + (double)j*j);

            if( r > radius )   

                 continue;           

            space_weight[maxk]=(float)std::exp( r*r*gauss_space_coeff);

            space_ofs[maxk++] = (int)(i*temp.step + j*cn);

        }

    }

    BilateralFilter_8u_Invoker body(dst, temp, radius, maxk, space_ofs, space_weight, color_weight);

    parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));

}

 

【注】上面bilateralFilter_8u函数计算了space_weight和color_weight，其中space_weight就是所谓的定义域核d(i,j,k,l)。

而color_weight并不是上文提到的值域核,而是按照差值大小形成的值域核。因此这个值会在BilateralFilter_8u_Invoker body中的操作符重载函数中继续计算。下面将列出操作符重载函数。

 

virtual void operator() (const Range& range) const  

{

    int i, j, cn = dest->channels(), k;

    size size = dest->size();

 

    for( i = range.start; i < range.end; i++ ) 

    {

        const uchar* sptr = temp->ptr(i+radius) + radius*cn;

        uchar* dptr = dest->ptr(i);

 

        if( cn == 1 )

        {

            for( j = 0; j < size.width; j++ )  

            {

                float sum = 0, wsum = 0;

                int val0 = sptr[j];

                k = 0;

                for( ; k < maxk; k++ )

                {

                    int val = sptr[j + space_ofs[k]];

                    float w = space_weight[k]*color_weight[std::abs(val - val0)];

                    sum += val*w;

                    wsum += w;

                }

                dptr[j] = (uchar)cvRound(sum/wsum);

            }

        }

        else

        {

            assert( cn == 3 );

            for( j = 0; j < size.width*3; j += 3 )  

            {

               float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;

               int b0=sptr[j],g0=sptr[j+1],r0=sptr[j+2];

               k = 0;

               for( ; k < maxk; k++ )

               {

                  const uchar* sptr_k = sptr + j + space_ofs[k];

                   int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];

                   float w = space_weight[k]*color_weight[std::abs(b - b0) +std::abs(g - g0) + std::abs(r - r0)];

                    sum_b += b*w; sum_g += g*w; sum_r += r*w;

                    wsum += w;

                }

                wsum = 1.f/wsum;

                b0 = cvRound(sum_b*wsum);

                g0 = cvRound(sum_g*wsum);

                r0 = cvRound(sum_r*wsum);

                dptr[j] = (uchar)b0;

               dptr[j+1] = (uchar)g0;

               dptr[j+2] = (uchar)r0;

            }

        }

    }

}

 

【注】上面重载函数中，注意一下代码：

int val0 = sptr[j];

int val = sptr[j + space_ofs[k]];

 float w = space_weight[k]* color_weight[std::abs(val - val0)];

 

其中w是上文的权重系数w(i,j,k,l)。

space_weight是在bilateralFilter_8u函数中计算的定义域核d(i,j,k,l)，而color_weight使用abs(val - val0)继续计算后，便是上文中的值域核。最后将每个空间上的对应系数和相应差值的对应系数相乘，便是权重系数。