全变差图像复原问题(Total Variation Image Restoration Problem)

本文主要讨论图像处理中的全变差图像复原问题，首先简要介绍一下数字图像复原问题：
令 x \in \mathcal{R}^n表示一个二维图像，令n=n_1 \cdot n_2表示该图像的所有像素，其中n_1 和n_2分别表示水平和竖直方向像素的数量。图像复原问题是从退化了的图像(记为x^0)复原出原图(记为\overline{x})，数学模型如下：
x^0 = K\overline{x}+\omega
其中，\omega \in \mathcal{R}^m表示加性噪声，K: \mathcal{R} \rightarrow \mathcal{R}^m表示一个线性变换(eg: 模糊核的卷积)。由于该模型通常是病态的，所以需要一定的正则化技术。其中最受欢迎的正则化技术是全变差正则，主要是因为它能够保持图像的边缘。令\partial_1 : \mathcal{R}^n \rightarrow \mathcal{R}^n和\partial_2 : \mathcal{R}^n \rightarrow \mathcal{R}^n分别表示水平和竖直方向上的有限差分算子。令\nabla := (\partial_1, \partial_2)表示梯度算子。考虑以下TV-l^2模型：
\min \||\nabla{x}|\|_1+\frac{\mu}{2}\|Kx-x^0\|^2
其中\mu>0是全变差正则化项和数据保真项之间的权衡常数，并且l^1-norm \||\cdot|\|_1定义于\mathcal{R}^n \times \mathcal{R}^n上
\||y|\|_1:=\|(|y|)\|_1, \quad \forall y = (y_1, y_2) \in \mathcal{R}^n \times \mathcal{R}^n
这里，|y| := sqrt(y_1^2+y_2^2) \in \mathcal{R}^n .
通过引入附加变量y \in \mathcal{R}^n \times \mathcal{R}^n，可重新把TV-l^2模型写为
\min \||y|\|_1+\frac{\mu}{2}\|Kx-x^0\|^2 \\ s.t. \quad \nabla x = y
这是模型
\min \{F(x)+G(y) | Ax+By=b, x \in \mathcal{X}, y \in \mathcal{Y}\}
中取F(x) = \frac{\mu}{2}\|Kx-x^0\|^2，G(y) = \||y|\|_1，b=0，A = \nabla 和 B = -I。该问题的拉格朗日增广函数为：
L_\beta(x, y, \lambda) = \||y|\|_1+\frac{\mu}{2}\|Kx-x^0\|^2-\lambda^T(\nabla x-y)+\frac{\beta}{2}\|\nabla x-y\|^2
其中\lambda为对偶变量，\beta表示惩罚参数。
则ADMM迭代如下：
\left \{ \begin{array}{l} x^{k+1} = \arg \min L_\beta(x, y^k, \lambda^k) \\ y^{k+1} = \arg \min L_\beta(x^{k+1}, y, \lambda^k) \\ \lambda^{k+1} = \lambda^k - \beta(\nabla x^{k+1} - y^{k+1}) \end{array} \right.
即
\left \{ \begin{array}{l} x^{k+1} = \arg \min \{\frac{\mu}{2}\|Kx-x^0\|^2-(\lambda^k)^T(\nabla x-y^k)+\frac{\beta}{2}\|\nabla x -y^k \|^2 \} \\ y^{k+1} = \arg \min \{\||y|\|_1 - (\lambda^k)^T(\nabla x^{k+1}-y) + \frac{\beta}{2}\|\nabla x^{k+1} - y\|^2\} \\ \lambda^{k+1} = \lambda^k-\beta(\nabla x^{k+1} - y^{k+1}) \end{array} \right.

• 其中x最小化子问题可解释为如下：
  \mu K^T(Kx^{k+1}-x^0) - \nabla^T(\lambda^k)+\beta \nabla^T(\nabla x^{k+1} - y^k) = 0
  即
  (\beta \nabla^T \nabla + \mu K^T K)x^{k+1} = \nabla^T(\beta y^k + \lambda^k)+\mu K^T x^0
  特别地，当K是空间不变卷积算子并且使用循环边界条件，算子K和\nabla可以通过离散傅里叶变换(DFT)对角化
  K = \mathcal{F}^{-1}D_K\mathcal{F} \quad 和 \quad \nabla = \mathcal{F}^{-1}D\mathcal{F}
  其中D_K和D是具有正对角项的对角矩阵，\mathcal{F}是傅里叶算子，因此，变量x^{k+1}子问题通过傅里叶变换和逆傅里叶变换容易求解。

• 变量y^{k+1}的更新如下：
  \begin{array}{l} y^{k+1} = \arg \min \{\||y|\|_1 - (\lambda^k)^T(\nabla x^{k+1}-y) + \frac{\beta}{2}\|\nabla x^{k+1} - y\|^2\} \\ y^{k+1} = shrink_{1/\beta}(\nabla x^{k+1}-\frac{\lambda^{k+1}}{\beta}) \end{array}
  定义如下算子
  shrink_{1/\beta}(t) = t-\min\{\frac{1}{\beta}, \|t\|\}\frac{t}{\|t\|}

• 变量\lambda^{k+1}的更新：
  \lambda^{k+1} = \lambda^k-\beta(\nabla x^{k+1} - y^{k+1})