高斯滤波，双边滤波，肤色检测

import numpy as npfrom scipy.misc import imreadimport matplotlib.pyplot as plt%matplotlib inlineplt.rcParams['figure.figsize'] = (12.0, 9.0) # set default size of plotsplt.rcParams['image.interpolation'] = 'nearest'plt.rcParams['image.cmap'] = 'gray'

def imshow_noax(img, normalize=True):    """ Tiny helper to show images as uint8 and remove axis labels """    if normalize:        img_max, img_min = np.max(img), np.min(img)        img = 255.0 * (img - img_min) / (img_max - img_min)    plt.imshow(img.astype('uint8'))    plt.gca().axis('off')

高斯滤波

在进行数学仿真或者误差评估时，往往认为传感器所引入的噪声服从正态分布（高斯白噪声）。为了消除这些噪声，引入了高斯滤波。其滤波的思路：对高斯函数进行离散化，以离散点上的高斯函数值为权值，对信号做一定范围邻域内的加权平均。对于图像（灰度矩阵），我们采用二维高斯滤波滤除去图像的高斯噪声。

G(x,y)=12πσ2e−(x−x0)2+(y−y0)22σ2

σ越大，高斯滤波的频带越宽，平滑程度高。
高斯滤波器的工作原理是用像素领域的加权值来代替该点的像素值。而每一领域像素点权值是随该点与中心点的距离单调增减的。

H(i,j)=12πσ2ei2+j22σ2,i,j∈[−D/2,D/2]

D为高斯核的维度。在opencv与sift中的源码中，D=1+2∗(int)(3.0∗σ)，因为在i(j)∈[i(j）−3σ,i(j)+3σ]的概率大于99%,所以高斯核的大小与sigma的取值有关。实际上，高斯核大小也可以手动选择，并没有严格的要求。在代码中，我们引入”truncate“作为高斯核大小的缩放因子,并将高斯核归一化。

H¯(i,j)=H(i,j)∑D∑DH(i0,j0)

def gaussian_kernel(sigma, truncate=3.0):    """    Return Gaussian that truncates at the given number of standard deviations.     """    sigma = float(sigma)    radius = int(truncate * sigma + 0.5)    x, y = np.mgrid[-radius:radius+1, -radius:radius+1]    sigma = sigma**2    G = np.exp(-0.5 * (np.power(x, 2) + np.power(y, 2)) / sigma)    G = G / np.sum(G)    return G

def convolve(src, F):    H, W = src.shape    HH, WW = F.shape    assert(H > HH)    assert(W > WW)    hh = HH / 2    ww = WW / 2    dst = np.zeros((H, W))    src_pad = np.pad(src, (hh, ww), 'symmetric')    for i in range(H):        for j in range(W):            overlap = src_pad[i : i + HH, j: j + WW]            average = np.sum(F * overlap)            dst[i, j] = average        return dst

二维高斯函数卷积的另一种方法：分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积。因此，高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长。

# 灰度图img = imread('test.jpg', flatten=True)f= gaussian_kernel(3)dst = convolve(img, f)plt.subplot(2, 2, 1)imshow_noax(img, normalize=False)plt.title('Original image')plt.subplot(2, 2, 2)imshow_noax(dst, normalize=False)plt.title('Gaussian image')

<matplotlib.text.Text at 0x7f0310613dd0>

双边滤波

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

ω(i,j,k,l)=exp(−(i−k)2+(j−l)22σ2d−||f(i,j)−f(k,l)||2σ2r)

def bflitGray(img, sigma_d, sigma_r):    H, W = img.shape    sigma_d = float(sigma_d)    radius = int(3.0 * sigma_d + 0.5)    R = 2 * radius + 1    dst = np.zeros((H, W))    img_pad = np.pad(img, (radius, radius), 'symmetric')    x, y = np.mgrid[-radius: radius+1, -radius: radius+1]    G = np.exp(-0.5 * (np.power(x, 2) + np.power(y, 2)) / sigma_d**2)    G = G / np.sum(G)    for i in range(H):        for j in range(W):            overlap = img_pad[i : i + R, j: j + R]            A = img[i, j]            H = np.exp(-0.5 * np.power((overlap - A), 2) / sigma_r**2)            F = np.dot(H, G)             F = F / np.sum(F)            dst[i, j] = np.sum(F * overlap)    return dst 

def bflitColor(img, sigma_d, sigma_r):    H, W, C = img.shape    sigma_d = float(sigma_d)    radius = int(3.0 * sigma_d + 0.5)    R = 2 * radius + 1    dst = np.zeros((H, W, C))    img_pad = np.pad(img, ((radius, radius), (radius, radius), (0, 0)), 'symmetric')    x, y = np.mgrid[-radius: radius+1, -radius: radius+1]    G = np.exp(-0.5 * (np.power(x, 2) + np.power(y, 2)) / sigma_d**2)    G = G / np.sum(G)    F = np.zeros(3)    for i in range(H):        for j in range(W):            overlap = img_pad[i : i + R, j: j + R, :]            A = img[i, j, :]            I0 = np.power(overlap[:, :, 0] - A[0], 2)            I1 = np.power(overlap[:, :, 1] - A[1], 2)            I2 = np.power(overlap[:, :, 2] - A[2], 2)            H = np.exp(-0.5 * (I0 + I1 + I2) / sigma_r**2)            F = np.dot(H, G)               F = F / np.sum(F)            dst[i, j, 0] = np.sum(F * overlap[:, :, 0])            dst[i, j, 1] = np.sum(F * overlap[:, :, 1])            dst[i, j, 2] = np.sum(F * overlap[:, :, 2])    return dst

imgColor = imread('test.jpg')dstColor = bflitColor(imgColor, 3, 0.3)plt.subplot(2, 2, 1)imshow_noax(imgColor, normalize=False)plt.title('Original image')plt.subplot(2, 2, 2)imshow_noax(dstColor, normalize=False)plt.title('bflitColor image')

Skin Detection

基于RGB颜色空间的简单阈值肤色识别

在RGB色域里，皮肤满足如下关系式
R>95 && G>40 && B>20 && R>G && R>B && Max(R,G,B)- Min(R,G,B)>15 Abs(R-G)>15 $$

另外还有其他基于YUV， YCbC等色彩空间的肤色检测，但都效果不佳（包括RGB， 例图：尽可能阐释算法的优点的图片）。

参考：Human skin color clustering for face detection

def skin_detect_with_rgb(image):    R = image[:, :, 0]    G = image[:, :, 1]    B = image[:, :, 2]    R = np.where(R > 95, R, 0)    G = np.where(G > 40, G, 0)    B = np.where(B > 20, B, 0)    R = np.where(R > G, R, 0)    R = np.where(R > B, R, 0)    G = np.where(R > 0, G, 0)    B = np.where(R > 0, B, 0)    RG = np.abs(R - G)    R = np.where(RG > 15, R, 0)    G = np.where(R > 0, G, 0)    B = np.where(R > 0, B, 0)    img_mask = np.zeros(image.shape)    img_mask[:, :, 0] = R    img_mask[:, :, 1] = G    img_mask[:, :, 2] = B    return img_mask

imgColor = imread('test.png')imgSkin = skin_detect_with_rgb(imgColor)imgBinary = np.where(imgSkin > 0, 255, 0)plt.subplot(2, 3, 1)imshow_noax(imgColor, normalize=False)plt.title('Original image')plt.subplot(2, 3, 2)imshow_noax(imgSkin, normalize=False)plt.title('Skin image')plt.subplot(2, 3, 3)imshow_noax(imgBinary, normalize=False)plt.title('Binary image')

<matplotlib.text.Text at 0x7f030c51e990>

磨皮是图像滤波的一个实际应用。比较常见滤波有中值滤波，均值滤波，高斯滤波，双边滤波。其中双边滤波对图像边缘保护较好，而且有相应的加速方法Real-Time O(1) Bilateral Filtering, 在美颜中有相当广泛的应用。当然例如导向滤波， on-Local以及BM3D等能保留边缘细节的滤波算法都可以作为磨皮的滤波算法。

保护好非皮肤区域，尤其是眼睛。肤色检测易受外界的影响，肤色检测结果往往会包含嘴唇、牙齿等细节部位，对整张图片质量影响较大。学好DL有助于图像处理，哈哈哈哈。

一个 Github 链接： SkinSmoothing
一个国内网站：微像素

注：图片来源与网络。