画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法
【課題】画像間対応パラメータを少ない計算量で安定且つ高速に高精度に同時推定できるようにした、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を提供する。
【解決手段】離散的な位置で得られた画像間のＮ次元類似度値がある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められたサブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、Ｎ個のＮ次元超平面の交点を求めるステップと、求められた交点をＮ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有する。
【解決手段】離散的な位置で得られた画像間のＮ次元類似度値がある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められたサブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、Ｎ個のＮ次元超平面の交点を求めるステップと、求められた交点をＮ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、例えば画像による位置計測、リモートセンシング、航空写真を利用した地図作製、ステレオビジョン、画像張り合わせ（モザイキング）、動画像位置合わせ、３次元モデリングなど多数の画像処理分野で使用できる、領域ベースマッチングを用いた推定方法に関し、特に、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法に関するものである。
【背景技術】
【０００２】
ステレオ画像処理をはじめとして、トラッキング、画像計測、リモートセンシング、画像レジストレーションなどの多くの画像処理分野では、画像間変位を得るために、領域ベースマッチングが基本的な処理として採用されている。
【０００３】
領域ベースマッチングは、注目領域のサイズや形状を任意に設定できること、注目画素に対して注目領域をオフセットできること、計算が直感的で直接的な実装が可能であることなどが特徴である。領域ベースマッチングで画像間の変位をサブピクセルで求めるときには、画素単位で離散的に得られる類似度評価値を利用してサブピクセル変位を推定する方法が一般的に利用されている。
【０００４】
画像間の対応を求めるためには、類似度と呼ぶ評価値を利用する。類似度とは、画像間の相対的な位置関係を入力、その位置関係のときの画像間の類似度を出力とする関数である。類似度は、画素単位の離散的な位置で値が決まるので、類似度だけをもとに画像間の対応位置関係を求めると、画素単位の分解能に制限される。そこで、類似度値をもとに補間することによって、サブピクセル分解能の画像間の対応を求める。
【０００５】
従来、位置分解能を拡大するために、類似度評価値をフィッティング内挿することでサブピクセル推定を行うことが多い。しかし、画像の平行移動量をサブピクセル分解能で求める際に、画像の水平方向と垂直方向を独立にフィッティング内挿していたため、十分な精度が得られなかったという問題点があった。
【０００６】
要するに、従来、デジタル画像を使った領域ベースのマッチングでは、変位の推定分解能を向上させるために、画像の水平／垂直方向について独立に類似度評価値を使ったサブピクセル推定を行っていた。
【０００７】
また、従来、濃度こう配法を利用して画像間の変位を求めると、はじめからサブピクセル変位量を得ることができる。濃度こう配法では、水平方向と垂直方向を合わせて扱っている。ただし、画像間変位が１画素以上のときには画像を縮小する必要がある。通常は、画像スケールの変更と同時に実装する。濃度こう配法の実装は繰り返し計算になるので、計算時間の見積もりや、ハードウェアへの実装が困難であるという欠点があった。
【０００８】
さらに、画像をＤＦＴしてから複素共役積をもとめて逆ＤＦＴする従来手法では、注目領域のサイズが２ｎでなければならず、各種テクニックが必要である。流体計測分野では多く用いられ、計算量を小さくできることが特徴である。しかし、ビジョン分野では、ほとんど利用されない。しかも、計測対象が平面であるという仮定を用いているので、凹凸がある計測対象には適用が困難であるという欠点があった。
【０００９】
従来は、以下のように詳細に説明されるように、画像の水平方向と垂直方向を独立に補間推定していた。しかし、従来のこのような方法によると、図１に示すように、推定誤差が発生するという問題点があった。
【００１０】
ここでは、画像を水平方向と垂直方向に独立と考え、変位のサブピクセル推定も独立に行う従来方法を「１次元推定方法」と称する。
【００１１】
従来の１次元推定方法の問題点を説明するために、まず、画像間の水平方向変位を求める問題を考える。真の画像間変位を(ｄｓ,ｄｔ)、異方性を持った類似度の回転角度をθｇとする。１次元推定方法では、離散化類似度として、Ｒ(−１,０)、Ｒ(０,０)、Ｒ(１,０)(図１の□)を使い、下記数１によって
を推定した。なお、図１は、従来の１次元サブピクセル推定方法を説明するための図である。図中の
は、ｓ＝−１,０,１での３つの類似度の値を用いた従来の１次元サブピクセル推定方法によって推定された位置で、ｄｔ≠０及びθｇ≠０の際に、この推定値は真のピーク位置に対して誤差を有する。
【数１】
【００１２】
仮に上記数１で直線ｔ＝０上の類似度最大位置を正しく推定できたとしても、図１から明らかなように、真の水平方向画像間変位ｄｓ(図１の▲の水平成分)に対して大きな推定誤差を含んでいる。すなわち、次の条件の全てが真のときには、水平方向推定誤差が発生して、
となる。
・垂直方向変位ｄｔ≠０。
・２次元類似度が異方性を持つ。
・異方性を持つ２次元類似度の回転角度θｇ≠０。
大部分の画像が上記条件に当てはまる。
【００１３】
また、垂直方向に関しても同様である。例えば、図２(ａ)に示す画像中のコーナー領域のＳＳＤ自己類似度を求めると、図２(ｂ)のようになっている。この自己類似度は異方性を持ち、θｇ≠０なので、従来の１次元推定方法では、サブピクセル推定誤差が発生する可能性があることを示している。また、図３(ａ)に示すテクスチャ画像を使ったサブピクセル推定においても、推定誤差が発生する可能性を示している。
【００１４】
また、コンピュータビジョンの分野によっては、エピポーラ拘束などの拘束条件を利用することで探索範囲を１次元に限定できるので、１次元信号の高精度なマッチングができれば十分な場合がある。しかしながら、２次元画像を使って、高精度な２次元変位を求める必要がある用途も多い。例えば、動き推定、動き推定による領域分割、ターゲットトラッキング、モザイキング、リモートセンシング、超解像などである。
【００１５】
画像の領域ベースマッチングにおいて、画像間変位が平行移動だけで表現できるときには、従来から多くの手法が提案され、実際に用いられてきた。たとえば、代表的な手法として、領域ベースマッチングと類似度補間手法を組み合わせた、サブピクセル推定手法は、次のようにして画像間のサブピクセル変位を推定する。
【００１６】
対応位置を求める画像を、ｆ(ｉ,ｊ)とｇ(ｉ,ｊ)としたとき、整数単位の変位(ｔｘ,ｔｙ)に対する画像間のＳＡＤを、次のように計算する。
【数２】
ただし、ＳＡＤは輝度差の絶対値の和(Sum of Absolute Differences)を表し、Ｗは対応を求める注目領域を表す。ＳＡＤは、画像が完全に一致したときに０になり、不一致が増えるに従って大きな値をとる。ＳＡＤは類似していない評価値を表すので、非類似度である。
【００１７】
ＳＡＤの代わりに、次のＳＳＤを用いることも多い。
【数３】
ただし、ＳＳＤは輝度差の絶対値の２乗和(Sum of Squared Differences)を表す。ＳＳＤも画像が完全に一致したときに０になり、不一致が増えるに従って大きな値をとる。ＳＳＤも類似していない評価値を表すので、非類似度である。
【００１８】
ＳＡＤやＳＳＤの代わりに、次の相関係数(ＺＮＣＣ:Zero-mean Normalized Cross-Correlation)を用いることもある。ＺＮＣＣは、注目領域内の画素濃度を統計データと見なし、注目領域間の統計的な相関値を計算している。
【数４】
ただし、
は、それぞれの領域内の濃度平均である。ＲＺＮＣＣは−１から１までの値をとり、画像が完全に一致したときに１になり、不一致が増えるに従って小さな値をとる。ＺＮＣＣの特徴は、画像間に濃度やコントラストの変化があってもほとんど同じ類似度が得られることである。
【００１９】
ＺＮＣＣは類似性の評価値を表すので類似度とよばれるが、ＳＡＤやＳＳＤは非類似度を表す。ここでは、説明を簡単にするために、以後はＳＡＤやＳＳＤ、ＺＮＣＣを「類似度」と統一して表記する。
【００２０】
類似度が最小値(ＳＡＤやＳＳＤのとき)または最大値(ＺＮＣＣのとき)になる整数単位の変位
を探索することで、画像間の整数単位変位を得ることができる。整数単位の変位を得た後で、次のようにサブピクセル変位を推定する。
【数５】
上記数５を用いて、最終的に得る画像間サブピクセル変位は、
となる。
【００２１】
このような従来手法では、画像間に回転やスケール変化があるときに正しい対応を得ることができない問題があった。画像間の回転角度がわずかでも、２つの画像が異なるものと判断され、対応位置が見つからないことや、誤った位置に対応を検出する問題があった。
【００２２】
一方、画像間の対応が平行移動だけでは表現できない場合、つまり、画像間に回転やスケール変化があるときには、回転やスケール変化も同時に推定する必要がある。従来、このときには、補間処理を用いて片方の画像(テンプレート画像)を平行移動と回転やスケール変化をさせてからマッチングを行い、類似度を計算しながら各パラメータを変化させ、類似度が最小値または最大値になるパラメータを探索する必要があった。次に、このパラメータ探索法を説明する。
【００２３】
ｆ(ｉ,ｊ)をテンプレート画像、ｇ(ｉ,ｊ)を入力画像とするとき、平行移動(ｔｘ,ｔｙ)および回転角度θとスケール変化ｓをパラメータとする類似度ＲＳＳＤを次のように計算する。
【数６】
ただし、
は、平行移動(ｔｘ,ｔｙ)および回転角度θとスケール変化ｓを与えたときの、位置(ｉ,ｊ)における画像ｇ(ｉ,ｊ)の補間画像を表す。補間画像を作るときには、双線形補間(Bi-Linear補間)や双３次補間(Bi-Cubic補間)などを利用する。また、ＳＳＤではなく、ＳＡＤやＺＮＣＣを利用してもよい。
【００２４】
何らかの方法で初期パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜０＞を求めてＲＳＳＤ(ｔｘ,ｔｙ,θ,ｓ)＜０＞を計算し、より小さなＲＳＳＤを与える更新パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜１＞を探索する。パラメータを更新しても結果に変化がなくなるまで、更新を続ける。このとき、Newton-Raphson法やSteepest (Gradent) Descent法(最急降下法)、Levenberg-Marquadt 法など数値的解法を用いる。
【００２５】
初期パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜０＞は、たとえば、画像間の対応が平行移動だけで表現できると仮定して通常の領域ベースマッチングによって推定した
を用いた
などを利用する。
【非特許文献１】清水雅夫・奥富正敏,「画像のマッチングにおける高精度なサブピクセル推定手法」,電子情報通信学会論文誌D-II,2001年7月,第J84-D-II巻,第7号,p.1409-1418
【非特許文献２】清水雅夫・奥富正敏,「画像のマッチングにおける高精度なサブピクセル推定の意味と性質」,電子情報通信学会論文誌D-II,2002年12月,第J85-D-II巻,第12号,p.1791-1800
【非特許文献３】マサオ・シミズ、マサトシ・オクトミ(Masao Shimizu and Masatoshi Okutomi),「プリサイス サブピクセル エスティメイション オン エリアベイセド マッチング(Precise Sub-Pixel Estimation on Area-Based Matching)」,プロク. 8th IEEE インターナショナル コンファレンス オン コンピュータ ビジョン (ICCV2001)（Proc. 8th IEEE International Conference on Computer Vision (ICCV2001)）,(カナダ,バンクーバー),2001年7月,p.90-97
【非特許文献４】マサオ・シミズ、マサトシ・オクトミ(Masao Shimizu and Masatoshi Okutomi),「アン アナリシス オフ サブピクセル エスティメイション エラー オン エリアベイセド イメージ マッチング(An Analysis of Sub-Pixel Estimation Error on Area-Based Image Matching)」,プロク. 14th インターナショナル コンファレンス オン デジタル シグナル プロセシング(DSP2002)（Proc. 14th International Conference on Digital Signal Processing (DSP2002)）,(ギリシア,サントリーニ),2002年7月,第II巻,p.1239-1242(W3B.4)
【非特許文献５】Ｃ.・クエンティン・デイビス、ゾハー・Ｚ．・カル、デニス・Ｍ.・フリーマン(C. Quentin Davis, Zoher Z. Karu and Dennis M. Freeman),「イクイバレンス オフ サブピクセル モーション エスティメイタス ベイセド オン オプティカル フロー アンド ブロック マッチング(Equivalence of subpixel motion estimators based on optical flow and block matching)」,IEEE インターナショナル シンポジウム フォー コンピュータ ビジョン1995（IEEE International Symposium for Computer Vision 1995）,(米国フロリダ州),1995年11月,p.7-12
【非特許文献６】サムソン・Ｊ.・ティモネー、デニス・Ｍ.・フリーマン(Samson J. Timoner and Dennis M. Freeman),「マルチイメージ グラジエントベイセド アルゴリズム フォー モーション エスティメイション(Multi-image gradient-based algorithms for motion estimation)」,オプティカル エンジニアリング（Optical Engineering）,(米国),2001年9月,第40巻,第9号,p.2003-2016
【非特許文献７】ジョナサン・Ｗ.・ブラント (Jonathan W. Brandt),「インプルーブド アキュラシ イン グラジエントベイセド オプティカル フロー エスティメイション(Improved Accuracy in gradient-based Optical Flow Estimation)」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International Journal of Computer Vision）,(米国),1997年,第25巻,第1号,p.5-22
【非特許文献８】Ｑ.・テン、Ｍ．・Ｎ．・ヒューニッス(Q. Tian and M. N. Huhns),「アルゴリズム フォー サブピクセル レジストレーション(Algorithms for Subpixel Registration)」,コンピュータ ビジョン,グラフィックス アンド イメージ プロセシング（Computer Vision, Graphics and Image Processing）,(米国),1986年,第35号,p.220-233
【非特許文献９】ショーン・ボーマン、 マーク・Ａ.・ロバートソン、ロバート・Ｌ.・スティブンソン(Sean Borman, Mark A. Robertson and Robert L. Stevenson),「ブロックマッチング サブピクセル モーション エスティメイション フロム ノイジー,アンダー―サンプルド フレームス---アン エンパイリカル パフォーマンス エバリュエーション(Block-Matching Sub-Pixel Motion Estimation from Noisy, Under-Sampled Frames---An Empirical Performance Evaluation)」,SPIE ビジュアル コミューニケイションズ アンド イメージ プロセシング 1999（SPIE Visual Communications and Image Processing 1999）,(米国カリフォルニア州)
【発明の開示】
【発明が解決しようとする課題】
【００２６】
このような多パラメータ最適化では、初期値が適切でないと正しい結果が得られない問題がある。また、画像パターンやノイズ、撮影レンズの光学歪みなどの影響で解が得られない問題もある。さらに、繰り返し計算を利用するので、計算時間を見積もることができないといった難点がある。このため、アルゴリズムを実装する上で、完成したシステムの総合応答時間を正確に見積もることができず、ハードウェア化することがきわめて困難であった。また、計算時間の面では、繰り返し計算の各段階で画像を補間する必要があるが、画像補間演算は演算量が多いため、多くの計算時間が必要になる。さらに、推定精度の面では、画像補間手法によって補間画像の性質が異なるので、採用する補間手法によって最終的なパラメータ推定精度や推定誤差の性質が異なる点も大きな問題であった。
【００２７】
本発明は上述のような事情よりなされたものであり、本発明の目的は、Ｎ次元類似度を考慮することで、上記従来手法の問題点を解決し、画像間対応パラメータを少ない計算量で安定且つ高速に高精度に同時推定できるようにした、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を提供することにある。
【課題を解決するための手段】
【００２８】
本発明は、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法に関し、本発明の上記目的は、Ｎ（Ｎは４以上の整数である）個の画像間対応パラメータを軸とするＮ次元類似度空間において、離散的な位置で得られた画像間のＮ次元類似度値を利用して、前記Ｎ個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法であって、前記画像間のＮ次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、前記Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、前記Ｎ個のＮ次元超平面の交点を求めるステップと、前記交点を前記Ｎ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有するようにすることによって効果的に達成される。
【００２９】
また、本発明は、画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法に関し、本発明の上記目的は、３個の画像間対応パラメータを軸とする３次元類似度空間において、離散的な位置で得られた画像間の３次元類似度値を利用して、前記３個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法であって、前記画像間の３次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する平面を求めるステップと、前記平面を各パラメータ軸に対して３個求めるステップと、前記３個の平面の交点を求めるステップと、前記交点を前記３次元類似度空間における３次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有するようにすることによって効果的に達成される。
【００３０】
更に、本発明は、画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法に関し、本発明の上記目的は、離散的に得られた画像間の２次元類似度値を利用して、連続領域における２次元類似度最大値または最小値の位置を高精度に推定する画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法であって、前記画像間の２次元類似度値が、水平軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（水平極値線ＨＥＬ）を求めるステップと、前記画像間の２次元類似度値が、垂直軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（垂直極値線ＶＥＬ）を求めるステップと、前記水平極値線ＨＥＬと前記垂直極値線ＶＥＬの交点を求めるステップと、前記交点を前記２次元類似度の最大値または最小値を与えるサブピクセル推定位置とするステップとを有するようにすることによって効果的に達成される。
【発明の効果】
【００３１】
本発明によれば、画像間の平行移動、回転、スケール変化を同時に高精度に推定することができる。即ち、本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法によれば、画像間の対応に、平行移動があるときだけでなく、回転やスケール変化があるときにも、繰り返し計算を利用することなく、高精度にこれらの画像間対応パラメータを同時に推定することができる。本発明では、繰り返し計算を利用しないので、従来のように計算ごとに画像を補間する必要が全くない。
【００３２】
また、本発明によれば、テンプレート画像を用いて少ない枚数の補間画像を作成しておき、これらの補間画像を用いた画像間の類似度値を利用して、平行移動と回転やスケール変化を同時に高精度に推定することができる。本発明は、一般的に利用されているようなパラメータ最適化手法と異なり、繰り返し計算を用いないため、計算が容易でしかも計算時間が予測可能で、繰り返し計算を行うときに問題となる収束性を保証する必要もなく、初期値による不安定性も生じない。つまり、本発明によれば、従来方法において生じる収束性の問題や初期値の問題を排除することができる。
【発明を実施するための最良の形態】
【００３３】
以下、本発明の好適な実施例を図面及び数式を参照しながら説明する。
【００３４】
実施例１：
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の実施例１では、画像のサンプリング周期で求まっている２次元類似度を使って、その２次元類似度の最大値又は最小値を与える２次元サブピクセル変位（つまり、水平及び垂直方向の画像間変位）を高精度に同時に推定するようにしている。
【００３５】
ここで、画像を水平方向と垂直方向に独立と考え、変位のサブピクセル推定も独立に行う従来方法を「１次元推定方法」と称するのに対して、本発明の実施例１に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法は、２次元画像から計算した離散位置における２次元類似度を使って、２次元サブピクセル推定を行なっているので、「画像のサブピクセルマッチングにおける高精度２次元推定方法」又は「２次元サブピクセル推定方法」と称し、更に、略して「２次元推定方法」と呼ぶようにしている。
【００３６】
後述する本発明の実施例１に係る２次元推定方法は、領域ベースマッチングを用いた高精度同時推定方法で、従来の「１次元推定方法」と比較して計算量がわずかに増えるだけで、圧倒的に高精度な２次元サブピクセル推定が可能である。
【００３７】
（１）２次元画像モデルと２次元類似度モデル
ここでは、何種類かの画像モデルをとりあげ、その２次元類似度を考える。異なる画像モデルから得られる２次元類似度が、同じ２次元類似度モデルで近似できることを示す。
【００３８】
領域ベースマッチングとは、２枚の画像の類似度を評価し、その最大または最小位置を２枚の画像間の変位として求めることである。この類似度は、通常は画像のサンプリング周期と一致して、とびとびの変位位置に対して値が得られる。画素単位でのマッチングは、これらの類似度の中から、最大値または最小値を探すことに対応する。サブピクセル推定は、これらの類似度を補間して、サブピクセルでの最大値または最小値を与えるサブピクセル変位位置を推定することである。以後の説明では、「最大または最小」を単に「最大」と表現する。ＳＡＤやＳＳＤ類似度を使うときには「最小」、相関類似度を使うときには「最大」を意味する。
【００３９】
（１―１）コーナーモデル
まず、領域ベースマッチングとサブピクセル推定が可能な２次元画像モデルとして、１次元画像モデル
【数７】
を単純に拡張して、下記数８で示すコーナー画像を考える。ここで、σimageは濃度エッジのシャープさである。このコーナー画像を図４(ａ)に示す。
【数８】
ただし、画像座標系をｕ−ｖ、類似度座標系（変位座標系）をｓ−ｔとする。
【００４０】
上記数８で表す画像を参照画像として、また、全く同じ画像を入力画像として２次元マッチングを行ったときのＳＳＤ類似度(本発明では、ＳＳＤ類似度を「自己類似度」と称することもある)は、下記数９で求めることができる。ＳＳＤ類似度は、値が小さいほど類似していることを示すので、この場合には暗い場所ほど類似性が高いことを示している。自己類似度は全く同じ画像を使った類似度なので、変位(ｓ,ｔ)＝(０,０)の位置の類似度が最も小さい。この類似度を図４(ｂ)に示す。
【数９】
ただし、総和範囲は、任意形状の注目範囲であるが、図４(ｂ)では１１×１１の矩形領域で計算した。
【００４１】
図４(ａ)では、コーナーに剪断成分が含まれない、つまり水平方向の濃度変化が垂直位置に依存しないので、２次元画像の水平方向と垂直方向を独立に扱うことができる。このとき、上記数９の類似度も、水平方向と垂直方向を独立に扱うことができる。従って、サブピクセル変位位置を推定するときにも、水平方向と垂直方向を独立に扱うことができる。
【００４２】
しかし、実際の画像では、必ずしも９０度のコーナーで構成されるわけではない。例えば、建築物を地上から撮影したステレオ画像を使って３次元情報を再構成するときには、建築物のコーナー領域を、他方の画像との対応を求めるが、コーナーが９０度に撮影されることはまれである（図２(ａ)参照）。
【００４３】
そこで、ｖ＝０に濃度エッジを持つ２枚の２次元画像
【数１０】
を、それぞれ左にθａ、右にθｂだけ回転した画像を考える。
【数１１】
【００４４】
このθａ、θｂを画像全体の回転角度θｇとコーナー角度θｃを使って、次のように定義する。
【数１２】
【００４５】
このとき、Ｉａ(ｕ,ｖ)とＩｂ(ｕ,ｖ)を使って表す２次元画像モデルＩｃ(ｕ,ｖ)を、下記数１３で定義する。
【数１３】
【００４６】
この２次元画像モデルは、図５に示すように、画像全体の回転角度θｇとコーナー角度θｃと２次元濃度エッジのシャープさσimageをパラメータに持つ。０≦θｇ≦π／４、０＜θｃ≦π／２の範囲を考慮する。これ以外の範囲は、画像の左右を反転することと、濃淡を反転することで表現することができる。θｇ＝π／４、θｃ＝π／２のときには、上記数８の単純な画像モデルと同一である。図６(ａ)、(ｂ)、(ｃ)、(ｇ)、(ｈ)、(ｉ)に、上記数１３で表す画像モデルの例を示す。図６(ｄ)、(ｅ)、(ｆ)、(ｊ)、(ｋ)、(ｌ)に、図６(ａ)、(ｂ)、(ｃ)、(ｇ)、(ｈ)、(ｉ)の２次元画像に対応した類似度を示す。
【００４７】
（１―２）繰り返しテクスチャモデル
類似度に方向依存性がないとき、つまり類似度が等方性なら、サブピクセル推定を水平方向と垂直方向に独立に行っても誤差が発生しない。しかし、類似度に方向依存性があるとき（本発明では「異方性」と呼ぶ）、つまり、方向によって微分値が異なるときには、サブピクセル推定を水平方向と垂直方向に独立に行うと誤差が発生する可能性がある。
上記数１３で示した２次元画像モデルＩｃ(ｕ,ｖ)では、θｃ≠π／２のときに自己類似度に異方性が現れた。しかし、直感的に理解できるように、これ以外にも自己類似度に異方性を生じさせるテクスチャが存在する。この例として、図３(ａ)に示す画像のＳＳＤ自己類似度を図３(ｂ)に示す。テクスチャに含まれるコーナーはθｃ＝π／２なのに、類似度に異方性が現れている。
【００４８】
θｃ＝π／２なのに自己類似度に異方性を生じさせる２次元画像モデルとして、下記数１４を考える。
【数１４】
【００４９】
上記数１４の２次元画像モデルＩｔ(ｕ,ｖ)は、図７に示すように、ｕ方向空間周波数ｆｕ、ｖ方向空間周波数ｆｖ、画像の回転角度θｇをパラメータに持つ。０≦θｇ≦π／２の範囲を考慮する。ｆｕ≠ｆｖのときには、自己類似度に異方性を生じる。自己類似度の異方性は、テクスチャにおける空間周波数の異方性に相当する。前節で導入した２次元画像モデルＩｃ(ｕ,ｖ)は、θｃ≠π／２のときに、空間周波数の異方性が発生していると考えることができる。
【００５０】
この２次元画像モデルとＳＳＤ自己類似度を図８に示す。
【００５１】
（１―３）ガウス関数モデル
自己類似度に異方性を生じさせる２次元画像モデルとして、下記数１５の２次元ガウス関数を考えることができる。
【数１５】
ただし、図９に示すように、σはガウス関数の標準偏差で、ｋは異方性係数(ｋ＞０)で、θｇは回転角度である。０≦θｇ≦π／２の範囲を考慮する。この２次元画像モデルとＳＳＤ自己類似度を図１０に示す。
【００５２】
（１―４）２次元類似度モデル
１次元サブピクセル推定と異なり、２次元サブピクセル推定では２次元画像モデルのパラメータ数が多いので、２次元画像モデルから得られるＳＳＤ類似度を直接検討することは得策ではない。ここで取り上げた画像モデルは、画像の性質の多くを含んでいるが、実際用途で扱う画像は、ここで取り上げた画像モデルを複数組み合わせたものや、より高次の画像が存在する。
【００５３】
そこで、以後の検討では、これまでに検討した何種類かの画像モデルを直接使わず、代わりに、これらの画像モデルから得られる２次元類似度を近似する２次元類似度モデルを検討する。２次元類似度モデルとして、下記数１６で表す２次元ガウス関数を採用する。
【数１６】
ただし、（ｄｓ,ｄｔ)は画像間の真の変位で、σはガウス関数の標準偏差で、ｋは異方性係数(ｋ＞０)で、θｇは回転角度である。０≦θｇ≦π／２の範囲を考慮する。
【００５４】
比較のため、ここで取り上げた２次元画像モデルから得られるＳＳＤ自己類似度の例と上記数１６の類似度の例を図１１に示す。ただし、(ｄｓ,ｄｔ)＝(０,０)、θｇ＝０、ｋ＝１としたときに、２次元画像モデルから計算した離散的な自己類似度を最もよく近似するように、上記数１６のパラメータを数値計算で求めた。
【００５５】
以上の説明は、連続領域で行ってきたが、実際の画像データから得られる類似度は、画素単位でのサンプリングになっている。この様子を図１２に示す。本発明の実施例１に係る２次元推定方法は、このように離散化単位で得られた２次元類似度を用いて、連続領域での類似度最大値を与える２次元変位を高精度に推定するものである。
【００５６】
（２）本発明の実施例１に係る２次元推定方法
（２―１）本発明の実施例１に係る２次元推定方法の原理＜連続領域での説明＞
ここでは、前述した２次元類似度モデルを使って、本発明の実施例１に係る２次元推定方法の原理を説明する。本発明は、離散化単位で得られた類似度値を使って、連続領域での類似度最大値位置を高精度に推定することを特徴としている。ここでは、本発明の実施例１の原理を連続領域で説明する。
【００５７】
図１３は２次元類似度モデルを等高線で表示した図である。この２次元類似度モデルのパラメータは、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１.０、k＝０.７、θｇ＝π／６である。等高線は、２次元類似度モデルの最大値に対して、１０％から９０％までレベルを示している。２次元類似度モデルの等高線は楕円になる。この楕円の中心が、２次元類似度モデルの最大値位置である。図１３(ａ)には、２次元類似度モデルの長軸を示してあるが、この角度が回転角度θｇに相当する。
【００５８】
数１６の２次元類似度モデルをｓで偏微分して０とおくことで、図１３(ｂ)の直線ｓ＝ａｔ＋ｂを表す関係を得ることができる。
【数１７】
【００５９】
考慮する条件がｋ＞０なので、上記数１７の各係数の分母≠０である。また、２次元類似度モデルが等方性のとき、つまりｋ＝１のとき、上記数１７はｓ＝ｄｓとなり、推定誤差が発生しない。本発明の実施例１では、この直線を水平極値線ＨＥＬ(Horizontal Extremal Line)と称する。ＨＥＬは図１３(ｂ)の等高線を表す楕円と水平直線との接点を通る。
【００６０】
一方、数１６の２次元類似度モデルをｔで偏微分して０とおくことで、図１３(ｂ)の直線ｔ＝Ａｓ＋Ｂを表す関係を得ることができる。
【数１８】
【００６１】
数１７と同様に、各係数の分母≠０である。また、ｋ＝１のとき、数１８は、ｔ＝ｄｔとなり、推定誤差が発生しない。本発明の実施例１では、この直線を垂直極値線ＶＥＬ(Vertical Extremal Line)と称する。ＶＥＬは図１３(ｂ)の等高線を表す楕円と垂直直線との接点を通る。
【００６２】
ＨＥＬとＶＥＬの交点は、数１６の２次元類似度モデルを異なる方向に偏微分したときに、どちらの方向にも０になる点であり、これは２次元類似度モデルの類似度値を最大にする位置に他ならない。実際に、この交点を計算すると、
【数１９】
となり、当然のことながら２次元類似度モデルの変位(ｄｓ,ｄｔ)を表している。
【００６３】
このときの分母≠０は、ＨＥＬとＶＥＬが交点を持つための条件である。分母を計算すると、次のようになる。
【数２０】
ｋの範囲はｋ＞０なので、この分母１−ａＡ≠０である。
【００６４】
まとめると、２次元類似度のｓ方向微分とｔ方向微分は、ＨＥＬとＶＥＬの２直線で表すことができ、この２直線の交点が２次元類似度の最大値位置になっている。従って、本発明の実施例１では、離散化単位で得られた類似度値を使ってＨＥＬとＶＥＬを求め、その交点を計算すれば、２次元サブピクセル推定を行うことができる。
【００６５】
（２―２）本発明の実施例１に係る２次元推定方法の実装方法＜離散領域への適用＞
以下では、離散的に算出された画像間の類似度値を利用して、連続領域における類似度最大値位置を推定する本発明の実施例１に係る２次元推定方法を具体的に説明する。
【００６６】
まず、離散的に得られている類似度値を使って、ＨＥＬを求める。ＨＥＬは水平方向について類似度最大値を通る直線なので、２本以上の水平ライン上でのサブピクセル類似度最大値位置が決定すれば、ＨＥＬを求めることができる。この様子を図１４(ａ)に示す。同図において、同心楕円は２次元類似度モデルの等高線を表す。同図中の□位置で離散的に得られた類似度値を使って、直線で示すＨＥＬを求める。
【００６７】
図１４(ａ)において、直線ｔ＝０上での最大類似度を与える位置を
変位位置(ｓ,ｔ)での類似度をＲ(ｓ,ｔ)とすると、
【数２１】
【００６８】
この推定結果には、前章で示したように推定誤差が含まれるが、このことについては次節で述べる。直線ｔ＝−１上と直線ｔ＝１上での最大類似度を与える位置をそれぞれ
とすると、
【数２２】
【数２３】
ただし、ｉｔ＝−1、ｉｔ＝１は、直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置オフセットである（後述する）。
【００６９】
これら３点の最大類似度位置を通る直線がＨＥＬである。実際には、画像パターンや画像に含まれるノイズや画像間の相違などのために、これら３点は完全には直線上には乗らない可能性があるので、最小二乗で近似直線を求める。この３点を最小二乗で近似する直線は、下記数２４で求めることができる。
【数２４】
【００７０】
次に、離散化単位で得られた類似度値を使って、ＶＥＬを求める。ＶＥＬは垂直方向について類似度最大値を通る直線なので、２本以上の垂直ライン上でのサブピクセル類似度最大値位置が決定すれば、ＶＥＬを求めることができる。この様子を、図１４(ｂ)に示す。同図において、同心楕円は、２次元類似度モデルの等高線を表す。同図中の□位置で離散化単位で得られた類似度値を使って、直線で示すＶＥＬを求める。
【００７１】
図１４(ｂ)において、直線ｓ＝０上での最大類似度を与える位置を
とすると、
【数２５】
【００７２】
直線ｓ＝−１上と直線ｓ＝１上での最大類似度を与える位置をそれぞれ
とすると、
【数２６】
【数２７】
ただし、ｉｓ＝−１、ｉｓ＝１は、直線ｓ＝−１上と直線ｓ＝１上で最大類似度を与える整数位置オフセットである。
【００７３】
これら３点の最大類似度位置を通る直線がＶＥＬである。この３点を最小二乗で近似する直線は、下記数２８で求めることができる。
【数２８】
【００７４】
ＨＥＬとＶＥＬの交点、即ち、数２４と数２８の交点が、連続領域における２次元類似度最大値のサブピクセル推定位置
になっている。
【数２９】
【００７５】
ところで、数２６、数２７における直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置オフセットｉｔ＝−１、ｉｔ＝１は、必ずしも０になる保証はない。例えば、図１５(ｂ)に示すのは、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、k＝０.３、θｇ＝π／８の２次元類似度モデルだが、このモデルに対するＨＥＬを求めるための、直線ｔ＝−１上の類似度最大値はｓ＝−２付近にある。従って、
を計算するときには、直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置を再探索する必要がある。垂直方向についても同様で
を計算するときには、直線ｓ＝−１上と直線ｓ＝１上で最大類似度を与える整数位置を再探索する必要がある。
【００７６】
このとき、図１５(ｂ)に示す±３の範囲の類似度を使い、−２≦ｉ≦＋２を考慮している。この探索範囲は、多数の現実的な画像をもとに決定した範囲である。もしも、２次元類似度モデルのパラメータ範囲に制限がないとすると、このときの再探索範囲は無限になってしまう。その理由は、(ｄｓ,ｄｔ)＝(０,０)のとき、数１７のＨＥＬについて、下記数３０のＤ(ｋ,θｇ)を考えると、
【数３０】
【００７７】
Ｄ(ｋ,θｇ)が最大値に近づくのは、θｇ→０、ｋ→０のときで、Ｄ(ｋ,θｇ)→cosθｇ／sinθｇ→∞だからである。ところで、ｋ→０のとき、２次元類似度の異方性が無限大になり、ＨＥＬとＶＥＬがほとんど一致する。もしもこのような２次元類似度を作り出す画像があったとしても、そのときには、ほとんど同一になったＨＥＬとＶＥＬ上で類似度最大値を検索すればよい。たとえば、数３０でＤ(ｋ,θｇ)→∞のときには、ＶＥＬ上の３点
における類似度値を求め、パラボラフィッティングを行えばよい。
【００７８】
従来の１次元推定では、図１５(ａ)に示す５個の類似度値を用いてサブピクセル位置を推定していた。このとき、直線ｔ＝０上の３点の類似度(□印)を用いてサブピクセル推定を行って水平方向サブピクセル変位(●印)とした。また、直線ｓ＝０上の３点の類似度(□印)を用いてサブピクセル推定を行って垂直方向サブピクセル変位(●印)とした。図１５(ａ)に示す２直線の交点が従来手法で求めた２次元サブピクセル推定値だが、大きな推定誤差を含むことがわかる。
【００７９】
これに対して本発明の実施例１に係る２次元推定方法は、図１５(ｂ)に示す２５個の類似度値を用いて２次元サブピクセル位置を推定する。ＨＥＬとＶＥＬの交点は、正確に２次元類似度のサブピクセル最大位置を通るので、推定は極めて高精度である。さらに、本発明の実施例１に係る２次元推定方法の計算量の増加は、従来の「１次元推定方法」と比較してわずかである。領域ベースマッチングでは、大部分の計算時間は類似度値を計算しているが、本発明の実施例１に係る２次元推定方法では既に得られている類似度値を利用するので、計算時間の増加は少ない。
【００８０】
（２―３）本発明の実施例１に係る２次元推定方法とサブピクセル推定誤差低減方法との組合せ
数２１から数２３、数２５から数２７において、３位置の類似度を使ってパラボラフィッティングによってサブピクセル位置を推定するときには、推定誤差が発生する。この推定誤差は、補間画像を利用することでキャンセルすることができる。
【００８１】
数２１から数２３、数２５から数２７において求めているサブピクセル位置は、従来どおりの１次元推定にすぎないので、非特許文献１に記載された、推定誤差を低減できるサブピクセル推定方法（以下、本発明の実施例１に係る２次元推定方法と区別するために、非特許文献１に記載されたサブピクセル推定方法を「サブピクセル推定誤差低減方法」と称する）をそのまま適用することができる。ここでは簡単に適用方法を説明する。
【００８２】
マッチングに使う２次元画像関数をＩ１(ｕ,ｖ)、Ｉ２(ｕ,ｖ)とするとき、水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、
【数３１】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０.５,０)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０.５,０)だけ変位していると考えることができるので、その変位を補正した推定結果は、下記数３２、数３３で求めることができる。
【数３２】
【数３３】
【００８３】
数３３も数２１と同様に推定誤差を含むが、その位相が逆になっているので、下記数３４によって推定誤差をキャンセルすることができる。
【数３４】
【００８４】
同様にして、水平ラインｔ＝−１、ｔ＝１について、
【数３５】
【数３６】
【００８５】
ただし、ｉｔ＝−１、ｉｔ＝１は、直線ｔ＝−１上と直線ｔ＝１上でＲｉｕ(ｓ,ｔ)の最大類似度を与える整数位置オフセットである。これらの値は、数２２、数２３とは異なる可能性がある。数３５と数３６を使って、下記数３７、数３８によって推定誤差をキャンセルすることができる。
【数３７】
【数３８】
【００８６】
垂直方向についても同様である。垂直方向補間画像Ｉ１ｉｖ(ｕ,ｖ)は、
【数３９】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。垂直方向補間画像Ｉ１ｉｖ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０,０.５)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０,０.５)だけ変位していると考えることができるので、その変位を補正した推定結果は、下記数４０、数４１で求めることができる。
【数４０】
【数４１】
【００８７】
数４１も数２５と同様に推定誤差を含むが、その位相が逆になっているので、下記数４２によって推定誤差をキャンセルすることができる。
【数４２】
【００８８】
同様にして、垂直ラインｓ＝−１、ｓ＝１について、
【数４３】
【数４４】
【００８９】
ただし、ｉｓ＝−１、ｉｓ＝１は、直線ｓ＝−１上と直線ｓ＝１上でＲｉｖ(ｓ,ｔ)の最大類似度を与える整数位置オフセットである。これらの値は、数２６、数２７とは異なる可能性がある。数４３と数４４を使って、下記数４５、数４６によって推定誤差をキャンセルすることができる。
【数４５】
【数４６】
【００９０】
数３４、数３７、数３８、数４２、数４５、数４６で計算した各ライン上のサブピクセル推定値は、従来の１次元推定方法によって推定されたサブピクセル推定値と比較して、推定誤差が低減されている。従って、サブピクセル推定誤差低減方法によって推定されたこれらの推定結果を使って、更に本発明の実施例１に係る２次元推定方法で２次元サブピクセル推定を行うことで、より高精度な２次元サブピクセル推定が可能になる。
【数４７】
補間画像は、水平方向と垂直方向にそれぞれ１回ずつ作ればよいので、計算量の増加は１次元推定と同程度である。
【００９１】
（３）本発明の実施例１に係る２次元推定方法による推定誤差と従来の推定方法による推定誤差との比較
以下では、２次元画像間の変位推定で、従来一般的に使用されている各種推定方法による推定誤差を本発明の実施例１に係る２次元推定方法による推定誤差と比較する。前述した２次元画像を使用して推定誤差を比較することで、同じ画像に対する比較を行う。
【００９２】
２次元サブピクセル推定誤差は、水平方向と垂直方向の推定誤差を考える必要がある。ところが、推定誤差に影響するのは、２次元画像が持つ３パラメータと画像間の２次元入力変位の合計５パラメータになるので、全てを網羅することは困難である。そこで、代表的な２次元画像を用いて、各方法による推定誤差を比較する。
【００９３】
（３―１）本発明の実施例１に係る２次元推定方法
代表的な２次元画像として、前述したコーナー画像を使用し、パラメータとしてθｇ＝０，π／８，π／４、θｃ＝π／４、σimage＝１.０を考える。０≦θｇ≦π／４以外の２次元画像は、画像の左右を反転することと、濃淡を反転することで表現することができるが、実際の検討は、さらに多くの２次元画像、回転角度、コーナー角度を用いて行っている。他の２次元画像でも、２次元類似度が同じような形になるときには、サブピクセル推定誤差は同じ傾向を示す。
【００９４】
サブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔは、数４７を使って、下記数４８で表すことができる。
【数４８】
【００９５】
サブピクセル推定誤差低減方法を適用しないときの本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔは、数２９を使って、下記数４９で表すことができる。
【数４９】
【００９６】
θｇ＝０,π／８,π／４の３画像に対する数４８、数４９の２次元サブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｃ)(ｄ)、(ｅ)(ｆ)に示す。(ｃ)と(ｄ)は、それぞれサブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。(ｅ)と(ｆ)は、それぞれサブピクセル推定誤差低減方法を適用しないときの本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【００９７】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じているにもかかわらず、サブピクセル推定誤差は極めて小さい。
【００９８】
（３―２）従来の１次元サブピクセル推定方法
従来の１次元推定方法は、数２１と数２５で求めることができる。また、非特許文献１に記載されたサブピクセル推定誤差低減方法を適用した推定値は、数３４と数４２で求めることができる。ここでは、数２１と数２５を使ったときの推定誤差を比較検討する。
【数５０】
【００９９】
θｇ＝０,π／８,π／４の３画像に対する上記数５０のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｇ)、(ｈ)に示す。(ｇ)と(ｈ)は、それぞれ従来の１次元推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１００】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、従来の１次元推定方法でも推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このときには大きな推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１０１】
（３―３）従来の濃度こう配法によるサブピクセル推定
濃度こう配法では、画像間における対応位置に濃度変化がないことを仮定して、水平方向移動量と垂直方向移動量の２つの未知数を含む拘束条件を画素ごとに得る。したがって、この２つの未知数はこのままでは求めることができないため、注目領域を設定し、注目領域中での拘束条件の２乗和を最小にするように、繰り返し計算によって未知数を求める。このようにして得られる変位量には画素単位という制約がなく、サブピクセル単位を得ることができる。
【０１０２】
一方、画像補間手法では、注目領域と周囲の探索領域を離散化単位よりも密に補間する。密に補間した単位で、類似度最大値を探索する。このときには、濃度こう配法と異なり、画像間の移動量に対する制限はない。
【０１０３】
この２種類のサブピクセル変位推定方法は、全く異なる実装に見え、異なる手法として扱われていたが、本質的に同一のものであることが示されている（非特許文献５を参照）。ここでは、補間画像を作るときに用いる補間関数の次数によってサブピクセル変位推定誤差が変化する様子を確認する。そこで、１次補間を利用した補間画像によるサブピクセル推定結果を、濃度こう配法による結果と見なす。さらに、高次の補間関数を利用した補間画像によるサブピクセル推定結果も検討する。
【０１０４】
（３―４）従来のバイリニア画像補間によるサブピクセル推定
サブピクセル変位を求める２枚の画像をＩ１(ｕ,ｖ)、Ｉ２(ｕ,ｖ)とする。これらの画像は離散化サンプリングされているとする。すなわち、ｕ,ｖは整数である。また、画像間の変位を(ｄｓ,ｄｔ)としたときに、０≦ｄｓ≦１、０≦ｄｔ≦１とする。画像間の真の変位がこの範囲にないときには、画素単位のマッチングによって画像全体を移動する。
【０１０５】
バイリニア補間(２方向１次補間)を利用したサブピクセル推定値
は、下記数５１で表すことができる。
【数５１】
ただし、総和範囲は任意形状の注目領域、argminは値を最小にするような
の組を求める演算を表す。
【０１０６】
θｇ＝０,π／８,π／４の３画像に対する上記数５１のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｉ)、(ｊ)に示す。(ｉ)と(ｊ)は、それぞれ水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１０７】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このとき推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１０８】
この推定誤差は、図１６、図１７、図１８の(ｃ)、(ｄ)の結果よりも大きい。すなわち、本発明の実施例１に係る２次元推定方法は、従来のバイリニア補間によるサブピクセル推定方法よりも高精度に変位を推定することができる。バイリニア補間によるサブピクセル推定方法は濃度こう配法と等価なので、本発明の実施例１に係る２次元推定方法は、従来の濃度こう配法よりも高精度にサブピクセル変位を推定することができるとも言える。
【０１０９】
また、濃度こう配法に相当する推定誤差は、図１６、図１７、図１８の(ｇ)、(ｈ)の１次元推定方法の結果よりもはるかに小さい。すなわち、従来広く利用されてきた、領域ベースマッチングを画素単位で行い、サブピクセル変位は類似度を利用して１次元推定する手法より、濃度こう配法ははるかに高精度に２次元サブピクセル変位を推定することができる。このため、領域ベースマッチングよりも濃度こう配法の方が高精度という認識があったが、本発明の実施例１に係る２次元推定方法を用いることで、濃度こう配法よりもさらに高精度にサブピクセル変位を推定することができる。
【０１１０】
（３―５）従来のバイキュービック画像補間によるサブピクセル推定
補間画像を作るときに、より高次の項まで考慮することで、高精度な補間画像を作ることができる。バイキュービック補間(２方向３次補間)は、補間したいサブピクセル画素位置の周囲４×４の画素値を使う補間方法である。バイキュービック補間を利用したサブピクセル推定値
は、下記数５２で表すことができる。
【０１１１】
【数５２】
ただし、総和範囲は任意形状の注目領域、argminは値を最小にするような
の組を求める演算を表す。このとき、重み関数ｈ(ｘ)は、sinc関数をもとに各種提案されていて、下記数５３が多く用いられている。
【０１１２】
【数５３】
ただし、ａは調整パラメータで、ａ＝−１やａ＝−０.５がよく用いられる。
【０１１３】
また、下記数５４も多く用いられる。
【数５４】
Ｂ、Ｃは調整パラメータで、Ｂ＝１、Ｃ＝０のときにはキュービックＢスプラインを近似し、Ｂ＝０、Ｃ＝０.５のときにはCatmull−Romキュービック特性になる。このように、重み係数の選び方によって補間特性が異なる。ここでは、数５３の重み関数において、パラメータａ＝−０.５５を採用した。ただし、画像処理の教科書等に最も多く紹介されているのはａ＝−１である。両者の特性の差については後述するが、ａ＝−０.５５はサブピクセル推定誤差を小さくするように数値的に決定した。
【０１１４】
θｇ＝０,π／８,π／４の３画像に対する数５２のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｋ)、(ｌ)に示す。(ｋ)と(ｌ)は、それぞれ水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１１５】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このとき推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１１６】
この推定誤差は、図１６、図１７、図１８の２次元推定(ｃ)、(ｄ)の結果と同程度である。すなわち、本発明の実施例１に係る２次元サブピクセル推定方法は、最適なパラメータを選択したときの高次補間画像を使ったサブピクセル推定と同程度の高精度であるということができる。
【０１１７】
ところで、数５３におけるパラメータａについて、具体的な計算例を示して説明する。数５３は、sinc関数を有限範囲で近似した重み係数と考えることができる。パラメータａは、その近似特性を調節するものである。図１９(ａ)、(ｂ)に、ａ＝−０.５５のときのｈ(ｘ)と、整数位置に値を持つ１次元関数ｆ(ｉ)の補間結果を示す。ｆ(ｉ)は、数７でσimage＝０.７としたものである。また、図１９(ｃ)、(ｄ)に、ａ＝−１.０のときのｈ(ｘ)と、整数位置に値を持つ１次元関数ｆ(ｉ)の補間結果を示す。
【０１１８】
ａ＝−１.０は、多くの画像処理の教科書などでバイキュービック補間として紹介されるパラメータである。sinc関数を良く近似しているが、少なくともここで検討した関数ｆ(ｉ)を良好に補間しているとは言い難い。ａ＝−０.５５は、sinc関数の近似としてはそれほど良くないが、関数ｆ(ｉ)を良好に補間している。
【０１１９】
図１６から図１８では、パラメータａ＝−０.５５を用いたが、この値はサブピクセル推定誤差を小さくするように数値的に求めた結果である。ａ＝−０.５を用いるとバイリニア補間の結果と同程度、ａ＝−１.０を用いるとバイリニア補間を用いたサブピクセル推定結果よりも誤差が大きくなる。
【０１２０】
（３―６）各方法による推定誤差の総合比較
以下では、２次元画像のパラメータをより広範囲に変更して、各手法のサブピクセル推定誤差を比較検討する。２次元画像は、前述したコーナー画像を使用する。考慮するパラメータとして、回転角度をθｇ＝０,π／１６，２π／１６，３π／１６，４π／１６、コーナー角度をθｃ＝π／４，π／３を考える。合計１０種類の２次元画像に対する推定誤差を検討する。また、σimage＝１.０を考える。
【０１２１】
画像間変位は、−０.５≦ｄｓ≦＋０.５、−０.５≦ｄｔ≦＋０.５の範囲で、０.１[画素]ごとに与える。このとき、２次元推定誤差の大きさ、
のＲＭＳ誤差を評価した。
【０１２２】
評価したサブピクセル推定方法は、前述した各方法の他に、従来の１次元推定方法に対するサブピクセル推定誤差低減方法を適用した推定結果も評価した。したがって、
・サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)
・サブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)
・サブピクセル推定誤差低減方法を用いた１次元推定(１Ｄ−ＥＥＣ)
・サブピクセル推定誤差低減方法を用いない１次元推定(１Ｄ)
・バイリニア補間を用いた２次元推定(Ｂｉ−Ｌｉｎｅａｒ)
・バイキュービック補間を用いた２次元推定(ａ＝−０.５５)(Ｂｉ−Ｃｕｂｉｃ)
の６種類の推定方法の推定誤差を評価した。評価した結果を表１に示す。
【０１２３】
【表１】
【０１２４】
表１から次のことが分かった。
まず、サブピクセル推定誤差低減方法を用いた１次元推定(１Ｄ−ＥＥＣ)は、回転角度θｇ＝０のときには、従来の１次元推定(１Ｄ)に対して推定誤差を低減できる。しかし、回転角度θｇ≠０のときには差が現れない。回転角度θｇは、異方性がある類似度の回転角度に対応する。
【０１２５】
次に、バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)やバイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定は、１次元推定(１Ｄ)に対して推定誤差が２桁程度小さい。この差は歴然である。バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)を用いた２次元推定は、濃度こう配法に対応する。濃度こう配法が高精度にサブピクセル推定ができると考えられてきた理由はここにある。
【０１２６】
それから、バイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定は、バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)を用いた２次元推定よりも高精度である。しかし、重み関数のパラメータによって結果は異なる。重み関数のパラメータを調整するということは、補間フィルタの特性を調整することに対応する。すなわち、画像に対して適切な補間フィルタを選択することで、画像補間を用いたサブピクセル推定手法では精度を向上させることができる。
【０１２７】
最後に、サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)は、バイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定よりも推定誤差が小さく、しかも、パラメータ調整の必要がない。サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)は、サブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)よりも計算コストが大きい。そこで、計算コストを重視して推定結果は濃度こう配法と同程度でよいときにはサブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)を利用し、さらに高精度な推定結果を得たいときにはサブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)を利用する、という使い分けも可能である。
【０１２８】
（４）本発明の実施例１に係る２次元推定方法を実装したコンピュータ・プログラムによる推定実験
本発明の実施例１に係る２次元推定方法を実装したコンピュータ・プログラムを、ＣＰＵを備えた情報処理端末（例えば、パソコン、ワークステーション等のコンピュータ）に実行させることによって、合成画像及び実画像を用いて、２次元サブピクセル推定実験を行った。なお、本発明でいうコンピュータ・プログラムは、便宜上プログラムと略して称することもある。
【０１２９】
（４―１）合成画像を使った推定実験
前述したコーナー画像を作り、２次元サブピクセル推定実験を行った。作成した合成画像は、８bitモノクロ画像で、画像サイズが６４×６４[画素]で、注目範囲が１１×１１[画素]である。画像にノイズは含まれない。画像のパラメータは、σimage＝１.０、コーナー角度θｃ＝π／４、回転角度θｇ＝０,π／８,π／４である。画像間変位は、−０.５≦ｄｓ≦＋０.５、−０.５≦ｄｔ≦＋０.５の範囲で、０.１[画素]ごとに与えた。
【０１３０】
図２０に従来の１次元推定と、サブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元推定についての結果を示す。図２０において、計算結果をメッシュで、実験結果を●で表す。実験結果と計算結果との相違の原因は、画像が８bit量子化階調（２５６階調）であることが考えられる。また、大きく段差がある部分は、整数での最大類似度を探索した結果が隣の画素位置に移動してしまったためと考えることができる。このように隣の位置に移動してしまうことは、従来の１次元推定方法の根本的な欠点である。
【０１３１】
（４―２）実画像を使った推定実験
マッチングに利用する画像から得られる２次元類似度の回転角度θｇ≒０のときには、従来の１次元推定方法でも推定誤差は発生しないが、θｇ≠０のときには、大きな推定誤差が生じる可能性がある。そこで、マッチングパターンとして２種類のパターンを用いたターゲットトラッキング実験を行った。
【０１３２】
使用したマッチングパターンは、円形パターンと傾いたエッジパターンである。これらのパターンを図２１(ａ)、(ｅ)に示す。円形パターンの大きさは、直径が約４１[画素]である。また、これらのパターンのＳＳＤ自己類似度を図２１(ｂ)、(ｆ)に示す。
【０１３３】
円形パターンの自己類似度は等方性なので、従来の１次元推定方法と本発明の実施例１に係る２次元推定方法のどちらでも、推定誤差が小さいことが予想できる。一方、傾いたエッジパターンに対応する自己類似度は異方性でかつ傾いているので、従来の１次元推定方法では推定誤差が大きく、本発明の実施例１に係る２次元推定方法では推定誤差が小さいことが予想できる。
【０１３４】
マッチングパターンを厚紙に印刷して、トラッキングターゲットとした。ターゲットは、ネジ送り式リニアステージ上を直線的に移動する。ターゲットの移動を約２５０枚の時系列画像で撮影し、最初の画像を参照パターンとして使い、以後の画像中からターゲットをトラッキングした。注目領域サイズは６０×６０[画素]（図２１(ａ)、(ｅ)中の黒い矩形領域）で、探索領域は移動予測位置に対して±８(１７×１７)である。ターゲットは画像上で右下方向にゆっくり移動する。マッチングが正しくでき、サブピクセル推定に誤差がなければ、計測軌跡は直線になる。従来の１次元推定方法と本発明の実施例１に係る２次元推定方法のどちらにもＳＳＤ自己類似度を使用し、サブピクセル推定誤差低減方法を適用したパラボラフィッティングによるサブピクセル推定を行った。
【０１３５】
図２１(ｃ)、(ｄ)に、トラッキングパターンとして円形パターン使ったときの計測結果を示す。図２１(ｃ)は従来の１次元推定方法、図２１(ｄ)は本発明の実施例１に係る２次元推定方法を使った結果である。どちらの場合も、サブピクセル推定誤差が小さく、軌跡は良好な直線になっている。
【０１３６】
図２１(ｇ)、(ｈ)に、トラッキングパターンとして傾いたエッジパターン使ったときの計測結果を示す。図２１(ｇ)は従来の１次元推定方法、(ｈ)は本発明の実施例１に係る２次元推定方法を使った結果である。従来の１次元推定方法では、非常に大きなサブピクセル推定誤差が発生していることが明らかである。
【０１３７】
本発明の実施例１に係る２次元推定方法では、トラッキングパターンによらず推定誤差が小さいが、円形パターンのときよりも計測軌跡が直線から外れている。この理由は、２次元類似度が２次元ガウスモデルと異なるため、ＨＥＬとＶＥＬを直線近似するときに誤差が含まれることが原因と考えられる。
【０１３８】
現実的には、このようなトラッキングパターンを利用することはまれである。しかし、この実験は、従来の１次元推定方法を使うと、マッチングパターンによっては、予想外に大きなサブピクセル推定誤差が発生する可能性があることを示している。この実験では、対象が直線上を移動していることがわかっているので、サブピクセル推定誤差を確認することができた。しかし、たとえば、異なる方向から撮影した複数の画像を使って、建築物の３次元情報を再構築するときなどには、傾いたエッジパターンの領域をマッチング領域として利用することもある。従来の１次元推定方法では、このようなときに大きな推定誤差が発生するという問題点があった。
【０１３９】
なお、本発明の実施例１では、水平極値線ＨＥＬと垂直極値線ＶＥＬを求める際、３点を最小二乗で近似する直線を求める代わりに、３点を通る２次曲線として求めることにより、直線近似による誤差を低減するように本発明の実施例１に係る２次元推定方法を拡張することもできる。
【０１４０】
実施例２及び実施例３：
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の実施例２及び実施例３を次のように説明する。本発明の実施例２及び実施例３の着眼点は以下の通りである。
【０１４１】
（Ａ）平行移動と回転、または、平行移動と回転とスケール変化を画像間対応パラメータ（本発明では、略してパラメータとも呼んでいる）にして、これらのパラメータを軸とする多次元空間（本発明では、パラメータ空間或いは多次元類似度空間とも称する）において離散的な位置で得られた類似度を利用する。類似度は離散的な位置で計算しておけばよいので、繰り返し計算の必要がない。
【０１４２】
（Ｂ）離散的な位置で得られた類似度は、多次元空間（つまり、パラメータ空間）における連続関数である類似度をサンプリングしたものであると考える。
【０１４３】
（Ｃ）現実的な仮定の下に類似度をモデル化したとき、多次元空間における各軸（つまり、パラメータ空間における各パラメータ軸）で類似度を偏微分して０とおくことで得られる超平面は、類似度の最大値または最小値を通る。例えば、類似性評価関数にＳＡＤやＳＳＤが採用された場合に、超平面は類似度の最小値を通る。また、類似度評価関数にＺＮＣＣが採用された場合に、超平面は類似度の最大値を通る。
【０１４４】
（Ｄ）このような超平面は、多次元空間の軸の数（つまり、パラメータ数）だけ求めることができる。これらの超平面の交点座標は、多次元空間における類似度の最大値または最小値を与えるパラメータと一致するため、これらの超平面の交点座標を求めることで、多次元空間における類似度の最大値または最小値を与えるパラメータを同時に求めることができる。
【０１４５】
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法は、サンプリンググリッドで離散的に得られた評価値を使って、高精度なサブサンプリンググリッド精度の評価値最大位置または最小位置を推定する方法と考えることができる。
【０１４６】
画像間の対応を平行移動のみに限定して考えるとき（つまり、本発明の実施例１の場合）には、このサンプリンググリッドは画像を構成する画素に相当する。このため、画素単位以上に高精度な平行移動量の表現として、サブピクセルという表現が適切だった。従って、前述したように、本発明の実施例１に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を「２次元サブピクセル推定方法」とも呼んでいる。
【０１４７】
しかし、平行移動と回転、または、平行移動と回転とスケール変化を画像間対応パラメータにする場合の本発明、つまり、本発明の実施例２及び実施例３に係る多パラメータ高精度同時推定方法では、画素やサブピクセルといった表現は一般性を欠いているので、以後ではそれぞれサンプリンググリッド、サブサンプリンググリッドという表現を採用する。
【０１４８】
＜１＞実施例２に係る３パラメータ同時推定方法
＜１−１＞３パラメータ同時推定方法の原理
ここでは、本発明の実施例２に係る３パラメータ同時推定方法の原理を説明する。本発明では、画像間の対応を探索するときには、何らかの類似度評価関数を利用して、類似度評価関数の最小または最大を探索する。このとき得られる類似度は、用いる類似度評価関数によっても異なるが、それ以上に入力として扱う画像によって異なる。
【０１４９】
ここでは、具体的な入力画像や類似度評価関数には触れず、画像間の類似度を次の３次元ガウス関数で近似して、これを３次元類似度モデルとする。
【数５５】
ただし、(ｓ１,ｓ２,ｓ３)は探索パラメータ（つまり、画像間対応パラメータ）で、σはガウス関数の標準偏差で、(ｋ２,ｋ３)は異方性係数(ｋ２,ｋ３＞０)である。また、(ｔ１,ｔ２,ｔ３)は、次のように(ｓ１,ｓ２,ｓ３)に対して回転と平行移動を行った結果である。
【０１５０】
【数５６】
ただし、(ｄ１,ｄ２,ｄ３)は探索パラメータ(ｓ１,ｓ２,ｓ３)の正解値で、(θ１,θ２,θ３)はそれぞれ(ｓ１,ｓ２,ｓ３)の各軸回りの回転角度で、０≦θ１,θ２,θ３≦π／２の範囲を考慮する。
【０１５１】
この３次元類似度モデルをｓ１,ｓ２,ｓ３で偏微分して０とおくと、次の３平面を得ることができる。
【数５７】
【数５８】
【数５９】
【０１５２】
Π１,Π２,Π３は、３次元空間内（つまり、ｓ１,ｓ２,ｓ３を軸とするパラメータ空間）の３平面を構成している。この３平面の交点は、明らかに(ｄ１,ｄ２,ｄ３)であり、正解パラメータと一致する。３平面の交点が求まらない条件は、３次元類似度モデルが３次元未満に縮退した場合である。このような条件では、ｋ２,ｋ３＝０または∞となり、３次元類似度モデルの条件と一致しない。
【０１５３】
＜１−２＞３パラメータ同時推定方法の実装
以下では、３パラメータ同時推定方法において、サブサンプリンググリッド推定値を求める具体的な方法を説明する。例えば、この３パラメータは、平行移動量(ｄｓ,ｄｔ)と回転角度θと考えることができる。
【０１５４】
前節『３パラメータ同時推定方法の原理』で示したように、パラメータ空間における類似度値を各パラメータ軸に関して微分して０とおくことで、各パラメータの最大値位置または最小値位置を通る３平面を作ることができる。従って、サンプリンググリッドにおいて得られた類似度値を用いて、この３平面を推定することができれば、各パラメータの類似度値の最大位置または最小位置をサブサンプリンググリッドにおいて求めることができる。
【０１５５】
まず、平面Π１の求め方を説明する。平面Π１上の点は、ｓ１軸に沿って類似度を微分した結果が０になる点だから、ｓ１軸に平行な直線上の類似度値を使ってサブサンプリング位置を推定することで、平面Π１上の点をいくつか求めることができる。図２２はパラメータ軸ｓ１に関する最大値を通る平面Π１を説明するための模式図である。図２２では、類似度値が得られているサンプリング位置を正方形で、ｓ１軸に平行な直線上で求めたサブサンプリング推定位置を黒丸で示している。類似度値のグリッド単位最大値位置は、(ｒ１,ｒ２,ｒ３)である。
【０１５６】
(ｓ１,ｓ２,ｓ３)における類似度値をＲ(ｓ１,ｓ２,ｓ３)とすると、Ｒ(ｒ１＋ｉ１,ｒ２＋ｉ２,ｒ３＋ｉ３)(ｉ１,ｉ２,ｉ３＝−１,０,１)の２７点の類似度値を使って、平面Π１上の９個の点(ｐ１(ｒ２＋ｉ２,ｒ３＋ｉ３),ｒ２＋ｉ２,ｒ３＋ｉ３)を推定することができる。推定には、次のパラボラフィッティングを用いることができる。
【数６０】
【０１５７】
また、後述するように、パラボラフィッティングによる推定誤差を低減する手法を適用することもできる。
【０１５８】
数６０のパラボラフィッティングを行うときには、Ｒ(ｒ１,ｒ２＋ｉ２,ｒ３＋ｉ３)の類似度値が最大であることを前提としているが、画像パターンによってはこの前提が成立しない場合もある。たとえば、画像に含まれるテクスチャのせん断変形が極端なときには、ｉ２＝ｉ３＝±１の位置でＲ(ｒ１,ｒ２＋ｉ２,ｒ３＋ｉ３)が最大値にならないことがある。このときにも(ｒ２＋ｉ２,ｒ３＋ｉ３)を通りパラメータ軸ｓ１に平行な直線上でサブサンプリンググリッド推定位置を求める必要がある。この場合には、この直線上のサンプリンググリッド上に求まっている類似度値の最大値位置を探索し、修正最大値位置ｒ１’に関して数６０のパラボラフィッティングを行えばよい。
【０１５９】
以後の説明では、数式表記を簡略化するために、ｓ１−ｒ１→ｓ１,ｓ２−ｒ２→ｓ２,ｓ３−ｒ３→ｓ３のように各パラメータ軸に沿って(ｒ１,ｒ２,ｒ３)だけ平行移動して(ｒ１,ｒ２,ｒ３)をパラメータ軸の原点位置と考える。このような平行移動は、グリッド単位での類似度最大位置を原点に移動することと等しく、本発明に関する一般性を失うものではない。本発明に係る多パラメータ高精度同時推定方法によって求めたサブサンプリンググリッド位置を最終的に(ｒ１,ｒ２,ｒ３)と加算することで、平行移動する前のサブサンプリンググリッド位置を簡単に計算することができる。
【０１６０】
平面Π１上にあると期待される９点のサブサンプリンググリッド位置は、実際の画像ではノイズや画像間の変形の影響で、完全には平面上にない可能性がある。このため、これらの９点を最小二乗法によって平面Π１にあてはめる。具体的には、９点と平面Π１とのｓ１軸の沿った距離の２乗和を最小にするように係数を決めればよい。平面Π１は、次のように求めることができる。
【０１６１】
【数６１】
となる。同様にして、平面Π２,Π３は、次のように求めることができる。
【０１６２】
【数６２】
【数６３】
【０１６３】
平面Π１,Π２,Π３の交点
が、サブサンプリング推定位置で、次のように求めることができる。
【数６４】
【数６５】
【０１６４】
＜２＞実施例３に係るＮパラメータ同時推定方法
＜２−１＞Ｎパラメータ同時推定方法の原理
ここでは、本発明の実施例３に係るＮパラメータ同時推定方法の原理を説明する。パラメータ数が３より多いときにも、本発明では同様の考え方でサブサンプリンググリッド高精度同時推定をすることができる。ｓ１からｓＮまでのＮ個のパラメータによって類似度値が決定するときに、類似度値を次のようにガウス関数のｎ個の乗算として表現する。
【０１６５】
【数６６】
ただし、
はＮ次元類似度空間における座標を表すＮ次元ベクトルで、
はＮ次元空間での回転と平行移動を表す行列で、
はガウス関数の異方性係数ベクトル(Ｎ−１次元)である。
【０１６６】
このとき、数６６で表されるＮ次元類似度モデルをｓｉ(ｉ＝１,２,…,Ｎ)で偏微分して０とおくことで、次のようなＮ個のＮ次元関係式を得ることができる。
【数６７】
【０１６７】
これらの関係式を連立させて解くことで、Ｎ次元類似度のサブサンプリンググリッド推定値を得ることができる。つまり、それらの関係式の解は、Ｎ次元類似度のサブサンプリンググリッド推定値になるわけである。なお、Ｎ次元類似度モデルが数６６のガウス関数乗算モデル近似からずれると、数６７は完全な超平面にはならないが、このときには、最小二乗法などで超平面に近似すればよい。
【０１６８】
＜２−２＞Ｎパラメータ同時推定方法の実装
以下では、Ｎパラメータのサブサンプリンググリッド同時推定を実際に行う具体的な方法を説明する。まず、サブサンプリンググリッド同時推定を行うときに、前提としてサンプリンググリッドでの類似度最大位置はわかっており、その(変位)位置を
とする。
【０１６９】
類似度値Ｒをパラメータ軸ｓｉ(ｉ＝１,２,…,Ｎ)について偏微分した結果が０になるＮ次元超平面Πｉ上には、ｓｉパラメータ軸と平行な直線上の類似度値
の３つの値をパラボラフィッティングによってサブサンプリンググリッド推定した３Ｎ−１個の位置、
が存在する。ただし、
【数６８】
である。
【０１７０】
ここで、
は、ｉ番目の要素だけが値ｋ(ｋ＝−１,０,＋１)であり、他の要素は全て値０であるようなＮ次元ベクトルである。たとえば、
である。また、
は、ｉ番目の要素が０で、他の各要素は、−１,０,１のいずれかの値をとるようなｊ番目(ｊ＝１,２,…,３Ｎ−１)のＮ次元ベクトルである。たとえば、
である。
【０１７１】
これらの３Ｎ−１個の位置、
は、Ｎ次元超平面Πｉ上に存在するか、少なくとも最小二乗の意味でＮ次元超平面Πｉにフィットする。従って、Ｎ次元超平面Πｉを
【数６９】
と表すと、３Ｎ−１個の位置、
に対して、
【数７０】
の関係が得られる。これを、次のように書き換える。
【０１７２】
【数７１】
ただし、
は、３Ｎ−１行Ｎ列の行列で、
は要素数Ｎのベクトルである。Ｎ次元超平面Πｉを表す係数
は、一般化逆行列を用いた最小二乗法によって、次のように求めることができる。
【０１７３】
【数７２】
このようにして求めたＮ個のＮ次元超平面Πｉ(ｉ＝１,２,…,Ｎ)の交点として、Ｎパラメータのサブサンプリンググリッド推定値を得ることができる。この交点は、次のように求めることができる。
【０１７４】
【数７３】
【数７４】
【０１７５】
以上の方法でＮパラメータの同時推定を行うことができるが、実際にはＮ次元超平面Πｉの式を、同じＮ個の未知数を含む次のような形式に表現してもよい。
【数７５】
【０１７６】
その理由は、Ｎ次元類似度関数がＮ次元対称になっても安定に超平面の交点が存在するようにするためである。
【０１７７】
また、Ｎ次元超平面Πｉ上の３Ｎ−１個の位置を求めることでＮ次元超平面の係数を求めたが、この個数を減らすこともできる。たとえば、Ｎ＝３のときには３Ｎ−１＝９だが、たとえばｊ＝３のときには、具体的に次のようになっている。
(−１,−１,０) (０,−１,０) (＋１,−１,０)
(−１,０,０) (０,０,０) (＋１,０,０)
(−１,＋１,０) (０,＋１,０) (＋１,＋１,０)
【０１７８】
この中から、次のような組を採用することで、個数を２(Ｎ−１)＋１個にすることができる。
(０,−１,０)
(−１,０,０) (０,０,０) (＋１,０,０)
(０,＋１,０)
【０１７９】
＜３＞本発明による実験結果と従来手法による実験結果との比較
図２３に示す時系列画像を人工的に作成した。この時系列画像は、濃淡で表現した２次元ガウス関数で構成されている。この実験で用いた２次元ガウス関数には方向性(異方性) があり、しかもその方向性が０または９０度以外の方向になるように設定されている。このため、この画像から計算する類似度にも異方性があり、かつ回転角度が０または９０度以外になる。
【０１８０】
図２３に示す画像は時系列画像の最初のフレームで、以後のフレームは最初のフレームに対して位置が変化せず回転角度だけがフレームごとに０.１度ずつ変化する。時系列画像は全部で３１フレームあり、最終フレームは先頭フレームに対して３度回転した画像になっている。これに対して画像の位置は回転中心に対して移動しない。
【０１８１】
図２４に、最初のフレームに対する以後のフレームの回転角度を計測した結果を示す。＋印で示す従来方法は、画像を６×６＝３６個の領域に分割し、各小領域ごとに平行移動量をマッチングによって計測し、得られた平行移動量を使って回転角度を推定した結果である。○印で示す本発明による方法は、従来方法で計測した回転角度を初期推定値として本発明を適用し、さらに高精度に回転角度を計測した結果である。フレームごとに０.１度ずつ回転角度が変化する状況が正確に計測できている。
【０１８２】
図２５に、最初のフレームに対する以後のフレームの位置を計測した結果を示す。＋印で示す従来方法は、画像を６×６＝３６個の領域に分割し、各小領域ごとに平行移動量をマッチングによって計測し、得られた平行移動量の平均として位置を推定した結果である。○印で示す本発明による方法は、従来方法で計測した位置を初期推定値として本発明を適用し、さらに高精度に位置を計測した結果である。位置が変化しない状況が正確に計測できている。
【０１８３】
＜４＞本発明の他の実施例
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の他の実施例を以下のように説明する。
【０１８４】
＜４−１＞サブピクセル推定誤差低減方法との組み合わせ
数６０でサブサンプリンググリッド推定位置を求めるときに、３位置の類似度を使ってパラボラフィッティングを用いているので、この推定値には固有の推定誤差を含んでいる。この推定誤差は、補間画像を利用することでキャンセルすることができる。ここでは、２パラメータのときに画像間の水平方向変位に関するサブピクセル位置を推定するときを例に説明する。
【０１８５】
マッチングに使う２次元画像関数をＩ１(ｕ,ｖ),Ｉ２(ｕ,ｖ)とするとき、ＳＳＤとパラボラフィッティングによって推定する水平方向サブピクセル位置は、次のように表すことができる。
【数７６】
【数７７】
【０１８６】
一方、画像Ｉ１(ｕ,ｖ)の水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、
【数７８】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。
【０１８７】
水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０.５,０)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０.５,０)だけ変位していると考えることができるので、その変位を補正した推定結果は、次式で求めることができる。
【数７９】
【数８０】
【０１８８】
数７８も数６０と同様に推定誤差を含むが、その位相が逆になっているので、次式によって推定誤差をキャンセルすることができる。
【数８１】
【０１８９】
＜４−２＞アフィンパラメータ推定
画像間の対応がアフィン変換で表現できるとき、下記数８２の変形モデル中の６個のパラメータを決定する必要がある。
【０１９０】
【数８２】
ただし、(ｕ１,ｖ１),(ｕ２,ｖ２)は、それぞれ画像Ｉ１,Ｉ２を表す座標で、画像Ｉ２を変形して画像Ｉ１に合わせている。
【０１９１】
離散的な位置で得られた類似度値を用いて、この６つのパラメータを同時に高精度に推定することができる。
【０１９２】
＜４−３＞射影パラメータ推定
画像間の対応が射影変換で表現できるとき、下記数８３の変形モデル中の８個のパラメータを決定する必要がある。
【０１９３】
【数８３】
ただし、(ｕ１,ｖ１),(ｕ２,ｖ２)は、それぞれ画像Ｉ１,Ｉ２を表す座標で、画像Ｉ２を変形して画像Ｉ１に合わせている。
【０１９４】
離散的な位置で得られた類似度値を用いて、この８つのパラメータを同時に高精度に推定することができる。
【０１９５】
＜４−４＞高速探索手法との組み合わせ
本発明によるサブサンプリンググリッド推定方法は、サンプリンググリッド単位での対応が正しく得られた後の高精度化処理と考えることもできる。したがって、サンプリンググリッド単位での探索には、すでに提案されている各種の高速化手法を用いることができる。
【０１９６】
例えば、画像の濃度情報を利用した高速探索手法としては、画像ピラミッドを用いた階層構造探索法、ＳＳＤやＳＡＤに上限値を設けて以後の積算を打ち切る手法などがある。また、画像の特徴を利用する高速探索方法としては、注目領域の濃度ヒストグラムを利用する方法などがある。
【０１９７】
なお、本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を実装したコンピュータ・プログラムを、ＣＰＵを備えた情報処理端末（例えば、パソコン、ワークステーション等のコンピュータ）に実行させることによって、合成画像或いは実画像を用いて、多パラメータを同時に高精度に推定することができる。
【産業上の利用可能性】
【０１９８】
以下のように、本発明を色々な画像処理分野に適用することができる。例えば、リモートセンシングや航空写真を使った地図作製では、複数の画像をつなぎ合わせること、つまりモザイキングを行うことが多い。このとき、画像間の対応位置を求めるために、領域ベースのマッチングを行う。このときには、対応が正しいことはもちろん、対応が画素単位より細かな分解能で求まることが重要になう。さらに、多数の画像を利用して対応を求めるので、計算の高速性も重要である。領域ベースマッチングを使って画素分解能以上の対応を求めるときには、画素単位のとびとびの位置で求まる類似度評価値を使って、フィッティング関数によるサブピクセル推定を行っているが、このときに本発明を適用することで、計算時間の増加がわずかで、はるかに高精度な対応位置を求めることができる。
【０１９９】
特に、画像パターンに斜め成分が含まれるとき、つまり、地図作製などで道路交差点などのパターンが画像座標系に対して傾いて撮影されているときに、本発明の効果が一層発揮される。
【０２００】
また、ステレオビジョンでは、複数枚の画像間の対応位置を求めることが基本的な処理であり、このときの対応位置精度が最終的な結果精度に大きく影響を及ぼす。画像間の対応位置を求めるときには、エピポーラ拘束などの拘束条件を利用するが、エピポーラ拘束条件自体を求めるときにも画像間の対応を調べる必要があり、このときにはエピポーラ拘束を利用できない。従って、本発明を適用することによって、高精度にエピポーラ拘束条件を求めることが可能になるため、結果として高精度なステレオ処理を行うことができるようになる。
【０２０１】
さらに、ＭＰＥＧ圧縮を行うときには、隣接フレーム間の対応位置を正確に求めることが、高品質で高圧縮率を達成するために必要である。隣接フレーム間の対応位置は、通常は１／２画素単位で求めるが、このときにはＳＡＤやＳＳＤ類似度評価値を利用した１次元の簡易手法を利用しているため、大きな誤差が含まれている可能性が在る。ところが、本発明は、従来の１次元推定方法に対してわずかな計算量の増加で、より確実で高精度な画像間の対応を求めることができるため、ＭＰＥＧ画像生成の圧縮率を向上することが可能になる。
【図面の簡単な説明】
【０２０２】
【図１】図１は、従来の１次元サブピクセル推定方法を説明するための図である。
【図２】図２は、３次元再構成アプリケーション用に一般的に使用される画像例（ａ）とその自己類似度（ｂ）を示す図である。自己類似度マップは、自己類似度マップがｋ≠１及びθｇ≠０の性質を有しているので、従来の１次元推定方法を用いることでサブピクセル推定誤差が生じることを示している。
【図３】図３は、θｃ＝π／２を有するテクスチャ画像例とその自己類似度を示す図である。
【図４】図４は、２次元画像モデルとその２次元類似度を示す図である。（ａ）はσimage＝０.７の画像モデルを示し、画像のサイズは１１×１１[画素]である。（ｂ）は画像モデルの自己類似度を示し、類似度の範囲は１１×１１である。（ｓ,ｔ）の暗い位置は高い類似度を有する２枚の画像の変位に対応する。最も暗い位置の（ｓ,ｔ）は、自己類似度の性質で（０,０）になる。
【図５】図５は、画像全体の回転角度θｇとコーナー角度θｃを有する２次元コーナー画像モデルを示す図である。
【図６】図６は、２次元コーナー画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、θｃ＝π／２、σimage＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、θｃ＝π／４、σimage＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。
【図７】図７は、画像全体の回転角度θｇと水平方向空間周波数ｆｕと垂直方向空間周波数ｆｖを有する２次元繰り返しテクスチャ画像モデルを示す図である。
【図８】図８は、２次元繰り返しテクスチャ画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、ｆｕ＝０.１、ｆｖ＝０.１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、ｆｕ＝０.０５、ｆｖ＝０.１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。自己類似度は一定の変位を有する２枚同様な画像から得られる。画像のサイズは１１×１１[画素]で、類似度の範囲は−５から５までである。
【図９】図９は、画像全体の回転角度θｇと水平標準偏差σと垂直標準偏差ｋσを有する２次元ガウス画像モデルを示す図である。
【図１０】図１０は、２次元ガウス画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、σ＝２、ｋ＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、σ＝２、ｋ＝０.５の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。自己類似度は一定の変位を有する２枚同様な画像から得られる。画像のサイズは１１×１１[画素]で、類似度の範囲は−５から５までである。
【図１１】図１１は、３種類異なった２次元画像モデルと、それらの２次元自己類似度と、２次元類似度モデルに比較して２次元自己類似度の部分図を示す図である。(ａ)はθｃ＝π／４、σimage＝１、θｇ＝０を有する２次元コーナー画像モデルを示す。(ｂ)はｆｕ＝０.０５１、ｆｖ＝０.１、θｇ＝０を有する２次元繰り返しテクスチャ画像モデルを示す。(ｃ)はσ＝２、ｋ＝０.５、θｇ＝０を有する２次元ガウス画像モデルを示す。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、整数位置(〇で表す)でサンプリングされた(ｄ)、(ｅ)、(ｆ)のｔ＝０の水平部分図と、２次元類似度モデル(点線で表す)を示す。
【図１２】図１２は、２次元連続的に類似度計算と離散的に類似度計算の例を示す図である。
【図１３】図１３は、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、ｋ＝０.７、θｇ＝π／６を有する２次元類似度モデルを示す図である。等高線は関数値の１０％から９０％までに対応している。(ａ)において、２次元類似度モデルの長軸がθｇ＝π／６を示す。(ｂ)において、ＨＥＬ（水平極値線）とＶＥＬ（垂直極値線）の両方が２次元類似度モデル(ｄｓ,ｄｔ)のピーク値位置を通るが、２次元類似度モデルの長軸或いは短軸と異なる。
【図１４】図１４は、ＨＥＬとＶＥＬの推定プロセスを説明するための図である。(ａ)において、●は３本の水平線ｔ＝−１、ｔ＝０、ｔ＝１上の類似度値（□）から得られた推定サブピクセル位置を示し、ＨＥＬは最小二乗で３つ●の位置から推定されることが可能である。(ｂ)において、ＶＥＬは同様なプロセスで推定されることが可能である。
【図１５】図１５は、従来の１次元推定方法及び本発明の実施例１に係る２次元推定方法にそれぞれ必要な類似度値の位置を示す図である。２次元類似度モデルは(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、ｋ＝０.３、θｇ＝π／８を有する。
【図１６】図１６は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝０、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１７】図１７は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝π／８、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１８】図１８は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝π／４、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１９】図１９は、重み関数とパラメータに依存する補間特性を説明するための図である。
【図２０】図２０は、合成画像を使った実験結果を示す図である。
【図２１】図２１は、実画像を使ったターゲットトラッキング実験結果を示す図である。
【図２２】図２２は、本発明の実施例２に係る３パラメータ同時推定方法において、パラメータｓ１に関する最大値を通る平面Π１を説明するための模式図である。
【図２３】図２３は、実験に用いた合成画像である。
【図２４】図２４は、画像の回転計測結果を示す図である。
【図２５】図２５は、画像の位置計測結果を示す図である。
【技術分野】
【０００１】
本発明は、例えば画像による位置計測、リモートセンシング、航空写真を利用した地図作製、ステレオビジョン、画像張り合わせ（モザイキング）、動画像位置合わせ、３次元モデリングなど多数の画像処理分野で使用できる、領域ベースマッチングを用いた推定方法に関し、特に、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法に関するものである。
【背景技術】
【０００２】
ステレオ画像処理をはじめとして、トラッキング、画像計測、リモートセンシング、画像レジストレーションなどの多くの画像処理分野では、画像間変位を得るために、領域ベースマッチングが基本的な処理として採用されている。
【０００３】
領域ベースマッチングは、注目領域のサイズや形状を任意に設定できること、注目画素に対して注目領域をオフセットできること、計算が直感的で直接的な実装が可能であることなどが特徴である。領域ベースマッチングで画像間の変位をサブピクセルで求めるときには、画素単位で離散的に得られる類似度評価値を利用してサブピクセル変位を推定する方法が一般的に利用されている。
【０００４】
画像間の対応を求めるためには、類似度と呼ぶ評価値を利用する。類似度とは、画像間の相対的な位置関係を入力、その位置関係のときの画像間の類似度を出力とする関数である。類似度は、画素単位の離散的な位置で値が決まるので、類似度だけをもとに画像間の対応位置関係を求めると、画素単位の分解能に制限される。そこで、類似度値をもとに補間することによって、サブピクセル分解能の画像間の対応を求める。
【０００５】
従来、位置分解能を拡大するために、類似度評価値をフィッティング内挿することでサブピクセル推定を行うことが多い。しかし、画像の平行移動量をサブピクセル分解能で求める際に、画像の水平方向と垂直方向を独立にフィッティング内挿していたため、十分な精度が得られなかったという問題点があった。
【０００６】
要するに、従来、デジタル画像を使った領域ベースのマッチングでは、変位の推定分解能を向上させるために、画像の水平／垂直方向について独立に類似度評価値を使ったサブピクセル推定を行っていた。
【０００７】
また、従来、濃度こう配法を利用して画像間の変位を求めると、はじめからサブピクセル変位量を得ることができる。濃度こう配法では、水平方向と垂直方向を合わせて扱っている。ただし、画像間変位が１画素以上のときには画像を縮小する必要がある。通常は、画像スケールの変更と同時に実装する。濃度こう配法の実装は繰り返し計算になるので、計算時間の見積もりや、ハードウェアへの実装が困難であるという欠点があった。
【０００８】
さらに、画像をＤＦＴしてから複素共役積をもとめて逆ＤＦＴする従来手法では、注目領域のサイズが２ｎでなければならず、各種テクニックが必要である。流体計測分野では多く用いられ、計算量を小さくできることが特徴である。しかし、ビジョン分野では、ほとんど利用されない。しかも、計測対象が平面であるという仮定を用いているので、凹凸がある計測対象には適用が困難であるという欠点があった。
【０００９】
従来は、以下のように詳細に説明されるように、画像の水平方向と垂直方向を独立に補間推定していた。しかし、従来のこのような方法によると、図１に示すように、推定誤差が発生するという問題点があった。
【００１０】
ここでは、画像を水平方向と垂直方向に独立と考え、変位のサブピクセル推定も独立に行う従来方法を「１次元推定方法」と称する。
【００１１】
従来の１次元推定方法の問題点を説明するために、まず、画像間の水平方向変位を求める問題を考える。真の画像間変位を(ｄｓ,ｄｔ)、異方性を持った類似度の回転角度をθｇとする。１次元推定方法では、離散化類似度として、Ｒ(−１,０)、Ｒ(０,０)、Ｒ(１,０)(図１の□)を使い、下記数１によって
を推定した。なお、図１は、従来の１次元サブピクセル推定方法を説明するための図である。図中の
は、ｓ＝−１,０,１での３つの類似度の値を用いた従来の１次元サブピクセル推定方法によって推定された位置で、ｄｔ≠０及びθｇ≠０の際に、この推定値は真のピーク位置に対して誤差を有する。
【数１】
【００１２】
仮に上記数１で直線ｔ＝０上の類似度最大位置を正しく推定できたとしても、図１から明らかなように、真の水平方向画像間変位ｄｓ(図１の▲の水平成分)に対して大きな推定誤差を含んでいる。すなわち、次の条件の全てが真のときには、水平方向推定誤差が発生して、
となる。
・垂直方向変位ｄｔ≠０。
・２次元類似度が異方性を持つ。
・異方性を持つ２次元類似度の回転角度θｇ≠０。
大部分の画像が上記条件に当てはまる。
【００１３】
また、垂直方向に関しても同様である。例えば、図２(ａ)に示す画像中のコーナー領域のＳＳＤ自己類似度を求めると、図２(ｂ)のようになっている。この自己類似度は異方性を持ち、θｇ≠０なので、従来の１次元推定方法では、サブピクセル推定誤差が発生する可能性があることを示している。また、図３(ａ)に示すテクスチャ画像を使ったサブピクセル推定においても、推定誤差が発生する可能性を示している。
【００１４】
また、コンピュータビジョンの分野によっては、エピポーラ拘束などの拘束条件を利用することで探索範囲を１次元に限定できるので、１次元信号の高精度なマッチングができれば十分な場合がある。しかしながら、２次元画像を使って、高精度な２次元変位を求める必要がある用途も多い。例えば、動き推定、動き推定による領域分割、ターゲットトラッキング、モザイキング、リモートセンシング、超解像などである。
【００１５】
画像の領域ベースマッチングにおいて、画像間変位が平行移動だけで表現できるときには、従来から多くの手法が提案され、実際に用いられてきた。たとえば、代表的な手法として、領域ベースマッチングと類似度補間手法を組み合わせた、サブピクセル推定手法は、次のようにして画像間のサブピクセル変位を推定する。
【００１６】
対応位置を求める画像を、ｆ(ｉ,ｊ)とｇ(ｉ,ｊ)としたとき、整数単位の変位(ｔｘ,ｔｙ)に対する画像間のＳＡＤを、次のように計算する。
【数２】
ただし、ＳＡＤは輝度差の絶対値の和(Sum of Absolute Differences)を表し、Ｗは対応を求める注目領域を表す。ＳＡＤは、画像が完全に一致したときに０になり、不一致が増えるに従って大きな値をとる。ＳＡＤは類似していない評価値を表すので、非類似度である。
【００１７】
ＳＡＤの代わりに、次のＳＳＤを用いることも多い。
【数３】
ただし、ＳＳＤは輝度差の絶対値の２乗和(Sum of Squared Differences)を表す。ＳＳＤも画像が完全に一致したときに０になり、不一致が増えるに従って大きな値をとる。ＳＳＤも類似していない評価値を表すので、非類似度である。
【００１８】
ＳＡＤやＳＳＤの代わりに、次の相関係数(ＺＮＣＣ:Zero-mean Normalized Cross-Correlation)を用いることもある。ＺＮＣＣは、注目領域内の画素濃度を統計データと見なし、注目領域間の統計的な相関値を計算している。
【数４】
ただし、
は、それぞれの領域内の濃度平均である。ＲＺＮＣＣは−１から１までの値をとり、画像が完全に一致したときに１になり、不一致が増えるに従って小さな値をとる。ＺＮＣＣの特徴は、画像間に濃度やコントラストの変化があってもほとんど同じ類似度が得られることである。
【００１９】
ＺＮＣＣは類似性の評価値を表すので類似度とよばれるが、ＳＡＤやＳＳＤは非類似度を表す。ここでは、説明を簡単にするために、以後はＳＡＤやＳＳＤ、ＺＮＣＣを「類似度」と統一して表記する。
【００２０】
類似度が最小値(ＳＡＤやＳＳＤのとき)または最大値(ＺＮＣＣのとき)になる整数単位の変位
を探索することで、画像間の整数単位変位を得ることができる。整数単位の変位を得た後で、次のようにサブピクセル変位を推定する。
【数５】
上記数５を用いて、最終的に得る画像間サブピクセル変位は、
となる。
【００２１】
このような従来手法では、画像間に回転やスケール変化があるときに正しい対応を得ることができない問題があった。画像間の回転角度がわずかでも、２つの画像が異なるものと判断され、対応位置が見つからないことや、誤った位置に対応を検出する問題があった。
【００２２】
一方、画像間の対応が平行移動だけでは表現できない場合、つまり、画像間に回転やスケール変化があるときには、回転やスケール変化も同時に推定する必要がある。従来、このときには、補間処理を用いて片方の画像(テンプレート画像)を平行移動と回転やスケール変化をさせてからマッチングを行い、類似度を計算しながら各パラメータを変化させ、類似度が最小値または最大値になるパラメータを探索する必要があった。次に、このパラメータ探索法を説明する。
【００２３】
ｆ(ｉ,ｊ)をテンプレート画像、ｇ(ｉ,ｊ)を入力画像とするとき、平行移動(ｔｘ,ｔｙ)および回転角度θとスケール変化ｓをパラメータとする類似度ＲＳＳＤを次のように計算する。
【数６】
ただし、
は、平行移動(ｔｘ,ｔｙ)および回転角度θとスケール変化ｓを与えたときの、位置(ｉ,ｊ)における画像ｇ(ｉ,ｊ)の補間画像を表す。補間画像を作るときには、双線形補間(Bi-Linear補間)や双３次補間(Bi-Cubic補間)などを利用する。また、ＳＳＤではなく、ＳＡＤやＺＮＣＣを利用してもよい。
【００２４】
何らかの方法で初期パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜０＞を求めてＲＳＳＤ(ｔｘ,ｔｙ,θ,ｓ)＜０＞を計算し、より小さなＲＳＳＤを与える更新パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜１＞を探索する。パラメータを更新しても結果に変化がなくなるまで、更新を続ける。このとき、Newton-Raphson法やSteepest (Gradent) Descent法(最急降下法)、Levenberg-Marquadt 法など数値的解法を用いる。
【００２５】
初期パラメータ(ｔｘ,ｔｙ,θ,ｓ)＜０＞は、たとえば、画像間の対応が平行移動だけで表現できると仮定して通常の領域ベースマッチングによって推定した
を用いた
などを利用する。
【非特許文献１】清水雅夫・奥富正敏,「画像のマッチングにおける高精度なサブピクセル推定手法」,電子情報通信学会論文誌D-II,2001年7月,第J84-D-II巻,第7号,p.1409-1418
【非特許文献２】清水雅夫・奥富正敏,「画像のマッチングにおける高精度なサブピクセル推定の意味と性質」,電子情報通信学会論文誌D-II,2002年12月,第J85-D-II巻,第12号,p.1791-1800
【非特許文献３】マサオ・シミズ、マサトシ・オクトミ(Masao Shimizu and Masatoshi Okutomi),「プリサイス サブピクセル エスティメイション オン エリアベイセド マッチング(Precise Sub-Pixel Estimation on Area-Based Matching)」,プロク. 8th IEEE インターナショナル コンファレンス オン コンピュータ ビジョン (ICCV2001)（Proc. 8th IEEE International Conference on Computer Vision (ICCV2001)）,(カナダ,バンクーバー),2001年7月,p.90-97
【非特許文献４】マサオ・シミズ、マサトシ・オクトミ(Masao Shimizu and Masatoshi Okutomi),「アン アナリシス オフ サブピクセル エスティメイション エラー オン エリアベイセド イメージ マッチング(An Analysis of Sub-Pixel Estimation Error on Area-Based Image Matching)」,プロク. 14th インターナショナル コンファレンス オン デジタル シグナル プロセシング(DSP2002)（Proc. 14th International Conference on Digital Signal Processing (DSP2002)）,(ギリシア,サントリーニ),2002年7月,第II巻,p.1239-1242(W3B.4)
【非特許文献５】Ｃ.・クエンティン・デイビス、ゾハー・Ｚ．・カル、デニス・Ｍ.・フリーマン(C. Quentin Davis, Zoher Z. Karu and Dennis M. Freeman),「イクイバレンス オフ サブピクセル モーション エスティメイタス ベイセド オン オプティカル フロー アンド ブロック マッチング(Equivalence of subpixel motion estimators based on optical flow and block matching)」,IEEE インターナショナル シンポジウム フォー コンピュータ ビジョン1995（IEEE International Symposium for Computer Vision 1995）,(米国フロリダ州),1995年11月,p.7-12
【非特許文献６】サムソン・Ｊ.・ティモネー、デニス・Ｍ.・フリーマン(Samson J. Timoner and Dennis M. Freeman),「マルチイメージ グラジエントベイセド アルゴリズム フォー モーション エスティメイション(Multi-image gradient-based algorithms for motion estimation)」,オプティカル エンジニアリング（Optical Engineering）,(米国),2001年9月,第40巻,第9号,p.2003-2016
【非特許文献７】ジョナサン・Ｗ.・ブラント (Jonathan W. Brandt),「インプルーブド アキュラシ イン グラジエントベイセド オプティカル フロー エスティメイション(Improved Accuracy in gradient-based Optical Flow Estimation)」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International Journal of Computer Vision）,(米国),1997年,第25巻,第1号,p.5-22
【非特許文献８】Ｑ.・テン、Ｍ．・Ｎ．・ヒューニッス(Q. Tian and M. N. Huhns),「アルゴリズム フォー サブピクセル レジストレーション(Algorithms for Subpixel Registration)」,コンピュータ ビジョン,グラフィックス アンド イメージ プロセシング（Computer Vision, Graphics and Image Processing）,(米国),1986年,第35号,p.220-233
【非特許文献９】ショーン・ボーマン、 マーク・Ａ.・ロバートソン、ロバート・Ｌ.・スティブンソン(Sean Borman, Mark A. Robertson and Robert L. Stevenson),「ブロックマッチング サブピクセル モーション エスティメイション フロム ノイジー,アンダー―サンプルド フレームス---アン エンパイリカル パフォーマンス エバリュエーション(Block-Matching Sub-Pixel Motion Estimation from Noisy, Under-Sampled Frames---An Empirical Performance Evaluation)」,SPIE ビジュアル コミューニケイションズ アンド イメージ プロセシング 1999（SPIE Visual Communications and Image Processing 1999）,(米国カリフォルニア州)
【発明の開示】
【発明が解決しようとする課題】
【００２６】
このような多パラメータ最適化では、初期値が適切でないと正しい結果が得られない問題がある。また、画像パターンやノイズ、撮影レンズの光学歪みなどの影響で解が得られない問題もある。さらに、繰り返し計算を利用するので、計算時間を見積もることができないといった難点がある。このため、アルゴリズムを実装する上で、完成したシステムの総合応答時間を正確に見積もることができず、ハードウェア化することがきわめて困難であった。また、計算時間の面では、繰り返し計算の各段階で画像を補間する必要があるが、画像補間演算は演算量が多いため、多くの計算時間が必要になる。さらに、推定精度の面では、画像補間手法によって補間画像の性質が異なるので、採用する補間手法によって最終的なパラメータ推定精度や推定誤差の性質が異なる点も大きな問題であった。
【００２７】
本発明は上述のような事情よりなされたものであり、本発明の目的は、Ｎ次元類似度を考慮することで、上記従来手法の問題点を解決し、画像間対応パラメータを少ない計算量で安定且つ高速に高精度に同時推定できるようにした、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を提供することにある。
【課題を解決するための手段】
【００２８】
本発明は、画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法に関し、本発明の上記目的は、Ｎ（Ｎは４以上の整数である）個の画像間対応パラメータを軸とするＮ次元類似度空間において、離散的な位置で得られた画像間のＮ次元類似度値を利用して、前記Ｎ個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法であって、前記画像間のＮ次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、前記Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、前記Ｎ個のＮ次元超平面の交点を求めるステップと、前記交点を前記Ｎ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有するようにすることによって効果的に達成される。
【００２９】
また、本発明は、画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法に関し、本発明の上記目的は、３個の画像間対応パラメータを軸とする３次元類似度空間において、離散的な位置で得られた画像間の３次元類似度値を利用して、前記３個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法であって、前記画像間の３次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する平面を求めるステップと、前記平面を各パラメータ軸に対して３個求めるステップと、前記３個の平面の交点を求めるステップと、前記交点を前記３次元類似度空間における３次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップとを有するようにすることによって効果的に達成される。
【００３０】
更に、本発明は、画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法に関し、本発明の上記目的は、離散的に得られた画像間の２次元類似度値を利用して、連続領域における２次元類似度最大値または最小値の位置を高精度に推定する画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法であって、前記画像間の２次元類似度値が、水平軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（水平極値線ＨＥＬ）を求めるステップと、前記画像間の２次元類似度値が、垂直軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（垂直極値線ＶＥＬ）を求めるステップと、前記水平極値線ＨＥＬと前記垂直極値線ＶＥＬの交点を求めるステップと、前記交点を前記２次元類似度の最大値または最小値を与えるサブピクセル推定位置とするステップとを有するようにすることによって効果的に達成される。
【発明の効果】
【００３１】
本発明によれば、画像間の平行移動、回転、スケール変化を同時に高精度に推定することができる。即ち、本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法によれば、画像間の対応に、平行移動があるときだけでなく、回転やスケール変化があるときにも、繰り返し計算を利用することなく、高精度にこれらの画像間対応パラメータを同時に推定することができる。本発明では、繰り返し計算を利用しないので、従来のように計算ごとに画像を補間する必要が全くない。
【００３２】
また、本発明によれば、テンプレート画像を用いて少ない枚数の補間画像を作成しておき、これらの補間画像を用いた画像間の類似度値を利用して、平行移動と回転やスケール変化を同時に高精度に推定することができる。本発明は、一般的に利用されているようなパラメータ最適化手法と異なり、繰り返し計算を用いないため、計算が容易でしかも計算時間が予測可能で、繰り返し計算を行うときに問題となる収束性を保証する必要もなく、初期値による不安定性も生じない。つまり、本発明によれば、従来方法において生じる収束性の問題や初期値の問題を排除することができる。
【発明を実施するための最良の形態】
【００３３】
以下、本発明の好適な実施例を図面及び数式を参照しながら説明する。
【００３４】
実施例１：
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の実施例１では、画像のサンプリング周期で求まっている２次元類似度を使って、その２次元類似度の最大値又は最小値を与える２次元サブピクセル変位（つまり、水平及び垂直方向の画像間変位）を高精度に同時に推定するようにしている。
【００３５】
ここで、画像を水平方向と垂直方向に独立と考え、変位のサブピクセル推定も独立に行う従来方法を「１次元推定方法」と称するのに対して、本発明の実施例１に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法は、２次元画像から計算した離散位置における２次元類似度を使って、２次元サブピクセル推定を行なっているので、「画像のサブピクセルマッチングにおける高精度２次元推定方法」又は「２次元サブピクセル推定方法」と称し、更に、略して「２次元推定方法」と呼ぶようにしている。
【００３６】
後述する本発明の実施例１に係る２次元推定方法は、領域ベースマッチングを用いた高精度同時推定方法で、従来の「１次元推定方法」と比較して計算量がわずかに増えるだけで、圧倒的に高精度な２次元サブピクセル推定が可能である。
【００３７】
（１）２次元画像モデルと２次元類似度モデル
ここでは、何種類かの画像モデルをとりあげ、その２次元類似度を考える。異なる画像モデルから得られる２次元類似度が、同じ２次元類似度モデルで近似できることを示す。
【００３８】
領域ベースマッチングとは、２枚の画像の類似度を評価し、その最大または最小位置を２枚の画像間の変位として求めることである。この類似度は、通常は画像のサンプリング周期と一致して、とびとびの変位位置に対して値が得られる。画素単位でのマッチングは、これらの類似度の中から、最大値または最小値を探すことに対応する。サブピクセル推定は、これらの類似度を補間して、サブピクセルでの最大値または最小値を与えるサブピクセル変位位置を推定することである。以後の説明では、「最大または最小」を単に「最大」と表現する。ＳＡＤやＳＳＤ類似度を使うときには「最小」、相関類似度を使うときには「最大」を意味する。
【００３９】
（１―１）コーナーモデル
まず、領域ベースマッチングとサブピクセル推定が可能な２次元画像モデルとして、１次元画像モデル
【数７】
を単純に拡張して、下記数８で示すコーナー画像を考える。ここで、σimageは濃度エッジのシャープさである。このコーナー画像を図４(ａ)に示す。
【数８】
ただし、画像座標系をｕ−ｖ、類似度座標系（変位座標系）をｓ−ｔとする。
【００４０】
上記数８で表す画像を参照画像として、また、全く同じ画像を入力画像として２次元マッチングを行ったときのＳＳＤ類似度(本発明では、ＳＳＤ類似度を「自己類似度」と称することもある)は、下記数９で求めることができる。ＳＳＤ類似度は、値が小さいほど類似していることを示すので、この場合には暗い場所ほど類似性が高いことを示している。自己類似度は全く同じ画像を使った類似度なので、変位(ｓ,ｔ)＝(０,０)の位置の類似度が最も小さい。この類似度を図４(ｂ)に示す。
【数９】
ただし、総和範囲は、任意形状の注目範囲であるが、図４(ｂ)では１１×１１の矩形領域で計算した。
【００４１】
図４(ａ)では、コーナーに剪断成分が含まれない、つまり水平方向の濃度変化が垂直位置に依存しないので、２次元画像の水平方向と垂直方向を独立に扱うことができる。このとき、上記数９の類似度も、水平方向と垂直方向を独立に扱うことができる。従って、サブピクセル変位位置を推定するときにも、水平方向と垂直方向を独立に扱うことができる。
【００４２】
しかし、実際の画像では、必ずしも９０度のコーナーで構成されるわけではない。例えば、建築物を地上から撮影したステレオ画像を使って３次元情報を再構成するときには、建築物のコーナー領域を、他方の画像との対応を求めるが、コーナーが９０度に撮影されることはまれである（図２(ａ)参照）。
【００４３】
そこで、ｖ＝０に濃度エッジを持つ２枚の２次元画像
【数１０】
を、それぞれ左にθａ、右にθｂだけ回転した画像を考える。
【数１１】
【００４４】
このθａ、θｂを画像全体の回転角度θｇとコーナー角度θｃを使って、次のように定義する。
【数１２】
【００４５】
このとき、Ｉａ(ｕ,ｖ)とＩｂ(ｕ,ｖ)を使って表す２次元画像モデルＩｃ(ｕ,ｖ)を、下記数１３で定義する。
【数１３】
【００４６】
この２次元画像モデルは、図５に示すように、画像全体の回転角度θｇとコーナー角度θｃと２次元濃度エッジのシャープさσimageをパラメータに持つ。０≦θｇ≦π／４、０＜θｃ≦π／２の範囲を考慮する。これ以外の範囲は、画像の左右を反転することと、濃淡を反転することで表現することができる。θｇ＝π／４、θｃ＝π／２のときには、上記数８の単純な画像モデルと同一である。図６(ａ)、(ｂ)、(ｃ)、(ｇ)、(ｈ)、(ｉ)に、上記数１３で表す画像モデルの例を示す。図６(ｄ)、(ｅ)、(ｆ)、(ｊ)、(ｋ)、(ｌ)に、図６(ａ)、(ｂ)、(ｃ)、(ｇ)、(ｈ)、(ｉ)の２次元画像に対応した類似度を示す。
【００４７】
（１―２）繰り返しテクスチャモデル
類似度に方向依存性がないとき、つまり類似度が等方性なら、サブピクセル推定を水平方向と垂直方向に独立に行っても誤差が発生しない。しかし、類似度に方向依存性があるとき（本発明では「異方性」と呼ぶ）、つまり、方向によって微分値が異なるときには、サブピクセル推定を水平方向と垂直方向に独立に行うと誤差が発生する可能性がある。
上記数１３で示した２次元画像モデルＩｃ(ｕ,ｖ)では、θｃ≠π／２のときに自己類似度に異方性が現れた。しかし、直感的に理解できるように、これ以外にも自己類似度に異方性を生じさせるテクスチャが存在する。この例として、図３(ａ)に示す画像のＳＳＤ自己類似度を図３(ｂ)に示す。テクスチャに含まれるコーナーはθｃ＝π／２なのに、類似度に異方性が現れている。
【００４８】
θｃ＝π／２なのに自己類似度に異方性を生じさせる２次元画像モデルとして、下記数１４を考える。
【数１４】
【００４９】
上記数１４の２次元画像モデルＩｔ(ｕ,ｖ)は、図７に示すように、ｕ方向空間周波数ｆｕ、ｖ方向空間周波数ｆｖ、画像の回転角度θｇをパラメータに持つ。０≦θｇ≦π／２の範囲を考慮する。ｆｕ≠ｆｖのときには、自己類似度に異方性を生じる。自己類似度の異方性は、テクスチャにおける空間周波数の異方性に相当する。前節で導入した２次元画像モデルＩｃ(ｕ,ｖ)は、θｃ≠π／２のときに、空間周波数の異方性が発生していると考えることができる。
【００５０】
この２次元画像モデルとＳＳＤ自己類似度を図８に示す。
【００５１】
（１―３）ガウス関数モデル
自己類似度に異方性を生じさせる２次元画像モデルとして、下記数１５の２次元ガウス関数を考えることができる。
【数１５】
ただし、図９に示すように、σはガウス関数の標準偏差で、ｋは異方性係数(ｋ＞０)で、θｇは回転角度である。０≦θｇ≦π／２の範囲を考慮する。この２次元画像モデルとＳＳＤ自己類似度を図１０に示す。
【００５２】
（１―４）２次元類似度モデル
１次元サブピクセル推定と異なり、２次元サブピクセル推定では２次元画像モデルのパラメータ数が多いので、２次元画像モデルから得られるＳＳＤ類似度を直接検討することは得策ではない。ここで取り上げた画像モデルは、画像の性質の多くを含んでいるが、実際用途で扱う画像は、ここで取り上げた画像モデルを複数組み合わせたものや、より高次の画像が存在する。
【００５３】
そこで、以後の検討では、これまでに検討した何種類かの画像モデルを直接使わず、代わりに、これらの画像モデルから得られる２次元類似度を近似する２次元類似度モデルを検討する。２次元類似度モデルとして、下記数１６で表す２次元ガウス関数を採用する。
【数１６】
ただし、（ｄｓ,ｄｔ)は画像間の真の変位で、σはガウス関数の標準偏差で、ｋは異方性係数(ｋ＞０)で、θｇは回転角度である。０≦θｇ≦π／２の範囲を考慮する。
【００５４】
比較のため、ここで取り上げた２次元画像モデルから得られるＳＳＤ自己類似度の例と上記数１６の類似度の例を図１１に示す。ただし、(ｄｓ,ｄｔ)＝(０,０)、θｇ＝０、ｋ＝１としたときに、２次元画像モデルから計算した離散的な自己類似度を最もよく近似するように、上記数１６のパラメータを数値計算で求めた。
【００５５】
以上の説明は、連続領域で行ってきたが、実際の画像データから得られる類似度は、画素単位でのサンプリングになっている。この様子を図１２に示す。本発明の実施例１に係る２次元推定方法は、このように離散化単位で得られた２次元類似度を用いて、連続領域での類似度最大値を与える２次元変位を高精度に推定するものである。
【００５６】
（２）本発明の実施例１に係る２次元推定方法
（２―１）本発明の実施例１に係る２次元推定方法の原理＜連続領域での説明＞
ここでは、前述した２次元類似度モデルを使って、本発明の実施例１に係る２次元推定方法の原理を説明する。本発明は、離散化単位で得られた類似度値を使って、連続領域での類似度最大値位置を高精度に推定することを特徴としている。ここでは、本発明の実施例１の原理を連続領域で説明する。
【００５７】
図１３は２次元類似度モデルを等高線で表示した図である。この２次元類似度モデルのパラメータは、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１.０、k＝０.７、θｇ＝π／６である。等高線は、２次元類似度モデルの最大値に対して、１０％から９０％までレベルを示している。２次元類似度モデルの等高線は楕円になる。この楕円の中心が、２次元類似度モデルの最大値位置である。図１３(ａ)には、２次元類似度モデルの長軸を示してあるが、この角度が回転角度θｇに相当する。
【００５８】
数１６の２次元類似度モデルをｓで偏微分して０とおくことで、図１３(ｂ)の直線ｓ＝ａｔ＋ｂを表す関係を得ることができる。
【数１７】
【００５９】
考慮する条件がｋ＞０なので、上記数１７の各係数の分母≠０である。また、２次元類似度モデルが等方性のとき、つまりｋ＝１のとき、上記数１７はｓ＝ｄｓとなり、推定誤差が発生しない。本発明の実施例１では、この直線を水平極値線ＨＥＬ(Horizontal Extremal Line)と称する。ＨＥＬは図１３(ｂ)の等高線を表す楕円と水平直線との接点を通る。
【００６０】
一方、数１６の２次元類似度モデルをｔで偏微分して０とおくことで、図１３(ｂ)の直線ｔ＝Ａｓ＋Ｂを表す関係を得ることができる。
【数１８】
【００６１】
数１７と同様に、各係数の分母≠０である。また、ｋ＝１のとき、数１８は、ｔ＝ｄｔとなり、推定誤差が発生しない。本発明の実施例１では、この直線を垂直極値線ＶＥＬ(Vertical Extremal Line)と称する。ＶＥＬは図１３(ｂ)の等高線を表す楕円と垂直直線との接点を通る。
【００６２】
ＨＥＬとＶＥＬの交点は、数１６の２次元類似度モデルを異なる方向に偏微分したときに、どちらの方向にも０になる点であり、これは２次元類似度モデルの類似度値を最大にする位置に他ならない。実際に、この交点を計算すると、
【数１９】
となり、当然のことながら２次元類似度モデルの変位(ｄｓ,ｄｔ)を表している。
【００６３】
このときの分母≠０は、ＨＥＬとＶＥＬが交点を持つための条件である。分母を計算すると、次のようになる。
【数２０】
ｋの範囲はｋ＞０なので、この分母１−ａＡ≠０である。
【００６４】
まとめると、２次元類似度のｓ方向微分とｔ方向微分は、ＨＥＬとＶＥＬの２直線で表すことができ、この２直線の交点が２次元類似度の最大値位置になっている。従って、本発明の実施例１では、離散化単位で得られた類似度値を使ってＨＥＬとＶＥＬを求め、その交点を計算すれば、２次元サブピクセル推定を行うことができる。
【００６５】
（２―２）本発明の実施例１に係る２次元推定方法の実装方法＜離散領域への適用＞
以下では、離散的に算出された画像間の類似度値を利用して、連続領域における類似度最大値位置を推定する本発明の実施例１に係る２次元推定方法を具体的に説明する。
【００６６】
まず、離散的に得られている類似度値を使って、ＨＥＬを求める。ＨＥＬは水平方向について類似度最大値を通る直線なので、２本以上の水平ライン上でのサブピクセル類似度最大値位置が決定すれば、ＨＥＬを求めることができる。この様子を図１４(ａ)に示す。同図において、同心楕円は２次元類似度モデルの等高線を表す。同図中の□位置で離散的に得られた類似度値を使って、直線で示すＨＥＬを求める。
【００６７】
図１４(ａ)において、直線ｔ＝０上での最大類似度を与える位置を
変位位置(ｓ,ｔ)での類似度をＲ(ｓ,ｔ)とすると、
【数２１】
【００６８】
この推定結果には、前章で示したように推定誤差が含まれるが、このことについては次節で述べる。直線ｔ＝−１上と直線ｔ＝１上での最大類似度を与える位置をそれぞれ
とすると、
【数２２】
【数２３】
ただし、ｉｔ＝−1、ｉｔ＝１は、直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置オフセットである（後述する）。
【００６９】
これら３点の最大類似度位置を通る直線がＨＥＬである。実際には、画像パターンや画像に含まれるノイズや画像間の相違などのために、これら３点は完全には直線上には乗らない可能性があるので、最小二乗で近似直線を求める。この３点を最小二乗で近似する直線は、下記数２４で求めることができる。
【数２４】
【００７０】
次に、離散化単位で得られた類似度値を使って、ＶＥＬを求める。ＶＥＬは垂直方向について類似度最大値を通る直線なので、２本以上の垂直ライン上でのサブピクセル類似度最大値位置が決定すれば、ＶＥＬを求めることができる。この様子を、図１４(ｂ)に示す。同図において、同心楕円は、２次元類似度モデルの等高線を表す。同図中の□位置で離散化単位で得られた類似度値を使って、直線で示すＶＥＬを求める。
【００７１】
図１４(ｂ)において、直線ｓ＝０上での最大類似度を与える位置を
とすると、
【数２５】
【００７２】
直線ｓ＝−１上と直線ｓ＝１上での最大類似度を与える位置をそれぞれ
とすると、
【数２６】
【数２７】
ただし、ｉｓ＝−１、ｉｓ＝１は、直線ｓ＝−１上と直線ｓ＝１上で最大類似度を与える整数位置オフセットである。
【００７３】
これら３点の最大類似度位置を通る直線がＶＥＬである。この３点を最小二乗で近似する直線は、下記数２８で求めることができる。
【数２８】
【００７４】
ＨＥＬとＶＥＬの交点、即ち、数２４と数２８の交点が、連続領域における２次元類似度最大値のサブピクセル推定位置
になっている。
【数２９】
【００７５】
ところで、数２６、数２７における直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置オフセットｉｔ＝−１、ｉｔ＝１は、必ずしも０になる保証はない。例えば、図１５(ｂ)に示すのは、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、k＝０.３、θｇ＝π／８の２次元類似度モデルだが、このモデルに対するＨＥＬを求めるための、直線ｔ＝−１上の類似度最大値はｓ＝−２付近にある。従って、
を計算するときには、直線ｔ＝−１上と直線ｔ＝１上で最大類似度を与える整数位置を再探索する必要がある。垂直方向についても同様で
を計算するときには、直線ｓ＝−１上と直線ｓ＝１上で最大類似度を与える整数位置を再探索する必要がある。
【００７６】
このとき、図１５(ｂ)に示す±３の範囲の類似度を使い、−２≦ｉ≦＋２を考慮している。この探索範囲は、多数の現実的な画像をもとに決定した範囲である。もしも、２次元類似度モデルのパラメータ範囲に制限がないとすると、このときの再探索範囲は無限になってしまう。その理由は、(ｄｓ,ｄｔ)＝(０,０)のとき、数１７のＨＥＬについて、下記数３０のＤ(ｋ,θｇ)を考えると、
【数３０】
【００７７】
Ｄ(ｋ,θｇ)が最大値に近づくのは、θｇ→０、ｋ→０のときで、Ｄ(ｋ,θｇ)→cosθｇ／sinθｇ→∞だからである。ところで、ｋ→０のとき、２次元類似度の異方性が無限大になり、ＨＥＬとＶＥＬがほとんど一致する。もしもこのような２次元類似度を作り出す画像があったとしても、そのときには、ほとんど同一になったＨＥＬとＶＥＬ上で類似度最大値を検索すればよい。たとえば、数３０でＤ(ｋ,θｇ)→∞のときには、ＶＥＬ上の３点
における類似度値を求め、パラボラフィッティングを行えばよい。
【００７８】
従来の１次元推定では、図１５(ａ)に示す５個の類似度値を用いてサブピクセル位置を推定していた。このとき、直線ｔ＝０上の３点の類似度(□印)を用いてサブピクセル推定を行って水平方向サブピクセル変位(●印)とした。また、直線ｓ＝０上の３点の類似度(□印)を用いてサブピクセル推定を行って垂直方向サブピクセル変位(●印)とした。図１５(ａ)に示す２直線の交点が従来手法で求めた２次元サブピクセル推定値だが、大きな推定誤差を含むことがわかる。
【００７９】
これに対して本発明の実施例１に係る２次元推定方法は、図１５(ｂ)に示す２５個の類似度値を用いて２次元サブピクセル位置を推定する。ＨＥＬとＶＥＬの交点は、正確に２次元類似度のサブピクセル最大位置を通るので、推定は極めて高精度である。さらに、本発明の実施例１に係る２次元推定方法の計算量の増加は、従来の「１次元推定方法」と比較してわずかである。領域ベースマッチングでは、大部分の計算時間は類似度値を計算しているが、本発明の実施例１に係る２次元推定方法では既に得られている類似度値を利用するので、計算時間の増加は少ない。
【００８０】
（２―３）本発明の実施例１に係る２次元推定方法とサブピクセル推定誤差低減方法との組合せ
数２１から数２３、数２５から数２７において、３位置の類似度を使ってパラボラフィッティングによってサブピクセル位置を推定するときには、推定誤差が発生する。この推定誤差は、補間画像を利用することでキャンセルすることができる。
【００８１】
数２１から数２３、数２５から数２７において求めているサブピクセル位置は、従来どおりの１次元推定にすぎないので、非特許文献１に記載された、推定誤差を低減できるサブピクセル推定方法（以下、本発明の実施例１に係る２次元推定方法と区別するために、非特許文献１に記載されたサブピクセル推定方法を「サブピクセル推定誤差低減方法」と称する）をそのまま適用することができる。ここでは簡単に適用方法を説明する。
【００８２】
マッチングに使う２次元画像関数をＩ１(ｕ,ｖ)、Ｉ２(ｕ,ｖ)とするとき、水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、
【数３１】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０.５,０)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０.５,０)だけ変位していると考えることができるので、その変位を補正した推定結果は、下記数３２、数３３で求めることができる。
【数３２】
【数３３】
【００８３】
数３３も数２１と同様に推定誤差を含むが、その位相が逆になっているので、下記数３４によって推定誤差をキャンセルすることができる。
【数３４】
【００８４】
同様にして、水平ラインｔ＝−１、ｔ＝１について、
【数３５】
【数３６】
【００８５】
ただし、ｉｔ＝−１、ｉｔ＝１は、直線ｔ＝−１上と直線ｔ＝１上でＲｉｕ(ｓ,ｔ)の最大類似度を与える整数位置オフセットである。これらの値は、数２２、数２３とは異なる可能性がある。数３５と数３６を使って、下記数３７、数３８によって推定誤差をキャンセルすることができる。
【数３７】
【数３８】
【００８６】
垂直方向についても同様である。垂直方向補間画像Ｉ１ｉｖ(ｕ,ｖ)は、
【数３９】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。垂直方向補間画像Ｉ１ｉｖ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０,０.５)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０,０.５)だけ変位していると考えることができるので、その変位を補正した推定結果は、下記数４０、数４１で求めることができる。
【数４０】
【数４１】
【００８７】
数４１も数２５と同様に推定誤差を含むが、その位相が逆になっているので、下記数４２によって推定誤差をキャンセルすることができる。
【数４２】
【００８８】
同様にして、垂直ラインｓ＝−１、ｓ＝１について、
【数４３】
【数４４】
【００８９】
ただし、ｉｓ＝−１、ｉｓ＝１は、直線ｓ＝−１上と直線ｓ＝１上でＲｉｖ(ｓ,ｔ)の最大類似度を与える整数位置オフセットである。これらの値は、数２６、数２７とは異なる可能性がある。数４３と数４４を使って、下記数４５、数４６によって推定誤差をキャンセルすることができる。
【数４５】
【数４６】
【００９０】
数３４、数３７、数３８、数４２、数４５、数４６で計算した各ライン上のサブピクセル推定値は、従来の１次元推定方法によって推定されたサブピクセル推定値と比較して、推定誤差が低減されている。従って、サブピクセル推定誤差低減方法によって推定されたこれらの推定結果を使って、更に本発明の実施例１に係る２次元推定方法で２次元サブピクセル推定を行うことで、より高精度な２次元サブピクセル推定が可能になる。
【数４７】
補間画像は、水平方向と垂直方向にそれぞれ１回ずつ作ればよいので、計算量の増加は１次元推定と同程度である。
【００９１】
（３）本発明の実施例１に係る２次元推定方法による推定誤差と従来の推定方法による推定誤差との比較
以下では、２次元画像間の変位推定で、従来一般的に使用されている各種推定方法による推定誤差を本発明の実施例１に係る２次元推定方法による推定誤差と比較する。前述した２次元画像を使用して推定誤差を比較することで、同じ画像に対する比較を行う。
【００９２】
２次元サブピクセル推定誤差は、水平方向と垂直方向の推定誤差を考える必要がある。ところが、推定誤差に影響するのは、２次元画像が持つ３パラメータと画像間の２次元入力変位の合計５パラメータになるので、全てを網羅することは困難である。そこで、代表的な２次元画像を用いて、各方法による推定誤差を比較する。
【００９３】
（３―１）本発明の実施例１に係る２次元推定方法
代表的な２次元画像として、前述したコーナー画像を使用し、パラメータとしてθｇ＝０，π／８，π／４、θｃ＝π／４、σimage＝１.０を考える。０≦θｇ≦π／４以外の２次元画像は、画像の左右を反転することと、濃淡を反転することで表現することができるが、実際の検討は、さらに多くの２次元画像、回転角度、コーナー角度を用いて行っている。他の２次元画像でも、２次元類似度が同じような形になるときには、サブピクセル推定誤差は同じ傾向を示す。
【００９４】
サブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔは、数４７を使って、下記数４８で表すことができる。
【数４８】
【００９５】
サブピクセル推定誤差低減方法を適用しないときの本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔは、数２９を使って、下記数４９で表すことができる。
【数４９】
【００９６】
θｇ＝０,π／８,π／４の３画像に対する数４８、数４９の２次元サブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｃ)(ｄ)、(ｅ)(ｆ)に示す。(ｃ)と(ｄ)は、それぞれサブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。(ｅ)と(ｆ)は、それぞれサブピクセル推定誤差低減方法を適用しないときの本発明の実施例１に係る２次元サブピクセル推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【００９７】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じているにもかかわらず、サブピクセル推定誤差は極めて小さい。
【００９８】
（３―２）従来の１次元サブピクセル推定方法
従来の１次元推定方法は、数２１と数２５で求めることができる。また、非特許文献１に記載されたサブピクセル推定誤差低減方法を適用した推定値は、数３４と数４２で求めることができる。ここでは、数２１と数２５を使ったときの推定誤差を比較検討する。
【数５０】
【００９９】
θｇ＝０,π／８,π／４の３画像に対する上記数５０のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｇ)、(ｈ)に示す。(ｇ)と(ｈ)は、それぞれ従来の１次元推定方法における水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１００】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、従来の１次元推定方法でも推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このときには大きな推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１０１】
（３―３）従来の濃度こう配法によるサブピクセル推定
濃度こう配法では、画像間における対応位置に濃度変化がないことを仮定して、水平方向移動量と垂直方向移動量の２つの未知数を含む拘束条件を画素ごとに得る。したがって、この２つの未知数はこのままでは求めることができないため、注目領域を設定し、注目領域中での拘束条件の２乗和を最小にするように、繰り返し計算によって未知数を求める。このようにして得られる変位量には画素単位という制約がなく、サブピクセル単位を得ることができる。
【０１０２】
一方、画像補間手法では、注目領域と周囲の探索領域を離散化単位よりも密に補間する。密に補間した単位で、類似度最大値を探索する。このときには、濃度こう配法と異なり、画像間の移動量に対する制限はない。
【０１０３】
この２種類のサブピクセル変位推定方法は、全く異なる実装に見え、異なる手法として扱われていたが、本質的に同一のものであることが示されている（非特許文献５を参照）。ここでは、補間画像を作るときに用いる補間関数の次数によってサブピクセル変位推定誤差が変化する様子を確認する。そこで、１次補間を利用した補間画像によるサブピクセル推定結果を、濃度こう配法による結果と見なす。さらに、高次の補間関数を利用した補間画像によるサブピクセル推定結果も検討する。
【０１０４】
（３―４）従来のバイリニア画像補間によるサブピクセル推定
サブピクセル変位を求める２枚の画像をＩ１(ｕ,ｖ)、Ｉ２(ｕ,ｖ)とする。これらの画像は離散化サンプリングされているとする。すなわち、ｕ,ｖは整数である。また、画像間の変位を(ｄｓ,ｄｔ)としたときに、０≦ｄｓ≦１、０≦ｄｔ≦１とする。画像間の真の変位がこの範囲にないときには、画素単位のマッチングによって画像全体を移動する。
【０１０５】
バイリニア補間(２方向１次補間)を利用したサブピクセル推定値
は、下記数５１で表すことができる。
【数５１】
ただし、総和範囲は任意形状の注目領域、argminは値を最小にするような
の組を求める演算を表す。
【０１０６】
θｇ＝０,π／８,π／４の３画像に対する上記数５１のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｉ)、(ｊ)に示す。(ｉ)と(ｊ)は、それぞれ水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１０７】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このとき推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１０８】
この推定誤差は、図１６、図１７、図１８の(ｃ)、(ｄ)の結果よりも大きい。すなわち、本発明の実施例１に係る２次元推定方法は、従来のバイリニア補間によるサブピクセル推定方法よりも高精度に変位を推定することができる。バイリニア補間によるサブピクセル推定方法は濃度こう配法と等価なので、本発明の実施例１に係る２次元推定方法は、従来の濃度こう配法よりも高精度にサブピクセル変位を推定することができるとも言える。
【０１０９】
また、濃度こう配法に相当する推定誤差は、図１６、図１７、図１８の(ｇ)、(ｈ)の１次元推定方法の結果よりもはるかに小さい。すなわち、従来広く利用されてきた、領域ベースマッチングを画素単位で行い、サブピクセル変位は類似度を利用して１次元推定する手法より、濃度こう配法ははるかに高精度に２次元サブピクセル変位を推定することができる。このため、領域ベースマッチングよりも濃度こう配法の方が高精度という認識があったが、本発明の実施例１に係る２次元推定方法を用いることで、濃度こう配法よりもさらに高精度にサブピクセル変位を推定することができる。
【０１１０】
（３―５）従来のバイキュービック画像補間によるサブピクセル推定
補間画像を作るときに、より高次の項まで考慮することで、高精度な補間画像を作ることができる。バイキュービック補間(２方向３次補間)は、補間したいサブピクセル画素位置の周囲４×４の画素値を使う補間方法である。バイキュービック補間を利用したサブピクセル推定値
は、下記数５２で表すことができる。
【０１１１】
【数５２】
ただし、総和範囲は任意形状の注目領域、argminは値を最小にするような
の組を求める演算を表す。このとき、重み関数ｈ(ｘ)は、sinc関数をもとに各種提案されていて、下記数５３が多く用いられている。
【０１１２】
【数５３】
ただし、ａは調整パラメータで、ａ＝−１やａ＝−０.５がよく用いられる。
【０１１３】
また、下記数５４も多く用いられる。
【数５４】
Ｂ、Ｃは調整パラメータで、Ｂ＝１、Ｃ＝０のときにはキュービックＢスプラインを近似し、Ｂ＝０、Ｃ＝０.５のときにはCatmull−Romキュービック特性になる。このように、重み係数の選び方によって補間特性が異なる。ここでは、数５３の重み関数において、パラメータａ＝−０.５５を採用した。ただし、画像処理の教科書等に最も多く紹介されているのはａ＝−１である。両者の特性の差については後述するが、ａ＝−０.５５はサブピクセル推定誤差を小さくするように数値的に決定した。
【０１１４】
θｇ＝０,π／８,π／４の３画像に対する数５２のサブピクセル推定誤差を、それぞれ図１６、図１７、図１８の(ｋ)、(ｌ)に示す。(ｋ)と(ｌ)は、それぞれ水平方向推定誤差errorｓと垂直方向推定誤差errorｔを示す。
【０１１５】
図１６は、θｇ＝０の画像を使った結果を示すが、このときには２次元類似度に異方性が生じない。このときには、推定誤差が生じない。図１７と図１８は、それぞれθｇ＝π／８とπ／４の画像を使った結果を示し、このときには２次元類似度に異方性が生じるが、このとき推定誤差が生じる。推定誤差の大きさは、異方性を持つ２次元類似度の回転角度と真の画像間変位に依存する。
【０１１６】
この推定誤差は、図１６、図１７、図１８の２次元推定(ｃ)、(ｄ)の結果と同程度である。すなわち、本発明の実施例１に係る２次元サブピクセル推定方法は、最適なパラメータを選択したときの高次補間画像を使ったサブピクセル推定と同程度の高精度であるということができる。
【０１１７】
ところで、数５３におけるパラメータａについて、具体的な計算例を示して説明する。数５３は、sinc関数を有限範囲で近似した重み係数と考えることができる。パラメータａは、その近似特性を調節するものである。図１９(ａ)、(ｂ)に、ａ＝−０.５５のときのｈ(ｘ)と、整数位置に値を持つ１次元関数ｆ(ｉ)の補間結果を示す。ｆ(ｉ)は、数７でσimage＝０.７としたものである。また、図１９(ｃ)、(ｄ)に、ａ＝−１.０のときのｈ(ｘ)と、整数位置に値を持つ１次元関数ｆ(ｉ)の補間結果を示す。
【０１１８】
ａ＝−１.０は、多くの画像処理の教科書などでバイキュービック補間として紹介されるパラメータである。sinc関数を良く近似しているが、少なくともここで検討した関数ｆ(ｉ)を良好に補間しているとは言い難い。ａ＝−０.５５は、sinc関数の近似としてはそれほど良くないが、関数ｆ(ｉ)を良好に補間している。
【０１１９】
図１６から図１８では、パラメータａ＝−０.５５を用いたが、この値はサブピクセル推定誤差を小さくするように数値的に求めた結果である。ａ＝−０.５を用いるとバイリニア補間の結果と同程度、ａ＝−１.０を用いるとバイリニア補間を用いたサブピクセル推定結果よりも誤差が大きくなる。
【０１２０】
（３―６）各方法による推定誤差の総合比較
以下では、２次元画像のパラメータをより広範囲に変更して、各手法のサブピクセル推定誤差を比較検討する。２次元画像は、前述したコーナー画像を使用する。考慮するパラメータとして、回転角度をθｇ＝０,π／１６，２π／１６，３π／１６，４π／１６、コーナー角度をθｃ＝π／４，π／３を考える。合計１０種類の２次元画像に対する推定誤差を検討する。また、σimage＝１.０を考える。
【０１２１】
画像間変位は、−０.５≦ｄｓ≦＋０.５、−０.５≦ｄｔ≦＋０.５の範囲で、０.１[画素]ごとに与える。このとき、２次元推定誤差の大きさ、
のＲＭＳ誤差を評価した。
【０１２２】
評価したサブピクセル推定方法は、前述した各方法の他に、従来の１次元推定方法に対するサブピクセル推定誤差低減方法を適用した推定結果も評価した。したがって、
・サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)
・サブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)
・サブピクセル推定誤差低減方法を用いた１次元推定(１Ｄ−ＥＥＣ)
・サブピクセル推定誤差低減方法を用いない１次元推定(１Ｄ)
・バイリニア補間を用いた２次元推定(Ｂｉ−Ｌｉｎｅａｒ)
・バイキュービック補間を用いた２次元推定(ａ＝−０.５５)(Ｂｉ−Ｃｕｂｉｃ)
の６種類の推定方法の推定誤差を評価した。評価した結果を表１に示す。
【０１２３】
【表１】
【０１２４】
表１から次のことが分かった。
まず、サブピクセル推定誤差低減方法を用いた１次元推定(１Ｄ−ＥＥＣ)は、回転角度θｇ＝０のときには、従来の１次元推定(１Ｄ)に対して推定誤差を低減できる。しかし、回転角度θｇ≠０のときには差が現れない。回転角度θｇは、異方性がある類似度の回転角度に対応する。
【０１２５】
次に、バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)やバイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定は、１次元推定(１Ｄ)に対して推定誤差が２桁程度小さい。この差は歴然である。バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)を用いた２次元推定は、濃度こう配法に対応する。濃度こう配法が高精度にサブピクセル推定ができると考えられてきた理由はここにある。
【０１２６】
それから、バイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定は、バイリニア補間(Ｂｉ−Ｌｉｎｅａｒ)を用いた２次元推定よりも高精度である。しかし、重み関数のパラメータによって結果は異なる。重み関数のパラメータを調整するということは、補間フィルタの特性を調整することに対応する。すなわち、画像に対して適切な補間フィルタを選択することで、画像補間を用いたサブピクセル推定手法では精度を向上させることができる。
【０１２７】
最後に、サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)は、バイキュービック補間(Ｂｉ−Ｃｕｂｉｃ)を用いた２次元推定よりも推定誤差が小さく、しかも、パラメータ調整の必要がない。サブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)は、サブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)よりも計算コストが大きい。そこで、計算コストを重視して推定結果は濃度こう配法と同程度でよいときにはサブピクセル推定誤差低減方法を用いない本発明の実施例１に係る２次元推定(２Ｄ)を利用し、さらに高精度な推定結果を得たいときにはサブピクセル推定誤差低減方法を用いた本発明の実施例１に係る２次元推定(２Ｄ−ＥＥＣ)を利用する、という使い分けも可能である。
【０１２８】
（４）本発明の実施例１に係る２次元推定方法を実装したコンピュータ・プログラムによる推定実験
本発明の実施例１に係る２次元推定方法を実装したコンピュータ・プログラムを、ＣＰＵを備えた情報処理端末（例えば、パソコン、ワークステーション等のコンピュータ）に実行させることによって、合成画像及び実画像を用いて、２次元サブピクセル推定実験を行った。なお、本発明でいうコンピュータ・プログラムは、便宜上プログラムと略して称することもある。
【０１２９】
（４―１）合成画像を使った推定実験
前述したコーナー画像を作り、２次元サブピクセル推定実験を行った。作成した合成画像は、８bitモノクロ画像で、画像サイズが６４×６４[画素]で、注目範囲が１１×１１[画素]である。画像にノイズは含まれない。画像のパラメータは、σimage＝１.０、コーナー角度θｃ＝π／４、回転角度θｇ＝０,π／８,π／４である。画像間変位は、−０.５≦ｄｓ≦＋０.５、−０.５≦ｄｔ≦＋０.５の範囲で、０.１[画素]ごとに与えた。
【０１３０】
図２０に従来の１次元推定と、サブピクセル推定誤差低減方法を適用した本発明の実施例１に係る２次元推定についての結果を示す。図２０において、計算結果をメッシュで、実験結果を●で表す。実験結果と計算結果との相違の原因は、画像が８bit量子化階調（２５６階調）であることが考えられる。また、大きく段差がある部分は、整数での最大類似度を探索した結果が隣の画素位置に移動してしまったためと考えることができる。このように隣の位置に移動してしまうことは、従来の１次元推定方法の根本的な欠点である。
【０１３１】
（４―２）実画像を使った推定実験
マッチングに利用する画像から得られる２次元類似度の回転角度θｇ≒０のときには、従来の１次元推定方法でも推定誤差は発生しないが、θｇ≠０のときには、大きな推定誤差が生じる可能性がある。そこで、マッチングパターンとして２種類のパターンを用いたターゲットトラッキング実験を行った。
【０１３２】
使用したマッチングパターンは、円形パターンと傾いたエッジパターンである。これらのパターンを図２１(ａ)、(ｅ)に示す。円形パターンの大きさは、直径が約４１[画素]である。また、これらのパターンのＳＳＤ自己類似度を図２１(ｂ)、(ｆ)に示す。
【０１３３】
円形パターンの自己類似度は等方性なので、従来の１次元推定方法と本発明の実施例１に係る２次元推定方法のどちらでも、推定誤差が小さいことが予想できる。一方、傾いたエッジパターンに対応する自己類似度は異方性でかつ傾いているので、従来の１次元推定方法では推定誤差が大きく、本発明の実施例１に係る２次元推定方法では推定誤差が小さいことが予想できる。
【０１３４】
マッチングパターンを厚紙に印刷して、トラッキングターゲットとした。ターゲットは、ネジ送り式リニアステージ上を直線的に移動する。ターゲットの移動を約２５０枚の時系列画像で撮影し、最初の画像を参照パターンとして使い、以後の画像中からターゲットをトラッキングした。注目領域サイズは６０×６０[画素]（図２１(ａ)、(ｅ)中の黒い矩形領域）で、探索領域は移動予測位置に対して±８(１７×１７)である。ターゲットは画像上で右下方向にゆっくり移動する。マッチングが正しくでき、サブピクセル推定に誤差がなければ、計測軌跡は直線になる。従来の１次元推定方法と本発明の実施例１に係る２次元推定方法のどちらにもＳＳＤ自己類似度を使用し、サブピクセル推定誤差低減方法を適用したパラボラフィッティングによるサブピクセル推定を行った。
【０１３５】
図２１(ｃ)、(ｄ)に、トラッキングパターンとして円形パターン使ったときの計測結果を示す。図２１(ｃ)は従来の１次元推定方法、図２１(ｄ)は本発明の実施例１に係る２次元推定方法を使った結果である。どちらの場合も、サブピクセル推定誤差が小さく、軌跡は良好な直線になっている。
【０１３６】
図２１(ｇ)、(ｈ)に、トラッキングパターンとして傾いたエッジパターン使ったときの計測結果を示す。図２１(ｇ)は従来の１次元推定方法、(ｈ)は本発明の実施例１に係る２次元推定方法を使った結果である。従来の１次元推定方法では、非常に大きなサブピクセル推定誤差が発生していることが明らかである。
【０１３７】
本発明の実施例１に係る２次元推定方法では、トラッキングパターンによらず推定誤差が小さいが、円形パターンのときよりも計測軌跡が直線から外れている。この理由は、２次元類似度が２次元ガウスモデルと異なるため、ＨＥＬとＶＥＬを直線近似するときに誤差が含まれることが原因と考えられる。
【０１３８】
現実的には、このようなトラッキングパターンを利用することはまれである。しかし、この実験は、従来の１次元推定方法を使うと、マッチングパターンによっては、予想外に大きなサブピクセル推定誤差が発生する可能性があることを示している。この実験では、対象が直線上を移動していることがわかっているので、サブピクセル推定誤差を確認することができた。しかし、たとえば、異なる方向から撮影した複数の画像を使って、建築物の３次元情報を再構築するときなどには、傾いたエッジパターンの領域をマッチング領域として利用することもある。従来の１次元推定方法では、このようなときに大きな推定誤差が発生するという問題点があった。
【０１３９】
なお、本発明の実施例１では、水平極値線ＨＥＬと垂直極値線ＶＥＬを求める際、３点を最小二乗で近似する直線を求める代わりに、３点を通る２次曲線として求めることにより、直線近似による誤差を低減するように本発明の実施例１に係る２次元推定方法を拡張することもできる。
【０１４０】
実施例２及び実施例３：
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の実施例２及び実施例３を次のように説明する。本発明の実施例２及び実施例３の着眼点は以下の通りである。
【０１４１】
（Ａ）平行移動と回転、または、平行移動と回転とスケール変化を画像間対応パラメータ（本発明では、略してパラメータとも呼んでいる）にして、これらのパラメータを軸とする多次元空間（本発明では、パラメータ空間或いは多次元類似度空間とも称する）において離散的な位置で得られた類似度を利用する。類似度は離散的な位置で計算しておけばよいので、繰り返し計算の必要がない。
【０１４２】
（Ｂ）離散的な位置で得られた類似度は、多次元空間（つまり、パラメータ空間）における連続関数である類似度をサンプリングしたものであると考える。
【０１４３】
（Ｃ）現実的な仮定の下に類似度をモデル化したとき、多次元空間における各軸（つまり、パラメータ空間における各パラメータ軸）で類似度を偏微分して０とおくことで得られる超平面は、類似度の最大値または最小値を通る。例えば、類似性評価関数にＳＡＤやＳＳＤが採用された場合に、超平面は類似度の最小値を通る。また、類似度評価関数にＺＮＣＣが採用された場合に、超平面は類似度の最大値を通る。
【０１４４】
（Ｄ）このような超平面は、多次元空間の軸の数（つまり、パラメータ数）だけ求めることができる。これらの超平面の交点座標は、多次元空間における類似度の最大値または最小値を与えるパラメータと一致するため、これらの超平面の交点座標を求めることで、多次元空間における類似度の最大値または最小値を与えるパラメータを同時に求めることができる。
【０１４５】
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法は、サンプリンググリッドで離散的に得られた評価値を使って、高精度なサブサンプリンググリッド精度の評価値最大位置または最小位置を推定する方法と考えることができる。
【０１４６】
画像間の対応を平行移動のみに限定して考えるとき（つまり、本発明の実施例１の場合）には、このサンプリンググリッドは画像を構成する画素に相当する。このため、画素単位以上に高精度な平行移動量の表現として、サブピクセルという表現が適切だった。従って、前述したように、本発明の実施例１に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を「２次元サブピクセル推定方法」とも呼んでいる。
【０１４７】
しかし、平行移動と回転、または、平行移動と回転とスケール変化を画像間対応パラメータにする場合の本発明、つまり、本発明の実施例２及び実施例３に係る多パラメータ高精度同時推定方法では、画素やサブピクセルといった表現は一般性を欠いているので、以後ではそれぞれサンプリンググリッド、サブサンプリンググリッドという表現を採用する。
【０１４８】
＜１＞実施例２に係る３パラメータ同時推定方法
＜１−１＞３パラメータ同時推定方法の原理
ここでは、本発明の実施例２に係る３パラメータ同時推定方法の原理を説明する。本発明では、画像間の対応を探索するときには、何らかの類似度評価関数を利用して、類似度評価関数の最小または最大を探索する。このとき得られる類似度は、用いる類似度評価関数によっても異なるが、それ以上に入力として扱う画像によって異なる。
【０１４９】
ここでは、具体的な入力画像や類似度評価関数には触れず、画像間の類似度を次の３次元ガウス関数で近似して、これを３次元類似度モデルとする。
【数５５】
ただし、(ｓ１,ｓ２,ｓ３)は探索パラメータ（つまり、画像間対応パラメータ）で、σはガウス関数の標準偏差で、(ｋ２,ｋ３)は異方性係数(ｋ２,ｋ３＞０)である。また、(ｔ１,ｔ２,ｔ３)は、次のように(ｓ１,ｓ２,ｓ３)に対して回転と平行移動を行った結果である。
【０１５０】
【数５６】
ただし、(ｄ１,ｄ２,ｄ３)は探索パラメータ(ｓ１,ｓ２,ｓ３)の正解値で、(θ１,θ２,θ３)はそれぞれ(ｓ１,ｓ２,ｓ３)の各軸回りの回転角度で、０≦θ１,θ２,θ３≦π／２の範囲を考慮する。
【０１５１】
この３次元類似度モデルをｓ１,ｓ２,ｓ３で偏微分して０とおくと、次の３平面を得ることができる。
【数５７】
【数５８】
【数５９】
【０１５２】
Π１,Π２,Π３は、３次元空間内（つまり、ｓ１,ｓ２,ｓ３を軸とするパラメータ空間）の３平面を構成している。この３平面の交点は、明らかに(ｄ１,ｄ２,ｄ３)であり、正解パラメータと一致する。３平面の交点が求まらない条件は、３次元類似度モデルが３次元未満に縮退した場合である。このような条件では、ｋ２,ｋ３＝０または∞となり、３次元類似度モデルの条件と一致しない。
【０１５３】
＜１−２＞３パラメータ同時推定方法の実装
以下では、３パラメータ同時推定方法において、サブサンプリンググリッド推定値を求める具体的な方法を説明する。例えば、この３パラメータは、平行移動量(ｄｓ,ｄｔ)と回転角度θと考えることができる。
【０１５４】
前節『３パラメータ同時推定方法の原理』で示したように、パラメータ空間における類似度値を各パラメータ軸に関して微分して０とおくことで、各パラメータの最大値位置または最小値位置を通る３平面を作ることができる。従って、サンプリンググリッドにおいて得られた類似度値を用いて、この３平面を推定することができれば、各パラメータの類似度値の最大位置または最小位置をサブサンプリンググリッドにおいて求めることができる。
【０１５５】
まず、平面Π１の求め方を説明する。平面Π１上の点は、ｓ１軸に沿って類似度を微分した結果が０になる点だから、ｓ１軸に平行な直線上の類似度値を使ってサブサンプリング位置を推定することで、平面Π１上の点をいくつか求めることができる。図２２はパラメータ軸ｓ１に関する最大値を通る平面Π１を説明するための模式図である。図２２では、類似度値が得られているサンプリング位置を正方形で、ｓ１軸に平行な直線上で求めたサブサンプリング推定位置を黒丸で示している。類似度値のグリッド単位最大値位置は、(ｒ１,ｒ２,ｒ３)である。
【０１５６】
(ｓ１,ｓ２,ｓ３)における類似度値をＲ(ｓ１,ｓ２,ｓ３)とすると、Ｒ(ｒ１＋ｉ１,ｒ２＋ｉ２,ｒ３＋ｉ３)(ｉ１,ｉ２,ｉ３＝−１,０,１)の２７点の類似度値を使って、平面Π１上の９個の点(ｐ１(ｒ２＋ｉ２,ｒ３＋ｉ３),ｒ２＋ｉ２,ｒ３＋ｉ３)を推定することができる。推定には、次のパラボラフィッティングを用いることができる。
【数６０】
【０１５７】
また、後述するように、パラボラフィッティングによる推定誤差を低減する手法を適用することもできる。
【０１５８】
数６０のパラボラフィッティングを行うときには、Ｒ(ｒ１,ｒ２＋ｉ２,ｒ３＋ｉ３)の類似度値が最大であることを前提としているが、画像パターンによってはこの前提が成立しない場合もある。たとえば、画像に含まれるテクスチャのせん断変形が極端なときには、ｉ２＝ｉ３＝±１の位置でＲ(ｒ１,ｒ２＋ｉ２,ｒ３＋ｉ３)が最大値にならないことがある。このときにも(ｒ２＋ｉ２,ｒ３＋ｉ３)を通りパラメータ軸ｓ１に平行な直線上でサブサンプリンググリッド推定位置を求める必要がある。この場合には、この直線上のサンプリンググリッド上に求まっている類似度値の最大値位置を探索し、修正最大値位置ｒ１’に関して数６０のパラボラフィッティングを行えばよい。
【０１５９】
以後の説明では、数式表記を簡略化するために、ｓ１−ｒ１→ｓ１,ｓ２−ｒ２→ｓ２,ｓ３−ｒ３→ｓ３のように各パラメータ軸に沿って(ｒ１,ｒ２,ｒ３)だけ平行移動して(ｒ１,ｒ２,ｒ３)をパラメータ軸の原点位置と考える。このような平行移動は、グリッド単位での類似度最大位置を原点に移動することと等しく、本発明に関する一般性を失うものではない。本発明に係る多パラメータ高精度同時推定方法によって求めたサブサンプリンググリッド位置を最終的に(ｒ１,ｒ２,ｒ３)と加算することで、平行移動する前のサブサンプリンググリッド位置を簡単に計算することができる。
【０１６０】
平面Π１上にあると期待される９点のサブサンプリンググリッド位置は、実際の画像ではノイズや画像間の変形の影響で、完全には平面上にない可能性がある。このため、これらの９点を最小二乗法によって平面Π１にあてはめる。具体的には、９点と平面Π１とのｓ１軸の沿った距離の２乗和を最小にするように係数を決めればよい。平面Π１は、次のように求めることができる。
【０１６１】
【数６１】
となる。同様にして、平面Π２,Π３は、次のように求めることができる。
【０１６２】
【数６２】
【数６３】
【０１６３】
平面Π１,Π２,Π３の交点
が、サブサンプリング推定位置で、次のように求めることができる。
【数６４】
【数６５】
【０１６４】
＜２＞実施例３に係るＮパラメータ同時推定方法
＜２−１＞Ｎパラメータ同時推定方法の原理
ここでは、本発明の実施例３に係るＮパラメータ同時推定方法の原理を説明する。パラメータ数が３より多いときにも、本発明では同様の考え方でサブサンプリンググリッド高精度同時推定をすることができる。ｓ１からｓＮまでのＮ個のパラメータによって類似度値が決定するときに、類似度値を次のようにガウス関数のｎ個の乗算として表現する。
【０１６５】
【数６６】
ただし、
はＮ次元類似度空間における座標を表すＮ次元ベクトルで、
はＮ次元空間での回転と平行移動を表す行列で、
はガウス関数の異方性係数ベクトル(Ｎ−１次元)である。
【０１６６】
このとき、数６６で表されるＮ次元類似度モデルをｓｉ(ｉ＝１,２,…,Ｎ)で偏微分して０とおくことで、次のようなＮ個のＮ次元関係式を得ることができる。
【数６７】
【０１６７】
これらの関係式を連立させて解くことで、Ｎ次元類似度のサブサンプリンググリッド推定値を得ることができる。つまり、それらの関係式の解は、Ｎ次元類似度のサブサンプリンググリッド推定値になるわけである。なお、Ｎ次元類似度モデルが数６６のガウス関数乗算モデル近似からずれると、数６７は完全な超平面にはならないが、このときには、最小二乗法などで超平面に近似すればよい。
【０１６８】
＜２−２＞Ｎパラメータ同時推定方法の実装
以下では、Ｎパラメータのサブサンプリンググリッド同時推定を実際に行う具体的な方法を説明する。まず、サブサンプリンググリッド同時推定を行うときに、前提としてサンプリンググリッドでの類似度最大位置はわかっており、その(変位)位置を
とする。
【０１６９】
類似度値Ｒをパラメータ軸ｓｉ(ｉ＝１,２,…,Ｎ)について偏微分した結果が０になるＮ次元超平面Πｉ上には、ｓｉパラメータ軸と平行な直線上の類似度値
の３つの値をパラボラフィッティングによってサブサンプリンググリッド推定した３Ｎ−１個の位置、
が存在する。ただし、
【数６８】
である。
【０１７０】
ここで、
は、ｉ番目の要素だけが値ｋ(ｋ＝−１,０,＋１)であり、他の要素は全て値０であるようなＮ次元ベクトルである。たとえば、
である。また、
は、ｉ番目の要素が０で、他の各要素は、−１,０,１のいずれかの値をとるようなｊ番目(ｊ＝１,２,…,３Ｎ−１)のＮ次元ベクトルである。たとえば、
である。
【０１７１】
これらの３Ｎ−１個の位置、
は、Ｎ次元超平面Πｉ上に存在するか、少なくとも最小二乗の意味でＮ次元超平面Πｉにフィットする。従って、Ｎ次元超平面Πｉを
【数６９】
と表すと、３Ｎ−１個の位置、
に対して、
【数７０】
の関係が得られる。これを、次のように書き換える。
【０１７２】
【数７１】
ただし、
は、３Ｎ−１行Ｎ列の行列で、
は要素数Ｎのベクトルである。Ｎ次元超平面Πｉを表す係数
は、一般化逆行列を用いた最小二乗法によって、次のように求めることができる。
【０１７３】
【数７２】
このようにして求めたＮ個のＮ次元超平面Πｉ(ｉ＝１,２,…,Ｎ)の交点として、Ｎパラメータのサブサンプリンググリッド推定値を得ることができる。この交点は、次のように求めることができる。
【０１７４】
【数７３】
【数７４】
【０１７５】
以上の方法でＮパラメータの同時推定を行うことができるが、実際にはＮ次元超平面Πｉの式を、同じＮ個の未知数を含む次のような形式に表現してもよい。
【数７５】
【０１７６】
その理由は、Ｎ次元類似度関数がＮ次元対称になっても安定に超平面の交点が存在するようにするためである。
【０１７７】
また、Ｎ次元超平面Πｉ上の３Ｎ−１個の位置を求めることでＮ次元超平面の係数を求めたが、この個数を減らすこともできる。たとえば、Ｎ＝３のときには３Ｎ−１＝９だが、たとえばｊ＝３のときには、具体的に次のようになっている。
(−１,−１,０) (０,−１,０) (＋１,−１,０)
(−１,０,０) (０,０,０) (＋１,０,０)
(−１,＋１,０) (０,＋１,０) (＋１,＋１,０)
【０１７８】
この中から、次のような組を採用することで、個数を２(Ｎ−１)＋１個にすることができる。
(０,−１,０)
(−１,０,０) (０,０,０) (＋１,０,０)
(０,＋１,０)
【０１７９】
＜３＞本発明による実験結果と従来手法による実験結果との比較
図２３に示す時系列画像を人工的に作成した。この時系列画像は、濃淡で表現した２次元ガウス関数で構成されている。この実験で用いた２次元ガウス関数には方向性(異方性) があり、しかもその方向性が０または９０度以外の方向になるように設定されている。このため、この画像から計算する類似度にも異方性があり、かつ回転角度が０または９０度以外になる。
【０１８０】
図２３に示す画像は時系列画像の最初のフレームで、以後のフレームは最初のフレームに対して位置が変化せず回転角度だけがフレームごとに０.１度ずつ変化する。時系列画像は全部で３１フレームあり、最終フレームは先頭フレームに対して３度回転した画像になっている。これに対して画像の位置は回転中心に対して移動しない。
【０１８１】
図２４に、最初のフレームに対する以後のフレームの回転角度を計測した結果を示す。＋印で示す従来方法は、画像を６×６＝３６個の領域に分割し、各小領域ごとに平行移動量をマッチングによって計測し、得られた平行移動量を使って回転角度を推定した結果である。○印で示す本発明による方法は、従来方法で計測した回転角度を初期推定値として本発明を適用し、さらに高精度に回転角度を計測した結果である。フレームごとに０.１度ずつ回転角度が変化する状況が正確に計測できている。
【０１８２】
図２５に、最初のフレームに対する以後のフレームの位置を計測した結果を示す。＋印で示す従来方法は、画像を６×６＝３６個の領域に分割し、各小領域ごとに平行移動量をマッチングによって計測し、得られた平行移動量の平均として位置を推定した結果である。○印で示す本発明による方法は、従来方法で計測した位置を初期推定値として本発明を適用し、さらに高精度に位置を計測した結果である。位置が変化しない状況が正確に計測できている。
【０１８３】
＜４＞本発明の他の実施例
本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法の他の実施例を以下のように説明する。
【０１８４】
＜４−１＞サブピクセル推定誤差低減方法との組み合わせ
数６０でサブサンプリンググリッド推定位置を求めるときに、３位置の類似度を使ってパラボラフィッティングを用いているので、この推定値には固有の推定誤差を含んでいる。この推定誤差は、補間画像を利用することでキャンセルすることができる。ここでは、２パラメータのときに画像間の水平方向変位に関するサブピクセル位置を推定するときを例に説明する。
【０１８５】
マッチングに使う２次元画像関数をＩ１(ｕ,ｖ),Ｉ２(ｕ,ｖ)とするとき、ＳＳＤとパラボラフィッティングによって推定する水平方向サブピクセル位置は、次のように表すことができる。
【数７６】
【数７７】
【０１８６】
一方、画像Ｉ１(ｕ,ｖ)の水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、
【数７８】
と表すことができる。ただし、このときの(ｕ,ｖ)は整数である。
【０１８７】
水平方向補間画像Ｉ１ｉｕ(ｕ,ｖ)は、Ｉ１(ｕ,ｖ)に対して(０.５,０)だけ変位していると考えることができる。この補間画像を使って、次に示すＳＳＤ類似度でサブピクセル推定を行うと、推定結果も(０.５,０)だけ変位していると考えることができるので、その変位を補正した推定結果は、次式で求めることができる。
【数７９】
【数８０】
【０１８８】
数７８も数６０と同様に推定誤差を含むが、その位相が逆になっているので、次式によって推定誤差をキャンセルすることができる。
【数８１】
【０１８９】
＜４−２＞アフィンパラメータ推定
画像間の対応がアフィン変換で表現できるとき、下記数８２の変形モデル中の６個のパラメータを決定する必要がある。
【０１９０】
【数８２】
ただし、(ｕ１,ｖ１),(ｕ２,ｖ２)は、それぞれ画像Ｉ１,Ｉ２を表す座標で、画像Ｉ２を変形して画像Ｉ１に合わせている。
【０１９１】
離散的な位置で得られた類似度値を用いて、この６つのパラメータを同時に高精度に推定することができる。
【０１９２】
＜４−３＞射影パラメータ推定
画像間の対応が射影変換で表現できるとき、下記数８３の変形モデル中の８個のパラメータを決定する必要がある。
【０１９３】
【数８３】
ただし、(ｕ１,ｖ１),(ｕ２,ｖ２)は、それぞれ画像Ｉ１,Ｉ２を表す座標で、画像Ｉ２を変形して画像Ｉ１に合わせている。
【０１９４】
離散的な位置で得られた類似度値を用いて、この８つのパラメータを同時に高精度に推定することができる。
【０１９５】
＜４−４＞高速探索手法との組み合わせ
本発明によるサブサンプリンググリッド推定方法は、サンプリンググリッド単位での対応が正しく得られた後の高精度化処理と考えることもできる。したがって、サンプリンググリッド単位での探索には、すでに提案されている各種の高速化手法を用いることができる。
【０１９６】
例えば、画像の濃度情報を利用した高速探索手法としては、画像ピラミッドを用いた階層構造探索法、ＳＳＤやＳＡＤに上限値を設けて以後の積算を打ち切る手法などがある。また、画像の特徴を利用する高速探索方法としては、注目領域の濃度ヒストグラムを利用する方法などがある。
【０１９７】
なお、本発明に係る画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法を実装したコンピュータ・プログラムを、ＣＰＵを備えた情報処理端末（例えば、パソコン、ワークステーション等のコンピュータ）に実行させることによって、合成画像或いは実画像を用いて、多パラメータを同時に高精度に推定することができる。
【産業上の利用可能性】
【０１９８】
以下のように、本発明を色々な画像処理分野に適用することができる。例えば、リモートセンシングや航空写真を使った地図作製では、複数の画像をつなぎ合わせること、つまりモザイキングを行うことが多い。このとき、画像間の対応位置を求めるために、領域ベースのマッチングを行う。このときには、対応が正しいことはもちろん、対応が画素単位より細かな分解能で求まることが重要になう。さらに、多数の画像を利用して対応を求めるので、計算の高速性も重要である。領域ベースマッチングを使って画素分解能以上の対応を求めるときには、画素単位のとびとびの位置で求まる類似度評価値を使って、フィッティング関数によるサブピクセル推定を行っているが、このときに本発明を適用することで、計算時間の増加がわずかで、はるかに高精度な対応位置を求めることができる。
【０１９９】
特に、画像パターンに斜め成分が含まれるとき、つまり、地図作製などで道路交差点などのパターンが画像座標系に対して傾いて撮影されているときに、本発明の効果が一層発揮される。
【０２００】
また、ステレオビジョンでは、複数枚の画像間の対応位置を求めることが基本的な処理であり、このときの対応位置精度が最終的な結果精度に大きく影響を及ぼす。画像間の対応位置を求めるときには、エピポーラ拘束などの拘束条件を利用するが、エピポーラ拘束条件自体を求めるときにも画像間の対応を調べる必要があり、このときにはエピポーラ拘束を利用できない。従って、本発明を適用することによって、高精度にエピポーラ拘束条件を求めることが可能になるため、結果として高精度なステレオ処理を行うことができるようになる。
【０２０１】
さらに、ＭＰＥＧ圧縮を行うときには、隣接フレーム間の対応位置を正確に求めることが、高品質で高圧縮率を達成するために必要である。隣接フレーム間の対応位置は、通常は１／２画素単位で求めるが、このときにはＳＡＤやＳＳＤ類似度評価値を利用した１次元の簡易手法を利用しているため、大きな誤差が含まれている可能性が在る。ところが、本発明は、従来の１次元推定方法に対してわずかな計算量の増加で、より確実で高精度な画像間の対応を求めることができるため、ＭＰＥＧ画像生成の圧縮率を向上することが可能になる。
【図面の簡単な説明】
【０２０２】
【図１】図１は、従来の１次元サブピクセル推定方法を説明するための図である。
【図２】図２は、３次元再構成アプリケーション用に一般的に使用される画像例（ａ）とその自己類似度（ｂ）を示す図である。自己類似度マップは、自己類似度マップがｋ≠１及びθｇ≠０の性質を有しているので、従来の１次元推定方法を用いることでサブピクセル推定誤差が生じることを示している。
【図３】図３は、θｃ＝π／２を有するテクスチャ画像例とその自己類似度を示す図である。
【図４】図４は、２次元画像モデルとその２次元類似度を示す図である。（ａ）はσimage＝０.７の画像モデルを示し、画像のサイズは１１×１１[画素]である。（ｂ）は画像モデルの自己類似度を示し、類似度の範囲は１１×１１である。（ｓ,ｔ）の暗い位置は高い類似度を有する２枚の画像の変位に対応する。最も暗い位置の（ｓ,ｔ）は、自己類似度の性質で（０,０）になる。
【図５】図５は、画像全体の回転角度θｇとコーナー角度θｃを有する２次元コーナー画像モデルを示す図である。
【図６】図６は、２次元コーナー画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、θｃ＝π／２、σimage＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、θｃ＝π／４、σimage＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。
【図７】図７は、画像全体の回転角度θｇと水平方向空間周波数ｆｕと垂直方向空間周波数ｆｖを有する２次元繰り返しテクスチャ画像モデルを示す図である。
【図８】図８は、２次元繰り返しテクスチャ画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、ｆｕ＝０.１、ｆｖ＝０.１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、ｆｕ＝０.０５、ｆｖ＝０.１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。自己類似度は一定の変位を有する２枚同様な画像から得られる。画像のサイズは１１×１１[画素]で、類似度の範囲は−５から５までである。
【図９】図９は、画像全体の回転角度θｇと水平標準偏差σと垂直標準偏差ｋσを有する２次元ガウス画像モデルを示す図である。
【図１０】図１０は、２次元ガウス画像モデルとその２次元自己類似度を示す図である。(ａ)、(ｂ)、(ｃ)は、σ＝２、ｋ＝１の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、σ＝２、ｋ＝０.５の画像モデルで、それぞれθｇ＝０，π／８，π／４である。(ｊ)、(ｋ)、(ｌ)は、それぞれ(ｇ)、(ｈ)、(ｉ)に対応する類似度である。自己類似度は一定の変位を有する２枚同様な画像から得られる。画像のサイズは１１×１１[画素]で、類似度の範囲は−５から５までである。
【図１１】図１１は、３種類異なった２次元画像モデルと、それらの２次元自己類似度と、２次元類似度モデルに比較して２次元自己類似度の部分図を示す図である。(ａ)はθｃ＝π／４、σimage＝１、θｇ＝０を有する２次元コーナー画像モデルを示す。(ｂ)はｆｕ＝０.０５１、ｆｖ＝０.１、θｇ＝０を有する２次元繰り返しテクスチャ画像モデルを示す。(ｃ)はσ＝２、ｋ＝０.５、θｇ＝０を有する２次元ガウス画像モデルを示す。(ｄ)、(ｅ)、(ｆ)は、それぞれ(ａ)、(ｂ)、(ｃ)に対応する類似度である。(ｇ)、(ｈ)、(ｉ)は、整数位置(〇で表す)でサンプリングされた(ｄ)、(ｅ)、(ｆ)のｔ＝０の水平部分図と、２次元類似度モデル(点線で表す)を示す。
【図１２】図１２は、２次元連続的に類似度計算と離散的に類似度計算の例を示す図である。
【図１３】図１３は、(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、ｋ＝０.７、θｇ＝π／６を有する２次元類似度モデルを示す図である。等高線は関数値の１０％から９０％までに対応している。(ａ)において、２次元類似度モデルの長軸がθｇ＝π／６を示す。(ｂ)において、ＨＥＬ（水平極値線）とＶＥＬ（垂直極値線）の両方が２次元類似度モデル(ｄｓ,ｄｔ)のピーク値位置を通るが、２次元類似度モデルの長軸或いは短軸と異なる。
【図１４】図１４は、ＨＥＬとＶＥＬの推定プロセスを説明するための図である。(ａ)において、●は３本の水平線ｔ＝−１、ｔ＝０、ｔ＝１上の類似度値（□）から得られた推定サブピクセル位置を示し、ＨＥＬは最小二乗で３つ●の位置から推定されることが可能である。(ｂ)において、ＶＥＬは同様なプロセスで推定されることが可能である。
【図１５】図１５は、従来の１次元推定方法及び本発明の実施例１に係る２次元推定方法にそれぞれ必要な類似度値の位置を示す図である。２次元類似度モデルは(ｄｓ,ｄｔ)＝(０.２,０.４)、σ＝１、ｋ＝０.３、θｇ＝π／８を有する。
【図１６】図１６は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝０、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１７】図１７は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝π／８、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１８】図１８は、推定誤差を比較するための図である。(ａ)は比較用に使用される画像モデルで、即ち、θｇ＝π／４、θｃ＝π／４、σimage＝１.０の２次元コーナーモデルである。(ｂ)は(ａ)の画像モデルの自己類似度を示す。
【図１９】図１９は、重み関数とパラメータに依存する補間特性を説明するための図である。
【図２０】図２０は、合成画像を使った実験結果を示す図である。
【図２１】図２１は、実画像を使ったターゲットトラッキング実験結果を示す図である。
【図２２】図２２は、本発明の実施例２に係る３パラメータ同時推定方法において、パラメータｓ１に関する最大値を通る平面Π１を説明するための模式図である。
【図２３】図２３は、実験に用いた合成画像である。
【図２４】図２４は、画像の回転計測結果を示す図である。
【図２５】図２５は、画像の位置計測結果を示す図である。
【特許請求の範囲】
【請求項１】
Ｎ（Ｎは４以上の整数である）個の画像間対応パラメータを軸とするＮ次元類似度空間において、離散的な位置で得られた画像間のＮ次元類似度値を利用して、前記Ｎ個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法であって、
前記画像間のＮ次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、
前記Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、
前記Ｎ個のＮ次元超平面の交点を求めるステップと、
前記交点を前記Ｎ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法。
【請求項２】
３個の画像間対応パラメータを軸とする３次元類似度空間において、離散的な位置で得られた画像間の３次元類似度値を利用して、前記３個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法であって、
前記画像間の３次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する平面を求めるステップと、
前記平面を各パラメータ軸に対して３個求めるステップと、
前記３個の平面の交点を求めるステップと、
前記交点を前記３次元類似度空間における３次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法。
【請求項３】
離散的に得られた画像間の２次元類似度値を利用して、連続領域における２次元類似度最大値または最小値の位置を高精度に推定する画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法であって、
前記画像間の２次元類似度値が、水平軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（水平極値線ＨＥＬ）を求めるステップと、
前記画像間の２次元類似度値が、垂直軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（垂直極値線ＶＥＬ）を求めるステップと、
前記水平極値線ＨＥＬと前記垂直極値線ＶＥＬの交点を求めるステップと、
前記交点を前記２次元類似度の最大値または最小値を与えるサブピクセル推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法。
【請求項１】
Ｎ（Ｎは４以上の整数である）個の画像間対応パラメータを軸とするＮ次元類似度空間において、離散的な位置で得られた画像間のＮ次元類似度値を利用して、前記Ｎ個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法であって、
前記画像間のＮ次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似するＮ次元超平面を求めるステップと、
前記Ｎ次元超平面を各パラメータ軸に対してＮ個求めるステップと、
前記Ｎ個のＮ次元超平面の交点を求めるステップと、
前記交点を前記Ｎ次元類似度空間におけるＮ次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける多パラメータ高精度同時推定方法。
【請求項２】
３個の画像間対応パラメータを軸とする３次元類似度空間において、離散的な位置で得られた画像間の３次元類似度値を利用して、前記３個の画像間対応パラメータを高精度に同時推定する画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法であって、
前記画像間の３次元類似度値が、ある１つのパラメータ軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する平面を求めるステップと、
前記平面を各パラメータ軸に対して３個求めるステップと、
前記３個の平面の交点を求めるステップと、
前記交点を前記３次元類似度空間における３次元類似度の最大値または最小値を与える前記画像間対応パラメータのサブサンプリンググリッド推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける３パラメータ高精度同時推定方法。
【請求項３】
離散的に得られた画像間の２次元類似度値を利用して、連続領域における２次元類似度最大値または最小値の位置を高精度に推定する画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法であって、
前記画像間の２次元類似度値が、水平軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（水平極値線ＨＥＬ）を求めるステップと、
前記画像間の２次元類似度値が、垂直軸に対して平行な直線上で最大または最小となるサブサンプリング位置を求め、求められた前記サブサンプリング位置を最も近似する直線（垂直極値線ＶＥＬ）を求めるステップと、
前記水平極値線ＨＥＬと前記垂直極値線ＶＥＬの交点を求めるステップと、
前記交点を前記２次元類似度の最大値または最小値を与えるサブピクセル推定位置とするステップと、
を有することを特徴とする画像のサブピクセルマッチングにおける２パラメータ高精度同時推定方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【公開番号】特開２００８−１１７４１６（Ｐ２００８−１１７４１６Ａ）
【公開日】平成２０年５月２２日（２００８．５．２２）
【国際特許分類】
【出願番号】特願２００７−３２３８６６（Ｐ２００７−３２３８６６）
【出願日】平成１９年１２月１４日（２００７．１２．１４）
【分割の表示】特願２００４−５２５６４０（Ｐ２００４−５２５６４０）の分割
【原出願日】平成１５年１０月２９日（２００３．１０．２９）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
【公開日】平成２０年５月２２日（２００８．５．２２）
【国際特許分類】
【出願日】平成１９年１２月１４日（２００７．１２．１４）
【分割の表示】特願２００４−５２５６４０（Ｐ２００４−５２５６４０）の分割
【原出願日】平成１５年１０月２９日（２００３．１０．２９）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
[ Back to top ]