この記事は
FFTってどうやるんだっけ? のレベルから始まるチョー簡単な位相限定相関法のお勉強ノートです。基本的には参考にしたページをぺたぺたしてセクションごとのまとめだけ書いていく感じで進めます。
そもそも何をしたいのか
- 位置固定で360度パンする動画の各フレームを合成して全周パノラマ画像を作る
- モデルは円筒で、貼りあわせ(モザイキング)は並進のみで行う(拡大と回転が無いと仮定するのでアフィン変換ですらない)
- つまりは2画像の位置合わせができればいいのだが、折角なので結構精度がいいらしいPOCを使ってみたい
要するにPOCベースで画像の位置合わせをしたいってこった!
離散フーリエ変換
POCは画像のフーリエ変換結果を使うので画像のフーリエ変換が出来ないとそもそも話にならない。でもフーリエ変換って具体的にどうやってやるんだ…? っていうことでお勉強しましょう。
実装の上で躓いたところとか重要なところを列挙すると
- RGB画像にDFTをかけるとスペクトル画像が出てくる
- DFTをかける時はRGB画像を複素数画像に移し替えてから
- 移し替える時はグレイスケール化して輝度値を実数部に格納。虚数部は0で初期化。
- スペクトル画像を振幅の8bit画像として可視化する時は
- 複素数は絶対値(実部と虚部のL2ノルムのsqrt)をとってlogをとる(振幅のlogを可視化)
- 最大255になるようにスケーリング
- logを取らないと直流成分しか見えなくなる(中央に白い点がある画像が出てくる)
- 2次元のDFTは
- 行ごとに1次元のDFTしてから列ごとに1次元のDFTをすればOK
- 実際は「行ごとにDFT→画像を転置→行ごとにDFT→画像を転置」
- ここで言う転置とは行列の転置とおなじような操作
- DFTをする前後で象限単位で点対称に画像を入れ替える
- 左上と右下を入れ替え、右上と左下を入れ替え
- DFTする画像の象限を入れ替えておく、DFTの結果も象限を入れ替える
- 入れ替えにより「画像中央が低周波数で外側が高周波数」になるようにする
- 画像左上を原点とすることによる弊害らしいが、よく分かっていない。
- DFTの式は間違えやすい(実際間違えまくった)
- 数え上げなので書き込み先はゼロクリア必須
- 実部/虚部どっちなのか、sin/cosどっちなのか、符号はどっちなのか?
- ただし、sin/cosに与える角度は共通
結局、画像に対するDFTの手順は
- 画像の象限を入れ替え
- 画像の行ごとにDFTをかける
- 画像を転置
- 画像の行ごとにDFTをかける
- 画像を転置
- 画像の象限を入れ替え
ということになる。
IDFTをする時は1次元のDFTの処理をIDFTにするだけで良い。上記の手順は全く同じ。
一次元のDFTは結局
数学の式で書くとこう
C++っぽいコードで書くとこう
const int N = xxx; //データサイズ const float s = 1.f; //IDFTの時は-1.f const float c = 1.f; //IDFTの時はN float xr[N], xi[N]; //入力データ.実部と虚部. float yr[N], yi[N]; //出力データ.実部と虚部. float t; //一時計算用.Theta. for(int i=0; i<size; ++i){ yr[i] = 0.f; yi[i] = 0.f; for(int j=0; j<size; ++j){ t = 2.f*M_PI*i*j/(float)N; yr[i] += xr[j]*cos(t) + s*xi[j]*sin(t); yi[i] += -s*xr[j]*sin(t) + xi[j]*cos(t); } yr[i] /= c; yi[i] /= c; }
位相画像を得る
さて、やっとDFT出来るようになったので、今度はPOCに使う位相画像(位相限定画像)を得る。
これはすっごく簡単で
- 8bit画像をDFTしてスペクトル画像化
- スペクトル画像を「ピクセルごとに振幅で正規化」する
参考にしたサンプルコードを見たら一発で分かるレベルの簡単さ。
可視化する時は実部だけを取り出してバイアスを加えると見やすいみたいです。
つまりは、各ピクセルの振幅情報を正規化によって消去するということなんですね。
サンプルコード
位相限定相関法
ようやっとメインディッシュです。
位相相関法の手順は
- 相関を取りたい2画像の位相画像(スペクトル)をそれぞれ生成
- 2枚の位相画像の相関をとる(まだスペクトル)
- 相関画像(スペクトル)をIDFT
- IDFTした相関画像中のピークを探索する
といった感じ。見つかったピークの位置がテンプレートの見つかった位置ということになる。
相関は実部と虚部でそれぞれ積をとるだけ良さ気なのだが、虚部の方は符号を反転させている。
でも、ピーク探索するときL1ノルムの大きいところを探索してるから符号関係無いですよね?
2015/08/06:追記
ピーク探索をするときはPOC画像の振幅の大きいピクセルを探せば良い模様。
つまり、実部と虚部のL2ノルムが大きいところを探せばいいってことですね。
参考にしたページのASのサンプルコードだと、実部と虚部の絶対値の和というノルムのようななにかで計算していますが、
いろいろ試した結果、L1ノルムかL2ノルムのどちらかで計算すれば良いようです。
むしろ、ノルムのような何かを使うと全然違う場所にピークがでたり、ピーク位置がズレたりで割とロクでもない結果になります。
実験結果
位置をずらしたドイツのトリを与えてみる。
入力画像1
入力画像2
入力画像1のDFT結果
POC結果
ピークが2個出てきてるんですが、これは符号の曖昧性ですかね?
反省
さて、完走した感想ですが(伝統芸能)
まず、DFTが死ぬほど遅い。本当にやってられないくらい遅い。
D言語マン曰く
次、今回は応用的に並進さえわかればいいし、精度もそれほど必要というわけでもないので、
POC画像をラスタスキャンしてピーク位置を探してますが、
サブピクセル精度の厳密な対応が欲しい場合はピークモデルを当てはめてより正確なピーク位置を推定するそうな。
次、細い所の理解がかなり怪しい。理解が浅いって勉強ノートとしてどうなんだ。
次、POCのピークを探索して得た位置ずれ量の符号の不定性を取り除く方法がわからない。今回のケースで言えば、点対称にピークが2つ出ているのだが、どっちを採用するべきなのかを判定する方法が不明。言い換えれば、左上方向にずれているのか、右下方向にずれているのかを判別する方法がわからない。
次、画像に窓関数をかけていない。普通はハニング窓をかけてからDFTするそうな。
次、画像サイズを無理やり2のべき乗数にしてFFTを適用できる。画像サイズが増えてでもそうしたほうが速い?
感想
かなり条件を制限しているとはいえ、POC自体はめちゃくちゃ簡単でびっくり。
むしろ、DFTの方が実装大変だった…。
パノラマ合成は
次の記事で!