【论文阅读】LIME：Low-light Image Enhancement via Illumination Map Estimation(笔记最全篇）

今天要整理的一篇论文是使用传统的方法解决低照度的问题，Guo等人于2016年发表在IEEE上，论文《LIME：Low-light Image Enhancement via Illumination Map Estimation》，DOI：10.1109/TIP.2016.2639450

Abstract

在低照度环境下拍摄的图像通常能见度都很低，这些图像除了在视觉效果上降低了美感以外，还让计算机视觉的显示效果降质了。为了解决这个问题，本文提出了一种简单有效的低照度图像增强算法（Low-light Image Enhancement，LIME）。具体来说，输入一张低照度图像，选取每个像素通道中的最大值初始化该图像光照图，然后通过强加入一种结构先验来细化这个初始光照图，最后根据Retinex理论合成增强图。实验表明，该算法较之前的方法在质量和性能上都取得不错的效果。

Introduction

能见度高的图像能反映出场景中清晰的细节，这对于一些采用视觉技术的工作有很大的影响，如物体检测和追踪等。但是，在低照度的环境下捕获的图像通常是低能见度，虽然目前有很多算法能提升图像亮度，但是过度提亮会有损图像的效果。直接放大低照度图像可能是调节暗区可见性最直观和最简单的方法，但此操作也会带来另一个问题，比如相对较亮的区域可能已经饱和，因此丢失了相应的细节。

直方图均衡化（HE） 策略是一种将原图像通过某种变换，得到一幅灰度直方图为均匀分布的新图像的方法，可以避免上面的问题。它的基本思想是对在图像中像素个数多的灰度级进行展宽，而对像素个数少的灰度级进行缩减（将像素灰度级范围规定在0-1），从而达到清晰图像的目的。直方图均衡化的核心是变换函数，而常用的变换函数有线性变换、分段线性变换、伽马校正，对数变换、幂变换等。在文中主要指出伽马校正的缺点：仅对每个像素进行操作，忽略了像素之间的相邻关系。

Retinex理论的理论依据依赖一个核心假设：（彩色）图像可以被分解成入射分量（光照）和反射分量两个主要成分，表达式如下：
$$\begin{array}{r}I=T\cdot R\end{array}\text{\hspace{1em}(1)}$$
其中$I$是原始低照度图，T是入射分量即光照，R是反射分量。最早基于Retinex理论的研究有SSR、MSR：它们都是把反射分量作为最后的增强结果，因此，只需要求解出入射分量即可。但是这样的处理方式，得到的结果通常是不自然的，且容易过度增强。往后的研究就偏向于优化入射分量，使其更自然，有通过融合多张图像来调整光照图，也有通过划分子块优化光照图等。最后来，研究的方向就偏向于加入权重模型，同时优化光照图和反射分量。

到这里，值得一提的方法是：Dong发现了反转低照度图像跟有雾图像的特征相似，因此提出了反转图像（作为伪雾图）应用去雾算法，同样能得到增强亮度的效果，这就将低照度增强领域扩展至图像去雾领域中去了，去雾算法可用在低照度增强中去。（对这个算法感兴趣的同学，可以去下载文章看看，这里提供论文DOI：10.1145/1836845.1836920）

（图）伪雾图与有雾图像

Contribution

本文提出的算法也是基于Retinex理论的，算法主要倾向于估计低照度图像的光照图来实现增强。这里值得注意的是，本文的算法跟以往的基于Retinex增强方法（直接将图像分解成L和R）不同，我们的算法是仅仅估计一种一个变量（光照），这大大地缩小解空间和减少计算量。本文的主要贡献是：

1. 提出了两种加速优化光照图的算法：增广拉格朗日乘数法（ALM）和权重策略。
2. 进行了大量的对比实验。

个人观点： 该文章最大的贡献其实是提出的第一种加速优化方法，他主要将优化光照的问题嵌入到最优化问题中，考虑图像的保真和结构性，自定义了一个优化的目标函数和相应的约束条件，然后采用传统的拉格朗日乘数法来解决这个优化问题。

Method

根据Retinex理论可知，我们理想的图像是R（入射分量）可以通过估计出光照图T和原始的低照度$I$得到：
$$\begin{array}{r}R=\frac{I}{T}\end{array}$$
为了简化计算，通常认为三通道（彩色）图像都是共享同一个光照图。首先初始化光照图$T$：
$$\begin{array}{r}{T}^c(x)=\underset{c\in \{R,G,B\}}{\max }{I}^c(x)\end{array}\text{\hspace{1em}(2)}$$
要保证这个的光照图不会使得增强的图像过于饱和，我们加入了一个非常小的参数$ϵ$来限制：
$$\begin{array}{r}R(x)=\frac{I(x)}{T(x)+ϵ}\end{array}\text{\hspace{1em}(3)}$$
$ϵ$是一个非常小的常量，主要是为了避免出现除零和R过度饱和的情况。文章明确地指出这项工作的目标是：非均匀地增强弱光图像的照明，而不是消除光源引起的色移。

在文中，作者还特意与Dong提出的反转图像去雾增强算法作对比。我觉得这里的工作整理得挺好的，这里帮助我更好地理解了基于大气物理散射模型的去雾算法和基于Retinex理论增强算法之间能融合及相似的地方。所以，我在整理这篇文章的笔记的时候，还是决定把这部分写上。

Dong使用反转低照度图像$(1-I)$作为伪雾图，使用大气散射模型（4）来实现增强：
$$\begin{array}{r}1-I=(1-R)\cdot T+\alpha (1-T)\end{array}\text{\hspace{1em}(4)}$$
在大气物理模型（4）理论中，T代表透射率，需要估计的，$\alpha$代表大气光强，是一个确定的常量。而（1-I）代表伪雾图，（1-R）代表去雾后的图像，即增强的图像。

式子（4）和（1）进行对比，其实有一定的相似性，它们之间存在一种关系。这需要引入著名的何凯明博士提出的暗通道先验理论(Dark Channel Prior，DCP)，即$I^{dark}={\min }_c{I}^c$，清晰无雾图的暗通道为0，使用暗通道先验理论来估计透射率，有：

$$\begin{array}{r}{T}^{\prime }(x)\leftarrow 1-\underset{c}{\min }\frac{1-I^c(x)}{\alpha }=1-\frac{1}{\alpha }+\underset{c}{\max }\frac{I^c(x)}{\alpha }\end{array}\text{\hspace{1em}(5)}$$
把（5）代入到（4）中，得到：
$$\begin{array}{r}R(x)=\frac{I(x)-1+\alpha }{(1-\frac{1}{\alpha }+{\max }_c\frac{I^c(x)}{\alpha }+ϵ)}+(1-\alpha )\end{array}\text{\hspace{1em}(6)}$$
当$\alpha =1$时，我们可以看到，(3)和(6)的结果是一样的；但是当$\alpha$远离1，模型(3)和(6)的等价关系就不存在了。

这个图展示了(3)和(6)两个模型的结果，虽然将大气光强设置成0.95，两个模型的效果仍然是非常明显的。本文主要是使用（3）的方式来获取增强图。经（2）获得初始化的光照图，然后使用两种不同的方法优化光照图，下面开始介绍文中的两种加速算法。

Speed-up Method(1)：ALM

ALM加速法，是将优化光照的问题归纳为最优化问题中，考虑图像的保真度和结构平滑度，自定义了一个优化的目标函数，使用F范数和L1正则分别衡量图像的保真度和结构平滑度：
$$\begin{array}{r}\underset{T}{\min }||T^{\prime }-T|{|}_F^2+k||W\cdot ▽T|{|}_1\end{array}\text{\hspace{1em}(7)}$$

这里的k是平衡F范数和L1正则的系数。(7)式的前项为F范数，考虑的是细化后的光照图L‘与初始化的光照图L之间的保真度，后项为L1正则，考虑的是平滑度，W为权重矩阵，$▽T$代表一阶导数滤波器，包括水平和垂直方向的。

传统上，可以通过交替方向乘子算法（ADMM）就能有效地解决问题（8）。 但从（8）中的目标可以看出，两个约束项都涉及T。这里引入了一个辅助变量G来代替$▽T$，以使问题可分离并因此易于解决。把$G=▽T$作为一项等式约束条件，最后我们的优化问题就等价于：
$$\begin{array}{r}\underset{T,G}{\min }||T^{\prime }-T|{|}_F^2+k||W\cdot G|{|}_1\\ s.t.▽T=G\end{array}\text{\hspace{1em}(8)}$$
根据优化问题的条件（关于优化问题，我也有整理，详情请看博文：梳理凸优化问题），这属于典型的带有等式约束优化，可使用增广拉格朗日函数（Augmented Lagrangain Method）来求解。根据定义，给约束项加上一个惩罚项，构造增广拉格朗日函数：
$$\begin{array}{r}L=||T^{\prime }-T|{|}_F^2+k||W\cdot G|{|}_1+\varphi (Z,▽T-G)\end{array}\text{\hspace{1em}(9)}$$
其中$\varphi (Z,▽T-G)=\frac{\mu }{2}||▽T-G|{|}_F^2+<Z,▽T-G>$，这里的<,>是指矩阵的内积，Z代表拉格朗日乘子。因为函数有三个变量：T、G、Z需要求解，这篇文章的作者采用分解子问题，分别使用ALM求解。

子问题 T

将函数（9）中有关T的部分抽离出来，单独优化：
$$\begin{array}{r}{T}^{(t+1)}\leftarrow \underset{T}{\mathrm{argmin}}||T^{\prime }-T|{|}_F^2+\varphi (Z^{(t)},▽T-G^{(t)})\end{array}\text{\hspace{1em}(10)}$$
式(10)是一个经典的最小二乘问题。

补充 最小二乘法问题的步骤：
（未知数>已知数）超定方程：$xb=y$
最小化残差平方和：$b^{\prime }=\mathrm{argmin}(||xb-y|{|}^2)$
对其进行位符求最值：$b^{\prime }=x^{T}y$
如果矩阵$x^{T}x$非奇异，则b有唯一解：$b^{\prime }=(x^{T}x{)}^{-1}{x}^{T}y$

因此T子问题的解决方案，直接对(10)求微分，令其为0，得到：
$$\begin{array}{rl}& 2(T-T^{\prime })+{\mu }^{(t)}{D}^T(DT-G^{(t)})+D^T{Z}^{(t)}=0\\ \Rightarrow & (2I+{\mu }^{(t)}{D}^{T}D)=2T^{\prime }+{\mu }^{(t)}{D}^T(G^{(t)}-\frac{Z^{(t)}}{{\mu }^{(t)}})\end{array}\text{\hspace{1em}(11)}$$
这里的$I$是指单位矩阵，$D$是具有前向差异的离散梯度算子的托普利兹矩阵(Toeplitz matrice)，包含了水平$D_h$和垂直$D_v$方向。为方便起见，操作$DX$和$D^{T}X$分别表示reshape($Dx$) 和reshape($D^{T}x$)，其中$x$是矢量化的$X$，reshape()表示将矢量整形为其矩阵形式的操作。 直接计算$2I+\mu (t)D^{T}D$的倒数是解决子问题T的直观方法。 但是，矩阵逆运算在计算上是庞大的，特别是对于像$D^{T}D$这样的大型矩阵而言。 幸运的是，可以通过假设循环边界条件，我们可以对上述问题应用二维FFT技术（快速傅里叶变换），这使我们能够快速计算解，有：
$$\begin{array}{r}{T}^{(t+1)}\leftarrow F^{-1}(\frac{F(2T^{\prime }+{\mu }^{(t)}{D}^T(G^{(t)}-\frac{Z^{(t)}}{{\mu }^{(t)}}))}{2+{\mu }^{(t)}{\sum }_{d\in \{h,v\}}\widehat{F(D_d)}\circ F(D_d)\}})\end{array}\text{\hspace{1em}(12)}$$

这里的$F()$代表二维快速傅里叶变换，$F^{-1}()$和$\widehat{F()}$分别代表二维反傅里叶变换和$F()$的复合共轭。这些运算都是对每个像素点的。

子问题 G

同样地，去掉(9)中与G无关的项，得到：
$$\begin{array}{r}{G}^{(t+1)}\leftarrow \underset{G}{\mathrm{argmin}}k||W\cdot G|{|}_1+\varphi (Z^{(t)},▽T^{(t)}-G)\end{array}\text{\hspace{1em}(13)}$$
通过执行收缩操作 （我也不太知道这个操作），可以轻松获得（13）的闭式解：
$$\begin{array}{r}{G}^{(t+1)}=S_{\frac{kW}{{\mu }^{(}t)}}[▽T^{(t+1)}+\frac{Z^{(t)}}{{\mu }^{(t)}}]\end{array}\text{\hspace{1em}(14)}$$
其中$S_{\epsilon >0}[]$代表收缩操作（原文：shrinkage operator），定义为：$S_{\epsilon }[x]=sgn(x)\max (|x|-\epsilon ,0)$，收缩算子对向量和矩阵的扩展是简单地逐个元素地处理数据，例如$S_A[X]$使用由A的相应项给定的阈值对X的元素执行收缩。

子问题 Z和$\mu$

更新Z和$\mu$可以通过下面的操作得到：
$$\begin{array}{rl}{Z}^{(t+1)}& \leftarrow Z^{(t)}+{\mu }^{(t+1)}(▽T^{(t+1)}-G^{(t+1)});\\ {\mu }^{(t+1)}& \leftarrow {\mu }^{(t)}\rho ,\rho >1\end{array}\text{\hspace{1em}(14)}$$
ALM算法中概述了问题(7）的精确求解器的整个过程。这个迭代更新操作的停止条件是$||▽T^{(t+1)}-G^{(t+1)}|{|}_F<\delta ||T^{\prime }|{|}_F$，这里的$\delta =10^{-5}$；或者设置的最大迭代次数已经到达。

到这里，已经介绍完ALM加速优化算法了，最后把算法的流程图放上来，就更清晰。

对于这一部分，原文中还分析了使用ALM算法的两个关键指标：收敛性和计算复杂度。在这里就不详细写了，有兴趣的可以翻看原文。下面介绍第二种加速算法。

Speed-up Method(2)：权重变量

尽管第一种优化方法的复杂度已经是相当低，但我们还是希望能进一步降低它。 让我们再回看一下问题（7）， 带来迭代过程的因素是稀疏加权梯度项，即$||W◦▽T|{|}_1$。 L1正则项与在T上求梯度的运算一起使得式子有些复杂，但是单独拆分出来分析，以下关系是成立：
$$\begin{array}{r}\underset{ϵ\to 0^{+}}{\lim }\sum _x\sum _{d\in \{h,v\}}\frac{W_d(x)({▽}_{d}T(x){)}^2}{|{▽}_d{T}^{\prime }(x)|+ϵ}\end{array}\text{\hspace{1em}(16)}$$
根据(16),我们可以使用${\sum }_x{\sum }_{d\in \{h,v\}}\frac{W_d(x)({▽}_{d}T(x){)}^2}{|{▽}_d{T}^{\prime }(x)|+ϵ}$近似等于$||W◦▽T|{|}_1$。所以，(7)可以被写成：
$$\begin{array}{r}\underset{T}{\min }||T^{\prime }-T|{|}_F^2+k\sum _x\sum _{d\in \{h,v\}}\frac{W_d(x)({▽}_{d}T(x){)}^2}{|{▽}_d{T}^{\prime }(x)|+ϵ}\end{array}\text{\hspace{1em}(17)}$$
尽管目标函数发生了变化，但与原始函数相比，从初始光照估计值T中提取照明结构的目标与理论是一致的。更具体地，当$|{▽}_{d}T(x)|$小时，$|{▽}_d{T}^{\prime }(x)|$将被抑制，值$\frac{W_d(x)({▽}_{d}T(x){)}^2}{|{▽}_d{T}^{\prime }(x)|+ϵ}$也将被抑制。 换句话说，目标T’被约束以避免在最初估计的照明图具有小幅度的梯度。 相反，如果$|{▽}_{d}T(x)|$强，则上述抑制得到缓解，因为该位置更可能在结构边界上而不是规则纹理上。

从上面的分析可以看出，问题（17）仅涉及二次项。 因此，可以通过解决以下问题直接解决此问题：
$$\begin{array}{r}(I+\sum _{d\in \{h,v\}}{D}_d^{T}Diag(w_d^{\prime })D_d)t=t^{\prime }\end{array}\text{\hspace{1em}(18)}$$
这里的$w_d^{\prime }$是$W{’}_d和W_d^{\prime }(x)\leftarrow \frac{W_d(x)}{|{▽}_d{T}^{\prime }(x)|+ϵ}$的矢量版本，操作$Diag(x)$是使用向量x构造对角矩阵。$(I+{\sum }_{d\in \{h,v\}}{D}_d^{T}Diag(w_d^{\prime })D_d)$是一个对称的正定拉普拉斯矩阵。这里，就只要优化算法消除迭代t就可以了，大大地减少了计算量。

对于初始光照图上的结构感知细化，关键是W的设计。在这一部分中，我们讨论以下三种可能的加权策略。

• 策略1： 将权重矩阵设置为：
  $$\begin{array}{r}{W}_h(x)\leftarrow 1;W_v(x)\leftarrow 1\end{array}\text{\hspace{1em}(20)}$$
  这使得问题（7）转变为经典的L2损失总方差最小化问题。
• 策略2： 将初始照明图的梯度用作权重是合理的，因此产生：
  $$\begin{array}{rl}{W}_h(x)& \leftarrow \frac{1}{|{▽}_h{T}^{\prime }(x)|+ϵ};\\ W_v(x)& \leftarrow \frac{1}{|{▽}_v{T}^{\prime }(x)|+ϵ}\end{array}\text{\hspace{1em}(21)}$$
• 策略3： 受相对总变异（RTV）的启发，对于每个位置，权重均通过以下方式设置：
  $$\begin{array}{rl}{W}_h(x)& \leftarrow \sum _{y\in \Omega (x)}\frac{G_{\sigma }(x,y)}{|{\sum }_{y\in \Omega (x)}{G}_{\sigma }(x,y){▽}_h{T}^{\prime }(x)|+ϵ};\\ W_v(x)& \leftarrow \sum _{y\in \Omega (x)}\frac{G_{\sigma }(x,y)}{|{\sum }_{y\in \Omega (x)}{G}_{\sigma }(x,y){▽}_v{T}^{\prime }(x)|+ϵ}\end{array}\text{\hspace{1em}(22)}$$
  这里的$G_{\sigma }(x,y)$是由具有标准偏差σ的高斯核产生。 形式上，$G_{\sigma }(x,y)$表示为：
  $$\begin{array}{r}{G}_{\sigma }(x,y)\propto \exp (-\frac{dist(x,y)}{2{\sigma }^2})\end{array}\text{\hspace{1em}(20)}$$
  其中，函数dist(x, y)用于测量位置x和y之间的空间欧几里得距离。 实际上，第二种加权策略就是其中一种。 当$\sigma \to 0^{+}$，这两种策略获得相同的权重矩阵。 我们注意到，与RTV不同，我们的权重矩阵是基于给定的T’构造的，而不是根据T进行迭代更新的，这意味着W只需要计算一次。

实验

最后，LIME算法的整体流程是：

放上一些实验结果：

conclusion

本文主要整理了低照度图像增强中的一种传统方法：LIME。呼~好累，不多说啦。最后按套路放上文章资源：LIME: Low-light Image Enhancement via Illumination Map Estimation
如果文章对您有帮助（呜呜呜，可怜一下我码了这么多字），点个赞再走呗！！！谢谢啦（试图出卖博主样貌，来收获欢心。。。。太难了）
尊重原创，转载请注明出处！