图像降噪：Image Smoothing via L0 Gradient Minimization

Image Smoothing via $L0$ Gradient Minimization

不同于 TV(Total Variation) 使用 $L_1$作为正则化项， WLS(Weighted Least Squares) 使用 $L_2$作为正则化项，本文采用 $L_0$作为正则化项来进行图像的重建，核心思想是通过限制图像非零梯度的数量来限制或增强图像中的高对比度的边缘。本方法的主要应用场景有降噪、边缘提取、数据压缩、HDR、分割、风格化等。

1D Smoothing

记一维输入信号为$g$，其经过平滑后的结果记为$f$，计算梯度发生变化的点的个数：
c(f)=#{p∣∣fp?fp+1∣≠0}(1)c(f)=\#\left\{p|{\quad}| f_{p}-f_{p+1} \mid \neq 0\right\}\tag{1} c(f)=#{ p∣∣fp??fp+1?∣??=0}(1)
上式中$p、p+1$表示相邻样本或像素，$|f_p?f_{p+1}|$表示梯度幅度(前向差分)，#{}\#\{\}#{ }表示计数算子，用于计算满足$|f_p?f_{p+1}|?=0$的点的个数，即梯度的$L_0$，$c(f)$不考虑信号变化的幅度，只考虑信号是否发生变化，$c(f)$的计算方式是本篇论文的核心

对于信号的平滑，或者对于信号主要结果的提取，应该要对信号的相似性上进行约束，即平滑后的信号与原信号间的结构是相似的，因此将此问题的目标函数做如下表示：
$$\underset{f}{\min }\sum _p{\left(f_p?g_p\right)}^2\text{ s.t. }c(f)=k\text{\hspace{1em}(2)}$$
其中$c(f)=k$表示平滑信号$f$当中梯度不为0的点的个数，下图给出了对一个含噪一维信号进行梯度$L_0$最小化的结果，及其与其他算法对该一维信号进行平滑结果的比较。

本文方法与其他算法的比较

$$?$$

本文算法平滑信号的结果

以上结果参数$k=6$，参数$k$决定了经过平滑后的信号所存在的梯度跳变的点数，对应上图中6次梯度变化，从这一结果可以看出，基于梯度$L_0$最小化的在平滑信号的过程中其结果不受信号跳变的影响，因此其平滑后的结果，不会像$a?d$所式结果一样，平滑后的结果依旧存在一定程度的波动，并且保证了平滑后的信号依旧保持锐利的边缘，而不是缓慢过度的边缘，即大的、决定图像结构的边缘不会被平滑所影响而变得模糊，即对于一维信号而言，基于梯度 $L_0$ 最小化的平滑就是在给定梯度跳变次数 $k$ 的前提下，寻找将发生梯度变化的放在何处，能够实现平滑后的信号与原信号的差异${\sum }_p{\left(f_p?g_p\right)}^2$的最小化。此算法相较于其他平滑算法而言，一个重要的性质是不论$k$取值多小，都不会模糊掉所剩下的边缘(当$k$过小时，可能没有边缘被检出)

将公式(2)当中的两项合并为更具一般性的框架，如下：
$$\underset{f}{\min }\sum _p{\left(f_p?g_p\right)}^2+\lambda ?c(f)\text{\hspace{1em}(3)}$$
其中参数$\lambda$用于控制信号平滑的程度，$\lambda$与$k$之间的关系可以用下图表示：

参数k与参数lambda之间的关系

可以看出，$k$随着$\frac{1}{\lambda }$的增加而单调递增，在一定程度上可以将$\frac{1}{\lambda } $看做是$k$，下图$g、h$给出了同一信号$\lambda$取不同值时，信号平滑结果的差异，可以看出$\lambda$的取值影响的是信号经过平滑之后依旧存在的峰值的个数，不会影响峰值的强度。

2D Formulation

记二维输入信号为$I$，其经过平滑后的结果记为$S$，并将图像的梯度表示为$?S_p={\left({?}_x{S}_p,{?}_y{S}_p\right)}^T$，计算图像中梯度发生变化的点的个数：
C(S)=#{p∣∣?xSp∣+∣?ySp∣≠0}(4)C(S)=\#\left\{p|| \partial_{x} S_{p}|+| \partial_{y} S_{p} \mid \neq 0\right\}\tag{4} C(S)=#{ p∣∣?x?Sp?∣+∣?y?Sp?∣??=0}(4)
同样给出求解$S$的目标函数：
$$\underset{S}{\min }\left\{\sum _p{\left(S_p?I_p\right)}^2+\lambda ?C(S)\right\}\text{\hspace{1em}(5)}$$

Wang‘s Algorithm

式(5)中包含计数函数和$L_0$，难以使用梯度下降等传统方法对以上优化问题进行求解，文章引用Wang在其论文《 A new alternating minimization algorithm for total variation image reconstruction》中提出的半二次分裂(Half Quadratic Splitting)方法进行优化求解，以下的内容节取自Wang的原文中对此方法的部分描述及证明。
首先，将信号做如下表示：
$$f=K{u}^0+w\text{\hspace{1em}(w.1)}$$
用于表示模糊和噪声对信号衰退的影响，其中$f$表示输入信号，$K$表示模糊算子，$w$表示噪声，$u_0$表示原始信号
将优化求解原始信号$u_0$的目标函数做如下表示：
$$\underset{u}{\min }\sum _{i=1}^{n^2}\Vert D_{i}u\Vert +\frac{\mu }{2}\Vert Ku?f{\Vert }_2^2\text{\hspace{1em}(w.2)}$$
其中$n$表示图像尺寸，$D_{i}u$表示图像梯度，$\sum ||D_{i}u||$为$u$的全变分，即是保边缘正则化项，不可微，为了处理这一不可微函数的优化问题，使用半二次分裂，引入一个辅助变量${\mathbf{w}}_i\in {\mathbb{R}}^2$，将$D_{i}u$从$||.||$中转移出来，并将$w_i$与$D_{i}u$之间的差值作为新的惩罚项加入到目标函数当中：
$$\underset{\mathbf{w},u}{\min }\sum _i{\Vert {\mathbf{w}}_i\Vert }_2+\frac{\beta }{2}\sum _i{\Vert {\mathbf{w}}_i?D_{i}u\Vert }_2^2+\frac{\mu }{2}\Vert Ku?f{\Vert }_2^2\text{\hspace{1em}(w.3)}$$
对引入的辅助变量的要求是尽可能接近$D_{i}u$，因此惩罚项的参数$\beta$应该尽可能大，当$\beta \to \infty$时式(w.3)收敛至(w.2)，上式(w.3)相对于参数$u、w$均是凸的，如果采用固定变量的方式求解将是比较简单的。添加这个辅助变量是半二次分裂的常见做法，具体可以参见论文《Analysis of Half-Quadratic Minimization Methods for Signal and Image Recovery》。
1.分离式(w.3)中所有包含$w_i$的项，将其中的$u$当做是常量，可得优化求解$w_i$的新的目标函数：
$$\underset{{\mathbf{w}}_i}{\min }\left(\Vert {\mathbf{w}}_i\Vert +\frac{\beta }{2}{\Vert {\mathbf{w}}_i?D_{i}u\Vert }^2\right)\text{\hspace{1em}(w.4)}$$
并解得$w_i$：
$${\mathbf{w}}_i=\max \left\{\Vert D_{i}u\Vert ?\frac{1}{\beta },0\right\}\frac{D_{i}u}{\Vert D_{i}u\Vert },i=1,?,n^2\text{\hspace{1em}(w.5)}$$
2.分离式(w.3)中所有包含$u$的项，将其中的$w_i$当做是常量，可以优化求解$u$的新的目标函数：
$$\underset{\mathbf{u}}{\min }\frac{\beta }{2}\sum _i{\Vert {\mathbf{w}}_i?D_{i}u\Vert }_2^2+\frac{\mu }{2}\Vert Ku?f{\Vert }_2^2\text{\hspace{1em}(w.6)}$$
对上式求导，并令导数的每一项为0，可得正规方程：
$$\left(\sum _i{D}_i^T{D}_i+\frac{\mu }{\beta }{K}^{T}K\right)u=\sum _i{D}_i^T{\mathbf{w}}_i+\frac{\mu }{\beta }{K}^{T}f\text{\hspace{1em}(w.7)}$$
即是：
$$\left({\left(D^{(1)}\right)}^T{D}^{(1)}+{\left(D^{(2)}\right)}^T{D}^{(2)}+\frac{\mu }{\beta }{K}^{T}K\right)u={\left(D^{(1)}\right)}^T{w}_1+{\left(D^{(2)}\right)}^T{w}_2+\frac{\mu }{\beta }{K}^{T}f\text{\hspace{1em}(w.8)}$$
并使用FFT加速，结果如下(这一步的推导并没有看懂)：
$$u={\mathcal{F}}^{?1}\left(\frac{\mathcal{F}{\left(D^{(1)}\right)}^{?}\mathcal{F}\left(w_1\right)+\mathcal{F}{\left(D^{(2)}\right)}^{?}\mathcal{F}\left(w_2\right)+(\mu /\beta )\mathcal{F}(K{)}^{?}\mathcal{F}(f)}{\mathcal{F}{\left(D^{(1)}\right)}^{?}\mathcal{F}\left(D^{(1)}\right)+\mathcal{F}{\left(D^{(2)}\right)}^{?}\mathcal{F}\left(D^{(2)}\right)+(\mu /\beta )\mathcal{F}(K{)}^{?}\mathcal{F}(K)}\right)\text{\hspace{1em}(w.9)}$$
其中$?$表示复共轭，$\mathcal{F}$表示傅里叶变换，式中的乘法、除法分别表示点乘、点除。

Solver

回到本篇论文本身优化问题的求解，本文作者直接套用了Wang的半二次分裂求解算法，即直接套用了式(w.5、w.9)来求解式(5)，首先在式(5)中引入辅助变量$h_p、v_p$，分别用于近似${?}_x{S}_p、{?}_y{S}_p$，并将$h_p、v_p$与${?}_x{S}_p、{?}_y{S}_p$之间的差值作为惩罚项加入到原目标函数，得到新的目标函数：
$$\underset{S,h,v}{\min }\left\{\sum _p{\left(S_p?I_p\right)}^2+\lambda C(h,v)+\beta \left({\left({?}_x{S}_p?h_p\right)}^2+{\left({?}_y{S}_p?v_p\right)}^2\right)\right\}\text{\hspace{1em}(6)}$$
其中C(h,v)=#{p∣∣hp∣+∣vp∣≠0}C(h, v)=\#\left\{p|| h_{p}|+| v_{p} \mid \neq 0\right\}C(h,v)=#{ p∣∣hp?∣+∣vp?∣??=0}
1.分离所有出带$S$的项，得到优化求解$S$的新的目标函数：
$$\underset{S}{\min }\left\{\sum _p{\left(S_p?I_p\right)}^2+\beta \left({\left({?}_x{S}_p?h_p\right)}^2+{\left({?}_y{S}_p?v_p\right)}^2\right)\right\}\text{\hspace{1em}(7)}$$
解得：
$$S={\mathcal{F}}^{?1}\left(\frac{\mathcal{F}(I)+\beta \left(\mathcal{F}{\left({?}_x\right)}^{?}\mathcal{F}(h)+\mathcal{F}{\left({?}_y\right)}^{?}\mathcal{F}(v)\right)}{\mathcal{F}(1)+\beta \left(\mathcal{F}{\left({?}_x\right)}^{?}\mathcal{F}\left({?}_x\right)+\mathcal{F}{\left({?}_y\right)}^{?}\mathcal{F}\left({?}_y\right)\right)}\right)\text{\hspace{1em}(8)}$$
其中$\mathcal{F}(1)$为$\delta $的傅里叶变换。2.分离出所有带$h、v$的项，得到优化求解$h、v$的新的目标函数：
$$\underset{h,v}{\min }\left\{\sum _p{\left({?}_x{S}_p?h_p\right)}^2+{\left({?}_y{S}_p?v_p\right)}^2\right)+\frac{\lambda }{\beta }C(h,v)\}\text{\hspace{1em}(9)}$$
$C(h,v)$表示$h、v$非零的个数，解得：
$$\left(h_p,v_p\right)=\{\begin{array}{ll}(0,0)& {\left({?}_x{S}_p\right)}^2+{\left({?}_y{S}_p\right)}^2\le \lambda /\beta \\ \left({?}_x{S}_p,{?}_y{S}_p\right)& \text{ otherwise }\end{array}\text{\hspace{1em}(10)}$$
除了式(8、10)的迭代，还需要确定参数$\lambda 、\beta$，在实际应用过程中先给定一个较小的参数${\beta }_0=2\lambda$，然后随着迭代次数$j$的增加而逐渐增大$\beta$，${\beta }_{j+1}=\kappa {\beta }_j$，且将$\beta$的最大值限制为$1E5$，这样做的目的是加快算法的收敛速度，考虑到需要平衡算法的收敛速度和效果，一般将$\kappa$设置为2，最后，算法一般需要迭代20-30次。下面将对比本篇论文的效果与其他类似算法的效果仿真结果的对比，可以看出本篇论文提出的梯度$L_0$最小化的方法是唯一能够做到保持边缘锐利的同时去除图像中的噪声。

不同算法的降噪效果对比