Structure-Revealing联合去噪增强模型论文解读

一、文章摘要概述

文章的题目是：
《Structure-Revealing Low-Light Image Enhancement Via Robust Retinex Model》
这是一篇2018年6月份的TIP(视觉顶刊)会议文章，文章针对弱光(低照度)图像存在密集噪声的问题，提出一种基于Rubost Retinex分解模型的
Structure-revealing弱光增强方法，考虑图像存在的噪声项，并且为了有效地解决模型优化问题，文章提供了一种基于拉格朗日乘子的交替方向最小化(ADM)算法代替对数变换。另外文章指出，该算法还可用于处理如用于水下或遥感的图像增强，以及在模糊或沙尘条件下的图像增强问题。与前几篇JED(联合增强和去噪)文章相比，这篇文章重点介绍了ADM在序列迭代过程中的具体计算过程。总结一下，文章在使用传统方法做联合图像增强和去噪时，做了3项创新：

• 在经典Retinex模型的基础上考虑噪声项的影响，提出一种Robust Retinex模型，并在模型基础上，首次尝试预测图像噪声，同时估计Structure-revealed 反射图和分段平滑的照明图。
• 提出一种基于增广Lagrange乘子的ADM算法代替对数变换，有效解决了模型优化问题。
• 该算法不仅可以用在含噪声弱光环境下，还可用在如水下图像增强、遥感图像增强、图像消除灰尘和灰尘的天气图像增强中。

二、背景(先验)知识

经典Retinex模型为：
$$\mathbf{I}=\mathbf{R}\circ \mathbf{L}$$式中，$\mathbf{I}$表示原始弱光图像，$\mathbf{R}$和$\mathbf{L}$分别代表分解后的反射图和照明图。一般为减少计算量会转换到对数域计算分量。
但是，传统Retinex一方面对原始图像的光照平滑度和对比度有较高要求，另一方面没有考虑到模型噪声。因此文章提出一个Robust Retinex模型，并证明模型的优越性，公式如下：
$$\mathbf{I}=\mathbf{R}\circ \mathbf{L}+\mathbf{N}$$式中，$\mathbf{N}$表示噪声项，文章认为噪声项较为均匀的分布在原始图像中，在前人的工作中也有考虑噪声的弱光增强方法，值得借鉴：

• Elad[PDF]使用两个双边滤波器，在对数域Retinex分解上抑制照明图和反射图噪声。
• Li[PDF]在反射图估计后使用边缘保持平滑法实现去噪。
• Yu[PDF]在反射图估计后使用导向滤波平滑法实现去噪。
  而文章不使用对数变换时，通过联合优化迭代方法分析模型和算法求解，如下。

二、Structure-Revealing弱光增强模型

联合去噪增强整体流程图：
文中给出了两套分解策略，作为对比，其中一套可以作为Baseline(不含噪声项)，另外一套就是文中所提的创新方法：

1、Baseline

Retinex分解模型如下所示：
$$\underset{\mathbf{R},\mathbf{L}}{\mathrm{argmin}}{\Vert \mathbf{R}\circ \mathbf{L}-\mathbf{I}\Vert }_F^2+\beta {\Vert ▽\mathbf{L}\Vert }_1+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2$$式中，$\beta$和$\omega$表示每一项系数，${\Vert \cdot \Vert }_F$和${\Vert \cdot \Vert }_1$分别表示F-norm和1-norm，$▽$表示一阶微分算子，$\circ$表示点乘操作，$\mathbf{G}$表示原始图像$\mathbf{I}$的调整梯度，以上各项分别代表：

• ${\Vert \mathbf{R}\circ \mathbf{L}-\mathbf{I}\Vert }_F^2$：原始图像$\mathbf{I}$和重构图像$\mathbf{R}\circ \mathbf{L}$之间的保真度。
• ${\Vert ▽\mathbf{L}\Vert }_1$：总体稀疏度，和$\mathbf{L}$的区域平滑度。
• ${\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2$：反射估计图$\mathbf{R}$和$\mathbf{G}$的距离，增强反射图的结构信息。

文章提到两个可解释性问题：
（1）、针对照明估计图$\mathbf{L}$的平滑度约束问题，在以往的工作中一般是使用L2-norm，在ML和DL中，L2-norm常常用在Loss函数中，防止模型过拟合提高模型的泛化能力，和L1-norm不同的是，L2-norm是对矩阵向量各元素平方和求平方根，所以其对矩阵中突变值更敏感，对应弱光图像中照度(光照)突变的区域，L2范数强制空间平滑照明会产生区域模糊现象，也叫边缘伪影现象，因此文中使用L1-norm范数来约束照明梯度，保持照明图像的整体结构，以取得更好的视觉效果。

（2）、反射图的梯度约束问题，考虑到低对比度通常表示较低的梯度范围，所以模型第三项通过调整反射率梯度来提高整体对比度，梯度矩阵$\mathbf{G}$是通过自适应调整因子和原始图像的梯度乘积得到的，公式如下：
$$\{\begin{array}{cc}\mathbf{G}=\mathbf{K}\circ ▽\mathbf{I}& \\ \mathbf{K}=1+\lambda e^{-|▽\mathbf{I}|/\sigma }& \end{array}$$式中，$\lambda$和$\sigma$分别控制放大程度和放大率，自适应调节因子$\mathbf{K}$与梯度$▽\mathbf{I}$成反比，使得调整梯度$\mathbf{G}$具有较为均匀的梯度变化，文章展示该项在整体增强效果中的作用，如图：

从肉眼来看，整体对比度信息差别不大，但是’白色方块’处局部对比度结构的确展示的要更加完整。

2、Robust Retinex噪声优化模型分析

考虑到自然弱光图像噪声不仅仅是加性(或乘性)高斯白噪声，通过某种分布很难评估噪声水平，因此文章从输入图像$\mathbf{S}$(HSV空间)直接估计一个噪声图，故通过Robust Retinex得到增强噪声优化模型如下所示：$$\underset{\mathbf{R},\mathbf{L},\mathbf{N}}{\mathrm{argmin}}{\Vert \mathbf{R}\circ \mathbf{L}+\mathbf{N}-\mathbf{I}\Vert }_F^2+\beta {\Vert ▽\mathbf{L}\Vert }_1+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2+\delta {\Vert \mathbf{N}\Vert }_F^2$$式中，$\mathbf{N}$表示噪声图，其他项如前所述，${\Vert \mathbf{N}\Vert }_F^2$项限制整体噪声，另外修改梯度矩阵$\mathbf{G}$，如下：
$$\{\begin{array}{cc}\mathbf{G}=\mathbf{K}\circ ▽\hat{\mathbf{I}}& \\ \mathbf{K}=1+\lambda e^{-|▽\hat{\mathbf{I}}|/\sigma }& \end{array}$$其中，$▽\hat{\mathbf{I}}$是一个绝对值函数，
$$▽\hat{\mathbf{I}}=\{\begin{array}{cc}0,& if\left|\hat{\mathbf{I}}\right|<\epsilon \\ \hat{\mathbf{I}},& otherwise\end{array}$$这样，小梯度(噪声水平)在放大前就被抑制了，避免后续增强过程中对极低照度下密集噪声的方法效果。

3、噪声优化模型求解

文章表示对于非凸优化问题，ADM是一个很好的工具，在此基础上，作者提出一个基于Lagrange乘子的ADM优化算法来求解该模型，总共分为五个步骤如下：
Step 1：用辅助变量$\mathbf{T}$替换$▽\mathbf{L}$，重写等式为：
$$\underset{\mathbf{R},\mathbf{L},\mathbf{N},\mathbf{T}}{\mathrm{argmin}}{\Vert \mathbf{R}\circ \mathbf{L}+\mathbf{N}-\mathbf{I}\Vert }_F^2+\beta {\Vert \mathbf{T}\Vert }_1+\delta {\Vert \mathbf{N}\Vert }_F^2+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2$$$$s.t.\mathbf{T}=▽\mathbf{L}$$Step 2：引入拉格朗日乘子$\mathbf{Z}$来移除等式约束，得到增广拉格朗日方程：
$$l(\mathbf{R},\mathbf{L},\mathbf{N},\mathbf{T},\mathbf{Z})={\Vert \mathbf{R}\circ \mathbf{L}+\mathbf{N}-\mathbf{I}\Vert }_F^2+\beta {\Vert \mathbf{T}\Vert }_1+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2+\delta {\Vert \mathbf{N}\Vert }_F^2+\Phi (\mathbf{Z},▽\mathbf{L}-\mathbf{T})$$其中：$$\Phi (\mathbf{Z},▽\mathbf{L}-\mathbf{T})=⟨\mathbf{Z},▽\mathbf{L}-\mathbf{T}⟩+(\mu /2){\Vert ▽\mathbf{L}-\mathbf{T}\Vert }_F^2$$$⟨\cdot ,\cdot ⟩$表示矩阵内积，$\mu$表示正系数。优化方程可以通过依次迭代更新每个变量来求解，同时将上一次迭代估计的其他变量作为常量。
Step 3：第k次迭代各变量求解

• $\mathbf{R}$求解：忽略与$\mathbf{R}$无关的变量，得到优化方程：$$\underset{\mathbf{R}}{\mathrm{argmin}}{\Vert \mathbf{R}\circ {\mathbf{L}}^{(\mathbf{k})}+{\mathbf{N}}^{(\mathbf{k})}-\mathbf{I}\Vert }_F^2+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2$$将第一项向量化，可得到一个最小二乘问题方程：
  $$\underset{\mathbf{R}}{\mathrm{argmin}}{\Vert \mathbf{r}{\tilde{\mathbf{I}}}^{(\mathbf{k})}+{\mathbf{n}}^{(\mathbf{k})}-\mathbf{i}\Vert }_F^2+\omega {\Vert ▽\mathbf{R}-\mathbf{G}\Vert }_F^2$$式中，$\mathbf{I}$是矩阵$\mathbf{L}$的向量化结果，$\tilde{\mathbf{I}}$是$\mathbf{I}$的对角矩阵，在下面有同样的向量化处理($\mathbf{r}$、$\mathbf{i}$、$\mathbf{n}$、$\mathbf{t}$、$\mathbf{g}$和$\mathbf{z}$分别对应于$\mathbf{R}$、$\mathbf{I}$、$\mathbf{N}$、$\mathbf{T}$、$\mathbf{G}$和$\mathbf{Z}$的向量化结果。)
  然后，对$\mathbf{R}$微分，并令微分方程为0，可得等式：
  $$2({\tilde{\mathbf{I}}}^{(k)}{)}^{\text{T}}({\tilde{\mathbf{I}}}^{(k)}\mathbf{r}+{\mathbf{n}}^{(k)}-\mathbf{i})+2\omega {\mathbf{D}}^{\text{T}}(\mathbf{D}\mathbf{r}-\mathbf{g})=0$$移项合并同类项可得：
  $$(f({\tilde{\mathbf{I}}}^{(k)})+\omega f(\mathbf{D}))\mathbf{r}={\tilde{\mathbf{I}}}^{(k)}(\mathbf{i}-{\mathbf{n}}^{(k)})+\omega {\mathbf{D}}^{\text{T}}\mathbf{g}$$式中，$\mathbf{D}$代表离散梯度算子，$f(\mathbf{x})={\mathbf{x}}^{\text{T}}\mathbf{x}$
  因此，可得到反射估计$\mathbf{r}$的第k+i次迭代结果：
  $${\mathbf{r}}^{(k+1)}=(f({\tilde{\mathbf{I}}}^{(k)})+\omega f(\mathbf{D}){)}^{-1}({\tilde{\mathbf{I}}}^{(k)}(\mathbf{i}-{\mathbf{n}}^{(k)})+\omega {\mathbf{D}}^{\text{T}}\mathbf{g})$$这里需要注意的是F-norm的求导计算，通过查阅相关资料，得F-norm的求导公式为：${▽}_{\mathbf{x}}{\Vert \mathbf{X}\Vert }_F^2=2\mathbf{X}$，其仿射变换的求导公式为：若$Y=AX+B$（设A，B均为常量矩阵），则：${▽}_{\mathbf{x}}f(A\mathbf{X}+B)=A^{\text{T}}{▽}_{\mathbf{Y}}f(Y)$。当然这只是一个简洁的记忆方法，对于一个搞图像的研究生就没必要推导了。
• $\mathbf{L}$求解：忽略与$\mathbf{L}$无关的变量，得到优化方程：
  $$\underset{\mathbf{L}}{\mathrm{argmin}}{\Vert {\mathbf{R}}^{(k+1)}\circ \mathbf{L}+{\mathbf{N}}^{(k)}-\mathbf{I}\Vert }_F^2+\Phi ({\mathbf{Z}}^{(k)},▽\mathbf{L}-{\mathbf{T}}^{(k)})$$即：$\underset{\mathbf{L}}{\mathrm{argmin}}{\Vert {\mathbf{R}}^{(k+1)}\circ \mathbf{L}+{\mathbf{N}}^{(k)}-\mathbf{I}\Vert }_F^2+⟨{\mathbf{Z}}^{(k)},▽\mathbf{L}-{\mathbf{T}}^{(k)}⟩+(\mu /2){\Vert ▽\mathbf{L}-{\mathbf{T}}^{(k)}\Vert }_F^2$
  将其向量化，可得到一个最小二乘问题方程：$$\underset{\mathbf{L}}{\mathrm{argmin}}{\Vert {\tilde{\mathbf{r}}}^{(k+1)}\circ \mathbf{I}+{\mathbf{n}}^{(k)}-\mathbf{i}\Vert }_F^2+⟨{\mathbf{z}}^{(k)},▽\mathbf{L}-{\mathbf{t}}^{(k)}⟩+(\mu /2){\Vert ▽\mathbf{L}-{\mathbf{t}}^{(k)}\Vert }_F^2$$
  然后，对$\mathbf{L}$微分，并令微分方程为0，可得等式：$$2({\tilde{\mathbf{r}}}^{(k+1)}{)}^{\text{T}}({\tilde{\mathbf{r}}}^{(k+1)}\mathbf{I}+{\mathbf{n}}^{(k)}-\mathbf{i})+{\mathbf{D}}^T{\mathbf{Z}}^{(k)}+\mu {\mathbf{D}}^T(\mathbf{D}\mathbf{I}-{\mathbf{t}}^{(k)})=0$$移项合并同类项可得：$$(2f({\tilde{\mathbf{r}}}^{(k+1)})+\mu f(\mathbf{D}))\mathbf{I}=2{\tilde{\mathbf{r}}}^{(k+1)}(\mathbf{i}-{\mathbf{n}}^{(k)})+\mu {\mathbf{D}}^{\text{T}}({\mathbf{t}}^{(k)}-\frac{{\mathbf{z}}^{(k)}}{\mu })$$式中，$\mathbf{D}$和$f(\mathbf{X})$和上式一样，注意内积形式$⟨\cdot ,\cdot ⟩$的导数推导，因此可得到照明估计$\mathbf{I}$的第k+1次迭代结果：$${\mathbf{I}}^{(k+1)}=(2f({\tilde{\mathbf{r}}}^{(k+1)})+\mu f(\mathbf{D}){)}^{-1}(2{\tilde{\mathbf{r}}}^{(k+1)}(\mathbf{i}-{\mathbf{n}}^{(k)})+\mu {\mathbf{D}}^{\text{T}}({\mathbf{t}}^{(k)}-\frac{{\mathbf{z}}^{(k)}}{\mu }))$$原文中，${\mathbf{n}}^{(k)}$写成了${\mathbf{n}}^{(k+1)}$应该是笔误。
• $\mathbf{N}$求解：忽略与$\mathbf{N}$无关的变量，得到优化方程：
  $$\underset{\mathbf{N}}{\mathrm{argmin}}{\Vert {\mathbf{R}}^{(k+1)}\circ {\mathbf{L}}^{(k+1)}+\mathbf{N}-\mathbf{I}\Vert }_F^2+\delta {\Vert \mathbf{N}\Vert }_F^2$$由于$\mathbf{N}$和$\mathbf{R}$、$\mathbf{L}$无乘积关系，因此可直接求出该二次问题的解析解：$${\mathbf{N}}^{(k+1)}=(\mathbf{I}-{\mathbf{R}}^{(k+1)}\circ {\mathbf{L}}^{(k+1)})/(1+\delta )$$文中这里是没有转换成最小二乘问题的，我也没看懂，WoC，为什么要加个’也’呢~
• $\mathbf{T}$求解：忽略与$\mathbf{T}$无关的变量，得到优化方程：
  $$\underset{\mathbf{T}}{\mathrm{argmin}}\beta {\Vert \mathbf{T}\Vert }_1+\Phi ({\mathbf{Z}}^{(k)},▽{\mathbf{L}}^{(k+1)}-\mathbf{T})$$文章介绍可通过Shrinkage操作直接得到，辅助变量$\mathbf{T}$的第k+1次迭代结果：
  $${\mathbf{T}}^{(k+1)}={\mathbf{S}}_{\frac{\beta }{{\mu }^{(k)}}}(▽{\mathbf{L}}^{(k+1)}+\frac{{\mathbf{Z}}^{(k)}}{{\mu }^{(k)}})$$式中，${\mathbf{S}}_{\epsilon }(\mathbf{x})=sign(\mathbf{x})max(\left|\mathbf{x}\right|-\epsilon ,0)$
• Lagrange乘子$\mathbf{Z}$和系数因子$\mu$的更新：
  $$\{\begin{array}{c}{\mathbf{Z}}^{(k+1)}={\mathbf{Z}}^{(k)}+{\mu }^{(k)}(▽{\mathbf{L}}^{(k+1)}-{\mathbf{T}}^{(k+1)})\\ {\mu }^{(k+1)}={\mu }^{(k)}\rho ,\rho >1\end{array}$$以上就是文中介绍的基于Lagrange乘子的ADM第k次序列迭代过程，当达到条件：
  （1）反射估计${\mathbf{R}}^{(k)}$到${\mathbf{R}}^{(k+1)}$(或者照明估计${\mathbf{L}}^{(k)}$到${\mathbf{L}}^{(k+1)})$变化量达到某一最低阈值时，
  （2）迭代达到某一指定次数时，
  以上条件(1)或(2)达到时，即停止更新，得到反射估计图$\mathbf{R}$、照明估计图$\mathbf{L}$和噪声图$\mathbf{N}$，整个算法的执行流程如下：
  

Step 4：照明分量调整以提高整体照度
文中通过Gamma矫正调整$\mathbf{L}$，公式如下：
$$\hat{\mathbf{L}}={\mathbf{L}}^{\frac{1}{\gamma }}$$根据经验，$\gamma$值设置为2.2。
Step 5：得到HSV空间V通道的去噪增强图像并转化到RGB空间，公式如下：
$$\hat{\mathbf{I}}=\mathbf{R}\circ \hat{\mathbf{L}}$$$$\hat{\mathbf{S}}\leftarrow \hat{\mathbf{I}}$$即可得到最终联合去噪与增强效果图像$\hat{\mathbf{S}}$。

三、实验效果及源码

待续…

源码;https://github.com/martinli0822/Low-light-image-enhancement