Boy Meets ML

機械工出身者が機械学習やその界隈の知識を考える

トップ > AOS > SfMについて #3.6：AOS（KAZEの非線形拡散フィルタの計算）

SfMについて #3.6：AOS（KAZEの非線形拡散フィルタの計算）

AOS KAZE SfM features 特徴量 特徴点

はじめに

前回記事ではKAZEの論文を確認し、どのように特徴量抽出を行っているのかを確認しました。 ただし、KAZEで利用している非線形拡散フィルタである偏微分方程式を解く方法AOS(Additive Operator Splitting)について、 引用されている文献を確認しないとよくわからなかったので、本記事はその補足を行うことを目的とします。

mlengineer.hatenablog.com

導入

KAZEの論文ではAOSに対してJ. Weickertらの論文1を引用しています。 1998年の論文であり、少し古いですが読み込んでKAZEを理解したいと思います。 実はKAZEの改良版であるA-KAZEではAOSは使われていませんが、A-KAZEでも非線形拡散フィルタを解くことに変わりはありませんので、 理解も進むかと思います。J. Weickertらの論文でも順を追ってAOSが必要になった理由が書かれていますので、 本記事もその流れに沿って概説していきます。

偏微分方程式を離散化する

画像をベクトル f \in \it{R}^mとして、このベクトルは\Omega = (0,a_1)\times\dots(0,a_M)というM次元空間から写像されたものだと考えます。 カラー画像はM=2の平面空間をm=3のベクトルに変換していると考えるようです。

非線形拡散フィルタは時間tに対して偏微分方程式によって表現されます。 これはCLMCフィルタと呼ばれるようで次式の形です。

\partial_t u = div(g(|\nabla u_{\sigma}|^2)|\nabla u)

これはu(x,t)という画像の時間的変化を表現した式であり、\nablaは勾配を計算するオペレータです。 この段階では細かいことは気にしないでください。 重要なことはあくまでも時間tで表現される連続値であることです。 コンピュータで扱うため、かつ時間tである空間を、画像処理で扱うスケール空間\sigmaに対応させていく必要があるのです。

連続値で表された微分方程式を離散化して解くにあたって、 時刻t=0のときの初期値は、u(x,t=0) = f(x)と定めます。 画像の端(論文では\partial \Omegaと書いてありましたが)に対しては境界条件として\nabla_{n} u=0とします。 これはいわゆるノイマン問題としての仮定だと思います。

拡散という現象は時間が経過するほど詳細な情報が薄れていき、均一に、つまり単純な情報を示すようになります。 画像処理に対して拡散フィルタを施すということは、時間を大きくするほど高周波数が減っていき単純な画像になっていくことになります。 一般的な物理現象としての拡散については以下の記事がわかりやすいです。

shimaphoto03.com

さて、非線形拡散フィルタを施すことで画像をボヤけさせていきますが、エッジは特徴的であると考えられるため、 エッジは他よりもボヤけさせたくありません。 なので、エッジの強さを表す |\nabla u_{\sigma}|が大きいほど、拡散の強さ gが低減するような形にしてあげます。

\nabla u_{\sigma}=\nabla (K_{\sigma}*u)

として、K_{\sigma}はガウシアンカーネルです。 そして、拡散の強さgを次のようにしてあげればいいわけです。

g(s) = \begin{cases} 1 & (s\leq0) \\ 1 - \exp{\left\{ -\frac{3.315}{\left(s/ \lambda \right)^4} \right\}} & (s>0) \end{cases}

KAZEの記事に書いたものと若干異っていますが、論文に従って記載しています。 これらの連続値である計算をexplicitな計算で離散化していきます。

整理すると、連続値バージョンでは、

\partial_t u = \partial_x (g(|\partial_x u_{\sigma}|^2) \partial_x u)

であった計算を離散化すると、

\frac{u^{k+1}_i - u^{k}_i }{\tau} = \sum_{j\in\it{N}(i)}^{}{\frac{g^k_j + g^k_i}{2h^2} (u^k_j - u^k_i)}

となるのです。 更に行列を駆使して整理してあげると、

\frac{u^{k+1} - u^{k} }{\tau} = A(u^k)u^k

と書けて、KAZE論文で出てきた式になります。 ただし、KAZE論文では謎であったものが、

A(u^k) = (a_{ij}(u^k)), \\ a_{ij}(u^k) = \begin{cases} \frac{g^k_i + g^k_j}{2h^2} & (j\in\it{N}(i)) \\ -\sum_{n\in\it{N}(i)}^{}{\frac{g^k_i + g^k_n}{2h^2} } & (j=i) \\ 0 & (else) \\ \end{cases}

であることが判明します。

結局のところexplicitになることが次式のように整理できるのでわかります。 explicitとは、陽的にかける、つまりk+1の状態をkのみで直接表現できて、他に方程式を解く必要がないことを意味します。

u^{k+1} = (I + \tau A(u^k))u^k

しかし、このexplicitな方程式を安定的に解くためには、いくつかの条件を満たす必要があり、 特に時間幅\tauを次のように設定しなくてはなりません。

\tau = \frac{1}{\max_{i}\sum_{j\neq i}^{}{a_{ij}(u^k)}}

これは文献によると画像処理においては、拡散率gの上限を1にするために、h=1となり、 これは0.5 \gt \tau となってしまいます。 このステップ幅はかなり厳しい制約条件と言えるため、離散化したときの正確さvs安定性という議論につながります。 そして、これを解決するためにsemi-explicitな解法としてAOSが求められた理由になります。

AOS

1次元の拡散フィルタ

まず1次元拡散フィルタの連続値バージョンである、

\partial_t u = \partial_x (g(|\partial_x u_{\sigma}|^2) \partial_x u)

という式を少し複雑に離散化する別の方法を次のように考えるようです。

\frac{u^{k+1} - u^{k} }{\tau} = A(u^k)u^{k+1}

これをいじると、

(I-\tau A(u^k))u^{k+1} = u^k

となりますが、この式ではu^{k+1}に対して陽的な表現にはなっていません。 なぜならば、これを解くためには左辺に存在する線形システムを解く必要があるためで、 このような表現をsemi-explicitといいます。

左辺の行列を

(I-\tau A(u^k))= B(u^k)

と考えます。 実はこの行列は優対角行列というものになります。 優対角行列とは、行列の要素を b_{ij}とすると、

|b_{ii}| \gt \sum_{j\neq i}^{}{|b_{ij}|}

となる特性をもつ行列を言います。 そして優対角行列は正則であるという特性も持ちます。 なので、

B^{-1}(u^k) = Q(u^k)

とできてしまい、semi-explicitにするとどのような時間幅を選んでも収束性を脅かすことなく、所望の正確さをもった離散化が行えてしまうのだ、そうです。

今、Bu=dという線形システムを考えて、Bが三重対角行列であるとします。 つまり次の形をしている行列であり、これはThomas algorithmというガウスの消去法の一種で解くことができます。

B = \begin{pmatrix} \alpha_1 & \beta_1 & & & \\ \gamma_1 & \alpha_2 & \beta_2 & & \\ & \ddots & \ddots & \ddots &\\ & & \gamma_{N-2}& \alpha_{N-1} &\beta_{N-1}\\ & & & \gamma_{N-1} &\alpha_{N}\\ \end{pmatrix}

Thomas algorithmは3ステップの処理を行う方法で、

  1. 行列BをLR分解する
    • B=LRという、下二重対角行列Lと上二重対角行列Rにします
    • すると、LRu=dという形になる
  2. 前進代入
    • Ly=dをyに関して解く
    • この時点におけるフィルタをcausal filterという、らしい
  3. 後退代入
    • Ru=yをuに関して解く
    • この時点におけるフィルタをanticausal filterという、らしい

と解きます。 この処理は計算量がかなり低く、\it{O}(N)で抑えることができるのです。 この辺りの計算は書き下してみると簡単にわかります。

m次元の拡散フィルタ

論文に沿って1次元の拡散フィルタをexpicitに解く場合と、semi-explicitに解く場合を見てきました。 m次元の場合も偏微分方程式のことを考えるとそのまま拡張すればOKです。

\partial_t u = \sum_{l=1}^{m}\partial_{x_l}(g(|\nabla u_{\sigma}|^2)\partial_{x_l}u)

explicitなときは、

u^{k+1} = (I+\tau \sum_{l=1}^{m}{A_l(u^k)})u^k

semi-explicitなときは、

u^{k+1} = (I - \tau \sum_{l=1}^{m}{A_l(u^k)})^{-1}u^k

となります。 m次元のexplicitなときは解の安定性条件が、

\tau \lt \frac{1}{2m}

となります。 m次元の場合は次元が大きくなるほどに時間幅がとても小さくする必要出てきてしまうわけです。

AOSではsemi-explicitなアプローチの利点をフル活用するために、

u^{k+1} = \frac{1}{m}\sum_{l=1}^{m}{(I - m \tau A_l(u^k))}^{-1}u^k

と細工を施す。 このように細工を施しても、元々の微分方程式を1次のテイラー展開してみると、implicitなもの、semi-implicitなものと一貫性を保っているらしい。

Splitting Operator(またはMethod)とは、微分方程式を解きやすくするテクニックで、1つの変数を複数の変数に分けていく手法を指すようです。 乗法(Multiplicative)ではなく加法(Additive)なアプローチであるため、 すべての次元方向に対しても同じ方法で分割(Split)することが可能です。 Multiplicative Splittingは可換性を保てれるとは限らないようです。

普通のガウシアンフィルタを使って拡散させていく方法は、g=1とした線形フィルタリングと等価です。 このとき時間とスケールの関係 t = \frac{\sigma ^2}{2}となりますので、 AOSでも同じ関係を使うことにするらしいです。

AOSは偏微分方程式を解くための工夫であり、数式を追っていくことに疲れて全体の流れを忘れてしまいそうになります。 これまでの処理フローは文献のV.ALGORITHMIC STRUCTUREに記載がありますので、 これを見たほうが全体を通したフローがわかりやすいと思います。

まとめ

今回はKAZEで使われているAOSに関して論文を確認しました。 熱力学や流体力学など機械系でも見覚えのある拡散フィルタでしたが、すっかり忘れてしまっていました。。 AOSと検索してもなかなか出てきませんが、拡散フィルタの離散化は画像処理とは別分野での資料が見つかる可能性が高そうなので、そちらで検索しつつ勉強すると、雰囲気がつかめると思います。

なお、上記の記事は間違いや誤認識が含まれる可能性もありますので、A-KAZEの論文も確認したら、適宜修正したいと思います。

文献

  1. J. Weickert, B. M. T. H. Romeny and M. A. Viergever, “Efficient and reliable schemes for nonlinear diffusion filtering,” in IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 398-410, March 1998.↩