DTAM論文を読んでみた - wacchoz’s note

どうせ暇だしDTAM論文を読んでみます。
結構むずかしいです・・・
画像処理は専門ではなくて趣味で読んでるだけなので、適当なことを書いてるかもしれないがそのつもりで。

Richard A. Newcombe, Steven J. Lovegrove and Andrew J. Davison.
DTAM: Dense Tracking and Mapping in Real-Time, ICCV 2011
https://www.robots.ox.ac.uk/~vgg/rg/papers/newcombe_davison__2011__dtam.pdf

著者によるデモ動画はこちら。
www.youtube.com

画面上のすべてのピクセルについて深度を計算する、しかもリアルタイムに、というのが凄いです。

f:id:wacchoz:20181106231101p:plain

フレームrでの画像のピクセルの輝度 $I_r(u)$ は、フレーム $m$ では、 $\displaystyle \boldsymbol{I}_m(\pi(KT_{mr}\pi^{-1}(\boldsymbol{u},d)))$ となります。
ここで $K$ はカメラの内部パラメータ行列、 $T_{mr}$ はカメラ座標 $x_r$ からカメラ座標 $x_m$ への変換行列、 $\pi$ は $(x,y,z)$ を $(x/z, y/z)$ に射影する関数。また $\pi^{-1}(u,d)$ は、射影後の点を、深度の逆数 $d$ としたときに、カメラ空間に戻す関数です。

同じ点は違う角度から見ても同じ輝度であると仮定します。

photometric errorを以下で定義します。
$\displaystyle \rho_r(\boldsymbol{I}_m,\boldsymbol{u}, d) =\boldsymbol{I}_r(\boldsymbol{u}) - \boldsymbol{I}_m(\pi(KT_{mr}\pi^{-1}(\boldsymbol{u},d)))$

すると以下の式を最小化するようなdを計算する問題になります。

$\displaystyle \boldsymbol{C}_r(\boldsymbol{u},d) = \frac{1}{|\mathcal{I} (r)|}\sum_{m\in\mathcal{I}(r)}||\rho_r(\boldsymbol{I}_m,\boldsymbol{u},d)||_1$

（右辺は隣接する数フレーム分の平均をとっています。）

しかしながら、これをそのまま計算しても、ノイズっぽくなってしまい残念な感じになります。
この計算自体は各ピクセル毎に独立に最小化すればいいわけですが、独立に行うことから滑らかさが一切ないことになります。

結果が滑らかになるように以下のようにregularisation term（第１項）を付け加えます。

$\displaystyle E_\boldsymbol{\xi}=\int_\Omega{g(\boldsymbol{u})||\nabla\boldsymbol{\xi}(\boldsymbol{u})||_\epsilon + \lambda \boldsymbol{C}(\boldsymbol{u},\boldsymbol{\xi}(\boldsymbol{u}) )}d\boldsymbol{u}$

$E_\boldsymbol{\xi}$ を最小化する問題として定式化します。
ここで $\boldsymbol{\xi}(\boldsymbol{u})$ ]は深度の逆数であり、求めたい関数です。
重み項 $g(\boldsymbol{u})$ は輝度の勾配から定義され、
$\displaystyle g(\boldsymbol{u}) = e^{-\alpha||\nabla\boldsymbol{I}_r(\boldsymbol{u})||_2^\beta}$

ここで||・|| $_\epsilon$ はHuberノルムで、
f:id:wacchoz:20181106233309p:plain
のように、 $L_2^2$ と $L_1$ の組み合わせです。
$L_2^2$ は深度を滑らかにするためであり、 $L_1$ は深度が不連続に変化することを許すそうです。

$g(\boldsymbol{u})$ の効果により、輝度の勾配が大きいところでは、深度が不連続になることを許しています。

$E_\boldsymbol{\xi}$ の第１項は凸関数で性質が良いのですが、第２項は凸関数ではなく、このままでは扱いづらい。
そこで以下の関数の最小化問題を解くことにします。

$\displaystyle E_{\boldsymbol{\xi},\boldsymbol{\alpha}}=\int_\Omega\{g(\boldsymbol{u})||\nabla\boldsymbol{\xi}(\boldsymbol{u})||_\epsilon + \frac{1}{2\theta}(\boldsymbol{\xi}(\boldsymbol{u})-\boldsymbol{\alpha}(\boldsymbol{u}))^2 + \lambda \boldsymbol{C}(\boldsymbol{u},\boldsymbol{\alpha}(\boldsymbol{u}) )\}d\boldsymbol{u}$

$\boldsymbol{\xi}$ と $\boldsymbol{\alpha}$ は両方とも未知なので、一見すると問題が複雑になっただけに見えます。
まずは $\theta\to 0$ とすると $\boldsymbol{\xi}=\boldsymbol{\alpha}$ となるのに注意。
そして $\boldsymbol{\alpha}$ を固定すると、
$\displaystyle E_{\boldsymbol{\xi}}=\int_\Omega\{g(\boldsymbol{u})||\nabla\boldsymbol{\xi}(\boldsymbol{u})||_\epsilon + \frac{1}{2\theta}(\boldsymbol{\xi}(\boldsymbol{u})-\boldsymbol{\alpha}(\boldsymbol{u}))^2 )\}d\boldsymbol{u}$
となりますが、この問題は $L_2^2$ -ROF denoising問題の類似であり、効率的な解法が知られています。
また $\boldsymbol{\xi}$ を固定すると、
$\displaystyle E_{\boldsymbol{\alpha}}=\int_\Omega\{ \frac{1}{2\theta}(\boldsymbol{\xi}(\boldsymbol{u})-\boldsymbol{\alpha}(\boldsymbol{u}))^2 + \lambda \boldsymbol{C}(\boldsymbol{u},\boldsymbol{\alpha}(\boldsymbol{u}) )\}d\boldsymbol{u}$
となり、これは各ピクセルごとに最小化することができます。

つまり、
(1)初期値：　 $\boldsymbol{\alpha}_0=\boldsymbol{\xi}_0=\min_d \boldsymbol{C}_r(\boldsymbol{u},d)$
(2) $\boldsymbol{\alpha}$ を固定し、 $E_{\boldsymbol{\xi},\boldsymbol{\alpha}}$ を最小化する。
(3) $\boldsymbol{\xi}$ を固定し、 $E_{\boldsymbol{\xi},\boldsymbol{\alpha}}$ を最小化する。
(4) $\theta$ を小さくする
を繰り返し収束させていけばよさげです。

ここでLegendre-Fenchel変換のお勉強を少し。
このあたりは以下の論文を読むとよいかもです。
Ankur Handa, Richard A. Newcombe, Adrien Angeli, Andrew J. Davison, Applications of Legendre-Fenchel transformation to computer vision problems

一般的な問題として、
$\displaystyle \min_{x\in X}\{ F(Kx)+G(X) \}$ という問題を考えます。
Legendre-Fenchel変換の定義から、
$\displaystyle F(Kx)=\max_{y\in Y}\{ \langle Kx, y \rangle - F^*(y) \}$
となります。
最初の最小化問題は、Dual Formとして
$\displaystyle \min_{x\in X} \max_{y\in Y}\{ \langle Kx, y\rangle - F^*(y) + G(x) \}$
と書けます。
$F^*$ は凸なので、 $G(x)$ が凸であれば、鞍点を求めるシンプルな問題になります。（多分）

さて元の問題に戻って、Legendre-Fenchelを適用してみます。その前に元の問題をベクトルの形で少し書き換えます。こんな感じで行列を１列に並べます。
f:id:wacchoz:20181107213357p:plain

$\boldsymbol{\xi}$ のベクトルバージョンとして $\boldsymbol{d}$ 、 $\boldsymbol{\alpha}$ のベクトルバージョンを $\boldsymbol{a}$ 、 $g(\boldsymbol{u})$ もベクトルバージョンを $\boldsymbol{g}$ と表記することにします。
ベクトル表記すると
$\displaystyle E=||AG\boldsymbol{d}||_\epsilon+\frac{1}{2\theta}||\boldsymbol{d}-\boldsymbol{a}||_2^2+\lambda \boldsymbol{C}(\boldsymbol{a})$
となります。
（とあるんだけど、ならないと思う。最初の積分の式で、 $g(\boldsymbol{u})$ を中に入れて、 $||\nabla \boldsymbol{g}\boldsymbol{\xi}||$ にしないとならない気がする。こうしても次頁のアルゴリズムとも整合性がとれないような…　結局のところ $||\boldsymbol{g}\nabla \boldsymbol{\xi}||$ ではなかろうか）

Legendre-Fenchel変換により
$\displaystyle ||AG\boldsymbol{d}||_\epsilon=\arg\max_{\boldsymbol{q},||\boldsymbol{q}||_2\le 1} \{ \langle AG\boldsymbol{d},\boldsymbol{q} \rangle -\delta_q(\boldsymbol{q}) - \frac{\epsilon}{2}||\boldsymbol{q}||_2^2\}$
となります。
ここで $A$ は、乗算でgradを計算できるような成分をもった行列です。
（Mingqiang Zhu, Fast Numerical Algorithms for Total Variation Based Image Restorationを参照）
また $G=\text{diag}(g)$ です。
f:id:wacchoz:20181107214509p:plain

これにより
$\displaystyle \max_{\boldsymbol{q},||\boldsymbol{q}||_2\le1}\{ \min_{\boldsymbol{d},\boldsymbol{a}}\boldsymbol{E}(\boldsymbol{d},\boldsymbol{a},\boldsymbol{q}) \}$

$\displaystyle \boldsymbol{E}(\boldsymbol{d},\boldsymbol{a},\boldsymbol{q}) = \{ \langle AG\boldsymbol{d},\boldsymbol{q}\rangle + \frac{1}{2\theta}||\boldsymbol{d}-\boldsymbol{a}||_2^2 + \lambda \boldsymbol{C}(\boldsymbol{a}) -\delta_q(\boldsymbol{q})-\frac{\epsilon}{2}||\boldsymbol{q}||_2^2 \}$
と書けます。

(1) $\displaystyle \frac{\partial E(\boldsymbol{d},\boldsymbol{a},\boldsymbol{q})}{\partial \boldsymbol{d}}$ , $\displaystyle \frac{\partial E(\boldsymbol{d},\boldsymbol{a},\boldsymbol{q})}{\partial \boldsymbol{q}}$ を計算する。
(2) $\boldsymbol{a}$ を固定して、 $\boldsymbol{q}$ に関してminを探しに山を下り、 $\boldsymbol{d}$ に関してmaxを探しに山を上る。
(3)次に $\boldsymbol{d}$ を固定して、 $\boldsymbol{a}$ に関してminを探す。これはピクセル毎に行える。
(4) $\theta$ を減らす
を繰り返せばよいようです。

論文中の(1)の $d^{n+1}$ の式は符号が１箇所間違っているかと思います。

やはり以下の論文の7.4あたりを読むといいかも。 Ankur Handa, Richard A. Newcombe, Adrien Angeli, Andrew J. Davison, Applications of Legendre-Fenchel transformation to computer vision problems

(2)のステップはGPU等で並列計算できます。多分(3)も出来ます。

PS. 論文の数式は何箇所か間違ってそうです。自分で式変形していったほうがいいかも。