Image Deconvolution with the Half-quadratic Splitting Method

在處理圖像重建或者逆問題的時候，我們經常會看到一種稱為 Half-quadratic Splitting（HQS）的方法，這是在優化領域里非常經典的一種方法，之前也斷斷續續地找了很多相關的資料，發現斯坦福大學的計算成像課里某一節 lecture note 把這個方法在圖像反卷積中的使用介紹地非常詳細。這篇文章也是基于這個 lecture note 的一個學習筆記。

Image Formation

一般來說，圖像的成像過程可以用下面的式子表示：

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

其中，$b$ 表示觀察到的模糊圖像，$c$ 表示一個卷積核，可以認為是光學系統的 PSF，$\eta$ 表示與信號無關的噪聲，$x$ 是希望重建出來的清晰圖像。

根據信號處理的理論，在空域的卷積相當于在頻域的乘積，所以上面的式子可以寫成：

$$b={\mathcal{F}}^{?1}\{\mathcal{F}(c)?\mathcal{F}(x)\}+\eta \text{\hspace{1em}(2)}$$

不過要注意的是，這個等價只有在卷積滿足 circular boundary conditions 時才成立。

為了后續的推導方便，這里將卷積運算轉換成矩陣運算，可以得到如下的式子：

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

Inverse Filtering and Wiener Deconvolution

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

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

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

$${\tilde{x}}_{if}={\mathcal{F}}^{?1}\{\frac{\mathcal{F}(b)}{\mathcal{F}(c)}\}\text{\hspace{1em}(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}$，具體的數學形式如下所示：

D x x ? x ? d x , d x = [ 0 0 0 0 ? 1 1 0 0 0 ] \mathbf{D}_x\mathbf{x} \Leftrightarrow \mathbf{x} \ast d_x, \quad d_x = \begin{bmatrix} 0 & 0 & 0\\ 0 & -1 & 1 \\ 0 & 0 & 0 \end{bmatrix} Dx?x?x?dx?,dx?= ?000?0?10?010? ?

D y x ? x ? d y , d y = [ 0 0 0 0 ? 1 0 0 1 0 ] \mathbf{D}_y\mathbf{x} \Leftrightarrow \mathbf{x} \ast d_y, \quad d_y = \begin{bmatrix} 0 & 0 & 0\\ 0 & -1 & 0 \\ 0 & 1 & 0 \end{bmatrix} Dy?x?x?dy?,dy?= ?000?0?11?000? ?

這個圖像的一階導數，既可以用矩陣直接計算，也可以用卷積的形式計算。

圖像的 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 矩陣。先定義如下的關系：

$$\begin{gathered} c?x={\mathcal{F}}^{?1}\{\mathcal{F}(c)?\mathcal{F}(x)\}=\mathbf{Cx} \\ {\mathcal{F}}^{?1}\{\mathcal{F}(c{)}^{?}?\mathcal{F}(x)\}={\mathbf{C}}^T\mathbf{x} \\ {\mathcal{F}}^{?1}\{\frac{\mathcal{F}(b)}{\mathcal{F}(c)}\}={\mathbf{C}}^{?1}\mathbf{x} \end{gathered}$$

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) 可以寫成如下形式：

$$({\mathbf{C}}^T\mathbf{C}+p{\mathbf{D}}^T\mathbf{D})?{\mathcal{F}}^{?1}\{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{c\}+p(\mathcal{F}\{d_x{\}}^{?}?\mathcal{F}\{d_x\}+\mathcal{F}\{d_y{\}}^{?}?\mathcal{F}\{d_y\})\}$$

$$({\mathbf{C}}^T\mathbf{b}+p{\mathbf{D}}^T\mathbf{z})?{\mathcal{F}}^{?1}\{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{b\}+p(\mathcal{F}\{d_x{\}}^{?}?\mathcal{F}\{z_1\}+\mathcal{F}\{d_y{\}}^{?}?\mathcal{F}\{z_2\})\}$$

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

$${\mathbf{prox}}_{f,p}(\mathbf{z})={\mathcal{F}}^{?1}\left(\frac{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{b\}+p(\mathcal{F}\{d_x{\}}^{?}?\mathcal{F}\{z_1\}+\mathcal{F}\{d_y{\}}^{?}?\mathcal{F}\{z_2\})}{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{c\}+p(\mathcal{F}\{d_x{\}}^{?}?\mathcal{F}\{d_x\}+\mathcal{F}\{d_y{\}}^{?}?\mathcal{F}\{d_y\})}\right)$$

Special Case of Denoising Reg

對于更為一般的正則項，$\mathbf{D}$ 可以認為是一個單位矩陣，公式 (20) 的求解將變得更為簡單：

$${\mathbf{prox}}_{f,p}(\mathbf{z})={\mathcal{F}}^{?1}\left(\frac{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{b\}+p\mathcal{F}\{z\}}{\mathcal{F}\{c{\}}^{?}?\mathcal{F}\{c\}+p}\right)$$

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}$$

那么，帶各向同性正則項的解卷積問題可以表示成：

min ? x 1 2 ∥ C x ? b ∥ 2 2 + λ ∑ i = 1 N ∥ [ z i z i + N ] ∥ 2 subject?to D x ? z = 0 \min_{\mathbf{x}} \frac{1}{2} \left \| \mathbf{C}\mathbf{x} - \mathbf{b} \right \|_2^{2} + \lambda \sum_{i=1}^{N} \left \| \begin{bmatrix} z_i\\ z_{i+N} \end{bmatrix} \right \|_{2} \\ \text{subject to} \quad \mathbf{D}\mathbf{x} - \mathbf{z} = 0 xmin?21?∥Cx?b∥22?+λi=1∑N? ?[zi?zi+N??] ?2?subject?toDx?z=0

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

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

z ← p r o x f , p ( v ) = arg?min ? z λ ∑ i = 1 N ∥ [ z i z i + N ] ∥ 2 + p 2 ∥ v ? z ∥ 2 2 v = D x \mathbf{z} \gets \mathbf{prox}_{f, p} (\mathbf{v}) = \argmin_{\mathbf{z}} \lambda \sum_{i=1}^{N} \left \| \begin{bmatrix} z_i\\ z_{i+N} \end{bmatrix} \right \|_{2} + \frac{p}{2} \left \| \mathbf{v} - \mathbf{z} \right \|_{2}^{2} \quad \mathbf{v} = \mathbf{D}\mathbf{x} z←proxf,p?(v)=zargmin?λi=1∑N? ?[zi?zi+N??] ?2?+2p?∥v?z∥22?v=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 方法的主要流程可以歸納如下：

總結

以上是生活随笔為你收集整理的Image Deconvolution with the Half-quadratic Splitting Method的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。