guided filter（导向滤波）导读

---

本文为笔者阅读Kaiming He于2013年发表的guided image filtering时的笔记和实验合集，作者在文中省略的一些步骤给予了补充，在此记录下来，希望对后来者有益。
本文分为如下两部分：
1.论文阅读及算法原理介绍
2.基于matlab的简单实现

---

1. 论文阅读及算法原理介绍
   先给出线性移-变滤波器的基本概念：给出输入图像p，指导图像I，得到输出图像，输入图像p和指导图像I可以是同一幅图像。滤波之后的图像可以用如下加权公式表示：
   $$q_i=\sum _j{W}_{ij}(I)p_j$$
   其中$i$为当前滤波点的索引，$j$为滤波器模版覆盖的所有像素点的索引，滤波器（卷积核）$W_{ij}$是指导图像$I$的函数，与输入图像$p$独立。

可以用如下图简单解释上述过程：

所谓线性移-变滤波器包含两个性质：（1）线性性：滤波器通过线性加权进行滤波；（2）移-变：滤波其系数$W_{ij}$随着图像中心的变化而变化，与之相对应的是移-不变滤波器，如传统的高斯滤波器/边缘检测算子等滤波器，一旦滤波器定下来不随着像素点的移动而改变。上述公式中的下标索引$i$为红色像素点所对应的空间位置（下标），下标索引$j$为所有蓝色和黄色像素点对应的空间位置。

一个简单的线性移-变滤波器就是联合双边滤波器，示意图如下图：

其核函数为：
$$W_{ij}^{bf}(I)=\frac{1}{K_i}\exp \left(-\frac{\Vert x_i-x_j{\Vert }^2}{{\sigma }_s^2}\right)\exp \left(-\frac{\Vert I_i-I_j{\Vert }^2}{{\sigma }_r^2}\right)$$
其中$K_i$为归一化系数，与传统的双边滤波器不同的是，这里在range方向上用指导图像$I$来衡量相似度，如果range方向上相似度用原始图像来作指导，那么联合滤波器就退化为传统的双边滤波器。
下面给出guided filter的建模及推导过程：guided filter是指导图像和输出图像在一个局部区域内的线性模型，即假定输出图像$q$是指导图像$I$在以$k$为中心的窗口${\omega }_k$线性变换：
$$q_i=a_k{I}_i+b_k,\forall i\in {\omega }_k$$
这里$(a_k,b_k)$在窗口${\omega }_k$中为常数。对公式两边求导可知$\nabla q=a\nabla I$，因此如果指导图像$I$中为边缘的地方，输出图像也有边缘。

通常的加性去噪模型假设输入图像$p$为输出图像$q$（恢复图像）加上某些分量$n$（这里的$n$可以理解成噪声或者某些细小的结构或振荡，比如纹理），因此有如下模型：
$$p_i=q_i+n_i$$
写成输出图像的形式：
$$q_i=p_i-n_i$$
guided filter建模的思想在于寻找这样$(a_k,b_k)$，使得$q_i=a_k{I}_i+b_k$指导滤波之后的输出图像$q_i$和输入图像尽可能相似（差异尽可能小），基于这样的想法，很自然就有了将${\omega }_k$窗口中滤波后的$q_i$和滤波之前的$p_i$的所有误差（差异）都累加起来，使得误差最小的$(a_k,b_k)$就是我们要找的guided filter，于是很自然就建立的下面的模型：
$$\min E(a_k,b_k)=\min \sum _{i\in {\omega }_k}(a_k{I}_i+b_k-p_i{)}^2$$
仔细观察这个模型会发现一个问题，我们寻找输出$q_i$是为了在结构上尽可能地和输入图像$p_i$相似，但是如果通过求解这个最小化模型给出的结果使得$q_i$和$p_i$相等，即输出图像和输入图像一模一样，那么做这个滤波就失去了意义，在该式中，令$b_k=0,a_k=p_i/I_i$即可得到描述的结果，此时$\min E(a_k,b_k)=0$，此乃我们在机器学习/模式识别/人工智能中常说的“过拟合”。于是需要考虑的是改造上面的模型，给它引入一个惩罚项，让输出图像$q_i$尽可能的和输入图像$p_i$尽可能相似，但又不能让它们完全相等，于是给模型加上一个$L^2$范数的正则项，并且给正则项加入一个$ϵ$来控制惩罚的强度，为什么$L^2$范数的正则项可以达到“泛化”和防止“过拟合”的目的，后面会给出一个详细的解释。于是上面的模型修正为：
$$\min E(a_k,b_k)=\min \sum _{i\in {\omega }_k}((a_k{I}_i+b_k-p_i{)}^2+ϵa_k^2)$$
这样才得到了作者在文中给出的模型，这个模型的求解很简单，不用去查所谓的回归分析等内容，这是一个无约束的二次最优话问题，显然是一个凸优化问题，极值点乃最优点，简单运用一下高等数学中的极值点的必要条件即可求解出$(a_k,b_k)$，这里送佛送到西，给出推导的详细过程：
使得$E(a_k,b_k)$取得极小值的必要条件是：
$$\begin{gathered} \frac{\partial E(a_k,b_k)}{\partial a_k}=0 \\ \frac{\partial E(a_k,b_k)}{\partial b_k}=0 \end{gathered}$$
对$E(a_k,b_k)$的具体模型有：
KaTeX parse error: Expected 'EOF', got '&' at position 43: …{\partial a_k} &̲=\sum_{i\in\ome…
先对$b_k$的偏导数进行化简：
$$\begin{array}{rl}& \sum _{i\in {\omega }_k}(a_k{I}_i+b_k-p_i)=0\\ & \sum _{i\in {\omega }_k}{b}_k=\sum _{i\in {\omega }_k}(p_i-a_k{I}_i)\\ & \sum _{i\in {\omega }_k}{b}_k=\sum _{i\in {\omega }_k}{p}_i-a_k\sum _{i\in {\omega }_k}(I_i)\end{array}$$
记${\omega }_k$中像素的个数有$|\omega |$个，那么
$$b_k|\omega |=\sum _{i\in {\omega }_k}{p}_i-a_k\sum _{i\in {\omega }_k}{I}_i$$
两边同时除以$|\omega |$
$$b_k=\frac{{\sum }_{i\in {\omega }_k}{p}_i-a_k{\sum }_{i\in {\omega }_k}{I}_i}{|\omega |}={\overline{p}}_k-a_k{\mu }_k$$
在对$a_k$的偏导数进行化简
$$\begin{array}{rl}& \sum _{i\in {\omega }_k}{I}_i(a_k{I}_i+b_k-p_i)+ϵa_k=0\\ & \sum _{i\in {\omega }_k}(a_k{I}_i^2+I_i{b}_k-I_i{p}_i+ϵa_k)=0\end{array}$$
将上面求出的$b_k$公式带入该式，有
$$\begin{array}{rl}& \sum _{i\in {\omega }_k}{a}_k{I}_i^2+I_i({\overline{p}}_k-a_k{\mu }_k)-I_i{p}_i+ϵa_k=0\\ & \sum _{i\in {\omega }_k}(I_i^2-I_i{\mu }_k+ϵ)a_k+I_i{p}_i=0\\ & a_k\sum _{i\in {\omega }_k}(I_i^2-I_i{\mu }_k+ϵ)=\sum _{i\in {\omega }_k}(I_i{p}_i-I_i{\overline{p}}_k)\end{array}$$

$$\begin{array}{rl}{a}_k& =\frac{{\sum }_{i\in {\omega }_k}(I_i{p}_i-I_i{\overline{p}}_k)}{{\sum }_{i\in {\omega }_k}(I_i^2-I_i{\mu }_k)+{\sum }_{i\in {\omega }_k}ϵ}\\ & =\frac{{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\sum }_{i\in {\omega }_k}{I}_i}{{\sum }_{i\in {\omega }_k}(I_i^2-I_i{\mu }_k)+|\omega |ϵ}\\ & =\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-\frac{{\overline{p}}_k{\sum }_{i\in {\omega }_k}{I}_i}{|\omega |}}{\frac{1}{\omega }{\sum }_{i\in {\omega }_k}(I_i^2-I_i{\mu }_k)+ϵ}\\ & =\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}(I_i^2-I_i{\mu }_k)+ϵ}\\ & =\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i^2-{\mu }_k\frac{{\sum }_{i\in {\omega }_k}{I}_i}{|\omega |}+ϵ}\\ & =\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i^2-{\mu }_k\frac{{\sum }_{i\in {\omega }_k}{I}_i}{|\omega |}+ϵ}\\ & =\frac{\frac{1}{|{\omega }_k|}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{\frac{1}{|{\omega }_k|}{\sum }_{i\in {\omega }_k}{I}_i^2-{\mu }_k{\mu }_k+ϵ}\\ & =\frac{\frac{1}{|\omega |}{\sum }_{i\in \omega }{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i^2-{\mu }_k^2+ϵ}\end{array}$$

由方差定义，
$$\begin{array}{rl}E(X-\overline{X}{)}^2& =E{X}^2-2EX\overline{X}+E{\overline{X}}^2\\ & =E{X}^2-2\overline{X}EX+E{\overline{X}}^2\\ & =E{X}^2-2\overline{X}EX+E{\overline{X}}^2\\ & =E{X}^2-2\overline{X}\overline{X}+{\overline{X}}^2\\ & =E{X}^2-{\overline{X}}^2\end{array}$$
于是有，
$$a_k=\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-{\overline{p}}_k{\mu }_k}{{\sigma }_k^2+ϵ}$$
至此，算法的建模和求解的理论推导过程完毕。
guided filter的计算算法的伪代码描述如下：

严格来说，至此算法部分已经结束，后面部分从理论和应用角度来理解guided filter在实际的滤波过程中和传统滤波器的关系及其相对于传统滤波器所体现出来的优点。
边缘保持特性：guided filter在用作平滑滤波器具有边缘保持特性，在此我么以$I=p$的特殊情况为例，从理论角度分析guided filter的边缘保持性质。
将$I=p$代入$(a_k,b_k)$的求解公式中，有
$$\begin{array}{rl}& a_k=\frac{{\sigma }_k^2}{{\sigma }_k^2+ϵ}\\ & b_k=(1-a_k){\mu }_k\end{array}$$
当$ϵ=0$时，有$a_k=1,b_k=0$,此时滤波后的图像为$q_i=a_k{I}_i+b_k=I_i=p_i$等价于完全没有滤波。
当$ϵ>0$时，有如下两种情况：
（1）指导图像像素值在${\omega }_k$的小窗口内有剧烈变化，即高方差区域，${\sigma }_k^2$远大于$ϵ$
$$a_k=1,b_k=0$$
滤波后图像为：
$$q_i=a_k{I}_i+b_k\approx I_i\approx p_i$$
也就是说在变化明显的区域，该区域滤波后几乎保持不变，从而具有很好的边缘保持特性。
（2）指导图像像素在${\omega }_k$的小窗口内变化不大，即平坦区域，${\sigma }_k^2$接近于0，${\sigma }_k^2$远小于$ϵ$
$$a_k\approx 0,b_k\approx {\mu }_k$$
滤波后图像为：
$$q_i=a_k{I}_i+b_k\approx {\mu }_k$$
也就是说在变化不怎么明显的平坦区域，该区域滤波后相当于对该区域做了均值滤波，从而对平坦区域就有很好的平滑作用。
将guided filter的参数和双边滤波器类比：$r\leftrightarrow {\sigma }_s,ϵ\leftrightarrow {\sigma }_r^2$参数支箭有一个广义上类似的作用。
在本部分的开头硕大guided filter时一个线性移-变滤波器，这里我们结合公式推导来给出guided filter的线性移-变滤波器形式，并给出其物理上的解释。移变滤波器公式为：
$$q_i=\sum _j{W}_{ij}(I)p_j$$
两边同时对$p_j$求偏导数，即可推出guided filter作为线性移-变滤波器的具体公式：
$$W_{ij}(I)=\frac{\partial q_i}{\partial p_j}$$
将$q_i=a_k{I}_i+b_k=\frac{1}{|\omega |}{\sum }_{k|iin{w}_k}{a}_k{I}_i+{\overline{p}}_k=\frac{1}{|\omega |}{\sum }_{k|iin{w}_k}{a}_k(I_i-{\mu }_k)+{\overline{p}}_k$以及$a_k$的具体公式代入上式，有
$$W_{ij}(I)=\frac{1}{\omega }\sum _k(\frac{\partial a_k}{\partial p_j}(I_i-{\mu }_k)+\frac{\partial {\overline{p}}_k}{\partial p_j})$$
$$\frac{\partial {\overline{p}}_k}{\partial p_j}=\frac{1}{|\omega |}\frac{\partial (p_1+\cdots +p_{|\omega |})}{\partial p_j}=\frac{1}{|\omega |}{\delta }_k\in {\omega }_k$$
$$\begin{array}{rlr}\frac{\partial a_k}{\partial p_j}& =\frac{1}{{\sigma }_k^2+ϵ}(\frac{1}{|\omega |}\sum _{i\in {\omega }_k}\frac{\partial p_i}{\partial p_j}{I}_i-\frac{\partial {\overline{q}}_k}{\partial p_j}{\mu }_k)& =\frac{1}{{\sigma }_k^2+ϵ}(\frac{1}{|\omega |}{I}_j-\frac{1}{|\omega |}{\mu }_k){\delta }_{k\in w_j}\end{array}$$
将二式代入移-变滤波器中，有
$$W_{ij}(I)=\frac{\partial q_i}{\partial p_j}=\frac{1}{|\omega {|}^2}\sum _{i\in {\omega }_i,j\in {\omega }_j}(1+\frac{(I_i-{\mu }_k)(I_j-{\mu }_k)}{{\sigma }_k^2+ϵ})$$

当$I_i-{\mu }_k$和$I_j-{\mu }_k$符号相同，即在以像素块${\omega }_k$的均值${\mu }_k$的度量下符号相同，即像素点相似，那么对应的移-变滤波器系数大，从而加权的权重更大；反之，符号相反，差异很大，权重系数小。从而，从滤波器的角度也解释了guided filter具有边缘保持和平滑滤波特性。如果把图像理解为随机场的话：
$$\frac{(I_i-{\mu }_k)(I_j-{\mu }_k)}{{\sigma }_k^2+ϵ}$$
可以理解为调制后的相关函数，相关函数越逼近1，说明相关性越大，从而两个像素对应的像素相似度越高，反之越逼近-1，说明相关性越小，图像相似度越小。

Guided filter的梯度保持特性：
设输入信号为$p$，它的的base部分为$q$,细节部分为$d=p-q$，在$I=p$时，$a_k=\frac{{\sigma }_k^2}{{\sigma }_k^2+ϵ}<1$
那么
$${\partial }_{x}d={\partial }_{x}p-{\partial }_{x}q=(1-a_k){\partial }_{x}p$$
这说明细节的梯度方向和原始信号的梯度方向一致，不会出现梯度方向反转。

上图中bilateral filter的细节部分出现了明显的梯度反转，而guided filter则保持梯度方向不变。这反应在图像上就是在边缘部分，一个有或多或少的振铃现象，另一个则能够避免这样的现象出现。

guide filter面对三位彩色图像的情形也很容易推广，在这里就不再赘述了。

这里补充一下为什么在建模的时候加入$L^2$范数正则化：
很多做过病态反问题的人给出的回答时，正则化的引入可以带来更稳定的数值解，而不至于在输入引入很小的扰动时，数值解很大地偏离了理论解。
做模式识别的人给出的回答是正则化的引入可以防止过拟合，给训练系统更好的泛化能力。在前面的建模部分，也分析了guided filter在指导图像和输入图像一致时，$L^2$范数正则化在一定程度上防止了过拟合情形的出现，下面给出$L^2$范数正则化防止过拟合的理论依据。
假定最优化模型已经建立，最优化目标函数为：$J(\omega ;X,y)$，给目标函数加入$L^2$范数正则项$||\omega ||=\frac{\alpha }{2}{\omega }^T\omega$，新的目标函数为：
$$J_{new}(\omega ;X,y)={\nabla }_{\omega }J(\omega ;X,y)+\alpha \omega$$
对其求梯度，有
$${\nabla }_{\omega }{J}_{new}(\omega ;X,y)={\nabla }_{\omega }J(\omega ;X,y)+\alpha \omega$$
对于${\omega }^{\star }=arg\underset{\omega }{\min }J(\omega ;X,y)$现在我们将多元函数$J(\omega ;X,y)$在${\omega }^{star}$做一个二阶的泰勒近似，于是有：
$$\tilde{J}(\omega )=J({\omega }^{\star })+\frac{1}{2}(\omega -{\omega }^{\star }{)}^{T}H(\omega -{\omega }^{\star })$$
这里H为J的Hessian矩阵，对其求偏导数，没有梯度项时因为极点值必要条件：梯度在极点为零，有
$${\nabla }_{\omega }\tilde{J}(\omega )=H(\omega -{\omega }^{\star })$$
在极点值处，显然上式为零，现在考虑加入$L^2$范数正则项，
$$\alpha \omega +H(\omega -{\omega }^{\star })=0$$
$$(H+\alpha I)\omega =H{\omega }^{\star }$$
$$\omega =(H+\alpha I{)}^{-1}{\omega }^{\star }$$
在$\alpha$很接近0时，此时的$\omega$和没有正则项的极值点很接近，那么$\alpha$逐渐增大呢？
由于Hessian矩阵时一个实的对称矩阵，因此可以对其做正交分解，即存在正交矩阵Q和对角矩阵$\Lambda$使得$H=Q\Lambda Q^T$,于是有
$$\begin{array}{rl}\omega & =(Q\Lambda Q^T+\alpha I{)}^{-1}Q\Lambda Q^T{\omega }^{\star }\\ & =[Q(\Lambda +\alpha I)Q^T{]}^{-1}Q\Lambda Q^T{\omega }^{\star }\\ & =Q(\Lambda +\alpha I{)}^{-1}{Q}^{T}Q\Lambda Q^T{\omega }^{\star }\\ & =Q(\Lambda +\alpha I{)}^{-1}\Lambda Q^T{\omega }^{\star }\end{array}$$
$$\begin{array}{rl}(\Lambda +\alpha I{)}^{-1}\Lambda & ={\left(\begin{array}{cccc}{\lambda }_1+\alpha & & & \\ & {\lambda }_2+\alpha & & \\ & & \ddots & \\ & & & {\lambda }_n+\alpha \end{array}\right)}^{-1}\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)\\ & =\left(\begin{array}{cccc}\frac{{\lambda }_1}{{\lambda }_1+\alpha }& & & \\ & \frac{{\lambda }_2}{{\lambda }_2+\alpha }& & \\ & & \ddots & \\ & & & \frac{{\lambda }_n}{{\lambda }_n+\alpha }\end{array}\right)\end{array}$$
这说明$\alpha$对特征值起到了衰减调制作用，远小于$\alpha$的特征值，在$\alpha$的调制下几乎不能对系统特征起到作用，而想法，在$\alpha$的调制下保留了那些相对于$\alpha$同一数量级甚至更高数量级的特征。这才是$L^2$正则化得以能够进行繁华的理论原因。这里的$\omega$就是guided filter滤波器中的$a_k$，$\alpha$即为guided filter中的$ϵ$。

到这里权重调制的作用已经很清晰了，可能有的读者线性代数功底不够，更进一步完善一下细节。
记$Q=[q_1,q_2,\cdots ,q_n]$，分量为列向量，那么$Q^T=[q_1,q_2,\cdots ,q_n{]}^T=$，代入$\omega$中，有

$$\begin{array}{rlr}\omega & =Q(\Lambda +\alpha I)Q^T{\omega }^{\ast }& \\ & =\left(\begin{array}{c}{q}_1,q_2,\cdots ,q_n\end{array}\right)\left(\begin{array}{cccc}\frac{{\lambda }_1}{{\lambda }_1+\alpha }& & & \\ & \frac{{\lambda }_2}{{\lambda }_2+\alpha }& & \\ & & \ddots & \\ & & & \frac{{\lambda }_n}{{\lambda }_n+\alpha }\end{array}\right)\left(\begin{array}{c}{q}_1^T\\ q_2^2\\ ⋮\\ q_n^2\end{array}\right){\omega }^{\ast }& =(\sum _{i=1}^n\frac{{\lambda }_i}{{\lambda }_i+\alpha }{q}_i{q}_i^T)\ast {\omega }^{\ast }\end{array}$$

到这里就非常清晰了，如果$\alpha$远大于某些特征值，那么这个值对应的权重为$\frac{1}{1+\infty }\approx 0$，对应的特征在正则项的调制下可以忽略不计，反之则该特征值对应的特征很强，对应的权重也很大，是主要特征。反映到guided filter中，就对应了$ϵ$对$a_k$的调制作用，指导图像的${\omega }_k$图像块的方差就对应了这里的特征值$\lambda$,$ϵ/{\sigma }_k^2$的相对大小，就决定了对应的权重调制到什么程度。

2.基于matlab的简单实现
该部分仅仅从guided filter的原理层面对其进行matlab简单实现，笔者一直认为matlab只是验证算法正确性的工具，代码尽可能贴近原理，算法性能优化及性能对比还是以c/cpp等工程应用编程语言为准，其O(N)的快速算法这里不涉及。matlab自带的imfilter本身经过加速，用for循环进行计算反而会降低其速度。
为了方便和原理对比，这里再贴一下算法的伪代码描述：

单通道灰度图像实现

function q = guided_filter(I, p, r, epsilon)
% demo realization of guided filter
% just used to illustrate the principle of guided filter
w = fspecial('average', 2*r+1);
mean_I = imfilter(I, w);
mean_p = imfilter(p, w);
corr_I = imfilter(I.*I, w);
corr_Ip = imfilter(I.*p, w);
var_I = corr_I - mean_I .* mean_I;
cov_Ip = corr_Ip - mean_I .* mean_p;
a = cov_Ip ./ (var_I + epsilon);
b = mean_p - a .* mean_p;
mean_a = imfilter(a, w);
mean_b = imfilter(b, w);
q = mean_a .* I + mean_b;
end

参考文献：

1. KaiMing He, Jian Sun and Xiaoou Tang, Guided Image Filtering, IEEE Trans on PAMI, Vol, No.6, June 2013, P1397-1409
2. Ian Goodfellow, Yoshua Bengio, Aaron Courville, Deep Learning, Chap7 Regularization for Deep Learning, P231-234.