2013-04-07

某女の子が歌ったり踊ったりするゲームの何か(2)

どうせ暇なので、過去のプログラムの在庫整理の続き。
シェーダで描画を改善してみる。
かなり見栄えが良くなりました。

シェーダの中身は某プログラムのパクリです。ごめんなさい。
いずれどうにかするかも。

ちゃんとフラグでシェーダを切り替えてないので、正しく描画されないモデルもあるけど、そのうち・・・
これは結構大変でいつになるやら・・・

「マルチバイト文字セットを使用する」でビルドして下さい。
ライブラリはCgを追加で使用しています。

ソースコード

github.com

(2013/5/5差し替え)
(2018/11 GitHubに移動)

こんな感じです。

f:id:wacchoz:20181110232247p:plain

2012-12-09

スキニングの話(3)

次はDual Quaternionを使ったスキニング手法。

Ladislav Kavan et al.,2008, Geometric skinning with approximate dual quaternion blending, ACM Trans. Graph. 27, 4, 105.

Quaternionは回転しか表せず、平行移動は表せないのがとっても不満です。
Dual Quaternionはその不満を解消するツールです。

Dual Quaternionは

$\displaystyle \hat{\mathbf{q}}=\mathbf{q}_0+\epsilon\mathbf{q}_\epsilon$

のように定義されます。
ここで $\epsilon$ は"dual unit"と呼ばれ、 $\epsilon^2=0$ という特殊な規則に従う数です。
$\mathbf{q}_0$ 、 $\mathbf{q}_\epsilon$ は普通のQuaternionであり、前者はnon-dual part、後者はdual partと呼ぶことにします。
Quaternionと $\epsilon$ の積は可換です。

さて、回転を表すdual quaternionはnon-dual partの $\mathbf{q}_\epsilon$ を0とし、

$\displaystyle \hat{\mathbf{q}}=\mathbf{q}_0$

とするだけ。つまり今までのquaternionと同じだ。

では今までと変わったのは何かというと、平行移動を表せることだ。
平行移動 $(t_0 ,t_1, t_2)$ を表すには

$\displaystyle \hat{\mathbf{t}}=1+\frac{\epsilon}{2}(t_0\mathbf{i}+t_1\mathbf{j}+t_2\mathbf{k})$

のようnon-dual partを1とし、dual partに平行移動成分を持たせたものとなる。
2で割るのは回転のquaternionのときと似ている。

平行移動と回転移動を合成するには普通に乗算すればよい。

さて、実際に頂点にdual quaternionを作用させるには、まず頂点座標 $\mathbf{v}$ を

$\displaystyle \hat{\mathbf{v}}=1+\epsilon(v_0\mathbf{i}+v_1\mathbf{j}+v_2\mathbf{k})$

のように変換し、

$\displaystyle \mathbf{v}'=\hat{\mathbf{q}}\hat{\mathbf{v}}\overline{\hat{\mathbf{q}}^*}$

とすればよい。

回転と平行移動を合成するというよりは、平行移動しながら回転するという螺旋的なイメージの方が合っている。
実際、

$\displaystyle \hat{\mathbf{q}}=\cos\frac{\hat{\theta}}{2}+\hat{\mathbf{s}}\sin\frac{\hat{\theta}}{2}$

というdual quaternionを考えてみる。 $\theta$ と $s$ は両方ともdual quaternionである。
これは螺旋軸に沿った回転と平行移動を表す。
$\theta_0$ は回転角を表し、 $\theta_\epsilon$ は螺旋軸に沿った平行移動量を表す。
また $s_0$ は螺旋軸の向きの単位ベクトルを表す。
$s_\epsilon$ だけわかりにくいのだが、回転中心を $r$ としたとき、 $s_\epsilon=r\times s_0$ となる。

さてこれをスキニングにどう使うのか。
まず回転を表すquaternionと平行移動を合成して、dual quaternionを作る。
これをボーンウェイトに従って線形補間する。
このdual quaternionを使って頂点を変換すればよい。
この変換について、論文中に高速化の方法が書かれているため、実装する際には参照するとよい。

2012-11-30

スキニングの話(2)

Linear Blend Skinningでしばしばおかしな変形を目にすることになる。
一例を図に示す。
（Ladislav Kavan et al.,Geometric Skinning with Approximate Dual Quaternion Blendingより引用）

f:id:wacchoz:20181110223349p:plain:w450

頂点 $\mathbf{v}$ は $j_1$ と $j_2$ から50%ずつ影響されているとする。
ボーンアニメーションにより、 $j_1$ は動かず、 $j_2$ が180°回転したとする。
このとき $j_1$ の影響で $C_{j1}\mathbf{v}$ に変換され、 $j_2$ の影響で $C_{j2}\mathbf{v}$ に変換されるため、50%ずつの重みにより平均した結果、 $\mathbf{v}'$ は $j_2$ に一致してしまい、完全にメッシュがつぶれてしまう。
"candy-wrapper"と呼ばれる現象である。

これは重みつき平均を計算する際、単純に線形に平均してしまったことに起因する。
$j_1$ は回転せず、 $j_2$ は180°回転しているため、本来であれば $\mathbf{v}'$ は90°回転となってほしいところである。
（+90°なのか-90°なのかという曖昧さは残っているが…）

この例を見ていると、回転をQuaternionでの球面線形補間とするという発想が生まれてくる。
この手法は以下の論文に論じられている。
KAVAN, L., AND ZARA, J. 2005. Spherical blend skinning: A　real-time deformation of articulated models. In 2005 ACM SIGGRAPH　Symposium on Interactive 3D Graphics and Games,　ACM Press, 9–16

回転を補間するのであるが、回転の中心点がどこであるかというのが論文中では詳しく書かれている。
上例のように隣接するボーン2つから影響される場合は単純である。回転中心は誰がどう見ても $j_2$ である。

さて $C=B(t)\times M^{-1}$ を回転部 $Q_k$ と移動部 $m_k$ に分解しよう。

$\displaystyle C_k= \begin{pmatrix} Q_k & \mathbf{m}_k \\\ \mathbf{0}^T & 1 \end{pmatrix}$

とし、さらに $Q_k$ をクォータニオン $q_k$ に変換しよう。

$C_{j1}$ と $C_{j2}$ の回転部のクォータニオンを適切に補間したものを $Q$ とすると、

$\displaystyle \mathbf{v}'=Q(\mathbf{v}-\mathbf{j}_2)+\mathbf{j}_2$ ]

とすればよい。
これがSpherical Blend Skinningの手法である。
クォータニオンの補間に球面線形補間を使うのが自然であるが、補間計算は頂点毎に行う必要があるため、球面線形補間を行うのは処理が重過ぎる。（arccosが1回にsinが2回）
単純に線形補間を行い正規化してもそれほど誤差がないそうである。

さて上例では隣接したボーンから影響される頂点について考えていた。
このため、 $j_2$ を回転中心に選ぶというシンプルな計算であった。

実際には3本以上のボーンから影響される頂点もあるし、隣接しないボーン２本から影響される場合もある。
このような場合にはどのように回転中心を選べばよいであろうか？

まず上例の場合では、回転中心 $r_c(=j_2)$ は

$\displaystyle C_{j_1}\mathbf{r}_c=C_{j_2}\mathbf{r}_c$

という関係が成り立っていたことに注目する。

一般に $n$ 個のボーンに影響されている場合、回転中心 $r_c$ は

$\displaystyle C_a\mathbf{r}_c=C_b\mathbf{r}_c,\ a,b\in \{j_1,\dots,j_n\}$

により定義することができる。

これは $n$ が3以上では優決定系であり、一般的には解が存在しないため、最小二乗法により解くことになる。
回転中心が決まれば、以下の式で変換後の頂点が計算できる。

$\displaystyle \mathbf{v}'=Q(\mathbf{v}-\mathbf{r}_c)+\sum_{i=1}^n w_i C_{j_i}\mathbf{r}_c$

しかしながらこの手法が3本以上のボーンなどの一般の場合に使われるかどうか微妙である。
何しろ回転中心の計算はボーンの組み合わせ毎に行わなければならず、面倒すぎるのだ。

2012-11-28

スキニングの話(1)

どうせ暇だしスキニングについて書いてみる。

スキニングとは3Dアニメーションで使われる手法で、ポリゴン内部にボーンを仕込んでおき、ボーンを動かすことでポリゴンを変形させる手法である。

骨格構造はジョイント（関節）の集合である。
そしてジョイントとは単なる点ではなく、ジョイントを原点とする座標系を意味すると考えたほうがよい。それぞれのジョイントに $xyz$ 座標系があるわけだ。
各ポリゴンはそれぞれのジョイント（で定義される座標系）からどれだけ影響されているかという重みを与えられている。
ではボーンは何かというと、親子関係のあるジョイント間を線で結んだだけのものである。人間にとってわかりやすくするための便宜上のものといえる。

f:id:wacchoz:20181108230226p:plain:w360

まず初期姿勢（rest pose）で考えよう。
両手を左右に伸ばして立っているあの姿勢である。

$M_{10}$ でジョイント $j_1$ を原点とする座標系からモデル空間へ変換行列を表すことにする。
たとえば $j_1$ を原点とする座標系の原点は、モデル空間での座標は、 $M_{10}\times (0,0,0,1)^T$ で表すことができる。

同様に、 $M_{21}$ でジョイント $j_2$ を原点とする座標系からジョイント $j_1$ を原点とする座標系に変換する行列としよう。

ジョイントは親から順に $1, 2, 3,\dots , k$ と繋がっているとすると、ジョイント $j_k$ を原点とする座標系をモデル空間に変換する行列は、

$\displaystyle M_{k0}=M_{k,k-1}\times\dots\times M_{21}\times M_{10}$

と書ける。またモデル空間からジョイント $j_k$ を原点とする座標系への変換は、この行列の逆行列であり、

$\displaystyle M_{0k}=M_{k0}{}^{-1}$

と書ける。

さてアニメーションとは、 $M_{k0}$ を時刻 $t$ によって刻々と変化させていくことに他ならない。

行列 $M$ は初期姿勢のときの行列として用いたため、アニメーションにより変化したときの行列は $B$ で表すことにしよう。

時刻 $t$ でのジョイント $j_k$ からモデル空間への変換行列を $B_{k0}(t)$ とする。

同様に

$\displaystyle B_{k0}(t)=B_{k,k-1}(t)\times\dots\times B_{21}(t)\times B_{10}(t)$

である。

ある初期姿勢でのメッシュ上の点 $v$ がジョイント $j_k$ にのみ影響されているとしよう。

このとき点 $\mathbf{v}$ の、ジョイント $j_k$ を原点とする座標系での座標は、 $M_{k0}{}^{-1} \mathbf{v}$ となる。

この点を時刻 $t$ でのモデル空間での座標に変換すると、

$\displaystyle v' = B_{k0}(t)\times M_{k0}{}^{-1}\times \mathbf{v}$

となる。これがボーンアニメーションの基本である。
しばしば $C_{k}(t) = B_{k0}(t) \times M_{k0}{}^{-1}$ とまとめて書かれることもある。

ここまでは点 $\mathbf{v}$ がジョイント $j_k$ のみに影響されるとしてきたが、複数のジョイントに影響される場合も同様であり、ジョイント $j_{k1}, j_{k2}, \dots, j_{kn}$ にそれぞれ重み $w1, w2, \dots , wn$ で影響されるとした場合、

$\displaystyle \mathbf{v'}=\sum_{i=1}^n w_i B_{ki,0}(t)M_{ki,0}{}^{-1}\mathbf{v}$

とする方法がLinear Blend Skinningである。この方法は広く用いられ、そしてしばしば問題を引き起こす。（つづく）

2012-11-04

Quaternionと回転の話

どうせ暇だしQuaternionについてまとめてみる。

Quaternion（クォータニオン）は３次元グラフィクス分野でよく出てくるツールである。
もともとはHamiltonによって考案された概念で、複素数の自然な拡張となっている。

単位複素数がGauss平面上での回転を表すことは高校数学で学んだ通り。
原点周りにθだけ回転させるには、複素数に

$\displaystyle u=e^{i\theta} = \cos\theta + i\sin\theta$

を乗算すればよいのであった。

さて３次元上での回転は、複素数を拡張したQuaternionで表すことができる。

以下のように定義される。
$\displaystyle \mathbf{u}=q_w+q_x\mathbf{i}+q_y\mathbf{j}+q_z\mathbf{k}$
ここで、i, j, kは
$\displaystyle \mathbf{i}^2=\mathbf{j}^2=\mathbf{k}^2=\mathbf{ijk}=-1$
と定義される。

この定義からすぐに従うが、
$\displaystyle \mathbf{ij}=-\mathbf{ji}$
となる。つまり乗法は非可換である。
３次元上の回転は回転の順序により結果が変わるというのは広く知られるところだが、Quaternionの非可換性はこれと対応するものだ。

Quaternionは複素数の拡張であると言ったが、類似する性質を見ていこう。
単位Quaternionを
$\displaystyle \mathbf{u}=x\mathbf{i}+y\mathbf{j}+z\mathbf{k}, \mbox{ where } x^2+y^2+z^2=1$
と書こう。
このとき複素数と同様に $\displaystyle \exp(\mathbf{u}\theta)=\cos\theta+\mathbf{u}\sin\theta$
が成り立つのである。

なお、一般に単位Quaternionは
$\displaystyle \mathbf{q}=\cos\theta+\mathbf{u}\sin\theta$
と書ける。

さて、Quaternionを用いてどのように点を回転させるのか見ていこう。
$\displaystyle \mathbf{q}=\cos\frac{\theta}{2}+v_x\sin\frac{\theta}{2}\mathbf{i}+v_y\sin\frac{\theta}{2}\mathbf{j}+v_z\sin\frac{\theta}{2}\mathbf{k}$
これが回転軸 $(v_x, v_y, v_z)$ の回りに回転角 $\theta$ での回転をあらわすQuaternionである。

複素数の場合とは違い、単純に頂点座標に乗算するわけではない。
頂点座標 $p=(p_x, p_y, p_z)$ を
$\displaystyle \mathbf{p}=p_x\mathbf{i}+p_y\mathbf{j}+p_z\mathbf{k}$
のようにQuaternionに変換しておく。
このとき、
$\displaystyle \mathbf{p}'=\mathbf{q}\mathbf{p}\mathbf{q}^{*}$
のように左右から乗算することで、p'の各成分が、頂点pを回転した座標となる。

回転の合成は回転行列の場合と同様で、Quaternionを単純に乗算すればよい。

Quaternionが重宝される理由の１つは回転の補間が容易にできることである。
３次元グラフィクスでは、回転 $q0$ と回転 $q1$ の間を最短距離で等速になるよう補間したいと思うことがしばしばある。
Quaternionを使うことにより、球面線形補間とよばれる式でこのような補間が可能となる。
以下の $\mathbf{q}'$ が $\mathbf{q}_0$ と $\mathbf{q}_1$ の球面線形補間である。
$\displaystyle \mathbf{q}'=\frac{\sin(\theta(1-t))}{\sin\theta}\mathbf{q_0}+\frac{\sin(\theta t)}{\sin\theta}\mathbf{q_1}$

ここで $\displaystyle \cos\theta=(\mathbf{q_0},\mathbf{q_1})$ である。
なお $\displaystyle (\mathbf{q_0},\mathbf{q_1})＜0$ のときは $\mathbf{q}_0$ を $-\mathbf{q}_0$ で置き換える必要があるとのこと。
（ $\mathbf{q}_0$ も $-\mathbf{q}_0$ も同じ回転を表す）

さて少し違う話題。

複素数に戻ろう。 $\displaystyle u=e^{i\theta} = \cos\theta + i\sin\theta$
について対数をとると $\displaystyle \log{u} = i\theta$
となる。（簡便のため、対数の主値のみを記載）
これを見ると、対数をとると回転角を取り出せることがわかる。

実は同様のアナロジーがQuaternionにも成り立つ。
$\displaystyle \mathbf{q}=\cos\frac{\theta}{2}+v_x\sin\frac{\theta}{2}\mathbf{i}+v_y\sin\frac{\theta}{2}\mathbf{j}+v_z\sin\frac{\theta}{2}\mathbf{k}$
について、
$\displaystyle \log{\mathbf{q}} = \frac{\theta}{2}(v_x\mathbf{i}+v_y\mathbf{j}+v_z\mathbf{k})$
が成り立つ。
つまり対数をとると、回転角と回転軸が取り出されるのである。

ネイピアが対数を発明したモチベーションは、乗算や除算を加減算で行うことで煩雑な計算を楽にすることであった。
Quaternionについて対数を用いて乗算を加算にすることができるだろうか。
そもそも対数で乗算を加算にできるのは、
$\displaystyle e^{x+y}=e^x e^y$
が成り立つためであった。
Quaternionの場合には乗算の非可換性のために、この式が成り立たない。
Baker-Campbell-Hausdorffの公式とよばれる $\displaystyle \exp(x) \exp(y)=\exp\left(x+y+\frac{1}{2}[x,y]+\frac{1}{12}[x,[x,y]])-\frac{1}{12}[y,[x,y]]-\dots\right)$
が成り立つ。
つまりQuaternionの場合には、対数をとり加算してから指数関数をとっても、乗算した場合と結果が一致しないのである。
しかしながら近似的には成り立ちそれなりに便利なため、たまに使われるようだ。

2012-10-28

某女の子が歌ったり踊ったりするゲームの何か

Standing on the shoulders of giants

２年以上前に書いてたプログラムの在庫整理。
とりあえずはモデルの表示まで。
かつてOpenGLバージョンとDirectXバージョンを作ったけど、今回はOpenGLバージョンを整理。
整理というほどは整理できてないんだけど。。。

ライブラリはGLUT、GLUI、libsquish、DirectX9を使用しています。
libsquishはVS2010以上の*.slnファイルが用意されていないので自分で作って下さい。
DirectX9は今回のうｐに限って言えばD3DXVECTOR3などの型を利用してるだけなので、ほぼ不要です。

GLUIは初めて使ってみたけど、回転操作が簡単に実装できるのが良いです。
正しい使い方をしているのか自信が無いけどあしからず。

テクスチャはDXT形式なので、DirectXならヘッダを付ければそのまま読めるのですが、OpenGLなのでlibsquishでRGBAに展開しています。

*.exeファイルのショートカットを作って、ショートカットの上にファイルをドラッグ＆ドロップすると便利です。

ソースコード

github.com

(2013/4/6 差替え／前にうｐしたものがコンパイル通ってませんでした)

(2018/11 GitHubに移動)

こんな感じです。

f:id:wacchoz:20181107222805p:plain

ワイヤフレーム表示

f:id:wacchoz:20181107222811p:plain

2012-10-08

DTAM論文を読んでみた

どうせ暇だし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. 論文の数式は何箇所か間違ってそうです。自分で式変形していったほうがいいかも。

wacchoz’s note

プログラミングとか数学について