论文阅读：Non-Local Image Dehazing

代码：https://github.com/danaberman/non-local-dehazing
本文中关于雾线的思想，和Meng（BCCR）方法中“推向论”的思想有异曲同工之妙。

1. 摘要

    本文提出了一种基于Non-Local先验的算法。 该算法基于以下假设：无雾图像的颜色可以很好地近似为数百种不同的颜色，这些颜色在RGB空间中形成紧密的簇。 通过观察发现，这些簇中包含的像素通常是非局部的，即他们分布在整个图像平面上，并且与相机的距离不同。在有雾的情况下，这些不同的距离转化为不同的传输系数。因此，清晰图像中每中颜色簇在雾图的RGB空间中为成为一条线（Meng：清晰像素在雾的影响下，被“推向”了大气光）。 使用这些雾线，本文提出的算法可以恢复传输率图和无雾图像。

2. 先验

    本文方法基于以下观察：图像中不同颜色的数量比像素的数量小几个数量级。使用K-means算法将一幅图像的RGB像素值聚类为500个聚类，并将图像中的每个像素替换为其所属聚类的中心像素。结果如下：

图1 量化前后对比

    可以看到，肉眼无法分辨其中的差别。
    图2(b)中展示了无雾图像中四个聚类，分别用四种颜色表示。其在原图中的像素用相同颜色的标出。可以看到证实了作者的观察：这些簇中包含的像素通常是非局部的，即他们分布在整个图像平面上，并且与相机的距离不同。 当给无雾图像渲染雾后，聚类中的点被延申成了一条线，即作者提出的雾线。不太理解的同学可以看我写的Meng(BCCR)方法的博客，里面介绍了“推向论”的思想。简单来说，就是无雾情况下，$I(x)\approx J(x)$，在雾的渲染下，$t(x)$减小，导致$I(x)=t(x)J(x)+(1-t(x))A$，$J(x)$的成分减少，$A$的成分增加。雾线的一端是$J(x)$（图2(d)中的彩色点），一端是$A$（图2(d)中的黑色点），所以像素点由$J(x)$被推向了$A$。因为聚类中的点的距离各不相同，雾的影响效果也各不相同（雾主要是通过深度来影响），所以雾线中点的分布也各不相同。

图2 雾线示例图

3. 去雾方法

算法一共四步，即1) 将像素聚类成雾线，2) 估计初始传输率，3) 正则化，4)去雾。算法过程如下：

3.1 Finding Haze-Lines

    本文使用的估计大气光的方法：“Air-light Estimation using Haze-Lines”, ICCP 2017。
    定义如下变量：
$${\mathbf{I}}_A(\mathbf{x})=\mathbf{I}(\mathbf{x})-\mathbf{A}\text{\hspace{1em}(1)}$$

    即将3D RGB 空间坐标系转换成原点为大气光的坐标系。结合大气散射模型，得到：
$${\mathbf{I}}_A(\mathbf{x})=t(\mathbf{x})\cdot [\mathbf{J}(\mathbf{x})-\mathbf{A}]\text{\hspace{1em}(2)}$$

    同时，可以将$I_A(x)$用球系坐标系表示：
$${\mathbf{I}}_A(\mathbf{x})=[r(\mathbf{x}),\theta (\mathbf{x}),\varphi (\mathbf{x})]\text{\hspace{1em}(3)}$$

    其中，$r$是像素点距离原点的距离（即$||\mathbf{I}-\mathbf{A}||$），$\theta$和$\varphi$分别代表经度和维度。该坐标系如图3所示。

图3 球系坐标系

    在图3中，每一个球面上的三角形表示一个无雾图像的聚类，球心为大气光。每一个三角形到球心都是一条雾线。在同一条雾线上的像素点，仅仅只是$r(\mathbf{x})$的不同。

3.2 Estimating Initial Transmission

    对于一条由$\mathbf{J}$和$\mathbf{A}$定义的雾线，$r(\mathbf{x})$与物体的深度有关：
$$r(\mathbf{x})=t(\mathbf{x})||\mathbf{J}(\mathbf{x})-\mathbf{A}||,0\le t(\mathbf{x})\le 1.\text{\hspace{1em}(4)}$$

    $||\mathbf{J}(\mathbf{x})-\mathbf{A}||$代表的是$\mathbf{J}(\mathbf{x})$和$\mathbf{A}$的距离，即由这两个点构成的雾线的长度。由于雾的存在，$\mathbf{I}(x)$被“推向”大气光$\mathbf{A}$，而“推向距离”与$t(\mathbf{x})$有关。

    当$t=1$时，最大半径坐标为：
$$r_{\max }\stackrel{\underset{def}{}}{=}||\mathbf{J}-\mathbf{A}||\text{\hspace{1em}(5)}$$

    结合$(4)(5)$，可以得到传输率：
$$t(\mathbf{x})=\frac{r(\mathbf{x})}{r_{\max }}\text{\hspace{1em}(6)}$$

    由于$J$是未知量，所以没办法通过$(5)$来求$r_{\max }$。作者提出了一个假设，即一条雾线中，离大气光最远的像素点为无雾的像素点（该假设不一定总是成立，例如天空区域像素点的聚类。后面有处理不成立情况下的方法）。所以定义${\hat{r}}_{\max }$为：
r ^ max ⁡ ( x ) = max ⁡ x ∈ H { r ( x ) } (7) \hat{r}_{\max}(\boldsymbol{x}) = \underset{x \in H}{\max}\{r(\boldsymbol{x})\} \tag{7} r^max​(x)=x∈Hmax​{r(x)}(7)

    结合$(6)(7)$，得到：
$$\tilde{t}(\mathbf{x})=\frac{r(\mathbf{x})}{{\hat{r}}_{\max }(\mathbf{x})}\text{\hspace{1em}(8)}$$

3.3 Regularization

    大气散射模型给出了传输率的下限(Lower Bound)：
t L B ( x ) = 1 − min ⁡ c ∈ { R , G , B } { I c ( x ) A c } (9) t_{LB}(\boldsymbol{x}) = 1 - \underset{c \in \{R,G,B\}}{\min} \{\frac{I_c(\boldsymbol{x})}{A_c}\} \tag{9} tLB​(x)=1−c∈{R,G,B}min​{Ac​Ic​(x)​}(9)

    该公式的推导过程：因为$J\ge 0$，所以有：
$$t(\mathbf{x})=\frac{\mathbf{A}-\mathbf{I}(\mathbf{x})}{\mathbf{A}-\mathbf{J}(\mathbf{x})}\ge \frac{\mathbf{A}-\mathbf{I}(\mathbf{x})}{\mathbf{A}}=1-\frac{\mathbf{I}(\mathbf{x})}{\mathbf{A}}$$

    而$t(x)$在上式中是一个三个通道的向量，实际上是一个数值，且目标是求$t(x)$的下限，于是得到了$(9)$的$t_{LB}(x)$。

    将$(9)$的传输率下限约束加入到$(8)$中，得到：
t ~ L B ( x ) = max ⁡ { t ~ ( x ) , t L B ( x ) } \tilde{t}_{LB}(\boldsymbol{x}) = \max\{\tilde{t}(\boldsymbol{x}),t_{LB}(\boldsymbol{x})\} t~LB​(x)=max{t~(x),tLB​(x)}

    最小化如下目标函数，得到最终的$\hat{t}(\mathbf{x})$：
$$\sum _x\frac{{\left[\hat{t}(\mathbf{x})-{\tilde{t}}_{LB}(\mathbf{x})\right]}^2}{{\sigma }^2(\mathbf{x})}+\lambda \sum _{\mathbf{x}}\sum _{\mathbf{y}\in N_x}\frac{[\hat{t}(\mathbf{x})-\hat{t}(\mathbf{y}){]}^2}{\Vert \mathbf{I}(\mathbf{x})-\mathbf{I}(\mathbf{y}){\Vert }^2}\text{\hspace{1em}(10)}$$

    分成两部分来分析上述目标函数：
    1) 数据项
$$\sum _x\frac{{\left[\hat{t}(\mathbf{x})-{\tilde{t}}_{LB}(\mathbf{x})\right]}^2}{{\sigma }^2(\mathbf{x})}$$

    $\sigma (\mathbf{x})$为${\tilde{t}}_{LB}$的标准差。$\sigma (\mathbf{x})$发挥了重要作用，由于$\sigma (\mathbf{x})$的限制，仅将估计值应用于假设成立的像素。当方差高时，初始估计的可靠性较低。$\sigma (\mathbf{x})$随着雾线中像素数的减少而增加。 当给定雾线中的半径分布较小时，雾线假设将不成立，因为没有观察到雾度不同的像素。 在这种情况下，$\sigma (\mathbf{x})$也增加。
    将数据项写成如下形式有助于理解：
$$\sum _{\mathbf{x}}\frac{1}{{\sigma }^2(\mathbf{x})}[\hat{t}(\mathbf{x})-{\tilde{t}}_{LB}(\mathbf{x}){]}^2$$

    以下两种情况假设不成立：1. 雾线像素点少；2. 雾线$r(\mathbf{x})$分布较小。在这两种情况下，${\tilde{t}}_{LB}(\mathbf{x})$不成立，所以不需要$\hat{t}(\mathbf{x})$和其相接近，即不需要最小化$[\hat{t}(\mathbf{x})-{\tilde{t}}_{LB}(\mathbf{x}){]}^2$该项。但总体还是需要最小化数据项，作者通过添加权重系数项来达到目的。在上述两种情况下，${\sigma }^2(\mathbf{x})$会增加，$1/{\sigma }^2(\mathbf{x})$会减小，整体数据项也就会最小化。
    这里是我搞不懂的地方，即为什么这两种情况下，${\sigma }^2(\mathbf{x})$会增加。${\sigma }^2(\mathbf{x})$是关于${\tilde{t}}_{LB}(\mathbf{x})$的方差，${\tilde{t}}_{LB}(\mathbf{x})$的表达式如下：
t ~ L B ( x ) = max ⁡ { r ( x ) r ^ max ⁡ ( x ) , 1 − min ⁡ c ∈ { R , G , B } { I c ( x ) A c } } \tilde{t}_{LB}(\boldsymbol{x}) = \max \{ \frac{r(\boldsymbol{x})}{\hat{r}_{\max}(\boldsymbol{x})}, 1 - \underset{c \in \{R,G,B\}}{\min} \{\frac{I_c(\boldsymbol{x})}{A_c}\} \} t~LB​(x)=max{r^max​(x)r(x)​,1−c∈{R,G,B}min​{Ac​Ic​(x)​}}

    第一，雾线像素点减少会导致${\sigma }^2(\mathbf{x})$增加。但是样本数量和方差并没有直接关系；第二，雾线$r(\mathbf{x})$分布减小会导致${\sigma }^2(\mathbf{x})$增加。但是分布减小，不应该是方差也减小吗？因为方差就是衡量样本的在期望周围的分散程度，分布减小，也就没那么分散，换句话说，就是更聚集了。这是还没搞懂的地方，以后搞懂了，会回来填上这个坑的。

    2) 平滑项
$$\lambda \sum _{\mathbf{x}}\sum _{\mathbf{y}\in N_x}\frac{[\hat{t}(\mathbf{x})-\hat{t}(\mathbf{y}){]}^2}{\Vert \mathbf{I}(\mathbf{x})-\mathbf{I}(\mathbf{y}){\Vert }^2}$$

    其中，$N_{\mathbf{x}}$为$\mathbf{x}$周围的像素点集合。当深度连续时，传输率应该相同，而一般相邻像素的深度是连续的，所以最小化相邻像素的传输率的差值，可以得到更平滑的传输率图。

3.4 Dehazing

    经过上述的求解，已经得到了传输率以及大气光值，最后一步直接使用大气散射模型得到去雾后的图像：
J ^ ( x ) = { I ( x ) − [ 1 − t ^ ( x ) ] A / t ^ ( x ) } (11) \hat{J}(\boldsymbol{x}) = \{ \boldsymbol{I}(\boldsymbol{x}) - [1 - \hat{t}(\boldsymbol{x})] \boldsymbol{A} / \hat{t}(\boldsymbol{x}) \} \tag{11} J^(x)={I(x)−[1−t^(x)]A/t^(x)}(11)

3.5 中间结果

    图4(a)和(b)分别为原图和去雾结果。(c)为每个像素点的$r(\mathbf{x})$，(d)为每个像素点所属雾线的${\hat{r}}_{\max }(\mathbf{x})$。(e)为每条雾线的${\hat{r}}_{\max }(\mathbf{x})$在图中的像素点。(f)滤除了雾线中具有最大半径，但$\sigma (\mathbf{x})＞\mathrm{２}$的像素，因为模型假设不适用于这些雾线。 由于RGB空间与空气的距离非常短，因此上述像素位于天空中。 因此，由于噪声，根据它们的角度将它们聚类是不可靠的。 在正则化步骤中，通过数据项权重$1/{\sigma }^2(\mathbf{x})$考虑到了这一事实，如(f)所示（暖色表示较高的值）。(g)是初步传输率估计图，(h)为最终传输率估计图。

图4 中间结果

4. 实验结果

图5 去雾结果

5. 总结

    本文介绍了一种用于单图像去雾方法，名为Non-Local方法。 该方法基于这样的假设，即可以仅用数百种不同的颜色来表示图像。 在RGB空间中，这对应于数百个紧密的颜色簇。在有雾图像中，这些紧密的颜色簇会由于雾而发生变化，并在RGB空间中形成穿过大气光坐标的雾线。本文提出了一种有效的算法来识别这些雾度线并基于它们估计每个像素的传输率。 同时在正则化过程中考虑了估计的方差，因此只有符合模型假设的像素才有助于结果。

6. 读后感

    和Meng的那篇文章一样，不知道第几次拿起这篇文章了。现在终于把这篇文章搞懂了。阅读完后发现，该方法和Meng的方法有着异曲同工之妙。一个是边界限制，一个是雾线。同样本文提出的方法，也是先求出一个粗糙的传输率图，再通过优化目标函数，得到最终的精细的传输率图。