Bilateral Filtering(双边滤波)算法研究

1. 简介

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

2. 原理

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

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

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

其中的s为基于像素间相似程度的高斯权重，同样用来对结果进行单位化。对两者进行结合即可以得到基于空间距离、相似程度综合考量的Bilateral滤波：

上式中的单位化分部综合了两种高斯权重于一起而得到，其中的c与s计算可以详细描述如下：

且有

上述给出的表达式均是在空间上的无限积分，而在像素化的图像中当然无法这么做，而且也没必要如此做，因而在使用前需要对其进行离散化。而且也不需要对于每个局部像素从整张图像上进行加权操作，距离超过一定程度的像素实际上对当前的目标像素影响很小，可以忽略的。限定局部子区域后的离散化公就可以简化为如下形式：

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

左图为原始的噪声图像；中间为高斯采样的权重；右图为Bilateral采样的权重。从图中可以看出Bilateral加入了相似程度分部以后可以将源图像左侧那些跟当前像素差值过大的点给滤去，这样就很好地保持了边缘。为了更加形象地观察两者间的区别，使用Matlab将该图在两种不同方式下的高度图3D绘制出来，如下：

上述三图从左到右依次为：双边滤波，原始图像，高斯滤波。从高度图中可以明显看出Bilateral和Gaussian两种方法的区别，前者较好地保持了边缘处的梯度，而在高斯滤波中，由于其在边缘处的变化是线性的，因而就使用连累的梯度呈现出渐变的状态，而这表现在图像中的话就是边界的丢失（图像的示例可见于后述）。

3. 代码实现

有了上述理论以后实现Bilateral Filter就比较简单了，其实它也与普通的Gaussian Blur没有太大的区别。这里主要包括3部分的操作: 基于空间距离的权重因子生成；基于相似度的权重因子的生成；最终filter颜色的计算。

3.1 Spatial Weight

这就是通常的Gaussian Blur中使用的计算高斯权重的方法，其主要通过两个pixel之间的距离并使用如下公式计算而来：

其中的就表示两个像素间的距离，比如当前像素与其右边紧邻的一个像素之间的距离我们就可以用来计算，也即两个二维向量{0
, 0}以及{0 , 1}之间的欧氏距离。直接计算一个区域上的高斯权重并单位化后就可以进行高斯模糊了。

3.2 Similarity Weight

与基于距离的高斯权重计算类似，只不过此处不再根据两个pixel之间的空间距离，而是根据其相似程度（或者两个pixel的值之间的距离）。

其中的表示两个像素值之间的距离，可以直接使用其灰度值之间的差值或者RGB向量之间的欧氏距离。

3.3 Color Filtering

有了上述两部分所必需的权重因子之后，那么具体的双边滤波的实现即与普通的高斯滤波无异。主要部分代码如下述：

[cpp] view
plain copy

UCHAR3 BBColor(int posX , int posY)
{
int centerItemIndex = posY * picWidth4 + posX * 3 , neighbourItemIndex;
int weightIndex;
double gsAccumWeight = 0;
double accumColor = 0;
// 计算各个采样点处的Gaussian权重，包括closeness，similarity
for(int i = -number ; i <= number ; ++i)
{
for(int j = -number ; j <= number ; ++j)
{
weightIndex = (i + number) * (number * 2 + 1) + (j + number);
neighbourItemIndex = min(noiseImageHeight - 1 , max(0 , posY + j * radius)) * picWidth4 +
min(noiseImageWidth - 1 , max(0 , posX + i * radius)) * 3;
pCSWeight[weightIndex] = LookupGSWeightTable(pSrcDataBuffer[neighbourItemIndex] , pSrcDataBuffer[centerItemIndex]);
pCSWeight[weightIndex] = pGSWeight[weightIndex] * pGCWeight[weightIndex];
gsAccumWeight += pCSWeight[weightIndex];
}
}
// 单位化权重因子
gsAccumWeight = 1 / gsAccumWeight;
for(int i = -number ; i <= number ; ++i)
{
for(int j = -number ; j <= number ; ++j)
{
weightIndex = (i + number) * (number * 2 + 1) + (j + number);
pCSWeight[weightIndex] *= gsAccumWeight;
}
}
// 计算最终的颜色并返回
for(int i = -number ; i <= number ; ++i)
{
for(int j = -number ; j <= number ; ++j)
{
weightIndex = (i + number) * (number * 2 + 1) + (j + number);
neighbourItemIndex = min(noiseImageHeight - 1 , max(0 , posY + j * radius)) * picWidth4 +
min(noiseImageWidth - 1 , max(0 , posX + i * radius)) * 3;
accumColor += pSrcDataBuffer[neighbourItemIndex + 0] * pCSWeight[weightIndex];
}
}
return UCHAR3(accumColor , accumColor , accumColor);
}