【论文阅读】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{aligned}I=T\cdot R \tag1\end{aligned}$I=T⋅R​(1)
其中$I$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$I得到：
$\begin{aligned}R=\frac IT \end{aligned}$R=TI​​
为了简化计算，通常认为三通道（彩色）图像都是共享同一个光照图。首先初始化光照图$T$T：
$\begin{aligned}T^c(x)=\max_{c\in\{R,G,B\}}I^c(x) \tag2\end{aligned}$Tc(x)=c∈{R,G,B}max​Ic(x)​(2)
要保证这个的光照图不会使得增强的图像过于饱和，我们加入了一个非常小的参数$\epsilon$ϵ来限制：
$\begin{aligned}R(x)=\frac {I(x)}{T(x)+\epsilon} \tag3\end{aligned}$R(x)=T(x)+ϵI(x)​​(3)
$\epsilon$ϵ是一个非常小的常量，主要是为了避免出现除零和R过度饱和的情况。文章明确地指出这项工作的目标是：非均匀地增强弱光图像的照明，而不是消除光源引起的色移。

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

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

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

$\begin{aligned} T'(x) \gets 1-\min_c\frac{1-I^c(x)}{\alpha}=1-\frac1\alpha+\max_c\frac{I^c(x)}{\alpha} \tag5\end{aligned}$T′(x)←1−cmin​α1−Ic(x)​=1−α1​+cmax​αIc(x)​​(5)
把（5）代入到（4）中，得到：
$\begin{aligned} R(x)=\frac{I(x)-1+\alpha}{(1-\frac1\alpha+\max_c\frac{I^c(x)}{\alpha}+\epsilon)}+(1-\alpha) \tag6\end{aligned}$R(x)=(1−α1​+maxc​αIc(x)​+ϵ)I(x)−1+α​+(1−α)​(6)
当$\alpha=1$α=1时，我们可以看到，(3)和(6)的结果是一样的；但是当$\alpha$α远离1，模型(3)和(6)的等价关系就不存在了。

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

Speed-up Method(1)：ALM

ALM加速法，是将优化光照的问题归纳为最优化问题中，考虑图像的保真度和结构平滑度，自定义了一个优化的目标函数，使用F范数和L1正则分别衡量图像的保真度和结构平滑度：
$\begin{aligned} \min_T||T'-T||^2_F+k||W\cdot \triangledown T||_1 \tag7\end{aligned}$Tmin​∣∣T′−T∣∣F2​+k∣∣W⋅▽T∣∣1​​(7)

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

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

子问题 T

将函数（9）中有关T的部分抽离出来，单独优化：
$\begin{aligned} T^{(t+1)}\gets \argmin_T||T'-T||^2_F+\phi(Z^{(t)},\triangledown T-G^{(t)}) \tag{10}\end{aligned}$T(t+1)←Targmin​∣∣T′−T∣∣F2​+ϕ(Z(t),▽T−G(t))​(10)
式(10)是一个经典的最小二乘问题。

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

因此T子问题的解决方案，直接对(10)求微分，令其为0，得到：
$\begin{aligned} &2(T-T')+\mu^{(t)}D^T(DT-G^{(t)})+D^TZ^{(t)}=0\\ \Rightarrow & (2I+\mu^{(t)}D^TD)=2T'+\mu^{(t)}D^T(G^{(t)}-\frac{Z^{(t)}}{\mu^{(t)}}) \tag{11}\end{aligned}$⇒​2(T−T′)+μ(t)DT(DT−G(t))+DTZ(t)=0(2I+μ(t)DTD)=2T′+μ(t)DT(G(t)−μ(t)Z(t)​)​(11)
这里的$I$I是指单位矩阵，$D$D是具有前向差异的离散梯度算子的托普利兹矩阵(Toeplitz matrice)，包含了水平$D_h$Dh​和垂直$D_v$Dv​方向。为方便起见，操作$DX$DX和$D^TX$DTX分别表示reshape($Dx$Dx) 和reshape($D^Tx$DTx)，其中$x$x是矢量化的$X$X，reshape()表示将矢量整形为其矩阵形式的操作。 直接计算$2I + \mu(t)D^TD$2I+μ(t)DTD的倒数是解决子问题T的直观方法。 但是，矩阵逆运算在计算上是庞大的，特别是对于像$D^TD$DTD这样的大型矩阵而言。 幸运的是，可以通过假设循环边界条件，我们可以对上述问题应用二维FFT技术（快速傅里叶变换），这使我们能够快速计算解，有：
$\begin{aligned} T^{(t+1)}\gets F^{-1}(\frac{F(2T'+\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)\}}) \tag{12}\end{aligned}$T(t+1)←F−1(2+μ(t)∑d∈{h,v}​F(Dd​)​∘F(Dd​)}F(2T′+μ(t)DT(G(t)−μ(t)Z(t)​))​)​(12)

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

子问题 G

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

子问题 Z和$\mu$μ

更新Z和$\mu$μ可以通过下面的操作得到：
$\begin{aligned} Z^{(t+1)}&\gets Z^{(t)}+\mu^{(t+1)}(\triangledown T^{(t+1)}-G^{(t+1)});\\ \mu^{(t+1)}&\gets\mu^{(t)}\rho,\rho\gt1 \tag{14}\end{aligned}$Z(t+1)μ(t+1)​←Z(t)+μ(t+1)(▽T(t+1)−G(t+1));←μ(t)ρ,ρ>1​(14)
ALM算法中概述了问题(7）的精确求解器的整个过程。这个迭代更新操作的停止条件是$||\triangledown T^{(t+1)}-G^{(t+1)}||_F\lt \delta||T'||_F$∣∣▽T(t+1)−G(t+1)∣∣F​<δ∣∣T′∣∣F​，这里的$\delta=10^{-5}$δ=10−5；或者设置的最大迭代次数已经到达。

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

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

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

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

从上面的分析可以看出，问题（17）仅涉及二次项。 因此，可以通过解决以下问题直接解决此问题：
$\begin{aligned} (I+\sum_{d\in\{h,v\}}D^T_dDiag(w'_d)D_d)t=t' \tag{18}\end{aligned}$(I+d∈{h,v}∑​DdT​Diag(wd′​)Dd​)t=t′​(18)
这里的$w'_d$wd′​是$W’_d和W'_d(x)\gets \frac{W_d(x)}{|\triangledown_dT'(x)|+\epsilon }$W’d​和Wd′​(x)←∣▽d​T′(x)∣+ϵWd​(x)​的矢量版本，操作$Diag(x)$Diag(x)是使用向量x构造对角矩阵。$(I+\sum_{d\in\{h,v\}}D^T_dDiag(w'_d)D_d)$(I+∑d∈{h,v}​DdT​Diag(wd′​)Dd​)是一个对称的正定拉普拉斯矩阵。这里，就只要优化算法消除迭代t就可以了，大大地减少了计算量。

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

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

实验

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

放上一些实验结果：

conclusion

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