暗通道与导引滤波
标签（空格分隔）： 论文学习笔记
1. 暗通道与图像去雾
雾霾是特定气候与人类活动相互作用的结果。高密度人口的经济生产及社会活动会排放大量细颗粒物，一旦排放量超过大气循环和承载能力，悬浮颗粒受静稳天气的影响持续积聚，极易出现大范围的雾霾。
现有的图像去雾（Image Dehazing）技术离不开一个简单的自然模型——大气散射模型（Atmospheric Scattering Model）。大气散射模型描述了，在雾霾和光照的共同作用下的成像机制：
阳光在物体表面形成反射光 J(x)，反射光在穿过雾霾的过程发生散射，只有部分能量 J(x)t(x) 能到达摄像头。与此同时，阳光也在悬浮颗粒表面散射形成大气光 α 被摄像头接收。因此，摄像头中的成像 I(x) 可由两部分组成，透射的物体亮度 J(x)t(x) 和散射的大气光照 α(1－t(x))：
I(x)=J(x)t(x)+A(1−t(x))
其中，t(x) 是媒介透射率（medium transmission），顾名思义表示能顺利透过雾霾到达摄像头的比率。因此，透射率跟物体与摄像头距离 d(x) 成反比，离摄像头越远的物体受雾霾影响更大。当距离 d(x) 趋于无穷大时，透射率 t(x) 趋于零，I(x) 趋近于 α，α=maxy∈{x|t(x)≤t0}I(y)。综上所述，去雾的核心是如何更精确地估计媒介透射率 t(x)。
1.1 暗通道理论
暗通道的定义：
在绝大多数非天空的局部区域里，某一些像素总会有至少一个颜色通道具有很低的值。换言之，该区域光强度的最小值是个很小的数。
我们给暗通道一个数学定义，对于任意的输入图像J，其暗通道可以用下式表达：
Jdark=miny∈Ω(x)(minc∈{r,g,b}Jc(y)).....(1)
其中
Ω(x) 表示以x为中心的图像窗口，
Jc(y) 表示rgb图像某一点y的某一通道。上述公式所表示的意义用代码表达也很简单，首先求出每个像素RGB分量中的最小值，存入一副和原始图像大小相同的灰度图中，然后再对这幅灰度图进行最小值滤波，滤波的半径由窗口大小决定，一般有WindowSize = 2 * Radius + 1;
暗通道先验的理论指出：
Jdark→0
实际生活中造成暗原色中地通道值主要有三个因素：
1)汽车、建筑物和城市中玻璃窗户的阴影，或者是树叶、树与岩石等自然景观的投影；
2)色彩鲜艳的物体或表面，在RGB的三个通道中有些通道的值很低（比如绿色的草地／树／植物，红色或黄色的花朵／叶子，或者蓝色的水面）；
3)颜色较暗的物体或者表面，例如灰暗色的树干和石头。总之，自然景物中到处都是阴影或者彩色，这些景物的图像的暗原色总是很灰暗的。
首先找一张没有雾的图做暗通道处理，（窗口大小为15*15，即最小值滤波的半径为7像素）
其暗通道如下：
再找一张有雾的图像进行暗通道处理:
暗通道处理结果如下：
通过对大量自然图像的暗通道处理统计,发现基本符合上述先验知识。
1.2雾的成像模型
在计算机视觉和计算机图形当中，雾图像的成像模型被广泛认为符合下式：
I(x)=J(x)t(x)+A(1−t(x)).......(2)
其中，I(x)是我们现有的带雾图像，J(x)是我们要恢复的无雾图像，A是全球大气光成分，t(x)为透射率。显然方程有无数解，那么要求得满足条件的解就需要先验条件进行约束。
下面来分析一下推导过程：
对(2)式进行变形
1.3暗通道理论去雾推导
Ic(x)Ac=t(x)Jc(x)Ac+1−t(x)......(3)
其中，c表示图像的rgb某通道。然后，假设在图像窗口内透射率t(x)为常数，表示为
tconst(x) ,然后对式(3)进行暗通道处理，得到下式：
miny∈Ω(x)(minc(Ic(y)Ac))=tconstminy∈Ω(x)(minc)Jc(y)Ac+1−tconst(x).....(4)
上式中，J是待求的无雾图像，根据暗通道先验理论有：
Jdark(x)=miny∈Ω(x)(mincJc(y))=0......(5)
则可得出公式(4)中等式的右边第一项为0，则公式(4)可简化为
miny∈Ω(x)(minc(Ic(y)Ac))=1−tconst(x).....(6)
通过上式我们可以求出透射率预估值。
在现实生活中，即使是晴天白云，空气中也存在着一些颗粒，因此，看远处的物体还是能感觉到雾的影响，另外，雾的存在让人类感到景深的存在，因此，有必要在去雾的时候保留一定程度的雾，这可以通过在式（6）中引入一个在[0,1] 之间的因子，则式（6）修正为：
tconst(x)=1−ωminy∈Ω(x)(minc(Ic(y)Ac)).....(7)
通常取
ω=0.8−0.95 中的某个值。
上述推论中都是假设全球达气光A值时已知的，在实际中，我们可以借助于暗通道图来从有雾图像中获取该值。具体步骤如下：
1） 从暗通道图中按照亮度的大小取前0.1%的像素。
2） 在这些位置中，在原始有雾图像I中寻找对应的具有最高亮度的点的值，作为A值。
到这一步，我们就可以进行无雾图像的恢复了。由式（1）可知： J = ( I – A)/t + A
现在I,A,t都已经求得了，因此，完全可以进行J的计算。
当投射图t 的值很小时，会导致J的值偏大，从而使淂图像整体向白场过度，因此一般可设置一阈值T0，当t值小于T0时，令t=T0，本文中所有效果图均以T0=0.1为标准计算。
因此，最终的恢复公式如下：
J(x)=I(x)−Amax(t0,t)+A.....(8)
单纯的利用上述公式进行图像去雾恢复，会出现比较严重的halo现象，在何凯明的论文中对恢复图像进行进一步的soft matting 处理，可达到较为理想的效果。
1.4 实验结果
1.5 该理论的缺陷
暗通道先验是一种统计的结果，是对大量户外无雾照片(outdoor haze-free images)的统计结果，如果目标场景内在的就和大气光类似，比如雪地、白色背景墙、大海等，则由于前提条件就不正确，因此一般无法获得满意的效果，而对于一般的风景照片这个算法能处理的不错。
2.导引滤波
导引滤波由何凯明博士在2013年PAMI期刊中发表文章Guided image filter首次提出，引起了不小的反向，导引滤波不仅能够在视觉效果上超越双边带滤波，更是在算法实现的复杂度达到了O(N)。接下来就来分析这篇文章。
2.1导引滤波概述
导引滤波相对于双边滤波有更好的边缘保持性，且不会出现梯度反转现象。在不同引导图像的引导下，可广泛应用于降噪、去雾、高动态范围压缩等。在导引滤波的定义中，用到了局部线性模型，至于该模型，可以用下图简单的理解：
该模型认为，某函数上一点与其邻近部分的点成线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需要求该函数上某一点的值时，只需计算所有包含该点的线性函数的值并做平均即可。这种模型，在表示非解析函数上，非常有用。
2.2数学推导
同理，我们可以认为图像是一个二维函数，而且没法写出解析表达式，因此我们假设该函数的输出与输入在一个二维窗口内满足线性关系，如下
qi=akIi+bk,∀i∈ωk.....(9)
其中，q是输出像素的值，I是输入图像的值，i和k是像素索引，a和b是当窗口中心位于k时该线性函数的系数。其实，输入图像不一定是待滤波的图像本身，也可以是其他图像即引导图像，这也是为何称为引导滤波的原因。对上式两边取梯度，可以得到
▽q=a▽I......(10)
即当输入图像I有梯度时，输出q也有类似的梯度，现在可以解释为什么引导滤波有边缘保持特性了。
下一步是求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值与真实值p之间的差距最小，也就是让下式最小
E(ak,bk)=∑i∈ωk((akIi+bk−pi)2+ϵa2k)......(11)
这里p只能是待滤波图像，并不像I那样可以是其他图像。同时，a之前的系数用于防止求得的a过大，也是调节滤波器滤波效果的重要参数。通过最小二乘法，我们可以得到:
ak=1|ω|∑i∈ωkIipi−μkpk¯¯¯δ2k+ϵ........(12)
bk=p¯k−akμk.....(13)
其中，
μk 是I在窗口
ωk 中的方差，
|w| 是窗口
ωk 中像素值的数量，
p¯k 表示待滤波图像p在窗口
ωk 中的均值。
在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下
qi=1|ω|∑k:i∈ωk(akIi+bk)=a¯iIi+b¯i.....(14)
当把引导滤波用作边缘保持滤波器时，往往有 I = p ，如果
ϵ =0，显然a=1, b=0是E(a,b)为最小值的解，从上式可以看出，这时的滤波器没有任何作用，将输入原封不动的输出。如果
ϵ >0，在像素强度变化小的区域（或单色区域），有a近似于（或等于）0，而b近似于（或等于）
p¯k ，即做了一个加权均值滤波；而在变化大的区域，a近似于1，b近似于0，对图像的滤波效果很弱，有助于保持边缘。而
ϵ 的作用就是界定什么是变化大，什么是变化小。在窗口大小不变的情况下，随着
ϵ 的增大，滤波效果越明显。
在滤波效果上，引导滤波和双边滤波差不多，在一些细节上，引导滤波较好。引导滤波最大的优势在于，可以写出时间复杂度与窗口大小无关的算法（打算在之后的文章中讨论），因此在使用大窗口处理图片时，其效率更高。
2.3实验效果
2.3.1去雾效果
下图为暗通道去雾和导引滤波去雾结果对比图，左图为暗通道去雾
2.3.2边缘保持效果
下图将输入图像同时作为导引图像，取得的实验结果为边缘保持的图像模糊效果，下图中左图为输入原图像，右图为滤波半径r = 8,矫正系数 ϵ=0.42 时的滤波结果。
下图为双边带滤波结果
2.4算法复杂度说明
何凯明博士在他的论文中给出了具体的算法步骤，如下图所示：
下面我们从数学推导公式的步骤一步一步分析代码实现:
ak=1|ω|∑i∈ωkIipi−μkpk¯¯¯δ2k+ϵ........(12)
bk=p¯k−akμk.....(13)
公式(12)(13)告诉我们怎么计算a和b的值，先计算
ak 的分子，Ip 在窗口
wk 中的和，再除以窗口中像素的个数，刚好就是盒式滤波(见
盒式滤波)，因此我们可以将输入的引导图像 I 和滤波图像 p 相乘，并对相乘后的图像做box filtering，即得第一项的结果。后面的
μk 和分别为 I 和 p 在窗口
wk 中均值，因此分别对 I 和 p 进行box filtering，再将box filtering之后的结果相乘即可。实际上，
ak 的分子就是 Ip 在窗口
wk 中的协方差。
接下来计算a_k的分母部分。是引导图 I 在窗口
wk 中的方差，学过概率论与数理统计的朋友应该知道，方差和期望（均值）之间是有关系的，如下式
DX=E(X2)−(EX)2
因此在计算 I 的方差时，我们可以先计算 I*I 的均值，再减去 I 均值的平方即
μk 的平方。在方法上，计算 I*I 的均值和计算 Ip 的均值是一样的。最后，对计算出来的方差图像，加上常量
ϵ （每个元素都加
ϵ ），分母就计算完了，自然，
ak 在所有窗口中的值也就得到了。
bk 的计算就更容易了。
输出图像的公式
qi=1|ω|∑k:i∈ωk(akIi+bk)=a¯iIi+b¯i.....(14)
输出值q又与两个均值有关，分别为a和b在窗口
wi 中的均值（不是
wk ），所以还是box filtering，我们将上一步得到两个图像都进行盒式滤波，得到两个新图：
ai 和
bi ，然后用
ai 乘以引导图像 I ，再加上
bi ，即得最终滤波之后的输出。
2.5加权导引滤波
加权导引滤波由Zhengguo Li博士等人在2015年的TIP上发表文章Weighted Guided Image Filtering上提出，该文章在何凯明博士的guided image filter 一文的基础上做了改进，在边缘保持上达到了更好的效果，并很大程度上消除了bilateral filter的halo现象。
读完文章，发现此篇文章建立的数学模型和何凯明博士建立的数学模型基本一致，但做了如下改进: 优化矫正系数 ϵ ，将图像像素的权重考虑到矫正系数当中。具体数学模型如下公式所示:
E=∑p∈Ωδ(p′)[(a′pG(p)+b′p−X(p))2+λΓG(p′)a′2p]
其中
ΓG(p′) 为图像的加权系数,文章中给了如下定义方式
ΓG(p′)=1N∑p=1Nδ2G,1(p′)+ϵδ2G,1(p)+ϵ
3.实际用途
文章中提到导引滤波多用来做边缘保持的图像模糊，图像去雾，高动态图像压缩等等。基于我们的实际用途，去除云阴影，还待进一步实验。
