Total Variation Minimization

LDPCはしばらくお休み(基本的に優先順位が低い).今日はTV(Total Variation)の実装について書く.

神奈川大の斉藤先生が国内では精力的に研究していらっしゃる.どうにも魅力的なツールに見えるのだが,意外と名前負けしている(すなわち,ネーミングセンスが抜群)感が拭えない.曰く,画像を骨格画像とテクスチャ画像に分離する手法である.なんとか実装が出来たので,公開しておく.ただし,まだ国際的な最前線はおろか,斉藤先生にも追いついていないことは明記しておく.

まず,TVとは何か,これを使うと何が出来るかについてざっくり話す.さらに,文献を挙げておく.最後に実装の方針とソースコードを載せる.

TVとは,不連続を許す空間におけるノルムのはかり方です(多分).不連続だから微分が出来ないところがあります.逆に,値が連続的でない場合も扱えます.実は,画像データはエッジが含まれているので,不連続を許す空間と相性が良さそうなのです.

実態としては,微分の絶対値和です.さっき微分できないと言いましたが,積分して微分は出来て,すなわち微分して積分は出来るのです(超おおざっぱ).詳しくは参考文献とプログラムを見てください.

さて,このTVと相対を為すものは何でしょうか?TVは単調増加や単調減少の時にノルムが小さく,激しく増減をしているとノルムが大きくなります.この逆,つまり相対を為すのは,変化が少ないとノルムが大きく,振動が大きいほどノルムが小さくなれば言いわけです.これをGノルムと名付けましょう.概念としてはこれで良いのですが,数学的に定義するのは大変そうです.これも参考文献を見てください.

ここで,TVノルムとGノルムの重み付け和を最小にする,という課題を考えます.すると,画像をTVノルムの属する空間とGノルムの属する空間に分離できます.すなわち,骨格画像とテクスチャ画像に分離できちゃうのです.それも,数学に裏打ちされた結果として.

最後に歴史的なところを適当にかいつまんで列挙しておきます.そもそもの発想はWeveletから来ています.Mayerさんによる,Waveletの高周波成分がテクスチャ画像にちがいない,というアレゲなところから始まっています.それから,Rudinさん,Osherさん,Fatemiさんが頑張りまして,ROFモデルというのを提案しました.さらに,VesaさんとOsherさんが頑張りました.そのころにはWaveletはなりを潜めています.この辺の順序はちょっと未確認です(発表文献と年数を追っていけばいいですけど).しばらくして,(Waveletを使わない)真っ当な数値計算法が提案されてきました.でもまだ,オイラーラグランジュ法を使っている「古典的」な手法でした.例えば,0除算を避けるイプシロンが必要でした(私はこのモデルを今年の春頃実装してました.一応完成はしたけど収束が遅かった).彗星のごとくChambolleさんが登場し,非常にスマートにこの計算方法を確立してしまいました.PDE(Partial Differential Equation)なる理論を用います(が,私にはさっぱり理解できません).どうも,凸空間の非線形問題で,解が収束したり,不動点定理を利用したりと,お祭り騒ぎのようです(が,これもよく分からない.私は動けばいい,使えればいい,ぐらいなんで・・・).満を持して,Aujolさん登場.Aujolさん,他3名によるA2BCモデルが提案され,Chambolleさんの解法でスマートに計算される,というのが現在までのところです.

必須文献はこれだ!

どれも大学関係者じゃないと入手は難しいかも.しかし,ここにも(ほとんど)同じ内容がある!

というわけで,これさえ読めばあっという間に追いつけるわけです.

次に,実装に移りましょう.2番目と4番目に詳しく書いてあるので,それで良いと思います.今日公開するのはROFモデルだけです.A2BCはまた日を改めて実装したいと思います.方針ですが,C言語で純関数言語的に実装してあります.確保したメモリは呼び出し側で解放してください.まぁ,メモリも潤沢ですので,まじめに解放してませんが.

array_r1_tは画像と同じ構造です.array_r2_tはベクトル画像の構造です.なるべく構造を隠蔽するために,el関数とvec1, vec2関数を使って要素にアクセスします(どちらもdefine文ですが).イテレーションはメイン関数の中で実行しています.そして,ここが縮小写像を計算しているところになります.回数は文献によると30回〜40回ぐらいらしいです.

画像データは0〜1に正規化しています.パラメータのτは1/4が最大にして最速で収束する値らしいです.これは固定.λは超くせ者です.0〜1に正規化しているので,これの数分の1が良いようです.とりあえず,16/256=1/16にしてあります.ざっくりいって,輝度が1/16ぐらいで振動しているのがテクスチャに分類されます.逆に言うと輝度が1/16以上の不連続な部分は,エッジとして保存されます.

プログラムを投下します.

とりあえず,大急ぎでTVとその実装を見てきました.手元で動かすといろいろ問題点が出てきますな(笑).続きはまた時間が出来たときに(1週間ぐらいこの話に費やしてしまったよ・・・パラメータが分からなくて).