General (Non-orthogonal) basis projection


x^{'}=\left(M^{T}M\right)^{-1}M^{T}x


M=
\left[\begin{matrix}
\vdots & \vdots & & \vdots \\
m_1 & m_2 & \cdots & m_n \\
\vdots & \vdots & & \vdots
\end{matrix}\right]

where m can be the non-orthogonal basis column vectors, x is a column vector in global coordinates, x^{'} is a column vector in non-orthogonal basis coordinates. Note the \left(M^{T}M\right)^{-1}M^{T} is the psudo inverse of M when M is not a square matrix.

derivation:


\begin{align*}
x &= Mx^{'} \\
M^{T}x &= M^{T}Mx^{'} \\
x^{'}&=\left(M^{T}M\right)^{-1}M^{T}x
\end{align*}

if M is an orthogonal basis, then


\begin{align*}
x^{'}&=\left(M^{T}M\right)^{-1}M^{T}x\\
&=I^{-1}M^{T}x\\
&=M^{T}x\\
&=M^{-1}x
\end{align*}

re-project:


\begin{align*}
Mx^{'} &=M\left(M^{T}M\right)^{-1}M^{T}x \\
&=MM^{-1}M^{-T}M^{T}x \\
&=x
\end{align*}

Two Point Interpolation Equations

1. Linear Interpolation

最簡單的interpolation,也稱做 linear alpha blending,後面的式子基本上都是這個變型而來,{\alpha}=[0,1]

x=\alpha{x_0}+\left(1-\alpha\right){x_1}

簡化版本,雖然不直觀,但計算上可以少一個乘法,適合程式最佳化使用

x=\alpha\left(x_0-x_1\right)+x_1

2. Cosine Interpolation

\alpha=0.5\left(1-cos(\pi\alpha)\right)

3. Smooth step Interpolation

可以用來近似cosine,比較好計算

\alpha={\alpha}^2(3-2\alpha)

Reference

  1. http://codeplea.com/simple-interpolation
  2. http://sol.gfxile.net/interpolation/

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"

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"