技术的真相 | 基于双边滤波的磨皮算法及优化

一．背景简介

现在视频类应用非常火热，直播、美妆、医美应用层出不穷。用户们在使用这类应用时都希望自己在屏幕上的样子美美的，皮肤细腻光滑。本文就介绍一种实现简单、效果很好的磨皮算法以及对它的优化思路。

二．磨皮算法

磨皮的基本方法

磨皮算法的核心就是让皮肤区域的像素变得平滑，过渡自然，这样皮肤上的瑕疵就不容易察觉了，皮肤看起来十分光滑。这里包括了两部分算法：

1. 图像滤波算法；

2. 皮肤区域检测算法。

整体的流程如下：

皮肤区域检测可参考技术的真相 | 从AR口红试妆了解人工智能试妆技术，此文不再赘述。

滤波算法一般有两类作用：一是让画面变得平滑甚至模糊，二是去除画面中的噪点。在磨皮算法里主要的作用是过滤掉画面中微小的瑕疵、让颜色过渡更加平滑。其主要的思想是让中心像素颜色和周围像素颜色取加权均值，然后更新中心像素的颜色值。滤波窗口中的各像素权重有很多种选取方式，由此产生了多种滤波算法。

常见的滤波算法

常用的滤波方法有盒型滤波（如均值滤波）、统计排序滤波（如中值滤波）、高斯滤波和保边滤波（如双边滤波）等。不同滤波算法产生的效果会有很大区别。其中一些对去除噪点效果很好，如中值滤波；一些产生平滑的画面效果，如均值滤波、高斯滤波、保边滤波等。

为什么要用保边的滤波算法

对比上面两张图片，可以看到非常明显的差异。

左图使用的是高斯滤波，权重符合二维高斯分布，中间权重高，越往外权重越低，在各个方向上距中心同一距离的位置上的权重是一致的。这样的卷积核简单粗暴，不管图像颜色信息如何，统一对待，所以在输出图像上可见山石边缘也变得非常模糊。

右图使用的是双边滤波，卷积核在高斯核的基础上增加了颜色差异的权重。对于颜色差异大的邻近点，它的权重变小，对中心点的平滑影响减弱，边缘可以被有效地保留下来。

人脸上有眉毛、眼睛、嘴等与皮肤颜色差异很大的区域，而且立体的五官上也会有明暗变化。这些边缘如果使用高斯滤波会模糊成一团。但是使用保边滤波，可以清晰地保留下五官和明暗区域，防止变平、失真。

让我们来看一下双边滤波和高斯滤波应用在真实人像上的效果对比（点击图片可查看大图）：

这是一张原始图片

应用高斯滤波算法（sigma=4.0），可以观察到人物脸部变得朦朦胧胧，眼镜边缘、鼻子边缘和嘴部边缘都变得不清晰了。

应用双边滤波算法（spaceSigma=4.0，colorSigma=18），可以观察到人物五官基本上仍然是清晰，眼镜边缘也没有什么变化。

可见在对人像磨皮的算法选择上，使用双边滤波更为合理。

三．双边滤波算法

双边滤波的原理

双边滤波算法由C. Tomasi和R. Manduchi在1998年提出，论文名称是Bilateral Filtering for Gray and Color Images。

双边滤波组合了空间上邻近程度的权重和像素颜色值相似程度的权重。公式如下：

其中，f是原始图像，h是滤波后图像，x是卷积核中心点，ξ是x的相邻点。卷积核包括了两部分：

空间权重c(ξ,x)衡量几何空间上的邻近程度（GeometricCloseness），颜色权重s(f(ξ),f(x))衡量色彩空间上的相似程度（Photometric Similarity）。

权重c和权重s都遵循高斯分布：离中心点的几何距离越远，权重c越小；和中心点的颜色差异越大，权重s越小。

上图从左到右依次为：输入图像（含有明显的颜色变化边界）；空间卷积核c；颜色卷积核s；叠加了空间卷积核和颜色卷积核的完整卷积核c x s；滤波后的图像（保留了颜色边界）

双边滤波的参数和效果

让我们对比一下高斯滤波和双边滤波的效果。这两种滤波器在OpenCV中都有实现，所以直接调用OpenCV的接口得到如下结果（点击图片可查看大图）：

注：实验图片使用经典的lena图像，512x512像素，tiff格式。

对比上面几幅图像，可以看到双边滤波可以明显地去处原图上的噪点，并且保留了相对清晰的物体轮廓。对比原图和spaceSigma=4，colorSigma=32的双边滤波结果，人物脸上和肩膀上的斑点都消失了，皮肤看起来光滑水润。随着参数改变，我们看到滤波结果发生了很大变化。

由于spaceSigma和colorSigma都是卷积核的方差。某一部分的方差变大，钟形曲线变得平缓，各个权重的差异变小，说明不强调这一部分的影响。

spaceSigma值越小，画面越清晰，值越大，画面越模糊，大到极限时，变为值域滤波（阈值化）。

colorSigma值越小，边缘越清晰，值越大，边缘越模糊，大到极限时，变为高斯滤波。

如果spaceSigma相对于colorSigma变小，说明卷积核更强调距离中心点近的点。比如colorSigma=32，spaceSigma由16变为4时，可以看到颜色接近的区域内部变得没那么模糊了。

如果colorSigma相对于spaceSigma变小，说明卷积核更强调与中心点颜色接近的点。比如spaceSigma=16，colorSigma由128变为32时，可以看到边缘没那么模糊了。

注1：

使用OpenCV时，如果不指定卷积核的尺寸，按照3通道sigma×6+1计算。因此计算此组图片时，当spaceSigma=4时，卷积核尺寸为25x25；当spaceSigma=16时，卷积核尺寸为97x97。

注2：

PSNR（Peak Signal to Noise Ratio），是一种最普遍，最广泛使用的评价图像的客观标准，单位是dB，其数值越大，失真越少。

双边滤波的缺点

双边滤波的卷积核是非线性的，因此计算复杂度高。不同位置的卷积核不同，不能预先计算或者执行FFT，计算起来比较费时。因此很多前辈做出了很多尝试，提出了一些优化方法可以近似双边滤波的效果。

双边滤波的实现和优化

0. 基本实现

最基本的实现方式就是按照上文中的公式，对每个像素进行二维卷积计算。这里使用OpenCV实现基本的双边滤波算法，以RGB3通道为例：

void BilateralBlur(const Mat &srcImage, Mat &dstImage, int kernelSize, double sigmaColor, double sigmaSpace) {
    int height = srcImage.rows;
    int width = srcImage.cols;
    int channel = srcImage.channels();
    if (dstImage.rows == 0 || dstImage.cols == 0)
        dstImage = Mat::zeros(srcImage.size(), srcImage.type());
    double ct = -0.5 / (sigmaColor * sigmaColor);
    double st = -0.5 / (sigmaSpace * sigmaSpace);
    int radius = (kernelSize - 1) / 2;
    double w;
    double sw;
    double cw;
    for (int row = 0; row < height; ++row) {
        for (int col = 0; col < width; ++col) {
            double sumr = 0, sumg = 0, sumb = 0, sumw = 0;
            Vec3b bgr0 = srcImage.at<Vec3b>(row, col);
            for (int n = -radius; n <= radius; n += step) {
                int y = clamp(row + n, 0, height - 1);
                for (int m = -radius; m <= radius; m += step) {
                    int x = clamp(col + m, 0, width - 1);
 
                    Vec3b bgr = srcImage.at<Vec3b>(y, x);
                    sw = exp((m * m + n * n) * st);
                    double c = abs(bgr0[0] - bgr[0]) + abs(bgr0[1] - bgr[1]) + abs(bgr0[2] - bgr[2]); // sum of difference of each channel
                    cw = exp(c * c * ct);
                    double w = sw * cw;
                    sumb += b * w;
                    sumg += g * w;
                    sumr += r * w;
                    sumw += w;
                }
            }
            dstImage.at<Vec3b>(row, col) = Vec3b(saturate_cast<uchar>(sumb / sumw),
                                                 saturate_cast<uchar>(sumg / sumw),
                                                 saturate_cast<uchar>(sumr / sumw));
        }
    }
}

1. 使用查找表减少计算

在计算高斯函数的时候，因为含有指数运算，运算量比较大。同时，自变量部分都都是整数，所以可以预先计算好权重查找表，在卷积时直接找到对应的权重值。权重表计算如下：

void BilateralBlur(const Mat &srcImage, Mat &dstImage, int kernelSize, double sigmaColor, double sigmaSpace) {
    int height = srcImage.rows;
    int width = srcImage.cols;
    int channel = srcImage.channels();
    if (dstImage.rows == 0 || dstImage.cols == 0)
        dstImage = Mat::zeros(srcImage.size(), srcImage.type());
    double ct = -0.5 / (sigmaColor * sigmaColor);
    double st = -0.5 / (sigmaSpace * sigmaSpace);
    int radius = (kernelSize - 1) / 2;
    double w;
    double sw;
    double cw;
    for (int row = 0; row < height; ++row) {
        for (int col = 0; col < width; ++col) {
            double sumr = 0, sumg = 0, sumb = 0, sumw = 0;
            Vec3b bgr0 = srcImage.at<Vec3b>(row, col);
            for (int n = -radius; n <= radius; n += step) {
                int y = clamp(row + n, 0, height - 1);
                for (int m = -radius; m <= radius; m += step) {
                    int x = clamp(col + m, 0, width - 1);


                    Vec3b bgr = srcImage.at<Vec3b>(y, x);
                    sw = exp((m * m + n * n) * st);
                    double c = abs(bgr0[0] - bgr[0]) + abs(bgr0[1] - bgr[1]) + abs(bgr0[2] - bgr[2]); // sum of difference of each channel
                    cw = exp(c * c * ct);
                    double w = sw * cw;
                    sumb += b * w;
                    sumg += g * w;
                    sumr += r * w;
                    sumw += w;
                }
            }
            dstImage.at<Vec3b>(row, col) = Vec3b(saturate_cast<uchar>(sumb / sumw),
                                                 saturate_cast<uchar>(sumg / sumw),
                                                 saturate_cast<uchar>(sumr / sumw));
        }
    }
}

2. 使用可分离的卷积近似计算

对于卷积核可以分解成一个行向量乘以一个列向量的形式时，例如：

就可以把一次二维卷积分解成两次一维卷积，单个像素卷积操作的计算复杂度从   下降为O(N)。对于较大的卷积核而言，优化幅度是相当可观的。实现时需要先对整幅图像做一次一维卷积，把卷积结果保存到临时图像上，再对临时图像做第二次一维卷积。

对于双边滤波的权重s部分，虽然分离后的结果不完全一样，但是误差不大，可以以这种方法做近似计算。

3. 使用递归方法近似计算

R. Deriche在1993年提出了一种递归高斯滤波方法，论文名称是Recursively implementating the Gaussian and its derivatives。在这种方法中，二维高斯卷积核被分解成两个一维高斯卷积核相乘，每一次卷积计算包含6次乘法和6次加法，与卷积核的尺寸无关。计算复杂度下降为O(1)。对于大小超过7x7的卷积核，性能都可以有所提升，卷积核越大提升越明显。

受到这篇论文的启发，很多学者都在探索如何迁移这种方法到双边滤波算法中。比如Q. Yang在2012年提出了递归双边滤波方法，论文名称是Recursive Bilateral Filtering。

本实验使用了https://github.com/ufoym/recursive-bf提供的代码。

以下是上述优化方法的结果对比（点击图片可查看大图）：

实验参数均为spaceSigma=4，colorSigma=32。

从PSNR数值来看，几种优化方法和标准算法区别不大。肉眼观察，可见一些细微的差异，但整体效果接近。

4. 使用SIMD指令集

在CPU上，还可以使用SIMD指令集做进一步的加速，加速比例一般可以做到接近于数据并行数目。一般来说，比如在64位宽的寄存器上使用8x8位数据进行运算，加速比可以达到6-8倍。受限于篇幅，本文不进行详细讨论。可以参考https://github.com/Fig1024/OP_RBF提供的代码。

三．参考文献

• Bilateral Filtering for Gray and Color Images. C. Tomasi, R. Manduchi

• Recursively implementating the Gaussian and its derivatives. R. Deriche

• Recursive Bilateral Filtering. Q. Yang