Image Deconvolution with the Half-quadratic Splitting Method

在处理图像重建或者逆问题的时候，我们经常会看到一种称为 Half-quadratic Splitting（HQS）的方法，这是在优化领域里非常经典的一种方法，之前也断断续续地找了很多相关的资料，发现斯坦福大学的计算成像课里某一节 lecture note 把这个方法在图像反卷积中的使用介绍地非常详细。这篇文章也是基于这个 lecture note 的一个学习笔记。

Image Formation

一般来说，图像的成像过程可以用下面的式子表示：

$$b=c\ast x+\eta \text{\hspace{1em}(1)}$$

其中，$b$ 表示观察到的模糊图像，$c$ 表示一个卷积核，可以认为是光学系统的 PSF，$\eta$ 表示与信号无关的噪声，$x$ 是希望重建出来的清晰图像。

根据信号处理的理论，在空域的卷积相当于在频域的乘积，所以上面的式子可以写成：

b = F − 1 { F ( c ) ⋅ F ( x ) } + η (2) b = \mathcal{F}^{-1} \{ \mathcal{F}(c) \cdot \mathcal{F}(x) \} + \eta \tag{2} b=F−1{F(c)⋅F(x)}+η(2)

不过要注意的是，这个等价只有在卷积满足 circular boundary conditions 时才成立。

为了后续的推导方便，这里将卷积运算转换成矩阵运算，可以得到如下的式子：

$$b=c\ast x\iff \mathbf{b}=\mathbf{Cx}\text{\hspace{1em}(3)}$$

Inverse Filtering and Wiener Deconvolution

在忽略噪声的情况下，式 (2) 可以表示成：

$$\mathcal{F}(b)=\mathcal{F}(c)\cdot \mathcal{F}(x)$$

进而可以直接算出 $x$:

x ~ i f = F − 1 { F ( b ) F ( c ) } (4) \tilde{x}_{if} = \mathcal{F}^{-1} \{ \frac{\mathcal{F}(b)}{\mathcal{F}(c)} \} \tag{4} x~if​=F−1{F(c)F(b)​}(4)

上面的 if 就是 inverse filtering 的缩写，如果已知卷积核的形式，就可以计算出该卷积核对应的傅里叶变换，然后除以观察图像的傅里叶变换，从而得到清晰图像在频域的值，最后在做一个反傅里叶变换，得到最终的清晰图像。

这种方法看起来直观高效，不过有一个问题，就是 $\mathcal{F}(c)$ 的值很小趋近于 0 的时候，将会放大观察图像中的噪声，维纳滤波通过加入一个阻尼系数，可以避免这个问题：

$${\tilde{x}}_{wf}={\mathcal{F}}^{-1}\left\{\frac{|\mathcal{F}(c){|}^2}{|\mathcal{F}(c){|}^2+\frac{1}{SNR}}\frac{\mathcal{F}(b)}{\mathcal{F}(c)}\right\}\text{\hspace{1em}(5)}$$

SNR 表示信噪比，如果观测图像中没有噪声，那么 SNR 会很高，这种情况下，维纳滤波就变成了 inverse filter。

维纳滤波相比于直接的逆滤波，可以得到更加合理的反卷积结果。不过，这类方法面临的主要问题就是无法利用自然图像的先验信息，接下来，我们将介绍如何将先验信息与优化方法结合。

Nature Image Priors

首先，我们来看什么是自然图像先验，自然图像先验一般是一种数学模型，告诉我们在自然图像中，像素的分布更倾向于什么形态，在求解病态逆问题的时候，可能存在无数种解都符合观测值，而自然图像先验，会帮助我们挑选那些看起来更合理的解。

自然图像先验可以对图像的分布进行建模，这里，我们用正则化来表示对图像的某些性质，正则化一般对应图像先验的负对数。通常的正则化，包括平滑性，稀疏性，稀疏梯度，自相似性等。比如，稀疏性，可以用一阶范数表示为：$|\mathbf{x}|=|\sum x_i|$，而平滑性，可以用梯度的二阶范数表示：$|\Delta \mathbf{x}{|}_2$，不过，在求解图像复原的逆问题中，应用最广泛的还是称为 total variation (TV) 的正则，TV 的建立，是基于下面的观察，大部分自然图像，都可以看到很多区域都是平滑的，只有在不同物体的交界处才会出现边界，在同一个物体的内部，可以认为像素值区域平滑或者相似，而在不同物体的边界处，会有一个陡然的变化。

为了对 TV 建模，需要计算图像的一阶导数，这里，将图像水平方向的一阶导数表示为 ${\mathbf{D}}_x\mathbf{x}$，而垂直方向的一阶导数表示为 ${\mathbf{D}}_y\mathbf{y}$，具体的数学形式如下所示：

$${\mathbf{D}}_x\mathbf{x}\iff \mathbf{x}\ast d_x,d_x=\left[\begin{array}{ccc}0& 0& 0\\ 0& -1& 1\\ 0& 0& 0\end{array}\right]$$

$${\mathbf{D}}_y\mathbf{x}\iff \mathbf{x}\ast d_y,d_y=\left[\begin{array}{ccc}0& 0& 0\\ 0& -1& 0\\ 0& 1& 0\end{array}\right]$$

这个图像的一阶导数，既可以用矩阵直接计算，也可以用卷积的形式计算。

图像的 TV 正则可以表示为图像梯度的一阶范数，TV 正则项有两种表达形式，一种称为各向异性 的 TV，另外一种称为各向同性的正则，分别如下所示：

• 各向异性的 TV 正则

$$T{V}_{anisotropic}(\mathbf{x})={\Vert {\mathbf{D}}_x\mathbf{x}\Vert }_1+{\Vert {\mathbf{D}}_y\mathbf{x}\Vert }_1=\sum _{i=1}^N\left|({\mathbf{D}}_x\mathbf{x}{)}_i\right|+\left|({\mathbf{D}}_y\mathbf{x}{)}_i\right|=\sum _{i=1}^N\sqrt{({\mathbf{D}}_x\mathbf{x}{)}_i^2}+\sqrt{({\mathbf{D}}_y\mathbf{x}{)}_i^2}$$

• 各向同性的 TV 正则

$$T{V}_{isotropic}(\mathbf{x})={\Vert \mathbf{Dx}\Vert }_{2,1}=\sum _{i=1}^N\sqrt{({\mathbf{D}}_x\mathbf{x}{)}_i^2+({\mathbf{D}}_y\mathbf{x}{)}_i^2}$$

这两种正则的差异，从表达式中可以看出，对于各向同性来说，每个像素 $i$ 的梯度约束是两个方向同时约束的，而对于各向异性来说，两个方向的梯度约束是分开约束的。

Regularized Deconvolution with Half-quadratic Splitting

The Half-quadratic Splitting Method

接下来，我们介绍 HQS 这种优化方法，在介绍用这种方法求解图像复原的逆问题之前，我们先探讨这种方法的一般形式，考虑如下的一个优化问题：

$$\mathbf{b}=\mathbf{Ax}+\eta$$

其中，$\mathbf{x}\in R^N$ 是一个待求解的向量，$\mathbf{b}\in R^M$ 表示一个观测向量，$\eta$ 表示噪声，$\mathbf{A}\in R^{M\times N}$ 表示转换矩阵，这个问题的一般求解形式如上所示：

$$\underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{Ax}-\mathbf{b}\Vert }_2^2+\lambda \Psi (\mathbf{x})$$

其中，前面的第一项一般称为数据项，第二项称为正则项，直接求解上面的式子，有的时候并不能很好地收敛，一种更为鲁棒的求解方式，应该是将上面的优化函数，改写成下面这种形式：

$$\begin{gathered} \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{Ax}-\mathbf{b}\Vert }_2^2+\lambda \Psi (\mathbf{z}) \\ \text{subject to}\mathbf{Dx}-\mathbf{z}=0 \end{gathered}$$

这个优化函数中，引入了一个中间变量 $\mathbf{z}$，这个额外的中间变量，可以将上面的优化函数拆成两部分，一部分是数据项，另外一部分是正则项，这两项依赖的变量是相互独立的，$\mathbf{x},\mathbf{z}$ 之间靠一个约束表达式联系，如果将这个约束项合入优化函数，则整个优化函数可以写成：

$$L_p(\mathbf{x},\mathbf{z})=f(\mathbf{x})+g(\mathbf{z})+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2$$

优化函数写成这种形式，可以通过相互迭代的方式求解，求解 $\mathbf{x}$ 与 求解 $\mathbf{z}$ 可以分开进行：

$$\mathbf{x}\leftarrow {\mathbf{prox}}_{f,p}(\mathbf{z})=\underset{\mathbf{x}}{\mathrm{argmin}}{L}_p(\mathbf{x},\mathbf{z})=\underset{\mathbf{x}}{\mathrm{argmin}}f(\mathbf{x})+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2\text{\hspace{1em}(13)}$$

$$\mathbf{z}\leftarrow {\mathbf{prox}}_{f,p}(\mathbf{Dx})=\underset{\mathbf{z}}{\mathrm{argmin}}{L}_p(\mathbf{x},\mathbf{z})=\underset{\mathbf{z}}{\mathrm{argmin}}g(\mathbf{z})+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2\text{\hspace{1em}(14)}$$

从上面的求解过程可以看出，当我们更新 $\mathbf{x}$ 的时候，只需要考虑 $f(\mathbf{x})$，而当我们更新 $\mathbf{z}$ 的时候，只需要考虑 $g(\mathbf{z})$，而不需要同时考虑这两个函数。这个性质，可以构建一个非常灵活的框架，让我们可以灵活地与各种正则函数相结合。这种方式也被称为 plug-and-play （即插即用）。

虽然 HQS 可以用于解决各种逆问题，不过我们这里还是讨论比较特殊的一种图像解卷积问题，我们讨论一种已知固定卷积核的情况，这样对应的矩阵是一个循环 Toeplitz 矩阵。先定义如下的关系：

c ∗ x = F − 1 { F ( c ) ⋅ F ( x ) } = C x F − 1 { F ( c ) ∗ ⋅ F ( x ) } = C T x F − 1 { F ( b ) F ( c ) } = C − 1 x c * x = \mathcal{F}^{-1} \{ \mathcal{F}(c) \cdot \mathcal{F}(x) \} = \mathbf{C}\mathbf{x} \\ \mathcal{F}^{-1} \{ \mathcal{F}(c)^{*} \cdot \mathcal{F}(x) \} = \mathbf{C}^{T}\mathbf{x} \\ \mathcal{F}^{-1} \{ \frac{\mathcal{F}(b)}{\mathcal{F}(c)} \} = \mathbf{C}^{-1}\mathbf{x} c∗x=F−1{F(c)⋅F(x)}=CxF−1{F(c)∗⋅F(x)}=CTxF−1{F(c)F(b)​}=C−1x

Standard Form of HQS with TV and Denoising Regularizers

接下来，我们考虑基于 TV 正则的 HQS 的优化方法，由上面的表达式，我们可以将带 TV 正则的优化函数写成如下形式：

$$\begin{gathered} \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{Cx}-\mathbf{b}\Vert }_2^2+\lambda \Psi (\mathbf{z}) \\ \text{subject to}\mathbf{Dx}-\mathbf{z}=0 \end{gathered}$$

其中，$\mathbf{D}=[{\mathbf{D}}_x^T,{\mathbf{D}}_y^T]\in R^{2N\times N}$ 表示 $x,y$ 方向的一阶导数，这里的 $\mathbf{z}\in R^{2N}$ 比 $\mathbf{x}\in R^N$ 要大一倍，因为每个像素，有 $x,y$ 两个方向的梯度。

对于更为一般的情况，我们可以使用一个简单的正则项，将待求解的图像投影到一个灵活的自然图像空间中，整个的 HQS 的形式可以写成如下所示：

$$\begin{gathered} \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{Cx}-\mathbf{b}\Vert }_2^2+\lambda \Psi (\mathbf{z}) \\ \text{subject to}\mathbf{x}-\mathbf{z}=0 \end{gathered}$$

Efficient Implementation of the x-Update using Inverse Filtering

前面介绍过，HQS 的方法，会交替迭代更新 $\mathbf{x},\mathbf{z}$，我们先来看 $\mathbf{x}$ 的更新，

$${\mathbf{prox}}_{f,p}(\mathbf{z})=\underset{\mathbf{x}}{\mathrm{argmin}}f(\mathbf{x})+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2=\underset{\mathbf{x}}{\mathrm{argmin}}\frac{1}{2}{\Vert \mathbf{Cx}-\mathbf{b}\Vert }_2^2+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2\text{\hspace{1em}(20)}$$

将上面的表达式展开，可以得到：

$$\begin{gathered} \frac{1}{2}{\Vert \mathbf{Cx}-\mathbf{b}\Vert }_2^2+\frac{p}{2}{\Vert \mathbf{Dx}-\mathbf{z}\Vert }_2^2 \\ =\frac{1}{2}(\mathbf{Cx}-\mathbf{b}{)}^T(\mathbf{Cx}-\mathbf{b})+\frac{p}{2}(\mathbf{Dx}-\mathbf{z})(\mathbf{Dx}-\mathbf{z}{)}^T \\ =\frac{1}{2}({\mathbf{x}}^T{\mathbf{C}}^T\mathbf{Cx}-2{\mathbf{x}}^T{\mathbf{C}}^T\mathbf{b}+{\mathbf{b}}^T\mathbf{b})+\frac{p}{2}({\mathbf{x}}^T{\mathbf{D}}^T\mathbf{Dx}-2{\mathbf{x}}^T{\mathbf{D}}^T\mathbf{z}+{\mathbf{z}}^T\mathbf{z}) \end{gathered}$$

将上面的表达式对 $\mathbf{x}$ 求导，可以得到：

$${\mathbf{C}}^T\mathbf{Cx}-{\mathbf{C}}^T\mathbf{b}+p{\mathbf{D}}^T\mathbf{Dx}-p{\mathbf{D}}^T\mathbf{z}$$

让导数为 0 ，进而可以求得 $\mathbf{x}$：

$$({\mathbf{C}}^T\mathbf{C}+p{\mathbf{D}}^T\mathbf{D}{)}^{-1}({\mathbf{C}}^T\mathbf{b}+p{\mathbf{D}}^T\mathbf{z})\text{\hspace{1em}(24)}$$

对于满足 circular boundary conditions 的 2D 图像解卷积问题，上面的式子，可以变换到傅里叶域，然后再进行求解，上面所有的矩阵相乘的形式，都可以找到对应的傅里叶域的形式。

Special Case of TV Regularizer

如果正则项是 TV 项，那么 $\mathbf{D}$ 就是有限差分算子，上面的公式 (24) 可以写成如下形式：

( C T C + p D T D ) ⇔ F − 1 { F { c } ∗ ⋅ F { c } + p ( F { d x } ∗ ⋅ F { d x } + F { d y } ∗ ⋅ F { d y } ) } ( \mathbf{C}^{T}\mathbf{C} + p\mathbf{D}^{T}\mathbf{D} ) \Leftrightarrow \mathcal{F}^{-1} \{ \mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{c\} + p (\mathcal{F}\{d_x\}^{*} \cdot \mathcal{F}\{d_x\} + \mathcal{F}\{d_y\}^{*} \cdot \mathcal{F}\{d_y\}) \} (CTC+pDTD)⇔F−1{F{c}∗⋅F{c}+p(F{dx​}∗⋅F{dx​}+F{dy​}∗⋅F{dy​})}

( C T b + p D T z ) ⇔ F − 1 { F { c } ∗ ⋅ F { b } + p ( F { d x } ∗ ⋅ F { z 1 } + F { d y } ∗ ⋅ F { z 2 } ) } (\mathbf{C}^{T} \mathbf{b} + p \mathbf{D}^{T} \mathbf{z} ) \Leftrightarrow \mathcal{F}^{-1} \{ \mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{b\} + p (\mathcal{F}\{d_x\}^{*} \cdot \mathcal{F}\{z_1\} + \mathcal{F}\{d_y\}^{*} \cdot \mathcal{F}\{z_2\}) \} (CTb+pDTz)⇔F−1{F{c}∗⋅F{b}+p(F{dx​}∗⋅F{z1​}+F{dy​}∗⋅F{z2​})}

由此，可以得到公式(20) 的解为：

p r o x f , p ( z ) = F − 1 ( F { c } ∗ ⋅ F { b } + p ( F { d x } ∗ ⋅ F { z 1 } + F { d y } ∗ ⋅ F { z 2 } ) F { c } ∗ ⋅ F { c } + p ( F { d x } ∗ ⋅ F { d x } + F { d y } ∗ ⋅ F { d y } ) ) \mathbf{prox}_{f, p} (\mathbf{z}) = \mathcal{F}^{-1} \left( \frac{\mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{b\} + p (\mathcal{F}\{d_x\}^{*} \cdot \mathcal{F}\{z_1\} + \mathcal{F}\{d_y\}^{*} \cdot \mathcal{F}\{z_2\})}{\mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{c\} + p (\mathcal{F}\{d_x\}^{*} \cdot \mathcal{F}\{d_x\} + \mathcal{F}\{d_y\}^{*} \cdot \mathcal{F}\{d_y\})} \right) proxf,p​(z)=F−1(F{c}∗⋅F{c}+p(F{dx​}∗⋅F{dx​}+F{dy​}∗⋅F{dy​})F{c}∗⋅F{b}+p(F{dx​}∗⋅F{z1​}+F{dy​}∗⋅F{z2​})​)

Special Case of Denoising Reg

对于更为一般的正则项，$\mathbf{D}$ 可以认为是一个单位矩阵，公式 (20) 的求解将变得更为简单：

p r o x f , p ( z ) = F − 1 ( F { c } ∗ ⋅ F { b } + p F { z } F { c } ∗ ⋅ F { c } + p ) \mathbf{prox}_{f, p} (\mathbf{z}) = \mathcal{F}^{-1} \left( \frac{\mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{b\} + p \mathcal{F}\{z\} } {\mathcal{F}\{c\}^{*} \cdot \mathcal{F}\{c\} + p } \right) proxf,p​(z)=F−1(F{c}∗⋅F{c}+pF{c}∗⋅F{b}+pF{z}​)

Updating z with the TV Regularizer

公式(14) 关于 $\mathbf{z}$ 的更新，可以表示成如下的形式：

$${\mathbf{prox}}_{g,p}(\mathbf{v})={\mathcal{S}}_{\lambda /p}(\mathbf{v})=\underset{\mathbf{z}}{\mathrm{argmin}}\lambda {\left|\mathbf{z}\right|}_1+\frac{p}{2}{\Vert \mathbf{v}-\mathbf{z}\Vert }_2^2$$

其中，$\mathbf{v}=\mathbf{Dx}$，$\mathcal{S}$ 是一个分段函数：

$${\mathcal{S}}_k(v)=\{\begin{array}{cc}v-k& v>k\\ 0& |v|<k\\ v+k& v<-k\end{array}$$

Isotropic TV Norm

对于各向同性的 TV 正则，正则项可以表示为：

$$g(\mathbf{z})=\lambda {\Vert \mathbf{Dx}\Vert }_{2,1}=\lambda \sum _{i=1}^N\sqrt{({\mathbf{D}}_x\mathbf{x}{)}_i^2+({\mathbf{D}}_y\mathbf{x}{)}_i^2}$$

那么，带各向同性正则项的解卷积问题可以表示成：

$$\begin{gathered} \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{Cx}-\mathbf{b}\Vert }_2^2+\lambda \sum _{i=1}^N{\Vert \left[\begin{array}{c}{z}_i\\ z_{i+N}\end{array}\right]\Vert }_2 \\ \text{subject to}\mathbf{Dx}-\mathbf{z}=0 \end{gathered}$$

$z_i$ 表示 $\mathbf{z}$ 的第 $i$ 个元素，其中，$1\le i\le N$ 表示水平方向的有限差分，$N+1\le i\le 2N$ 表示垂直方向的有限差分。

$\mathbf{z}$ 的更新可以表示成：

$$\mathbf{z}\leftarrow {\mathbf{prox}}_{f,p}(\mathbf{v})=\underset{\mathbf{z}}{\mathrm{argmin}}\lambda \sum _{i=1}^N{\Vert \left[\begin{array}{c}{z}_i\\ z_{i+N}\end{array}\right]\Vert }_2+\frac{p}{2}{\Vert \mathbf{v}-\mathbf{z}\Vert }_2^2\mathbf{v}=\mathbf{Dx}$$

最终 $\mathbf{z}$ 的更新可以表示为：

$$\left[\begin{array}{c}{z}_i\\ z_{i+N}\end{array}\right]\leftarrow {\mathcal{S}}_{\lambda /p}\left(\left[\begin{array}{c}{v}_i\\ v_{i+N}\end{array}\right]\right),1\le i\le N$$

Updating z with DnCNN or any Gaussian Denoiser as the Regularizer

如果我们进一步审视公式 (14) ，不考虑矩阵 $\mathbf{D}$ 的情况下，

$$\begin{gathered} \underset{\mathbf{z}}{\mathrm{argmin}}g(\mathbf{z})+\frac{p}{2}{\Vert \mathbf{x}-\mathbf{z}\Vert }_2^2 \\ =\underset{\mathbf{z}}{\mathrm{argmin}}\lambda \Psi (\mathbf{z})+\frac{p}{2}{\Vert \mathbf{x}-\mathbf{z}\Vert }_2^2 \\ =\underset{\mathbf{z}}{\mathrm{argmin}}\Psi (\mathbf{z})+\frac{p}{2\lambda }{\Vert \mathbf{x}-\mathbf{z}\Vert }_2^2\text{\hspace{1em}(36)} \end{gathered}$$

公式 (36) 可以看成是一个降噪问题，可以等价成如下的表达式：

$$\underset{\mathbf{x}}{\mathrm{argmin}}\Psi (\mathbf{x})+\frac{1}{2{\sigma }^2}{\Vert \mathbf{v}-\mathbf{x}\Vert }_2^2$$

其中，$\mathbf{v}$ 可以看成是观测到的有噪图像，$\mathbf{x}$ 表示我们需要恢复的无噪图像，基于这个观测，对于降噪的正则项，那么对 $\mathbf{z}$ 的更新可以简单地变成一个降噪过程，这个降噪可以用传统的降噪方法，也可以用基于CNN 的降噪方法，

$${\mathbf{prox}}_{\mathcal{D},p}(\mathbf{x})=\mathcal{D}\left(\mathbf{x},{\sigma }^2=\frac{\lambda }{p}\right)$$

最后，做一个总结，基于 TV 正则和基于降噪正则的 HQS 方法的主要流程可以归纳如下：