位相限定相関法で位置合わせ - ぬうぱんの備忘録

この記事は

　FFTってどうやるんだっけ？　のレベルから始まるチョー簡単な位相限定相関法のお勉強ノートです。基本的には参考にしたページをぺたぺたしてセクションごとのまとめだけ書いていく感じで進めます。

そもそも何をしたいのか

位置固定で３６０度パンする動画の各フレームを合成して全周パノラマ画像を作る
モデルは円筒で、貼りあわせ（モザイキング）は並進のみで行う（拡大と回転が無いと仮定するのでアフィン変換ですらない）
つまりは２画像の位置合わせができればいいのだが、折角なので結構精度がいいらしいPOCを使ってみたい

要するにPOCベースで画像の位置合わせをしたいってこった！

何を勉強するのか（ヤクの毛刈り）

必要な知識を列挙するとこんな感じ。

位相限定相関法(Phase Only Correlation, POC)
画像のフーリエ変換
離散フーリエ変換（主にFFT）

別にそんなに多くなかった。
じゃけんボトムアップで調べて行きましょうね。

離散フーリエ変換

　POCは画像のフーリエ変換結果を使うので画像のフーリエ変換が出来ないとそもそも話にならない。でもフーリエ変換って具体的にどうやってやるんだ…？　っていうことでお勉強しましょう。

フーリエ級数展開から画像のフーリエ変換までざっくり解説

フーリエ変換と画像圧縮の仕組み from yuichi takeda

www.slideshare.net

サンプルコード(C)

サンプルコード(AS)

実装の上で躓いたところとか重要なところを列挙すると

RGB画像にDFTをかけるとスペクトル画像が出てくる
DFTをかける時はRGB画像を複素数画像に移し替えてから
- 移し替える時はグレイスケール化して輝度値を実数部に格納。虚数部は０で初期化。
スペクトル画像を振幅の8bit画像として可視化する時は
- 複素数は絶対値(実部と虚部のL2ノルムのsqrt)をとってlogをとる(振幅のlogを可視化)
- 最大255になるようにスケーリング
- logを取らないと直流成分しか見えなくなる（中央に白い点がある画像が出てくる）
２次元のDFTは
- 行ごとに１次元のDFTしてから列ごとに１次元のDFTをすればOK
- 実際は「行ごとにDFT→画像を転置→行ごとにDFT→画像を転置」
- ここで言う転置とは行列の転置とおなじような操作
DFTをする前後で象限単位で点対称に画像を入れ替える
- 左上と右下を入れ替え、右上と左下を入れ替え
- DFTする画像の象限を入れ替えておく、DFTの結果も象限を入れ替える
- 入れ替えにより「画像中央が低周波数で外側が高周波数」になるようにする
- 画像左上を原点とすることによる弊害らしいが、よく分かっていない。
DFTの式は間違えやすい（実際間違えまくった）
- 数え上げなので書き込み先はゼロクリア必須
- 実部/虚部どっちなのか、sin/cosどっちなのか、符号はどっちなのか？
- ただし、sin/cosに与える角度は共通

結局、画像に対するDFTの手順は

画像の象限を入れ替え
画像の行ごとにDFTをかける
画像を転置
画像の行ごとにDFTをかける
画像を転置
画像の象限を入れ替え

ということになる。
IDFTをする時は１次元の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)}$ でそれを２次元でやるんだから正方形の画像だとすると ${\mathcal{O}(n^3)}$ ですかね？　重い。
で、高速化したかったらFFTとかで実装しろってことになるみたいです。でもなんか難しそうな事やってるし今回はスルーで。
ちなみにfftwというFFTのライブラリがあって、これは非常に高速な実装らしいです。なので、速度が大事なら素直にライブラリ使いましょう。

位相画像を得る

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

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

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

サンプルコード

位相限定相関法

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

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

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