双边滤波

高斯滤波是最常用的图像去噪方法之一，它能很好地滤除掉图像中随机出现的高斯噪声，但是在之前的博客中提到过，高斯滤波是一种低通滤波（有兴趣的点击这里，查看之前的博客），它在滤除图像中噪声信号的同时，也会对图像中的边缘信息进行平滑，表现出来的结果就是图像变得模糊，如下图所示；

高斯滤波之所以会导致图像变得模糊，是因为它在滤波过程中只关注了位置信息；例如，以 qq q q 为中心的窗口中，某一点 pp p p 在高斯滤波过程中的权重的计算方法如下式；
G(p)=12πσ2e−∥p−q∥22σ2G(p)=12πσ2e−∥p−q∥22σ2 G(p)=\frac{1}{2\pi \sigma ^{2}} e^{-\frac{\left \| p-q \right \|^{2}}{2\sigma ^{2}}} G(p)=2πσ21e−2σ2∥p−q∥2
即在滤波窗口内，距离中心点越近的点的权重越大；这种只关注距离的思想在某些情况下是可行的，例如在平坦的区域，距离越近的区域其像素分布也越相近，自然地，这些点的像素值对滤波中心点的像素值更有参考价值。但是在像素值出现跃变的边缘区域，这种方法会适得其反，损失掉有用的边缘信息。此时就出现了一类算法——边缘保护滤波方法，双边滤波就是最常用的边缘保护滤波方法（另一种常用来与双边滤波对比的边缘保护滤波方法——引导滤波）。

双边滤波的思想很简单，在高斯滤波的基础上加入了像素值权重项，也就是说既要考虑距离因素，也要考虑像素值差异的影响，像素值越相近，权重越大。将像素值权重表示为 GrGr G_{r} Gr ，空间距离权重表示为 GsGs G_{s} Gs。

Gs=exp(−∥p−q∥22σ2s)Gs=exp(−∥p−q∥22σs2) G_{s}=exp(-\frac{\left \| p-q \right \|^{2}}{2\sigma _{s}^{2}}) Gs=exp(−2σs2∥p−q∥2)
Gr=exp(−∥Ip−Iq∥22σ2r)Gr=exp(−∥Ip−Iq∥22σr2) G_{r}=exp(-\frac{\left \| I_{p}-I_{q} \right \|^{2}}{2\sigma _{r}^{2}}) Gr=exp(−2σr2∥Ip−Iq∥2)

那么整个滤波器可以表示为BFBF BF BF，那么滤波结果为；
BF=1Wq∑p∈SGs(p)Gr(p)∗Ip=1Wq∑p∈Sexp(−∥p−q∥22σ2s)exp(−∥Ip−Iq∥22σ2r)∗IpBF=1Wq∑p∈SGs(p)Gr(p)∗Ip=1Wq∑p∈Sexp(−∥p−q∥22σs2)exp(−∥Ip−Iq∥22σr2)∗Ip \begin{aligned} BF&=\frac{1}{W_{q}}\sum_{p\in S}G_{s}(p)G_{r}(p)*I_{p}\\&=\frac{1}{W_{q}}\sum_{p\in S}exp(-\frac{\left \| p-q \right \|^{2}}{2\sigma _{s}^{2}})exp(-\frac{\left \| I_{p}-I_{q} \right \|^{2}}{2\sigma _{r}^{2}})*I_{p}\end{aligned} BF=Wq1p∈S∑Gs(p)Gr(p)∗Ip=Wq1p∈S∑exp(−2σs2∥p−q∥2)exp(−2σr2∥Ip−Iq∥2)∗Ip
其中WqWq {W_{q}} Wq 为滤波窗口内每个像素值的权重和，用于权重的归一化；
Wq=∑p∈SGs(p)Gr(p)=∑p∈Sexp(−∥p−q∥22σ2s)exp(−∥Ip−Iq∥22σ2r)Wq=∑p∈SGs(p)Gr(p)=∑p∈Sexp(−∥p−q∥22σs2)exp(−∥Ip−Iq∥22σr2) {W_{q}}=\sum_{p\in S}G_{s}(p)G_{r}(p)=\sum_{p\in S}exp(-\frac{\left \| p-q \right \|^{2}}{2\sigma _{s}^{2}})exp(-\frac{\left \| I_{p}-I_{q} \right \|^{2}}{2\sigma _{r}^{2}}) Wq=p∈S∑Gs(p)Gr(p)=p∈S∑exp(−2σs2∥p−q∥2)exp(−2σr2∥Ip−Iq∥2)
在平坦区域，滤波器中每个像素点的GrGr G_{r} Gr 值相近，空间距离权重GsGs G_{s} Gs主导滤波效果。在边缘区域，边缘同侧的GrGr G_{r} Gr 值相近，且远大于边缘另一侧的GrGr G_{r} Gr 值，此时另一侧的像素点的权重对滤波结果几乎不影响，边缘信息得到保护。表现出了一定的自适应性。

在平坦区域出现噪声点时，噪声点周围的信号的权值都很小，在进行归一化后，这些权值得到了提升，因此对噪声点也有滤除效果。

双边滤波OpenCV实现

在OpenCV中已经有编辑好的双边滤波的函数；

void bilateralFilter( InputArray src, OutputArray dst, int d,double sigmaColor, double sigmaSpace,int borderType = BORDER_DEFAULT );

函数中的参数依次表示src：输入图像，dst: 输出图像，d：滤波窗口的直径（函数注释中使用的是Diameter，那么很可能函数中选取的窗口是圆形窗口），sigmaColor：像素值域方差，sigmaSpace：空间域方差，以及边缘处理方式。

根据算法原理，我编写了自己的双边滤波函数；其中的核心部分如下；

//-------------------------------------------------------------
//作者：不用先生，2018.11.26
//自实现的图像双边滤波算法
//bilateral.cpp
//-------------------------------------------------------------
for (int i = 0; i < row; i++)    //对每一个点进行处理
{for (int j = 0; j < col; j++){double weightSum = 0;double filterValue = 0;for (int row_d = -(d / 2); row_d <= (d / 2); row_d++)   //以图像中的一点为中心，d为边长的方形区域内进行计算{for (int col_d = -(d / 2); col_d <= (d / 2); col_d++){double distance_Square = row_d*row_d + col_d*col_d;double value_Square = pow((scr.at<uchar>(i, j) - copyBorder_dst.at<uchar>(i + (d / 2) + row_d, j + (d / 2) + col_d)), 2);double weight = exp(-1 * (distance_Square / (2 * sigmaSpace*sigmaSpace) + value_Square / (2 * sigmaColor*sigmaColor)));weightSum += weight;               //求滤波窗口内的权重和，用于归一化；filterValue += (weight*copyBorder_dst.at<uchar>(i + (d / 2) + row_d, j + (d / 2) + col_d));}}dst.at<uchar>(i, j) = filterValue / weightSum;}
}

从代码中可以看出，我选择的是一个方形区域而不是圆形区域进行滤波计算，大家如果觉得欠妥，在自己编写代码时可以修改成圆形区域，也很简单。完整代码，可以到这里下载，或者直接在评论区留下自己的邮箱，我看到留言后会尽快发的您的邮箱。

测试结果

用自己编写的代码、OpenCV自带的双边滤波函数以及高斯滤波函数分别对一张彩色图像今进行处理，并对比了结果。

这里出现了很奇怪的现象，所有输入参数相同时（我选取参数的是d=21；sigmaColor=10；sigmaSpace = 10；），我的函数计算的结果比OpenCV自带的 bilateralFilter 函数的处理结果更加平滑。当把 bilateralFilter 函数的sigmaColor和sigmaSpace参数调整为两到三倍的时候， bilateralFilter 函数的结果与我的计算结果相近。这种差异性的结果不应该是窗口选取方式不同导致的。具体原因，我觉得可能与 bilateralFilter 函数中方差的应用有关。但是总的来说，两个函数都对图像有一定的平滑效果，并很好地保护了边缘信息。

本博客中只讲述了双边滤波最基础的原理以及实现过程，这种计算方法比较耗时，后续的人们做了一系列的工作，用于提升双边滤波的计算效率。大家感兴趣的可以去了解一下。

参考

Tomasi C, Manduchi R. Bilateral filtering for gray and color images[C]//Computer Vision, 1998. Sixth International Conference on. IEEE, 1998: 839-846.
He K, Sun J, Tang X. Guided image filtering[J]. IEEE transactions on pattern analysis & machine intelligence, 2013 (6): 1397-1409.

已完。。