Intelligence Architecture けんきうノート

◆ GitBookに移行します

TV最小化(ROF版)ノイズ除去

[Math rendered by jsMath]

- 概要

こっちのほうがオリジナル。 TV最小化アルゴリズムはChambolle版(途中の導出のしかたはかなり異なるけど)です。

Color-TVにしたりとか、制約条件変えてみたりとかしたい場合、TV最小化アルゴリズムと同じようにはできないかもーな気がするので、復習です。

- 問題定式化

平均がゼロ、分散が既知のノイズが乗った観測画像 \(u_0\) から、元の画像 \(u\) を以下のように推定します。

\[\eqalign{ u & = \arg\min\int|\nabla u|dx \\ \text{ s.t. } & \int u \ dx = \int u_0 dx \\ & \int \frac12(u-u_0)^2 dx = \sigma^2 }\]

- 変分法適用

2つの制約条件に対するラングランジュ乗数を \(\alpha,\lambda\) とすると、 ラグランジアンは \[ L(u,\alpha,\lambda) = \int|\nabla u|dx + \alpha \int(u-u_0)dx + \lambda \left( \int \frac12(u-u_0)^2 dx - \sigma^2 \right) \] となります。

あとは \(L\) を \(u\)(と\(\alpha,\lambda\))で微分してゼロとすれば最適性条件が出るわけですが、\(L\) は関数 \(u\) を変関数とする汎関数なので、やっぱり微分ではなく変分を使います。

Euler-Lagrange方程式を使ってもいいのですが自分でやってみましょう!

\(u\rightarrow u+\delta u\) としたときの \(L\) の第一変分 \(\delta L\) を計算します。

\[\eqalign{ \delta L &= L(u+\delta u) - L(u) \\ &= \int|\nabla u+\nabla\delta u|dx + \alpha \int(\delta u+u-u_0)dx + \lambda \left( \int \frac12(\delta u+u-u_0)^2 dx - \sigma^2 \right) - L(u) \\ &= \int\sqrt{(\nabla u)^2+2\nabla u\cdot\nabla\delta u+O(\delta u^2)}dx - \int|\nabla u|dx + \alpha \int \delta u \ dx + \lambda \int (\delta u(u-u_0)+O(\delta u^2)) dx \\ &= \int|\nabla u|\sqrt{1+2{\nabla u\cdot\nabla\delta u \over (\nabla u)^2}}dx - \int|\nabla u|dx + \alpha \int \delta u \ dx + \lambda \int \delta u(u-u_0) dx \\ &= \int|\nabla u|\left(1+{\nabla u\cdot\nabla\delta u \over (\nabla u)^2}+O(\delta u^2)\right)dx - \int|\nabla u|dx + \alpha \int \delta u \ dx + \lambda \int \delta u(u-u_0) dx \\ &= \int{\nabla u \over |\nabla u|}\cdot\nabla\delta u \ dx + \alpha \int \delta u \ dx + \lambda \int \delta u(u-u_0) dx \\ &= -\int\left(\nabla\cdot{\nabla u \over |\nabla u|}\right)\delta u \ dx + \alpha \int \delta u \ dx + \lambda \int \delta u(u-u_0) dx \\ &= -\int\delta u\left( \nabla\cdot{\nabla u \over |\nabla u|} -\alpha -\lambda(u-u_0) \right) \ dx }\]

途中、自明でなさそうなとこを以下に列挙&注記します:

  • 4つ目と6つ目の等号:第一変分のみを考えればいいので \(O(\delta u^2) \rightarrow 0\) としています。
  • 5つ目の等号:ルートの部分は \(\sqrt{1+x}\) のマクローリン展開をしています。
  • 7つ目の等号:部分積分とグリーンの公式を使います(TVについての補足参照)。 ここで、定義域の境界上で \(\nabla u\) と境界線とが平行でなければならないという条件が出ます。

したがって、任意の \(\delta u\) に対して \(\delta L=0\) となるためには、 \[ \nabla\cdot{\nabla u \over |\nabla u|} -\alpha -\lambda(u-u_0)=0 \] が必要条件であり、これがすなわち \(u\) の最適性条件となります。

- 最急降下法

TV最小化アルゴリズムでの最急降下法と同様に考えて、 \[ u \leftarrow u - \Delta t \ \delta L|_{\delta u=\delta(x)} \] と更新します。 時間発展を \(u(n)=u(x; t=n\Delta t), \ (n=0,1,2\ldots)\) のように表すことにして、 \[ u(n+1) = u(n) + \Delta t \left( \nabla\cdot{\nabla u(n) \over |\nabla u(n)|} -\lambda(n)(u(n)-u_0) \right) \] とします。 ここで \(\alpha\) の項を勝手に落としていますが、それに対応する制約 \(\int u \ dx = \int u_0 dx \) は結果的に守られることを以下のように確認できます。

まず更新式を両辺積分します。 \[ \int u(n+1)dx = \int u(n)dx + \Delta t \left( \int \nabla\cdot{\nabla u(n) \over |\nabla u(n)|}dx -\lambda(n)\int(u(n)-u_0)dx \right) \] ここで \(\int \nabla\cdot{\nabla u(n) \over |\nabla u(n)|}dx = \oint {\nabla u(n) \over |\nabla u(n)|}dn\) で \(\nabla u\) は \(dn\) に直交するので結局0になります。 よって \[ \int u(n+1)dx = \int u(n)dx - \Delta t \ \lambda(n) \int(u(n)-u_0)dx \] となります。 あとは \(u(0)=u_0\) を前提とすれば帰納的に、

  • 前提より \(\int u(0)dx=\int u_0dx\) なので \(n=0\) のとき成立。
  • \(n=k\) のとき成立する(\(\int u(k)dx=\int u_0dx\))ならば、\(\int u(k+1)dx = \int u(k)dx - \Delta t \ \lambda(k) \int(u(k)-u_0)dx = \int u_0dx \) なので \(n=k+1\) でも成立。
  • よって全ての \(n\) で成立。

とできます。

最後に \(\lambda\) は、最適性条件の式に \((u-u_0)\) をかけてから積分し、 \[ \int \nabla\cdot{\nabla u \over |\nabla u|} \ (u-u_0) dx -\alpha\int(u-u_0)dx -\lambda\int(u-u_0)^2dx=0 \] 2つの制約条件を代入して、 \[ \int \nabla\cdot{\nabla u \over |\nabla u|} \ (u-u_0) dx -2\lambda\sigma^2=0 \] 第1項は部分積分&グリーンの公式を使って、 \[ -\int {\nabla u \over |\nabla u|} \cdot \nabla(u-u_0) dx -2\lambda\sigma^2=0 \] したがって \[ \lambda = -\frac1{2\sigma^2} \int \left( |\nabla u| - {\nabla u\cdot\nabla u_0 \over |\nabla u|} \right) dx \] を得ます。

まとめると \[\left\{\eqalign{ u(0) &= u_0 \\ \lambda(n) &= -\frac1{2\sigma^2} \int \left( |\nabla u(n)| - {\nabla u(n)\cdot\nabla u_0 \over |\nabla u(n)|} \right) dx \\ u(n+1) &= u(n) + \Delta t \left( \nabla\cdot{\nabla u(n) \over |\nabla u(n)|} -\lambda(n)(u(n)-u_0) \right) }\right.\] です。 あとは \(x\) を適切に離散化すればプログラムも書けます。

- サンプルコード

ここにあります。 \(x\) 空間の離散化は色々試しましたが、結局は中心差分に落ち着きました。 ただし \[ \nabla\cdot{\nabla u \over |\nabla u|} = {u_{xx}u_y^2+u_{yy}u_x^2-2u_xu_yu_{xy} \over |\nabla u|^3} \] と変形してから差分化しています。

収束性は残念ながらあまりよくないようです。 この点ではChambolle版(TV最小化アルゴリズム)のほうがいいですね。