読者です 読者をやめる 読者になる 読者になる

ぬうぱんの備忘録

技術系のメモとかいろいろ

作った曲一覧はこちら

位相限定相関法で位置合わせ

CV C++

この記事は

 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の手順は
  1. 画像の象限を入れ替え
  2. 画像の行ごとにDFTをかける
  3. 画像を転置
  4. 画像の行ごとにDFTをかける
  5. 画像を転置
  6. 画像の象限を入れ替え

ということになる。
IDFTをする時は1次元のDFTの処理をIDFTにするだけで良い。上記の手順は全く同じ。

一次元のDFTは結局

数学の式で書くとこう
{
X(k)
=
\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn}
}

{
e^{-j\frac{2\pi}{N}kn}
=
\cos(\frac{2\pi}{N}kn)
 -j \sin(\frac{2\pi}{N}kn)
}

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は超重い。参考にしたサンプルコードかナイーブな書き方をしてるからとかではなくて、アルゴリズム的にそもそも重たい。
一次元ですでに{\mathcal{O}(n^2)}でそれを2次元でやるんだから正方形の画像だとすると{\mathcal{O}(n^3)}ですかね? 重い。
で、高速化したかったらFFTとかで実装しろってことになるみたいです。でもなんか難しそうな事やってるし今回はスルーで。
ちなみにfftwというFFTのライブラリがあって、これは非常に高速な実装らしいです。なので、速度が大事なら素直にライブラリ使いましょう。

位相画像を得る

さて、やっとDFT出来るようになったので、今度はPOCに使う位相画像(位相限定画像)を得る。
これはすっごく簡単で

  1. 8bit画像をDFTしてスペクトル画像化
  2. スペクトル画像を「ピクセルごとに振幅で正規化」する

参考にしたサンプルコードを見たら一発で分かるレベルの簡単さ。
可視化する時は実部だけを取り出してバイアスを加えると見やすいみたいです。
つまりは、各ピクセルの振幅情報を正規化によって消去するということなんですね。

位相限定相関法

ようやっとメインディッシュです。
位相相関法の手順は

  1. 相関を取りたい2画像の位相画像(スペクトル)をそれぞれ生成
  2. 2枚の位相画像の相関をとる(まだスペクトル)
  3. 相関画像(スペクトル)をIDFT
  4. IDFTした相関画像中のピークを探索する

といった感じ。見つかったピークの位置がテンプレートの見つかった位置ということになる。
相関は実部と虚部でそれぞれ積をとるだけ良さ気なのだが、虚部の方は符号を反転させている。
でも、ピーク探索するときL1ノルムの大きいところを探索してるから符号関係無いですよね?
2015/08/06:追記
ピーク探索をするときはPOC画像の振幅の大きいピクセルを探せば良い模様。
つまり、実部と虚部のL2ノルムが大きいところを探せばいいってことですね。
参考にしたページのASのサンプルコードだと、実部と虚部の絶対値の和というノルムのようななにかで計算していますが、
いろいろ試した結果、L1ノルムかL2ノルムのどちらかで計算すれば良いようです。
むしろ、ノルムのような何かを使うと全然違う場所にピークがでたり、ピーク位置がズレたりで割とロクでもない結果になります。

実験結果

位置をずらしたドイツのトリを与えてみる。


入力画像1
f:id:NU_Pan:20150727015502p:plain
入力画像2
f:id:NU_Pan:20150727015534p:plain
入力画像1のDFT結果
f:id:NU_Pan:20150727015549p:plain
POC結果
f:id:NU_Pan:20150727015609p:plain
POC画像を探索した結果、画像1と画像2のズレは(32, 18)ピクセルと出てきた。ぴったり合ってる。
ピークが2個出てきてるんですが、これは符号の曖昧性ですかね?

反省

さて、完走した感想ですが(伝統芸能)
まず、DFTが死ぬほど遅い。本当にやってられないくらい遅い。
D言語マン曰く



らしいので素直にfftw使いましょう。

次、今回は応用的に並進さえわかればいいし、精度もそれほど必要というわけでもないので、
POC画像をラスタスキャンしてピーク位置を探してますが、
サブピクセル精度の厳密な対応が欲しい場合はピークモデルを当てはめてより正確なピーク位置を推定するそうな。

次、細い所の理解がかなり怪しい。理解が浅いって勉強ノートとしてどうなんだ。

次、POCのピークを探索して得た位置ずれ量の符号の不定性を取り除く方法がわからない。今回のケースで言えば、点対称にピークが2つ出ているのだが、どっちを採用するべきなのかを判定する方法が不明。言い換えれば、左上方向にずれているのか、右下方向にずれているのかを判別する方法がわからない。

次、画像に窓関数をかけていない。普通はハニング窓をかけてからDFTするそうな。

次、画像サイズを無理やり2のべき乗数にしてFFTを適用できる。画像サイズが増えてでもそうしたほうが速い?

感想

かなり条件を制限しているとはいえ、POC自体はめちゃくちゃ簡単でびっくり。
むしろ、DFTの方が実装大変だった…。

パノラマ合成は

次の記事で!