图像增强（3）

图像增强（3）

这是图像去雾方向的一个里程碑式方法【1】，简单有效，是何凯明成名作（CVPR09 best paper）
                                                           
                                                                                     图1. 左：原图；右：暗通道去雾结果

算法原理

一般用式(1)对有雾图像建模
I ( x ) = J ( x ) t ( x ) + A ( 1 − t ( x ) ) (1) I(x)=J(x)t(x)+A(1-t(x)) tag{1} I(x)=J(x)t(x)+A(1−t(x))(1)
其中 I I I是观察到的有雾的图像
       J J J是理想图像
       t t t是大气透射率
      A A A是环境光的估计
式（1）的 J ( x ) t ( x ) J(x)t(x) J(x)t(x)被称为 direct attenuation，描述场景辐射的衰减
                 A ( 1 − t ( x ) ) A(1-t(x)) A(1−t(x))被称为 airlight，描述因为空气颗粒导致的光的散射（会造成场景颜色的偏移）

去雾算法的目标就是根据已有的 I I I和一些先验知识估计 A A A和 t t t，最终恢复 J J J

而暗通道方法最核心的一块就是找到了一个具有普遍性质的先验知识：对非天空区域，至少有一个颜色通道的值会非常低

接下来分析如何利用这个先验得到理想图像 J J J:

对图像上一个小窗口，可以认为大气透射率 t t t是一个常量，不妨记为 t ~ ( x ) widetilde{t}(x) t (x)，于是(1)式可写为
min ⁡ y ∈ Ω ( x ) I c ( y ) = t ~ ( x ) min ⁡ y ∈ Ω ( x ) ( J c ( y ) + ( 1 + t ~ ( x ) ) ) A c (2) min_{yin Omega(x)}I^c(y) = widetilde{t}(x)min_{yinOmega(x)}(J^c(y)+(1+widetilde{t}(x)))A^c tag{2} y∈Ω(x)min​Ic(y)=t (x)y∈Ω(x)min​(Jc(y)+(1+t (x)))Ac(2)
其中 min ⁡ y ∈ Ω ( x ) min_{yinOmega(x)} miny∈Ω(x)​表示在 x x x的邻域 Ω ( x ) Omega(x) Ω(x)中取最小值
       上标 c c c表示分别对 r g b rgb rgb三个通道应用该式

对式(2)稍做变形并对等号两边同取颜色通道最小值，有
min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( I c ( y ) A c ) ) = t ~ ( x ) min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( J c ( y ) A c ) ) + ( 1 − t ~ ( x ) ) (3) min_{c}(min_{yinOmega(x)}(frac{I^c(y)}{A^c}))=widetilde{t}(x)min_{c}(min_{yinOmega(x)}(frac{J^c(y)}{A^c}))+(1-widetilde{t}(x)) tag{3} cmin​(y∈Ω(x)min​(AcIc(y)​))=t (x)cmin​(y∈Ω(x)min​(AcJc(y)​))+(1−t (x))(3)

根据 暗通道先验，有
J d a r k ( x ) = min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( J c ( y ) ) ) = 0 (4) J^{dark}(x)=min_{c}(min_{yinOmega(x)}(J^c(y)))=0 tag{4} Jdark(x)=cmin​(y∈Ω(x)min​(Jc(y)))=0(4)

而 A c > 0 A^c>0 Ac>0，因此可以推得
t ~ ( x ) = 1 − min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( I c ( y ) A c ) ) (5) widetilde{t}(x)=1-min_{c}(min_{yinOmega(x)}(frac{I^c(y)}{A^c})) tag{5} t (x)=1−cmin​(y∈Ω(x)min​(AcIc(y)​))(5)

注意到暗通道先验针对的是**“非天空区域”**。对于天空区域， I c ( y ) I^c(y) Ic(y)与环境光 A c A^c Ac非常接近，也就是说
min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( I c ( y ) A c ) ) → 0 ⇒ t ~ ( x ) → 0 (6) min_{c}(min_{yinOmega(x)}(frac{I^c(y)}{A^c}))to0\ Rightarrow widetilde{t}(x)to 0 tag{6} cmin​(y∈Ω(x)min​(AcIc(y)​))→0⇒t (x)→0                (6)
将式(6)代入式(3)，同样可以推得式(5)

因此式(5)对大气投射率的估计无论对于天空区域还是非天空区域，均成立

在实际场景中其实并不存在完全无雾( I ( x ) I(x) I(x))的情况，事实上人眼已经学会根据雾的分布来估计目标的景深。
因此在式(5)加权重保留一定的雾更贴合人眼的观察
t ~ ( x ) = 1 − ω min ⁡ c ( min ⁡ y ∈ Ω ( x ) ( I c ( y ) A c ) ) (6) widetilde{t}(x)=1-omegamin_{c}(min_{yinOmega(x)}(frac{I^c(y)}{A^c})) tag{6} t (x)=1−ωcmin​(y∈Ω(x)min​(AcIc(y)​))(6)
其中 ω ∈ [ 0 , 1 ] omegain[0,1] ω∈[0,1]，文献【1】建议取 ω = 0.95 omega=0.95 ω=0.95

环境光 A A A的估计比较简单：
1）将暗通道图中像素按亮度从大到小排列，选择前0.1%的像素；
2）在原始有雾图像 I I I中查找对应位置像素中最高亮度值作为 A A A，或者对应所有像素亮度平均作为 A A A。

至此，透射率 t ( x ) t(x) t(x)和环境光 A A A均得到，根据式(1)很容易得到去雾图：
J ( x ) = I ( x ) − A m a x ( t ( x ) , t 0 ) + A (7) J(x)=frac{I(x)-A}{max(t(x), t_0)}+A tag{7} J(x)=max(t(x),t0​)I(x)−A​+A(7)
其中 m a x ( t ( x ) , t 0 ) max(t(x),t_0) max(t(x),t0​)是为了避免 t ( x ) t(x) t(x)过小导致图像过饱和
                                                 
                                                          图1. 左上：原图；右上：暗通道图；左下：透射图；右下：去雾图

从图1可以看到，应用(7)式可以去掉原图的雾，但前景的边缘会有“白边”
这是由于窗口操作 Ω ( x ) Omega(x) Ω(x)带来的膨胀效应（图1窗口半径 r = 9 r=9 r=9），可以预想窗口半径变小会改善“白边”问题
                 
                                                                        图2. 左1列：原图；2列：透射图；右3列：去雾图
图2依然应用(7)式，但窗口半径 r = 0 r=0 r=0，
可以看到透射图的膨胀效应没有了，细节更清楚，但去雾图与图1比明显偏暗，特别是图2右下角的火车场景最为明显

原因出现在式(7)的 I ( x ) − A m a x ( t ( x ) , t 0 ) frac{I(x)-A}{max(t(x), t_0)} max(t(x),t0​)I(x)−A​：
窗口变小后透射图 t ( x ) t(x) t(x)中前景（高亮度部分）收缩，上式分母变小；
注意到环境光 A A A是亮度最大的部分，因此 I ( x ) − A I(x)-A I(x)−A中前景部分是负数；
因此窗口变小将导致式(7)亮度降低

经过以上讨论，为了改善“膨胀效应”，同时避免亮度显著降低，需要对透射图使用“保边滤波器”代替图1的“均值/高斯滤波器”

soft matting

作者在文献【1】中提出的方案是 soft matting
matting 研究的是如何将一幅图像的前景信息和背景信息分离的问题（抠图）
I i = α F i + ( 1 − α ) B i (8) I_i=alpha F_i+(1-alpha)B_i tag{8} Ii​=αFi​+(1−α)Bi​(8)
其中 α alpha α是透明度， F F F是前景像素， B B B是背景像素，原图 I I I是两者的线性加权

对比式(1)和式(8)，可以发现两者形式上是一样的，因此可以用该领域的东西来优化透射图 t ( x ) t(x) t(x)

guided filter

soft matting 可以改善图1提出的问题，但计算复杂度和空间复杂度高，几乎没有实用价值

guided filter 是作者在文章【2】中提出的一种滤波器，可以用来解决“膨胀效应”：
q i = ∑ j ∈ ω i W i j ( I ) p j (9) q_i=sum_{jinomega_i}W_{ij}(I)p_j tag{9} qi​=j∈ωi​∑​Wij​(I)pj​(9)
式(9)对原图 p p p进行滤波得到 q q q，滤波权重 W i j W_{ij} Wij​由另外的一张引导图像 I I I获得
当引导图像 I I I就是原图 p p p本身的时候，导向滤波就变成了一个保边滤波器

导向滤波器假设在小窗口 ω k omega_k ωk​内引导图像和滤波输出图像呈线性关系：
q i = a k I i + b k , ∀ i ∈ ω k (10) q_i=a_kI_i+b_k,forall i in omega_k tag{10} qi​=ak​Ii​+bk​,∀i∈ωk​(10)
其中 ( a k , b k ) (a_k,b_k) (ak​,bk​)是窗口 w k w_k wk​的加权参数
线性关系的假设可以确保在该窗口内平滑区域和纹理区域的对比关系得以保留，也就是不丢失纹理信息

由此可以建立优化目标函数
a r g min ⁡ ∑ i ∈ ω k ( a k I i + b k − p i ) 2 (11) argminsum_{iinomega_k}(a_kI_i+b_k-p_i)^2 tag{11} argmini∈ωk​∑​(ak​Ii​+bk​−pi​)2(11)
上式可以理解为：在确保小窗口内由引导图像 I I I体现的纹理尽量不丢失的前提下，让输出的图像 q q q尽量与原图 p p p接近

求解式(11)可以参考最小二乘法的推导
1）首先在式(11)中加入正则化项
E ( a k , b k ) = ∑ i ∈ ω k ( ( a k I i + b k − p i ) 2 + ϵ a k 2 ) (12) E(a_k,b_k)=sum_{iinomega_k}((a_kI_i+b_k-p_i)^2+epsilon a_k^2) tag{12} E(ak​,bk​)=i∈ωk​∑​((ak​Ii​+bk​−pi​)2+ϵak2​)(12)
2）式(12)分别对 a k a_k ak​和 b k b_k bk​求偏导
a k = 1 ∣ ω ∣ ∑ i ∈ ω k I i p i − μ k p ‾ k σ k 2 + ϵ b k = p ‾ k − a k μ k (13) a_k=frac{frac{1}{|omega|}sum_{iinomega_k}{I_ip_i-mu_koverline{p}_k}}{sigma_k^2+epsilon}\ b_k=overline{p}_k-a_kmu_k tag{13} ak​=σk2​+ϵ∣ω∣1​∑i∈ωk​​Ii​pi​−μk​p​k​​bk​=p​k​−ak​μk​                    (13)
      其中 μ k mu_k μk​和 σ k 2 sigma_k^2 σk2​分别是引导图 I I I在小窗口 ω k omega_k ωk​内的均值和方差
              p ‾ k overline{p}_k p​k​是原图 p p p在小窗口 ω k omega_k ωk​内的均值
             ∣ ω ∣ |omega| ∣ω∣是小窗口 ω k omega_k ωk​的像素个数
3）第2步为每个小窗口 ω omega ω求得一组 ( a k , b k ) (a_k,b_k) (ak​,bk​)，但每个像素会被包含在相邻的不同小窗口中（8-邻域），那到底用哪一组参数？
      论文【2】将相邻的参数求平均，作为当前点的参数
q i = 1 ∣ ω ∣ ∑ k : i ∈ ω k ( a k I i + b k ) = a ‾ i I i + b ‾ i (13) q_i=frac{1}{|omega|}sum_{k:iinomega_k}(a_kI_i+b_k)\ =overline{a}_iI_i+overline{b}_i tag{13} qi​=∣ω∣1​k:i∈ωk​∑​(ak​Ii​+bk​)=ai​Ii​+bi​                (13)
      其中 a ‾ i = 1 ∣ ω ∣ ∑ k ∈ ω i a k overline{a}_i=frac{1}{|omega|}sum_{kinomega_i}a_k ai​=∣ω∣1​∑k∈ωi​​ak​， b ‾ i = 1 ∣ ω ∣ ∑ k ∈ ω i b k overline{b}_i=frac{1}{|omega|}sum_{kinomega_i}b_k bi​=∣ω∣1​∑k∈ωi​​bk​

值得一提的是，OpenCV3.x 已经实现了导向滤波
                      
                                                                                        图3. 应用了导向滤波的去雾结果
对比图1、图2和图3，可以看到，应用了导向滤波后既改善了图1的**“膨胀效应”，又有效改善了图2的“亮度降低”**问题

代码

这个三方版本的去雾算法写的比较清晰，但是没加入soft matting和导向滤波
导向滤波作者开源了matlab代码

文献

【1】 He K, Sun J, Tang X. Single image haze removal using dark channel prior[J]. IEEE transactions on pattern analysis and machine intelligence, 2010, 33(12): 2341-2353.
【2】He K, Sun J, Tang X. Guided image filtering[C]//European conference on computer vision. Springer, Berlin, Heidelberg, 2010: 1-14.