保边滤波算法及python实现

三大保边滤波算法：引导滤波、双边滤波、加权最小二乘平滑滤波算法

双边滤波

各向异性滤波，基于小波变换的滤波（Edge-avoiding wavelets and their applications），基于域变换的滤波（Domain Transform for Edge-Aware Image and Video Processing），基于geodesic的滤波(Geodesic Image and Video Editing)
双边滤波器：非线性滤波器，输出图像每个像素都是其邻域各像素的加权平均，权值随空间距离和数值差异的增加而减小。

$$BLF(g{)}_p=\frac{1}{k_p}\sum _q{G}_{{\sigma }_s}(||p-q)G_{{\sigma }_r}(||g_p-g_q)$$
$$k_p=\sum _q{G}_{{\sigma }_s}(||p-q||)G_{{\sigma }_r}(||g_p-g_q||)$$
其中，g为输入图像，p和q表示像素的空间位置，核函数$G_{{\sigma }_s}$、$G_{{\sigma }_r}$通常高斯核，${\sigma }_s$决定空间支撑、${\sigma }_r$决定值域支撑调节边缘敏感度
双边滤波器性能定性分析：
1.给定像素p，其原始像素值为$g_p$ （a particular pixel p, whose unfiltered value is $g_p$ ）。随${\sigma }_s$的增大，更多的邻域像素q参与像素p的平滑，其中与$g_p$相近的$g_q$ 贡献最大，导致平滑后p点像素值与g_p差值很小（as ${\sigma }_s$ is increased, more and more distant pixels q , whose value $g_q$ is close to $g_p$ are being averaged together, and as a result the filtered value does not stray too far from $g_p$ ）。当${\sigma }_s\to \infty$时，双边滤波器等同于值域滤波器（in the limit (${\sigma }_s\to \infty$), the bilateral filter becomes a range filter）。因此，为实现大尺度平滑，不仅需要增大${\sigma }_s$，${\sigma }_r$也需要相应增加（more aggressive smoothing cannot be achieved only by increasing ${\sigma }_s$ , and the range support ${\sigma }_r$ must be increased as well）。
2.增大${\sigma }_r$会降低双边滤波器的保边性能，导致部分边缘模糊（increasing σ_r reduces the ability of the bilateral filter to preserve edges, and some of them become blurred）。当${\sigma }_r\to \infty$，双边滤波器等同于高斯滤波器（in the limit (${\sigma }_r\to \infty$ ) the bilateral filter becomes a linear Gaussian filter）。

imgfilter = cv2.bilateralFilter(src, d, sigmaColor, sigmaSpace,borderType)
# 使用方法
imgfilter = cv2.bilateralFilter(img, 11, 25, 25)

第一个参数src,表示输入图像
第二个参数d，表示在滤波过程中每个像素邻域的直径
第三个参数sigmaColor，颜色空间滤波器的sigma值。这个值越大，表示该像素领域内有更宽的颜色被混合到一起
第四个参数sigmaSpace，坐标空间中滤波器的sigma值，值越大，表示更大的区域足够相似的颜色获取相同的颜色。当d＞0时，无论sigmaSpace取何值，d指定邻域大小。否则d正比与sigmaSpace。
第五个参数borderType，处理边界模式，可直接选默认值

导数滤波

文章链接：链接1
有博主做了翻译：链接2
引导滤波的定义中，用到了局部线性模型，该模型认为某函数上一点与其邻近部分的点成线性关系，一个复杂的函数可以用多个局部的线性函数来表示，当需求该函数上某一点的值时，只需计算所有包含该点的线性函数的值做平均，引导滤波函数的输入输出在二维窗口内满足线性关系：
$$q_i=a_k{I}_i+b_k,\forall i\in w_k$$
p是输入图像，I是引导图像，q是滤波输出图像，$a_k$和$b_k$是系数，$w_k$是窗口
对等式两边求导$\nabla q=a\nabla I$
当引导图像$I$有梯度时，输出q也有类似的梯度
最小化输出图q和输入图p的差别，定义代价函数：
$$E(a_k,b_k)=\sum _{i\in w_k}((a_k{I}_i+b_k-p_i{)}^2+ϵa_k^2)$$
这里，ϵ是正则化参数，对代价函数求最小值

$p_k$和${\mu }_k$是p和I窗口内的均值。求$a_k$


当作为滤波器时，通常I=p，
如果$ϵ=0$，从上式中可以看出，滤波器没有任何作用，输出等于输入
$ϵ>0$，在像素强度变化小的区域，做了一个加权均值滤波，在变化大的区域，对图像的滤波效果很弱，有助于保持边缘
$ϵ$的作用就是界定变化的大小，在窗口不变的情况下，随着ϵ的增大，滤波效果更加明显。
指导滤波算法流程：

以彩色图像做指导滤波
上面的结果，图像对应位置相乘，都是单通道的情况，指导图为单通道，如果要用彩色图像为指导图，则略微复杂一点

import cv2
import numpy as np
from skimage.metrics import peak_signal_noise_ratio


def guide_filter(I, p, r, eps):
    """
    :param I: 引导图像
    :param p: 输入图像
    :param r: 窗口半径
    :param eps:
    :return:
    """
    I = I / 255  # 引导滤波图像必须归一化
    p = p / 255
    winSize = (2 * r + 1, 2 * r + 1)
    # I的均值
    mean_I = cv2.blur(I, winSize)
    # p的均值
    mean_p = cv2.blur(p, winSize)

    # I*I和I*p的均值平滑
    mean_II = cv2.blur(I*I, winSize)
    mean_Ip = cv2.blur(I*p, winSize)

    # 方差
    var_I = mean_II - mean_I * mean_I

    # 协方差
    cov_Ip = mean_Ip - mean_I * mean_p

    a = cov_Ip / (var_I + eps)
    b = mean_p - a*mean_I

    # 对a、b进行均值平滑
    mean_a = cv2.blur(a, winSize)
    mean_b = cv2.blur(b, winSize)

    q = mean_a * I + mean_b
    q = q * 255
    q[q > 255] = 255

    return q.astype(np.uint8)

if __name__ == '__main__':

    img = cv2.imread('1.jpg', 0)

    # 引导滤波
    imgguidefilter = guide_filter(img, img, 16, 1e-2)
 
    cv2.imshow('img', img)
    cv2.imshow('guide_filter', imgguidefilter)
    cv2.waitKey()

加权最小二乘法滤波原理

链接1
基于双边滤波的算法无法在多尺度提取到很好的细节信息，并可能出现伪影。加权最小二乘滤波目的即是使得结果图像u与原图像p经过平滑后尽量相似，但是在边缘部分尽量维持原状
     ${\sum }_p((u_p-g_p{)}^2+\lambda (a_{x,p}(g)(\frac{{\partial }_u}{{\partial }_x}{)}_p^2+a_{y,p}(g)(\frac{{\partial }_u}{{\partial }_y}{)}_p^2)$     式（1）
其中$a_x$，$a_y$为权重系数。目标函数第一项代表输入图像和输出图像越相似越好；第二项是正则项，通过最小化的偏导，使得输出图像越平滑越好，上式可以改写为矩阵形式：
    $(u-g{)}^T(u-g)+\lambda (u^T{D}_x^T{A}_x{D}_{x}u+u^T{D}_y^T{A}_y{D}_{y}u)$  式（2）
其中，$A_x$,$A_y$为以$a_x$，$a_y$为对角元素的对角矩阵，$D_x$,$D_y$为前向差分矩阵，和是后向差分算子，使得上式最小值，u满足：
     $(I+\lambda L_g)u=g$                     式（3）
（
注：式3推导，式2对u求导，导数为0。$u^{T}u=u^2$
$$\frac{\partial ((u-g{)}^2+\lambda (u^2{D}_x^T{A}_x{D}_x+u^2{D}_y^T{A}_y{D}_y))}{\partial u}=0$$
$$2(u-g)+2\lambda u(D_x^T{A}_x{D}_x+D_y^T{A}_y{D}_y)=0$$
$$u(I+\lambda (D_x^T{A}_x{D}_x+D_y^T{A}_y{D}_y))=g$$
)

其中，$L_g=D_x^T{A}_x{D}_x+D_y^T{A}_y{D}_y$，作者取的平滑权重系数为
$$a_{x,p}(g)=(|\frac{\partial l}{\partial x}(p){|}^{\alpha }+\epsilon {)}^{-1},a_{y,p}(g)=(|\frac{\partial l}{\partial y}(p){|}^{\alpha }+\epsilon {)}^{-1}$$
其中I表示Log，ε一般取0.0001，式3求出u为：
$$u=F_{\lambda }(g)=(I+\lambda L_g{)}^{-1}g$$
当所选区域连续，平滑权重系数可以近似为a_x≈a_y≈a，那么
$$F_{\lambda }(g)\approx (I+\lambda L_g{)}^{-1}g$$
L为拉普拉斯齐次矩阵

注：拉普拉斯矩阵

拉普拉斯矩阵(Laplacian matrix) 也叫做导纳矩阵、基尔霍夫矩阵或离散拉普拉斯算子，是图论中用于表示图的一种重要矩阵。
定义：给定一个具有n个顶点的简单图G=(V,E)，V为顶点集合，E为边集合，其拉普拉斯矩阵可定义为：
$$L=D-A$$
其中A∈R为邻接矩阵，$D\in R^{n\times n}$为度矩阵
（邻接矩阵是表示顶点之间相邻关系的矩阵。

$$A_{i}j=\{\begin{array}{ll}{w}_{ij}& ，(v_i,v_j)或<v_i,v_j>\in E(G)\\ 0或\infty & ，(v_i,v_j)或<v_i,v_j>\notin E(G)\end{array}$$

$w_{i}j$表示边上的权重，$E（G）$边的集合）
注意A中的元素仅仅可能是0和1，且其对角元素全为0。D是一个对角矩阵，其对角线上的元素计算方式$D_{ii}={\sum }_j{A}_{ij}$ 。L的每个元素值的具体计算方式如下：
$$L_{ij}=\{\begin{array}{ll}{D}_{ij}& ifi=j\\ -1& ifi\ne jand{v}_{i}isadjacentto{v}_j\\ 0& otherwise\end{array}$$
拉普拉斯矩阵都是对称的