综述：图像滤波常用算法实现及原理解析（2）

双边滤波的代码

opencv中提供了bilateralFilter()函数来实现双边滤波操作，其原型如下：

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

• InputArray src: 输入图像，可以是Mat类型，图像必须是8位整型或浮点型单通道、三通道的图像。
• OutputArray dst: 输出图像，和原图像有相同的尺寸和类型。
• int d: 表示在过滤过程中每个像素邻域的直径范围。如果这个值是非正数，则函数会从第五个参数sigmaSpace计算该值。
• double sigmaColor: 颜色空间过滤器的值，这个参数的值月大，表明该像素邻域内有越宽广的颜色会被混合到一起，产生较大的半相等颜色区域。 （这个参数可以理解为值域核的 和 )
• double sigmaSpace: 坐标空间中滤波器的sigma值，如果该值较大，则意味着越远的像素将相互影响，从而使更大的区域中足够相似的颜色获取相同的颜色。当d>0时，d指定了邻域大小且与sigmaSpace无关，否则d正比于sigmaSpace. （这个参数可以理解为空间域核的 和 ）
• int borderType=BORDER_DEFAULT: 用于推断图像外部像素的某种边界模式，有默认值BORDER_DEFAULT.

具体代码如下：

#include <iostream>#include <opencv2/opencv.hpp>
using namespace std;using namespace cv;
//定义全局变量const int g_ndMaxValue = 100;const int g_nsigmaColorMaxValue = 200;const int g_nsigmaSpaceMaxValue = 200;int g_ndValue;int g_nsigmaColorValue;int g_nsigmaSpaceValue;
Mat g_srcImage;Mat g_dstImage;
//定义回调函数void on_bilateralFilterTrackbar(int, void*);
int main(){    g_srcImage = imread("lena.jpg");
    //判断图像是否加载成功    if(g_srcImage.empty())    {        cout << "图像加载失败!" << endl;        return -1;    }    else        cout << "图像加载成功!" << endl << endl;
    namedWindow("原图像", WINDOW_AUTOSIZE);    imshow("原图像", g_srcImage);
    //定义输出图像窗口属性和轨迹条属性    namedWindow("双边滤波图像", WINDOW_AUTOSIZE);    g_ndValue = 10;    g_nsigmaColorValue = 10;    g_nsigmaSpaceValue = 10;
    char dName[20];    sprintf(dName, "邻域直径 %d", g_ndMaxValue);
    char sigmaColorName[20];    sprintf(sigmaColorName, "sigmaColor %d", g_nsigmaColorMaxValue);
    char sigmaSpaceName[20];    sprintf(sigmaSpaceName, "sigmaSpace %d", g_nsigmaSpaceMaxValue);
    //创建轨迹条    createTrackbar(dName, "双边滤波图像", &g_ndValue, g_ndMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_ndValue, 0);
    createTrackbar(sigmaColorName, "双边滤波图像", &g_nsigmaColorValue,                     g_nsigmaColorMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_nsigmaColorValue, 0);
    createTrackbar(sigmaSpaceName, "双边滤波图像", &g_nsigmaSpaceValue,                    g_nsigmaSpaceMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_nsigmaSpaceValue, 0);
    waitKey(0);
    return 0;}
void on_bilateralFilterTrackbar(int, void*){    bilateralFilter(g_srcImage, g_dstImage, g_ndValue, g_nsigmaColorValue, g_nsigmaSpaceValue);    imshow("双边滤波图像", g_dstImage);}

需要有高斯滤波和双边滤波的相关知识背景才能更好的理解导向滤波。在导向滤波中，首先利用了局部线性模型。这个模型认为某函数上一点与其近邻部分的点成线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需要求该函数上某一点的值时，只需要计算所有包含该点的线性函数的值并取平均值即可。这种模型，在表示非解析函数上，非常有用。

同理，我们可以认为图像是一个二维函数，并且假设该函数的输出与输入在一个二维窗口内满足线性关系，如下：

其中，是输出像素的值，是输入图像的值，和是像素索引，和是当窗口中心位于k时该线性函数的系数。其实，输入图像不一定是待滤波的图像本身，也可以是其他图像即引导图像，这也是为何称为引导滤波的原因。对上式两边取梯度，可以得到:

即当输入图像有梯度时，输出也有类似的梯度，现在可以解释为什么引导滤波有边缘保持特性了。下一步是求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值与真实值之间的差距最小，也就是让下式最小:

这里只能是待滤波图像，并不像那样可以是其他图像。同时，之前的系数用于防止求得的过大，也是调节滤波器滤波效果的重要参数(相当于L2正则化的权重惩罚)。接下来利用最小二乘法的原理令  和 得到2个二元一次方程，求解得到：

其中 是 在窗口的平均值， 是 在窗口 的方差， 是窗口 中的像素个数，  是待滤波图像在窗口 中的均值。在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下:

这里， 是所有包含像素 的窗口， 是其中心位置。

当把引导滤波用作边缘保持滤波器时，往往有  ，如果 ，显然是为最小值的解，从上式可以看出，这时的滤波器没有任何作用，将输入原封不动的输出。如果 ，在像素强度变化小的区域（或单色区域），有近似于（或等于0，而近似于（或等于） ，即做了一个加权均值滤波；而在变化大的区域，近似于1，近似于0，对图像的滤波效果很弱，有助于保持边缘。而 的作用就是界定什么是变化大，什么是变化小。在窗口大小不变的情况下，随着的增大，滤波效果越明显。

在滤波效果上，引导滤波和双边滤波差不多，然后在一些细节上，引导滤波较好(在PS的磨皮美白中，经过实践，效果更好)。引导滤波最大的优势在于，可以写出时间复杂度与窗口大小无关的算法，因此在使用大窗口处理图片时，其效率更高。

同样，OpenCV中也有导向滤波的接口。具体代码如下：

void cv::ximgproc::guidedFilter ( InputArray  guide,
InputArray  src,
OutputArray  dst,
int  radius,
double  eps,
int  dDepth = -1 
) 

这边有个基于scipy实现的python代码，可以参考一下：

import numpy as npimport scipy as spimport scipy.ndimage

def box(img, r):    """ O(1) box filter        img - >= 2d image        r   - radius of box filter    """    (rows, cols) = img.shape[:2]    imDst = np.zeros_like(img)

    tile = [1] * img.ndim    tile[0] = r    imCum = np.cumsum(img, 0)    imDst[0:r+1, :, ...] = imCum[r:2*r+1, :, ...]    imDst[r+1:rows-r, :, ...] = imCum[2*r+1:rows, :, ...] - imCum[0:rows-2*r-1, :, ...]    imDst[rows-r:rows, :, ...] = np.tile(imCum[rows-1:rows, :, ...], tile) - imCum[rows-2*r-1:rows-r-1, :, ...]
    tile = [1] * img.ndim    tile[1] = r    imCum = np.cumsum(imDst, 1)    imDst[:, 0:r+1, ...] = imCum[:, r:2*r+1, ...]    imDst[:, r+1:cols-r, ...] = imCum[:, 2*r+1 : cols, ...] - imCum[:, 0 : cols-2*r-1, ...]    imDst[:, cols-r: cols, ...] = np.tile(imCum[:, cols-1:cols, ...], tile) - imCum[:, cols-2*r-1 : cols-r-1, ...]
    return imDst
def _gf_color(I, p, r, eps, s=None):    """ Color guided filter    I - guide image (rgb)    p - filtering input (single channel)    r - window radius    eps - regularization (roughly, variance of non-edge noise)    s - subsampling factor for fast guided filter    """    fullI = I    fullP = p    if s is not None:        I = sp.ndimage.zoom(fullI, [1/s, 1/s, 1], order=1)        p = sp.ndimage.zoom(fullP, [1/s, 1/s], order=1)        r = round(r / s)
    h, w = p.shape[:2]    N = box(np.ones((h, w)), r)
    mI_r = box(I[:,:,0], r) / N    mI_g = box(I[:,:,1], r) / N    mI_b = box(I[:,:,2], r) / N
    mP = box(p, r) / N
    # mean of I * p    mIp_r = box(I[:,:,0]*p, r) / N    mIp_g = box(I[:,:,1]*p, r) / N    mIp_b = box(I[:,:,2]*p, r) / N
    # per-patch covariance of (I, p)    covIp_r = mIp_r - mI_r * mP    covIp_g = mIp_g - mI_g * mP    covIp_b = mIp_b - mI_b * mP
    # symmetric covariance matrix of I in each patch:    #       rr rg rb    #       rg gg gb    #       rb gb bb    var_I_rr = box(I[:,:,0] * I[:,:,0], r) / N - mI_r * mI_r;    var_I_rg = box(I[:,:,0] * I[:,:,1], r) / N - mI_r * mI_g;    var_I_rb = box(I[:,:,0] * I[:,:,2], r) / N - mI_r * mI_b;
    var_I_gg = box(I[:,:,1] * I[:,:,1], r) / N - mI_g * mI_g;    var_I_gb = box(I[:,:,1] * I[:,:,2], r) / N - mI_g * mI_b;
    var_I_bb = box(I[:,:,2] * I[:,:,2], r) / N - mI_b * mI_b;
    a = np.zeros((h, w, 3))    for i in range(h):        for j in range(w):            sig = np.array([                [var_I_rr[i,j], var_I_rg[i,j], var_I_rb[i,j]],                [var_I_rg[i,j], var_I_gg[i,j], var_I_gb[i,j]],                [var_I_rb[i,j], var_I_gb[i,j], var_I_bb[i,j]]            ])            covIp = np.array([covIp_r[i,j], covIp_g[i,j], covIp_b[i,j]])            a[i,j,:] = np.linalg.solve(sig + eps * np.eye(3), covIp)
    b = mP - a[:,:,0] * mI_r - a[:,:,1] * mI_g - a[:,:,2] * mI_b
    meanA = box(a, r) / N[...,np.newaxis]    meanB = box(b, r) / N
    if s is not None:        meanA = sp.ndimage.zoom(meanA, [s, s, 1], order=1)        meanB = sp.ndimage.zoom(meanB, [s, s], order=1)
    q = np.sum(meanA * fullI, axis=2) + meanB
    return q

def _gf_gray(I, p, r, eps, s=None):    """ grayscale (fast) guided filter        I - guide image (1 channel)        p - filter input (1 channel)        r - window raidus        eps - regularization (roughly, allowable variance of non-edge noise)        s - subsampling factor for fast guided filter    """    if s is not None:        Isub = sp.ndimage.zoom(I, 1/s, order=1)        Psub = sp.ndimage.zoom(p, 1/s, order=1)        r = round(r / s)    else:        Isub = I        Psub = p

    (rows, cols) = Isub.shape
    N = box(np.ones([rows, cols]), r)
    meanI = box(Isub, r) / N    meanP = box(Psub, r) / N    corrI = box(Isub * Isub, r) / N    corrIp = box(Isub * Psub, r) / N    varI = corrI - meanI * meanI    covIp = corrIp - meanI * meanP

    a = covIp / (varI + eps)    b = meanP - a * meanI
    meanA = box(a, r) / N    meanB = box(b, r) / N
    if s is not None:        meanA = sp.ndimage.zoom(meanA, s, order=1)        meanB = sp.ndimage.zoom(meanB, s, order=1)
    q = meanA * I + meanB    return q

def _gf_colorgray(I, p, r, eps, s=None):    """ automatically choose color or gray guided filter based on I's shape """    if I.ndim == 2 or I.shape[2] == 1:        return _gf_gray(I, p, r, eps, s)    elif I.ndim == 3 and I.shape[2] == 3:        return _gf_color(I, p, r, eps, s)    else:        print("Invalid guide dimensions:", I.shape)

def guided_filter(I, p, r, eps, s=None):    """ run a guided filter per-channel on filtering input p        I - guide image (1 or 3 channel)        p - filter input (n channel)        r - window raidus        eps - regularization (roughly, allowable variance of non-edge noise)        s - subsampling factor for fast guided filter    """    if p.ndim == 2:        p3 = p[:,:,np.newaxis]
    out = np.zeros_like(p3)    for ch in range(p3.shape[2]):        out[:,:,ch] = _gf_colorgray(I, p3[:,:,ch], r, eps, s)    return np.squeeze(out) if p.ndim == 2 else out

def test_gf():    import imageio    cat = imageio.imread('cat.bmp').astype(np.float32) / 255    tulips = imageio.imread('tulips.bmp').astype(np.float32) / 255
    r = 8    eps = 0.05
    cat_smoothed = guided_filter(cat, cat, r, eps)    cat_detail = cat / cat_smoothed    print(cat_detail.shape)    cat_smoothed_s4 = guided_filter(cat, cat, r, eps, s=4)
    imageio.imwrite('cat_smoothed.png', cat_smoothed)    imageio.imwrite('cat_smoothed_s4.png', cat_smoothed_s4)    imageio.imwrite('cat_smoothed_detailed.png',cat_detail)
    tulips_smoothed4s = np.zeros_like(tulips)    tulips_detailed = np.zeros_like(tulips)    for i in range(3):        tulips_smoothed4s[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps, s=4)
    tulips_detailed = tulips / tulips_smoothed4s    imageio.imwrite('tulips_detailed.png',tulips_detailed)    imageio.imwrite('tulips_smoothed4s.png', tulips_smoothed4s)
    tulips_smoothed = np.zeros_like(tulips)    for i in range(3):        tulips_smoothed[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps)    imageio.imwrite('tulips_smoothed.png', tulips_smoothed)
if __name__ == '__main__':    test_gf()

一副图像，经过mask是图像本身的导向滤波之后，得到一张细节图和一张滤波图。下面从左到右分别是原图，细节图和滤波图。其实这是现在很多low-level领域的预处理步骤。拿到细节图之后可以用卷积神经网络做下面的处理。

这里还推荐一个很好的轮子，C++实现的。

https://github.com/atilimcetin/guided-filter

结语 

以上就是常见的四种滤波算法的介绍。