位相限定相関法における周波数領域の誤差に着目した
画像の高精度な平行移動量推定
High Accuracy Subpixel Displacement Estimation for Images by Minimizing
the Frequency-Domain Error of Phase-Only Correlation
岩田駿人
阿部正英
川又政征
東北大学大学院工学研究科電子工学専攻
Hayato IWATA
Masahide ABE
Masayuki KAWAMATA
Department of Electronic Engineering, Graduate School of Engineering, Tohoku University
アブストラクト 本稿では位相限定相関法における周波 数領域の誤差に着目し，勾配法を用いてこれを最小化す ることにより 2 つの画像の平行移動量を高精度に推定す る手法を提案する．提案法では，まず位相限定相関関数 のピーク座標を検出することでピクセル単位の平行移動 量を推定し，その値を起点として正規化クロスパワース ペクトルの理論値との誤差を用いて最急降下法を実行し， 誤差が最も小さくなるような推定値に収束させる．特に 雑音の多い環境において，提案法は従来の手法よりも高 い精度での推定が可能であること実験により示す． 1 はじめに 画像レジストレーション [1] とは 2 つの画像の位置合わ せを行う技術であり，画像処理・映像処理の様々な分野 において基礎となる重要な技術である．2 つの画像の対 応関係を推定するための手法として，近年位相限定相関 法 (Phase-Only Correlation, POC)[2] に注目が集まって いる． POC は高いロバスト性を持つと同時に，後述する相関 ピークへのフィッティングを行うことにより，画像のピク セル間隔よりも細かい精度 (サブピクセル精度) での推定 が可能という特徴がある．これらの特徴を生かし，動画 の動き推定や 3 次元計測など幅広い分野で応用される． 筆者らの研究グループでは，以前より劣化した映像の 位置ずれを修復するための技術として POC を用いた画像 の平行移動量推定に着目してきた．劣化した映像は様々 なノイズを含むことから高いロバスト性が要求され，ま た最終的に推定値を積算して処理を行うため高い精度で の推定が必要となる．先行研究において POC 関数を用い ることより高精度な位置ずれ修復が可能であることが示 されている [3]． しかし近年 POC 関数の統計的な性質について研究が行 われた結果，雑音のある環境での POC 関数の挙動が明ら かになってきた [4][5]．その結果を踏まえると POC 関数 の理論的な形状に基づく既存の手法では，雑音のある環 境において必ずしも正確な推定を行うことができないと 考えられる． そこで本論文では新たに周波数領域での理論値との誤 差に着目し，勾配法を用いてこれを最小化することによっ て 2 つの画像の平行移動量を推定する手法を提案する．提 案法ではまず POC 関数のピークを検出することによりピ クセル精度での平行移動量を推定し，その値を起点とし て最急降下法により周波数領域での誤差を最小とするよ う値に収束させることで平行移動量を推定する．特に雑 音の多い環境において，提案法は従来の手法よりも高い 精度での推定が可能であること実験により示す． 2 位相限定相関法を用いた画像の平行移動量推定 本章では位相限定相関法を用いて 2 つの画像の平行移 動量を推定する手法について解説する．2.1 節で位相限定 相関法の基礎について述べ，2.2 節および 2.3 節，2.4 節 において位相限定相関法を用いてサブピクセル精度での 平行移動量推定を行う既存の各手法について述べる． 2.1 位相限定相関法 位相限定相関法 [2] とは，入力された 2 つの画像に対し て周波数領域で振幅の正規化を行った上で相互相関を取 る手法である．2 つの画像 a および b の周波数スペクト ルをそれぞれ A(m, n), B(m, n) としたとき，2 つの画像 の正規化クロスパワースペクトルは式 (1) で表され，これ を逆フーリエ変換することにより POC 関数が得られる． R(m, n) = B(m, n)A ∗(m, n) |B(m, n)A(m, n)| (1) r(x, y) =F−1[R(m, n)] (2) 第31回信号処理シンポジウム 2016年11月8日∼10日（関西大学）
Image a Image b 図 1: POC による画像の平行移動量推定 ここで m, n は横方向縦方向それぞれの周波数インデック ス，A∗は A の複素共役，F−1は逆フーリエ変換をそれ ぞれ表している． 2 つの画像が循環シフトの関係にある時，POC 関数はシ フト量に相当する位置にピークを持つデルタ関数となる． r(x, y) = { 1, (x = x0, y = y0) 0, (otherwise) (3) ここで x0, y0は横方向縦方向それぞれの平行移動量である． 一般的に画像の平行移動は循環シフトとは異なるが，図 1 に示すように一般的な平行移動の関係においても POC 関 数はデルタ関数状の鋭いピークをもつ．よってこのピー ク座標を検出することによりピクセル精度で 2 つの画像 の平行移動量を推定することが可能である． 2.2 相関ピークの理論式を用いた平行移動量推定 2 つの画像が小数単位の平行移動量を持つとき，POC 関数は sinc 関数で近似できるような形状を持つことが知 られている [6]． r(x, y)≃ sinc(x − x0)sinc(y− y0) (4) この性質を利用し，POC 関数の系列に対して sinc 関数を フィッティングすることにより，x = x0, y = y0となる本 来のピーク座標を求め，サブピクセル精度で画像の平行移 動量を推定する手法 [7] が知られている．また，sinc 関数の 性質を利用したピーク評価式 (Peak Evaluation Formula, PEF) を用いることで，ピーク近傍の数点のサンプルから 高い精度で平行移動量を推定する手法 [8] が提案されてお り，高速な処理を可能としている．
しかし PEF は POC 関数が理想的に sinc 関数である 時のみに成立するものである．雑音のある環境において POC 関数はピークが減衰しサイドローブが上昇するとい う変動をする [5] ため，理想的な相関ピークの形状を前提 とした手法では正しく推定ができないと考えられる． 2.3 二次関数フィッティングを用いた平行移動量推定 PEF を用いた手法より雑音に強い方法として，POC 関 数を補間によって拡大し二次関数でフィッティングして ピーク座標を求める手法 (POC-QCFM)[3] がある．この 手法ではピーク近傍の 3 点のみを用いて平行移動量を推 定するため，PEF を用いた手法より雑音の影響を受けに くいことが明らかにされている [9]． しかし POC-QCFM の推定精度は POC 関数の拡大率 に依存し，平行移動量の値によっては近似による誤差が 生じる．また，空間領域でのフィッティングに基づいてい ることは変わりないため，雑音によって POC 関数のピー クに変動が生じた場合においては正しく推定できないこ とが予想される． 2.4 固有値分解を用いた平行移動量推定 上記 2 つの手法とは異なるアプローチで平行移動量を 推定する手法として，正規化クロスパワースペクトルか ら直接平行移動量を推定する手法 [10] がある．具体的に は特異値分解を用いて正規化クロスパワースペクトルを ランク 1 の行列で近似し，縦横それぞれの特異ベクトル の位相から最小二乗法によって平行移動量を推定する． しかし正規化クロスパワースペクトルから得られる位 相情報には 2π の任意性があるため，unwrap の処理が必 須である．雑音のある環境においてはこの unwrap の処理 がうまくいかない場合があることから，推定精度が著し く低下するという問題がある． 3 周波数領域の誤差に着目した新しい手法 本章では平行移動量を推定する新しいアプローチとし て周波数領域での理論値との誤差に着目し，勾配法を用い てこれを最小化することにより平行移動量を高精度に推 定する手法を提案する．3.1 節において推定値における理 想的な空間シフトと正規化クロスパワースペクトルとの 平均 2 乗誤差 (Mean Square Error, MSE) を導出し，3.2 節において最急降下法を用いてこれを最小化することで 平行移動量を推定する手法について述べる． 3.1 空間シフトの理論値との誤差 正規化クロスパワースペクトルは 2.1 節で述べたよう に式 (1) で表される．2 つの画像が循環シフトの関係にあ るとき，それらの正規化クロスパワースペクトルは周波 数領域における空間シフトと等しい． R(m, n) = e−j2π(x0Mm+y0Nn) (5) ここで x0, y0は横方向縦方向それぞれの平行移動量，M, N は横方向縦方向それぞれの画像サイズを表している．
Angle 正規化クロスパワースペクトル 理想的な空間シフトの系列 図 2: 正規化クロスパワースペクトルと空間シフト 図 2 に平行移動がある場合の正規化クロスパワースペク トルと理想的な空間シフトそれぞれの偏角を示す．画像端 の影響により位相がわずかに変動しているが，R(m, n) は おおよそ式 (5) に近い形となることが分かる．そこでこの R(m, n) と理想的な空間シフトの差に着目する．R(m, n) と推定値 (x0, y0) における理想的な空間シフトの MSE ε(x0, y0) を以下のように定義する． ε(x0, y0) = 1 M N ∑ m,n R(m, n) − e−j2π(x0 Mm+y0Nn) 2 (6) ここで ∑m,n は∑(Mm=−1)/2−(M−1)/2∑(Nn=−1)/2−(N−1)/2を表して いる．そしてこの ε(x0, y0) を最小とするような推定値 (˜x0, ˜y0) が正しい推定値であると考えられる． (˜x0, ˜y0) = arg min x0,y0 [ε(x0, y0)] (7) 誤差を最小化する手法を考えるために，まずこの誤差 の性質を解析する．ここでは正規化クロスパワースペク トルが理想的に空間シフトと等しいと仮定し，横方向お よび縦方向における正確な平行移動量と推定値の差をそ れぞれ ∆x = ˜x0− x0, ∆y = ˜y0− y0とする．正規化クロ スパワースペクトルと推定値における理想的な空間シフ トとの MSE は以下のように解くことができる． ε(∆x, ∆y) = 1 M N ∑ m,n 1 − e−j2π(∆x Mm+ ∆y Nn) 2 = 2 [ 1− 1 M N ∑ m,n cos { 2π ( ∆x M m + ∆y N n )}] = 2 1 −M N1 sin(π∆x) sin(π∆y) sin(π∆xM )sin ( π∆y N ) (8) なお，式 (8) は M および N が奇数の場合の式である．M および N が偶数の場合はナイキスト周波数である m = 図 3: R(m, n) と理想的な空間シフトとの MSE M/2, n = N/2 の部分が不連続点となるため、この部分を 除外し MSE は以下の式のようになる． ε(∆x, ∆y) = 2 [ 1− 1 (M− 1)(N − 1)
sin{(1−M1)π∆x}sin{(1−N1)π∆y}
sin(π∆x M ) sin ( π∆y N ) (9) M, N が奇数である場合の ε(∆x, ∆y) を図 3 に示す． M, N が偶数である場合についても同様の形状となる． 図より−1 ≤ ∆x ≤ 1, −1 ≤ ∆y ≤ 1 の範囲におい て ε(∆x, ∆y) は単峰性の形状をしており，その頂点が (∆x, ∆y) = (0, 0)，すなわち正しい推定値であることが 分かる． 3.2 最急降下法を用いた誤差の最小化 誤差を最小化する問題を考えるとき，一般的には線形 最小 2 乗法などのアルゴリズムが理想的である．しかしこ の問題においては，前節において求めたように MSE が非 線形の式によって表されるため，R(m, n) から直接計算に よって平行移動量の推定値 ˜x0, ˜y0を求めることは難しい．
正規化 クロスパワー スペクトル IDFT 誤差の計算 最急 降下法 ピーク 検出 推定値 振幅正規化 振幅正規化 DFT DFT a b 図 4: 最急降下法による画像の平行移動量推定 よって勾配法などの繰り返しの演算によって正しい推 定値を探索する手法を考える．前節において−1 ≤ ∆x ≤ 1,−1 ≤ ∆y ≤ 1 の範囲において MSE が単峰性の形状を 持つことが明らかになったので，正しい推定値が±1[pixel] の範囲に存在することが分かれば勾配法によって収束さ せることが可能である．そこで POC 関数のピーク座標検 出を併用することにより平行移動量の推定を行う． 図 4 に提案法の概要を示す．提案法ではまず正規化ク ロスパワースペクトルを IDFT することにより POC 関数 を求め，そのピーク座標を検出することによりピクセル精 度の平行移動量推定値を求める．次に求めたピクセル精度 の推定値を起点として，R(m, n) と推定値 (˜x0(k), ˜y0(k)) における理想的な空間シフトとの誤差を用いて最急降下 法を実行し，誤差が最も小さくなるような推定値に収束 させる．提案法の最急降下法における推定値の更新式は 以下の式のようになる． ( ˜ x0(k + 1) ˜ y0(k + 1) ) = ( ˜ x0(k) ˜ y0(k) ) − µ∇ε(˜x0(k), ˜y0(k)) (10) ここで ˜x0(k), ˜y0(k) は横方向縦方向の k 回目の更新にお ける推定値を表しており，µ は 1 回あたりの更新量を調整 するステップサイズパラメーターである． 4 実験と比較 本章では，前章で提案した最急降下法を用いる画像の 平行移動量推定の手法について画像を用いた実験を行う． 4.1 節において実際に最急降下法を用いることで正しい推 定値に収束するかどうかを実験し，4.2 節において画像に 雑音がある環境において既存法および提案法の推定精度 を比較する．
(a) Lenna (512×512) (b) Peppers (512×512)
図 5: 実験に使用した画像 図 6: 最急降下法による推定値の収束 4.1 提案法における推定値の収束 提案法を用いて画像の平行移動量を推定を行うために， まず人工的に位置ずれを発生させた画像を用意した．画像 は図 5 に示す 512 × 512[pixel] の Lenna を用い，フーリ エ変換を用いて横方向に x0= 1.3[pixel]，縦方向に y0 = 0.6[pixel] の平行移動を発生させた．また，画像端の影響 を低減させるために Hanning 窓による窓かけを行った． 初めに 2 つの画像から POC 関数を求めたところ， (x, y) = (1, 1) の位置にピークが立つ系列が得られた．次 に得られたピーク座標を元に推定値の初期値を ˜x0(0) = 1, ˜y0(0) = 1 とし，式 (10) に基づいてステップサイズ µ = 0.1 として最急降下法を実行したところ図 6 に示 すように正しい推定値 ˜x0 = 1.3, ˜y0 = 0.6 に収束す ることが確認された．15 回目の更新において推定値は ˜ x0(15) = 1.2997, ˜y0(15) = 0.5998 となり，雑音のない 環境においては 1/100[pixel] を下回る非常に高い精度で の推定が可能であることが分かった． 4.2 推定精度の比較 雑音のある環境での提案法および既存の手法について， 推定精度を比較する実験を行った．実験には図 5 に示す 2 種類の画像を用いた．まず平行移動量あたりの推定精度 を比較するために，x0 =−2 ∼ 2, y0=−2 ∼ 2 の範囲で
表 1: シミュレーション諸元 PEF based フィッティング点数 p = 6 + 6 スペクトル重み幅 U1= U2= 255 POC-QCFM 拡大率 q = 4 Proposed 更新回数 K = 15 ステップサイズ µ = 0.1∼ 0.6 図 7: シフト量あたりの推定誤差 人工的に平行移動を発生させた画像を用意し，片方の画 像に分散 σ2 η= 0.001 の白色ガウス雑音を付加した．なお， 水平方向の平行移動を発生させたときは垂直方向の移動 量は y0= 0，垂直方向の平行移動を発生させた時は水平 方向の移動量は x0= 0 とした．生成した二つの画像につ いて，表 1 に示した諸元に基づき平行移動量の推定を行 い，正確な値との誤差を算出した．付加雑音を変えて各 1000 回のシミュレーションを行い，誤差の平均と標準偏 差を算出した．画像 Lenna を用いた場合の結果を図 7 に 示す．画像 Peppers を用いた場合もほぼ同様の傾向が見 られたため結果は省略する． PEF および QCFM を用いた従来の手法における推定 値は平行移動量の値に対して一定の傾向を持っているが， 提案法における推定値は正確な値を中心とした分布をし ていることが分かる．提案法は周波数領域の誤差を最小 表 2: 実験に使用したステップサイズ σ2 η ≤ 10−4 ≤ 10−3 2× 10−3 5× 10−3 10−2 µ 0.1 0.2 0.3 0.4 0.6 (a) 画像Lennaを用いた場合 (b) 画像Peppersを用いた場合 図 8: 雑音分散あたりの推定誤差 化することで，元の画像に付加された雑音の影響を最小 化できているためだと考えられる． 次に，より一般的な環境を想定し平行移動を x0=−2 ∼ 2, y0 =−2 ∼ 2 の範囲でそれぞれ一様分布でランダムに 発生させ，生成した画像のうち片方に σ2 η = 10−5 ∼ 10−2 の白色ガウス雑音を付加した．シミュレーションをそれ ぞれの雑音分散の値に対して 10000 回行い，推定値の平 均 2 乗誤差の平方根 (Root Mean Square Error, RMSE) を算出した．なおステップサイズ µ は，事前の実験にお いて正しく収束することを確認した表 2 に示すものを用 いた。結果を図 8 に示す． 提案法は従来法と比較して雑音を付加した場合の推定 精度が総合的に向上していることが分かる．特に雑音分 散 σ2 η = 10−2 という厳しい環境においても 1/100[pixel]
程度の精度での推定が可能であり，劣化した映像の修復 において効果的に高い精度の推定を行うことが可能であ ると考えられる． 5 まとめ 本稿では周波数領域における正規化クロスパワースペ クトルと理論値との MSE に着目し，最急降下法を用いて これを最小化することでサブピクセル精度で画像の平行 移動量を推定する手法を提案した．また実験によって特 に雑音の多い環境下において従来の手法よりも高い精度 での推定が可能であることを示した． 今後の検討課題としては，まず最急降下法の探索にお ける終了条件の検討があげられる．今回の実験では一定 回数の更新で収束したものと見なしたが，実際の応用に おいては勾配の値などの条件によって探索を終了させる 手法を考える必要がある．また，雑音のある環境におい ては MSE の底となる部分が上昇することから全体として MSE の勾配が緩くなるため，収束に時間がかかってしま うことが分かっている．これを解決するために，POC 関 数のピークの値などのパラメーター利用してステップサ イズ µ を適応的に決めてやる必要があると考えられる． 参考文献
[1] B. Zitova and J. Flusser, “Image registration meth-ods: A survey,” Image and vision computing, vol.21, no.11, pp.977–1000, Oct. 2003.
[2] C. Kuglin, “The phase correlation image alignment method,” Proc. International Conference on Cyber-netics and Society, 1975, pp.163–165, Jan. 1975.
[3] M. Hagiwara, M. Abe, and M. Kawamata, “Estima-tion method of frame displacement for old films us-ing phase-only correlation,” Journal of Signal Pro-cessing, vol.8, no.5, pp.421–429, Sep. 2004.
[4] S. Yamaki, J. Odagiri, M. Abe, and M. Kawamata, “Effects of stochastic phase spectrum differences on phase-only correlation functions part I: Statis-tically constant phase spectrum differences for fre-quency indices,” 2012 3rd IEEE International ference on Network Infrastructure and Digital Con-tent, pp.360–364, Sept. 2012. [5] 福井一弘，八巻俊輔，阿部正英，川又政征，“白色ガウ ス雑音に起因する位相差の変動を持つ実信号の位相 限定相関関数の統計的解析，” 2015 年電子情報通信学 会基礎・境界ソサイエティ大会講演論文集，no.A-4-2， p.61，9 月 2015．
[6] H. Foroosh, J.B. Zerubia, and M. Berthod, “Ex-tension of phase correlation to subpixel registra-tion,” IEEE Trans. on Image Processing, vol.11, no.3, pp.188–200, Mar. 2002.
[7] K. Takita, Y. Sasaki, T. Higuchi, and K. Kobayashi, “High-accuracy subpixel image registration based on phase-only correlation,” IEICE Trans. on Funda-mentals of Electronics, Communications and Com-puter Sciences, vol.86, no.8, pp.1925–1934, Aug. 2003.
[8] S. Nagashima, T. Aoki, T. Higuchi, and K. Kobayashi, “A subpixel image matching technique using phase-only correlation,” 2006 International Symposium on Intelligent Signal Processing and Communications, pp.701–704, Dec. 2006. [9] 松本圭右，八巻俊輔，阿部正英，川又政征，“位相限 定相関を用いた平行移動量推定における加法性雑音 の影響の評価，” 2015 年電子情報通信学会基礎・境界 ソサイエティ大会講演論文集，no.A-4-22，p.89，3 月 2013．
[10] W.S. Hoge, “A subspace identification extension to the phase correlation method,” IEEE Trans. on Medical Imaging, vol.22, no.2, pp.277–280, Feb. 2003.
