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次)。

Continue reading "Solving the Discrete Poisson Equation using Conjugate Gradient Method"

Gradient domain high dynamic range compression with local Poisson solver

flowsolver本篇參考論文 A real-time implementation of gradient domain high dynamic range compression using a local Poisson solver

 

1. Model

前情提要,在前一篇中我們討論到了gradient domain的Poisson reconstruction,給定已知的divergence \operatorname{div}{G},嘗試重建I


\begin{align}
\nabla^2{I}&=\nabla_{x} \cdot \Phi\cdot\nabla_{x} f +\nabla_{y} \cdot \Phi\cdot\nabla_{y} f \\
&=\operatorname{div}{G} \label{eq:model}
\end{align}

2. Integral Problem

基本上這可看為一個積分問題,用一維的gradient思考,\nabla{f}套用衰減函數\Phi後由任何一邊開始積分就可以返回I,如下式可以積回任何第n點的I_n,其中C為一個常量

I_n=\sum_{i=0}^{n}{\nabla{f}_i{\Phi_i}}+C

一個直白的說法就是給你數個位置的個坡度,然後以平滑的方式重建回斜坡

但在二為空間就不是如此,從任何一邊積分的結果都會有所不同,所以在無法兼顧的情況下,我們僅能將其轉換為一個optimization問題來看,使其符合目標式\eqref{eq:model}從任何單點來看都需要整個畫面的divergence才能完整的重建

 

Continue reading "Gradient domain high dynamic range compression with local Poisson solver"

Camera response curve calibration

公式參考引用數驚人的[1]"Paul E. Debevec and Jitendra Malik. Recovering High Dynamic Range Radiance Maps from Photographs. In SIGGRAPH 97, August 1997."

1. Introduction

在影像處理中,有幾個主題如HDR, Lens shading correction, shape from shading等,都需精準的由intensity判斷一些重要訊息,然而在真實的相機中,input irradiance E和輸出的intensity Z未必是與曝光時間\Delta{t}呈線性關係,而是一個類似gamma curve的非線性的對應曲線。這會造成計算上的誤差,甚至是許多artifact的產生,因此校正這個關係使得整張影像的intensity Z與input irradiance E呈線性變成一個研究主題。

  • Remark: 這部分是我自己的碎碎念

Continue reading "Camera response curve calibration"

From implicit filtering to explicit filtering

在影像處理中,我們常常對影像進行許多filtering,如跑下面這個3x3的box blur kernel W

W=\begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix}

每個output影像pixel O_i,可視為在kernel內的input影像pixel I_j的線性加總

O_{i}=\sum_{j{\in}W}{w_{j}I_{j}}

我們令這個動作被應用在全影像成一個線性系統,這邊把OI都展開成向量

A\vec{O}=\vec{I}

\vec{O}=A^{-1}\vec{I}

這邊的A^{-1}作用等同於對影像做一個W kernel的filtering,在這裡求解\vec{O}=A^{-1}\vec{I}系統最佳化我們稱之為implicit filtering,其實等同於做一個explicit的box filter,而後者明顯比較容易實現! 也更適合硬體實作。

不幸的是,我們定義問題時通常以一個quadratic form的線性系統Ax=b來做最佳化;令這些式子難以被實現成硬體加速(通常是memory trade-off或iteration問題)。於是將implicit filtering轉換成explicit filtering變成一個特別的研究方向。比方說L-1 regularization的最佳化問題可以被轉換成一個iterative median filtering [1]。

[1]Y. Li and S. Osher, “A New Median Formula with Applications to PDE Based Denoising,” Comm. Math. Sciences, vol. 7, pp. 741-753, 2009.

Jacobi method

講到iterative method,得從最簡單的方法開始─Jacobi method

假設有一方程組

Ax=b

這邊我們舉例

\begin{bmatrix}{A}_{11}&{A}_{12}&{A}_{13}\\{A}_{21}&{A}_{22}&{A}_{23}\\{A}_{31}&{A}_{32}&{A}_{33}\end{bmatrix}\begin{bmatrix}{x}_{1}\\{x}_{2}\\{x}_{3}\end{bmatrix}=\begin{bmatrix}{b}_{1}\\{b}_{2}\\{b}_{3}\end{bmatrix}

Continue reading "Jacobi method"

Linear quadratic form in optimization problem Ax=b

談到最佳化(optimization)就不能逃避這個問題,為什麼大多的論文最佳化都喜歡來個

 \underset{x}{min}{\left \|Ax-b\right \|^2}

然後解

Ax=b

因為最好微分也最好解方程式! 但便宜的代價就是結果未必符合期望...解出來的x會有平均化residual的特性,在大多某些實際問題上不適用,以機率角度來說就是Gaussian prior,通常在實驗結果是個結果悽慘的對照組,但仍是一個非常適合入門的最佳化式子,因為我們需要一個夠爛的基準解,只要比它稍好就談得上學術貢獻拉!

本篇提供一個簡易的推導,故事是這樣來的,先從最簡單的二次式開始吧

Continue reading "Linear quadratic form in optimization problem Ax=b"