技術(shù)的真相 | 基于雙邊濾波的磨皮算法及優(yōu)化

一．背景簡(jiǎn)介

現(xiàn)在視頻類應(yīng)用非常火熱，直播、美妝、醫(yī)美應(yīng)用層出不窮。用戶們?cè)谑褂眠@類應(yīng)用時(shí)都希望自己在屏幕上的樣子美美的，皮膚細(xì)膩光滑。本文就介紹一種實(shí)現(xiàn)簡(jiǎn)單、效果很好的磨皮算法以及對(duì)它的優(yōu)化思路。

二．磨皮算法

磨皮的基本方法

磨皮算法的核心就是讓皮膚區(qū)域的像素變得平滑，過渡自然，這樣皮膚上的瑕疵就不容易察覺了，皮膚看起來(lái)十分光滑。這里包括了兩部分算法：

1. 圖像濾波算法；

2. 皮膚區(qū)域檢測(cè)算法。

整體的流程如下：

皮膚區(qū)域檢測(cè)可參考技術(shù)的真相 | 從AR口紅試妝了解人工智能試妝技術(shù)，此文不再贅述。

濾波算法一般有兩類作用：一是讓畫面變得平滑甚至模糊，二是去除畫面中的噪點(diǎn)。在磨皮算法里主要的作用是過濾掉畫面中微小的瑕疵、讓顏色過渡更加平滑。其主要的思想是讓中心像素顏色和周圍像素顏色取加權(quán)均值，然后更新中心像素的顏色值。濾波窗口中的各像素權(quán)重有很多種選取方式，由此產(chǎn)生了多種濾波算法。

常見的濾波算法

常用的濾波方法有盒型濾波（如均值濾波）、統(tǒng)計(jì)排序?yàn)V波（如中值濾波）、高斯濾波和保邊濾波（如雙邊濾波）等。不同濾波算法產(chǎn)生的效果會(huì)有很大區(qū)別。其中一些對(duì)去除噪點(diǎn)效果很好，如中值濾波；一些產(chǎn)生平滑的畫面效果，如均值濾波、高斯濾波、保邊濾波等。

為什么要用保邊的濾波算法

對(duì)比上面兩張圖片，可以看到非常明顯的差異。

左圖使用的是高斯濾波，權(quán)重符合二維高斯分布，中間權(quán)重高，越往外權(quán)重越低，在各個(gè)方向上距中心同一距離的位置上的權(quán)重是一致的。這樣的卷積核簡(jiǎn)單粗暴，不管圖像顏色信息如何，統(tǒng)一對(duì)待，所以在輸出圖像上可見山石邊緣也變得非常模糊。

右圖使用的是雙邊濾波，卷積核在高斯核的基礎(chǔ)上增加了顏色差異的權(quán)重。對(duì)于顏色差異大的鄰近點(diǎn)，它的權(quán)重變小，對(duì)中心點(diǎn)的平滑影響減弱，邊緣可以被有效地保留下來(lái)。

人臉上有眉毛、眼睛、嘴等與皮膚顏色差異很大的區(qū)域，而且立體的五官上也會(huì)有明暗變化。這些邊緣如果使用高斯濾波會(huì)模糊成一團(tuán)。但是使用保邊濾波，可以清晰地保留下五官和明暗區(qū)域，防止變平、失真。

讓我們來(lái)看一下雙邊濾波和高斯濾波應(yīng)用在真實(shí)人像上的效果對(duì)比（點(diǎn)擊圖片可查看大圖）：

這是一張?jiān)紙D片

應(yīng)用高斯濾波算法（sigma=4.0），可以觀察到人物臉部變得朦朦朧朧，眼鏡邊緣、鼻子邊緣和嘴部邊緣都變得不清晰了。

應(yīng)用雙邊濾波算法（spaceSigma=4.0，colorSigma=18），可以觀察到人物五官基本上仍然是清晰，眼鏡邊緣也沒有什么變化。

可見在對(duì)人像磨皮的算法選擇上，使用雙邊濾波更為合理。

三．雙邊濾波算法

雙邊濾波的原理雙邊濾波算法由C. Tomasi和R. Manduchi在1998年提出，論文名稱是Bilateral Filtering for Gray and Color Images。

雙邊濾波組合了空間上鄰近程度的權(quán)重和像素顏色值相似程度的權(quán)重。公式如下：

其中，f是原始圖像，h是濾波后圖像，x是卷積核中心點(diǎn)，ξ是x的相鄰點(diǎn)。卷積核包括了兩部分：

空間權(quán)重c(ξ,x)衡量幾何空間上的鄰近程度（GeometricCloseness），顏色權(quán)重s(f(ξ),f(x))衡量色彩空間上的相似程度（Photometric Similarity）。

權(quán)重c和權(quán)重s都遵循高斯分布：離中心點(diǎn)的幾何距離越遠(yuǎn)，權(quán)重c越??；和中心點(diǎn)的顏色差異越大，權(quán)重s越小。

上圖從左到右依次為：輸入圖像（含有明顯的顏色變化邊界）；空間卷積核c；顏色卷積核s；疊加了空間卷積核和顏色卷積核的完整卷積核c x s；濾波后的圖像（保留了顏色邊界）

雙邊濾波的參數(shù)和效果

讓我們對(duì)比一下高斯濾波和雙邊濾波的效果。這兩種濾波器在OpenCV中都有實(shí)現(xiàn)，所以直接調(diào)用OpenCV的接口得到如下結(jié)果（點(diǎn)擊圖片可查看大圖）：

注：實(shí)驗(yàn)圖片使用經(jīng)典的lena圖像，512x512像素，tiff格式。

對(duì)比上面幾幅圖像，可以看到雙邊濾波可以明顯地去處原圖上的噪點(diǎn)，并且保留了相對(duì)清晰的物體輪廓。對(duì)比原圖和spaceSigma=4，colorSigma=32的雙邊濾波結(jié)果，人物臉上和肩膀上的斑點(diǎn)都消失了，皮膚看起來(lái)光滑水潤(rùn)。隨著參數(shù)改變，我們看到濾波結(jié)果發(fā)生了很大變化。

由于spaceSigma和colorSigma都是卷積核的方差。某一部分的方差變大，鐘形曲線變得平緩，各個(gè)權(quán)重的差異變小，說(shuō)明不強(qiáng)調(diào)這一部分的影響。

spaceSigma值越小，畫面越清晰，值越大，畫面越模糊，大到極限時(shí)，變?yōu)橹涤驗(yàn)V波（閾值化）。

colorSigma值越小，邊緣越清晰，值越大，邊緣越模糊，大到極限時(shí)，變?yōu)楦咚篂V波。

如果spaceSigma相對(duì)于colorSigma變小，說(shuō)明卷積核更強(qiáng)調(diào)距離中心點(diǎn)近的點(diǎn)。比如colorSigma=32，spaceSigma由16變?yōu)?時(shí)，可以看到顏色接近的區(qū)域內(nèi)部變得沒那么模糊了。

如果colorSigma相對(duì)于spaceSigma變小，說(shuō)明卷積核更強(qiáng)調(diào)與中心點(diǎn)顏色接近的點(diǎn)。比如spaceSigma=16，colorSigma由128變?yōu)?2時(shí)，可以看到邊緣沒那么模糊了。

注1：使用OpenCV時(shí)，如果不指定卷積核的尺寸，按照3通道sigma×6+1計(jì)算。因此計(jì)算此組圖片時(shí)，當(dāng)spaceSigma=4時(shí)，卷積核尺寸為25x25；當(dāng)spaceSigma=16時(shí)，卷積核尺寸為97x97。注2：PSNR（Peak Signal to Noise Ratio），是一種最普遍，最廣泛使用的評(píng)價(jià)圖像的客觀標(biāo)準(zhǔn)，單位是dB，其數(shù)值越大，失真越少。

雙邊濾波的缺點(diǎn)雙邊濾波的卷積核是非線性的，因此計(jì)算復(fù)雜度高。不同位置的卷積核不同，不能預(yù)先計(jì)算或者執(zhí)行FFT，計(jì)算起來(lái)比較費(fèi)時(shí)。因此很多前輩做出了很多嘗試，提出了一些優(yōu)化方法可以近似雙邊濾波的效果。

雙邊濾波的實(shí)現(xiàn)和優(yōu)化

0. 基本實(shí)現(xiàn)

最基本的實(shí)現(xiàn)方式就是按照上文中的公式，對(duì)每個(gè)像素進(jìn)行二維卷積計(jì)算。這里使用OpenCV實(shí)現(xiàn)基本的雙邊濾波算法，以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. 使用查找表減少計(jì)算在計(jì)算高斯函數(shù)的時(shí)候，因?yàn)楹兄笖?shù)運(yùn)算，運(yùn)算量比較大。同時(shí)，自變量部分都都是整數(shù)，所以可以預(yù)先計(jì)算好權(quán)重查找表，在卷積時(shí)直接找到對(duì)應(yīng)的權(quán)重值。權(quán)重表計(jì)算如下：

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. 使用可分離的卷積近似計(jì)算對(duì)于卷積核可以分解成一個(gè)行向量乘以一個(gè)列向量的形式時(shí)，例如：

就可以把一次二維卷積分解成兩次一維卷積，單個(gè)像素卷積操作的計(jì)算復(fù)雜度從  下降為O(N)。對(duì)于較大的卷積核而言，優(yōu)化幅度是相當(dāng)可觀的。實(shí)現(xiàn)時(shí)需要先對(duì)整幅圖像做一次一維卷積，把卷積結(jié)果保存到臨時(shí)圖像上，再對(duì)臨時(shí)圖像做第二次一維卷積。

對(duì)于雙邊濾波的權(quán)重s部分，雖然分離后的結(jié)果不完全一樣，但是誤差不大，可以以這種方法做近似計(jì)算。

3. 使用遞歸方法近似計(jì)算

R. Deriche在1993年提出了一種遞歸高斯濾波方法，論文名稱是Recursively implementating the Gaussian and its derivatives。在這種方法中，二維高斯卷積核被分解成兩個(gè)一維高斯卷積核相乘，每一次卷積計(jì)算包含6次乘法和6次加法，與卷積核的尺寸無(wú)關(guān)。計(jì)算復(fù)雜度下降為O(1)。對(duì)于大小超過7x7的卷積核，性能都可以有所提升，卷積核越大提升越明顯。

受到這篇論文的啟發(fā)，很多學(xué)者都在探索如何遷移這種方法到雙邊濾波算法中。比如Q. Yang在2012年提出了遞歸雙邊濾波方法，論文名稱是Recursive Bilateral Filtering。本實(shí)驗(yàn)使用了https://github.com/ufoym/recursive-bf提供的代碼。

以下是上述優(yōu)化方法的結(jié)果對(duì)比（點(diǎn)擊圖片可查看大圖）：

實(shí)驗(yàn)參數(shù)均為spaceSigma=4，colorSigma=32。

從PSNR數(shù)值來(lái)看，幾種優(yōu)化方法和標(biāo)準(zhǔn)算法區(qū)別不大。肉眼觀察，可見一些細(xì)微的差異，但整體效果接近。

4. 使用SIMD指令集在CPU上，還可以使用SIMD指令集做進(jìn)一步的加速，加速比例一般可以做到接近于數(shù)據(jù)并行數(shù)目。一般來(lái)說(shuō)，比如在64位寬的寄存器上使用8x8位數(shù)據(jù)進(jìn)行運(yùn)算，加速比可以達(dá)到6-8倍。受限于篇幅，本文不進(jìn)行詳細(xì)討論?？梢詤⒖糷ttps://github.com/Fig1024/OP_RBF提供的代碼。

三．參考文獻(xiàn)

• Bilateral Filtering for Gray and Color Images. C. Tomasi, R. ManduchiRecursively implementating the Gaussian and its derivatives. R. DericheRecursive Bilateral Filtering. Q. Yang