Solving the Discrete Poisson Equation using Conjugate Gradient Method

在 Tone mapping: A quick survey. (Tone mapping 原理) 中我們介紹了 gradient domain 的解法，本篇暨 FFT(DST) 解法後，提供一個 spatial domain 的 optimization 方式。

與原始論文使用 Multigrid 不完全相同，但Conjugate gradient 的作法可被 Multigrid 利用來取代 Gauss–Seidel method 作為 smooth 的方法。與 Multigrid 相比直接求解可能需要更多次的 iteration (大約2000次)。

等號右邊 \(\operatorname{div}{G}\) 已知，求一個 \(I\) 做 discrete Laplace operator 後滿足該式
\[
\begin{align*}
\nabla^2{I}&=\nabla_{x} \cdot \Phi\cdot\nabla_{x} f +\nabla_{y} \cdot \Phi\cdot\nabla_{y} f \\
&=\operatorname{div}{G}
\end{align*}
\]

這邊我們將上式改寫為 \(Ax=b\) 的形式

\[\nabla^2{x}=Ax\]

\[\operatorname{div}{G}=b\]

我們不需要真的建立一個巨大的 sparse matrix \(A\) (雖然看到真的有人這麼做)，事情可以更簡單易懂，而 \(Ax\) 就相當於將 \(x\) 應用一個 3×3 的 Laplace operator ，求 \(x\) 與 \(L\) 的 convolution，當然 \(L\) 也就只是一個普通的 linear kernel 罷了，如同 image blur 一樣

Note: \(x\)是目標影像 reshape 的 1x\((M\)x\(N)\) 的向量，\(A\)則是一個\((M\)x\(N)\)x\((M\)x\(N)\)的矩陣，其中每個 row 中僅有5個非零項。這邊\(M\)和\(N\)分別表影像寬跟高。

\[Ax=x\otimes{L}\]

\[L=\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}\]

注意 \(L\) 的正負號與實作時一次微分的方向有關(左減右或右減左)，不妨都試試看結果；應會有正片或負片效果XD

由於這邊的 convolution 會有 1 pixel 的 boundary ，這邊注意邊界條件 ( boundary condition ) 很多種都可以試試，常見的幾種如下 ( 絕對超級簡單的解釋! 詳見opencv copyMakeBorder)

1. Dirichlet boundary condition: zero padding
2. Neumann boundary condition: repeat border pixel

當然 Dirichlet boundary 比較簡單，但是數值容易在邊緣衰減很快； Neumann boundary 則是較能維持邊緣數值。

這邊由於 sparse matrix 乘法已經被轉換成一個 linear convolution ，所以很 \(Ax\) 容易理解，直接上Wiki的 conjugate gradient 就可以解出來。

\[
\begin{align}
& \mathbf{r}_0 := \mathbf{b} – \mathbf{A x}_0 \\
& \mathbf{p}_0 := \mathbf{r}_0 \\
& k := 0 \\
& \hbox{repeat} \\
& \qquad \alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k} \\
& \qquad \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \\
& \qquad \mathbf{r}_{k+1} := \mathbf{r}_k – \alpha_k \mathbf{A p}_k \\
& \qquad \hbox{if } r_{k+1} \hbox{ is sufficiently small then exit loop} \\
& \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k} \\
& \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\
& \qquad k := k + 1 \\
& \hbox{end repeat} \\
& \hbox{The result is } \mathbf{x}_{k+1}
\end{align}
\]

Matlab Code，記得把 \(Ax\) 換成 imfilter