(57)【要約】【課題】より高精度な画像マッチングを行うこと。【解決手段】画像マッチング装置において、参照画像を分割した領域ごとにラベルを割当てる動きクラスタリング処理部１０９と、動きモデルの第１の動きパラメータを推定する動きパラメータ推定部１１０と、第１の動きパラメータによって定まる対象画像上の点と参照画像上の点の間のポテンシャルエネルギーが最小になるように第１の動きパラメータを修正する動きパラメータ修正部１１３と、第１の動きパラメータから第２の動きパラメータを決定してたラベルに対応する領域内の各第２の格子点に割当てる動きパラメータ割り当て部１１１と、第２の動きパラメータと第１の運動方程式の解法処理結果とから補正項を算出する補正項計算部１１２と、第１の運動方程式に補正項を加えた第２の運動方程式の解法処理を行う第２の解法処理部１０８とを備えた。【選択図】        図１
2007257518.tif 000002
【特許請求の範囲】
【請求項１】
対象画像と参照画像との間の対応関係を求める画像マッチング装置であって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算部と、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力を計算する弾性エネルギー力計算部と、
前記各第２の格子点に作用する摩擦力を計算する摩擦力計算部と、
前記画像エネルギー力と前記第１の弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する第１の運動方程式を数値解析により解法処理を行う第１の解法処理部と、
前記第１の解法処理部による解法処理の結果から、前記参照画像を複数の領域に分割して、分割した領域ごとに領域を識別するラベルを割り当てる動きクラスタリング処理部と、
前記ラベルが割り当てられた領域における前記第２の格子点と前記第１の格子点の写像関係を示す動きモデルを規定するパラメータである第１の動きパラメータを数値解析により推定する動きパラメータ推定部と、
推定された前記第１の動きパラメータによって定まる前記対象画像上の点および前記参照画像上の点の間の前記ポテンシャルエネルギーが最小になるように前記第１の動きパラメータを修正する動きパラメータ修正部と、
修正された前記第１の動きパラメータから、前記ラベルに対応する領域内の各第２の格子点に対して最適な前記第１の動きパラメータである第２の動きパラメータを決定し、決定された前記第２の動きパラメータを、前記ラベルに対応する領域内の各第２の格子点に対して割り当てる動きパラメータ割り当て部と、
割り当てられた前記第２の動きパラメータと前記第１の解法処理部による解法処理の結果とに基づいて、前記第１の運動方程式を補正する補正項を算出する補正項計算部と、
前記第１の運動方程式に前記補正項を加えた前記各第２の格子点に関する第２の運動方程式を数値解析により解法処理を行う第２の解法処理部と
前記第２の解法処理部による解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理部と、
を備えたことを特徴とする画像マッチング装置。
【請求項２】
前記動きパラメータ推定部は、前記動きパラメータを、前記ラベルに基づいて最小二乗法による数値解析処理により推定することを特徴とする請求項１に記載の画像マッチング装置。
【請求項３】
前記動きパラメータ推定部は、前記動きパラメータを、前記ラベルに基づいてＭ推定による数値解析により推定することを特徴とする請求項１に記載の画像マッチング装置。
【請求項４】
前記動きパラメータ修正部は、最急降下法を用いて前記ポテンシャルエネルギーが最小となる点を求めることにより、前記動きパラメータを修正することを特徴とする請求項１に記載の画像マッチング装置。
【請求項５】
前記弾性エネルギー力計算部は、さらに、前記第２の格子点と前記第２の動きパラメータによって定められる前記第２の格子点の動き後の第３の格子点との間の弾性エネルギーから受ける第２の弾性エネルギー力を計算し、
前記第２の解法処理部は、前記画像エネルギー力と前記第２の弾性エネルギー力と前記摩擦力と前記補正項とによる前記第２の運動方程式を数値解析により解法処理を行うことを特徴とする請求項１に記載の画像マッチング装置。
【請求項６】
前記弾性エネルギー力計算部は、前記ラベルが同一の領域における前記第２の格子点と前記第３の格子点に関しては前記第２の弾性エネルギー力を計算し、前記ラベルが異なる領域における前記第２の格子点と前記第３の格子点に関しては前記第２の弾性エネルギー力が前記ラベルが同一の領域の場合より小さい値になるように計算することを特徴とする請求項５に記載の画像マッチング装置。
【請求項７】
前記補正項計算部は、前記第２の動きパラメータに基づいた動きベクトルと、前記第２の格子点の動きベクトルの差に基づいて前記補正項を算出することを特徴とする請求項１に記載の画像マッチング装置。
【請求項８】
前記補正項計算部は、前記第２の動きパラメータに基づいた動きベクトルと、前記第２の格子点の動きベクトルの差を入力とした線形関数の出力を前記補正項として算出することを特徴とする請求項７に記載の画像マッチング装置。
【請求項９】
前記補正項計算部は、前記第２の動きパラメータに基づいた動きベクトルと、前記第２の格子点の動きベクトルの差が所定値以上となる場合には、前記差を入力とした非線形関数の出力を前記補正項として算出することを特徴とする請求項７に記載の画像マッチング装置。
【請求項１０】
前記動きパラメータ割り当て部は、前記第２の格子点と前記第１の格子点の画素値の差分値に基づいたコスト関数により前記ラベルを規定し、前記コスト関数が最小値となる前記ラベルを前記第２の動きパラメータとして前記第２の格子点に割り当てることを特徴とする請求項１に記載の画像マッチング装置。
【請求項１１】
前記動きパラメータ割り当て部は、前記ラベルの連続性の拘束を加えた前記コスト関数により前記ラベルを規定し、前記コスト関数が最小値となる前記ラベルを前記第２の動きパラメータとして前記第２の格子点に割り当てることを特徴とする請求項１０に記載の画像マッチング装置。
【請求項１２】
前記動きパラメータ割り当て部は、最適であると決定された前記ラベルの前記コスト関数の値と当該ラベルの次に第２の前記ラベルの前記コスト関数の値とが第１の範囲内にある場合には、前記第２の動きパラメータを前記ラベルに対応する領域内の各第２の格子点に対して割り当てないことを特徴とする請求項１０に記載の画像マッチング装置。
【請求項１３】
前記動きパラメータ割り当て部は、さらに、最適であると決定された前記ラベルの前記コスト関数の値と当該ラベルの次に最適な前記ラベルの前記コスト関数の値とが第２の範囲内にある場合には、前記第２の動きパラメータを前記ラベルに対応する領域内の各第２の格子点に対して割り当てることを特徴とする請求項１２に記載の画像マッチング装置。
【請求項１４】
前記対象画像に存在するが前記参照画像に存在しない領域であるオクルージョン領域を検出するオクルージョン検出部と、
前記動きパラメータ割り当て部は、前記第２の格子点が前記オクルージョン領域に含まれる場合に、前記オクルージョン領域に含まない前記第２の格子点に対してのみ前記第２の動きパラメータを割り当てることを特徴とする請求項１に記載の画像マッチング装置。
【請求項１５】
前記オクルージョン検出部は、前記第２の格子点の周囲の格子点から形成される面積と前記参照画像における前記第２の格子点に対応する前記第１の格子点の周囲の格子点から形成される面積との比に基づいて前記オクルージョン領域を検出することを特徴とする請求項１４に記載の画像マッチング装置。
【請求項１６】
対象画像と参照画像との間の対応関係を求める画像マッチング方法であって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算ステップと、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力を計算する弾性エネルギー力計算ステップと、
前記各第２の格子点に作用する摩擦力を計算する摩擦力計算ステップと、
前記画像エネルギー力と前記第１の弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する第１の運動方程式を数値解析により解法処理を行う第１の解法処理ステップと、
前記第１の解法処理ステップによる解法処理の結果から、前記参照画像を複数の領域に分割して、分割した領域ごとに領域を識別するラベルを割り当てる動きクラスタリング処理ステップと、
前記ラベルが割り当てられた領域における前記第２の格子点と前記第１の格子点の写像関係を示す動きモデルを規定するパラメータである第１の動きパラメータを数値解析により推定する動きパラメータ推定ステップと、
推定された前記第１の動きパラメータによって定まる前記対象画像上の点および前記参照画像上の点の間の前記ポテンシャルエネルギーが最小になるように前記第１の動きパラメータを修正する動きパラメータ修正ステップと、
前記修正された前記第１の動きパラメータから、前記ラベルに対応する領域内の各第２の格子点に対して最適な前記第１の動きパラメータである第２の動きパラメータを決定し、決定された前記第２の動きパラメータを、前記ラベルに対応する領域内の各第２の格子点に対して割り当てる動きパラメータ割り当てステップと、
割り当てられた前記第２の動きパラメータと前記第１の解法処理ステップによる解法処理の結果とに基づいて、前記第１の運動方程式を補正する補正項を算出する補正項計算ステップと、
前記第１の運動方程式に前記補正項を加えた前記各第２の格子点に関する第２の運動方程式を数値解析により解法処理を行う第２の解法処理ステップと
前記第２の解法処理ステップによる解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理ステップと、
を含むことを特徴とする画像マッチング方法。
【請求項１７】
対象画像と参照画像との間の対応関係を求める画像マッチングプログラムであって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算ステップと、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力を計算する弾性エネルギー力計算ステップと、
前記各第２の格子点に作用する摩擦力を計算する摩擦力計算ステップと、
前記画像エネルギー力と前記第１の弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する第１の運動方程式を数値解析により解法処理を行う第１の解法処理ステップと、
前記第１の解法処理ステップによる解法処理の結果から、前記参照画像を複数の領域に分割して、分割した領域ごとに領域を識別するラベルを割り当てる動きクラスタリング処理ステップと、
前記ラベルが割り当てられた領域における前記第２の格子点と前記第１の格子点の写像関係を示す動きモデルを規定するパラメータである第１の動きパラメータを数値解析により推定する動きパラメータ推定ステップと、
推定された前記第１の動きパラメータによって定まる前記対象画像上の点および前記参照画像上の点の間の前記ポテンシャルエネルギーが最小になるように前記第１の動きパラメータを修正する動きパラメータ修正ステップと、
前記修正された前記第１の動きパラメータから、前記ラベルに対応する領域内の各第２の格子点に対して最適な前記第１の動きパラメータである前記第２の動きパラメータを決定し、決定された前記第２の動きパラメータを、前記ラベルに対応する領域内の各第２の格子点に対して割り当てる動きパラメータ割り当てステップと、
割り当てられた前記第２の動きパラメータと前記第１の解法処理ステップによる解法処理の結果とに基づいて、前記第１の運動方程式を補正する補正項を算出する補正項計算ステップと、
前記第１の運動方程式に前記補正項を加えた前記各第２の格子点に関する第２の運動方程式を数値解析により解法処理を行う第２の解法処理ステップと
前記第２の解法処理ステップによる解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理ステップと、
をコンピュータに実行させる画像マッチングプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象画像と参照画像の２つの画像からそれぞれの対応点を検出して対応付ける画像マッチング装置、画像マッチング方法および画像マッチングプログラムに関する。
【背景技術】
【０００２】
一つの対象画像から他方の参照画像の各画素それぞれの対応点を検出して対応関係を求める画像マッチング技術は、例えば、動画像における動き検出、ステレオマッチング、画像モーフィング、画像認識、動画像符号化等の画像処理技術の分野において利用されている。
【０００３】
このような画像マッチング技術としては、オプティカルフロー手法、ブロックベース手法、勾配法、ベイジアンメソッドという主として４種類に分類することができる（例えば、非特許文献１参照）。
【０００４】
オプティカルフロー手法は、「輝度の変化は一定である」というオプティカルフロー式を導出しそのオプティカルフロー式を拘束条件としてフローを求めるものである。ブロックベース手法は、画像を所定のブロックに分割し、ブロック毎のテンプレートマッチングによって動きを求める手法である。勾配法は、画像の輝度勾配が減少する方向にマッチングをおこなう手法である。ベイジアンメソッドは、確率的に尤もらしいマッチングを求める手法である。
【０００５】
また、従来の画像マッチング技術として、複数の多重解像度フィルタを用い、複数の多重解像度画像ピラミッドを生成し、生成された画像ピラミッドを上位の階層から下位の階層に順にマッチング処理を行うことによって、大きな動きから小さな動きまで画像の対応付けを行えるロバスト性の高い画像マッチング技術がある（例えば、特許文献１参照）。
【０００６】
【非特許文献１】A.MuratTekalp,“DigitalVideoProcessing”,PrenticeHall,1995
【特許文献１】特許第２９２７３５０号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、このような従来の画像マッチング技術には次のような問題がある。オプティカルフロー手法では、ノイズに敏感で速い動きに対応することが本質的に困難であるという問題があり、ブロックベースの手法では、画像の中のブロック毎に画像マッチングを行っているので、画像中のオブジェクトの平行移動などの動きには高い信頼性を有するが、画像中のオブジェクトが変形したり、回転する等の動きに対して本質的に対応することが難しいという問題がある。勾配法では、画像の輝度勾配が減少する方向にマッチングを行うため、安定してオブジェクトの動きを探索することが困難であるという問題がある。ベイジアンメソッドでは、大域的最適点の求め方が困難であるという問題がある。
【０００８】
一方、特許文献１に開示された技術では、原理的に複数の多重解像度フィルタを用い、多重解像度画像ピラミッドの最上位階層から最下位階層までのマッチング処理を行っているが、基本的には局所的な最適化処理であり大域的に最適になるとは限らない。
【０００９】
本発明は、上記に鑑みてなされたものであって、より高精度な画像マッチングを行うことができる画像マッチング装置、画像マッチング方法および画像マッチングプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
上述した課題を解決し、目的を達成するために、本発明は、対象画像と参照画像との間の対応関係を求める画像マッチング装置であって、前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算部と、前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力を計算する弾性エネルギー力計算部と、前記各第２の格子点に作用する摩擦力を計算する摩擦力計算部と、前記画像エネルギー力と前記第１の弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する第１の運動方程式を数値解析により解法処理を行う第１の解法処理部と、前記第１の解法処理部による解法処理の結果から、前記参照画像を複数の領域に分割して、分割した領域ごとに領域を識別するラベルを割り当てる動きクラスタリング処理部と、前記ラベルが割り当てられた領域における前記第２の格子点と前記第１の格子点の写像関係を示す動きモデルを規定するパラメータである第１の動きパラメータを数値解析により推定する動きパラメータ推定部と、推定された前記第１の動きパラメータによって定まる前記対象画像上の点および前記参照画像上の点の間の前記ポテンシャルエネルギーが最小になるように前記第１の動きパラメータを修正する動きパラメータ修正部と、修正された前記第１の動きパラメータから、前記ラベルに対応する領域内の各第２の格子点に対して最適な前記第１の動きパラメータである第２の動きパラメータを決定し、決定された前記第２の動きパラメータを、前記ラベルに対応する領域内の各第２の格子点に対して割り当てる動きパラメータ割り当て部と、割り当てられた前記第２の動きパラメータと前記第１の解法処理部による解法処理の結果とに基づいて、前記第１の運動方程式を補正する補正項を算出する補正項計算部と、前記第１の運動方程式に前記補正項を加えた前記各第２の格子点に関する第２の運動方程式を数値解析により解法処理を行う第２の解法処理部と、前記第２の解法処理部による解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理部と、を備えたことを特徴とする。
【００１１】
また、本発明は、上記画像マッチング装置で実行することができる画像マッチング方法および画像マッチングプログラムである。
【発明の効果】
【００１２】
本発明によれば、参照画像を複数の領域に分割して領域ごとに、第２の格子点の動きを規定した動きモデルを分類するための第１の動きパラメータの最適な動きパラメータである第２の動きパラメータを求めて、動きモデルを示す運動方程式の数値解析による解法処理を行って、対象画像と前記参照画像との対応関係を求めているので、領域境界の不連続性と領域内の一様性をロバストに表現できるようになり、より高精度な画像マッチングを行うことができ、その結果、より高画質な補間フレームを生成することができるという効果を奏する。
【発明を実施するための最良の形態】
【００１３】
以下に添付図面を参照して、この発明にかかる画像マッチング装置、画像マッチング方法および画像マッチングプログラムの最良な実施の形態を詳細に説明する。
【００１４】
なお、以下の実施の形態にかかる画像マッチング装置は、動画像の動き補償を行って補償画像を生成する動き補償装置に適用する他、ステレオマッチング、画像モーフィング、
画像認識、動画像符号化等の画像処理技術に適用することもできる。
【００１５】
（実施の形態１）
図１は、実施の形態１にかかる画像マッチング装置の構成を示すブロック図である。本実施の形態にかかる画像マッチング装置１００は、図１に示すように、第１解法処理部１０１と、画像エネルギー力計算部１０３と、弾性エネルギー力計算部１０４と、摩擦力計算部１０５と、第２解法処理部１０８と、動きクラスタリング処理部１０９と、動きパラメータ推定部１１０と、動きパラメータ修正部１１３と、動きパラメー割り当て部１１１と、補正項計算部１１２と、マッピング処理部１０７と、フレームメモリ１０６とを主に備えている。
【００１６】
実施の形態１にかかる画像マッチング装置１００では、対象画像および参照画像のそれぞれを所定の格子状に分割した各格子点を力学上の質点として扱い、対象画像の第１の格子点と第１の格子点に１対１に対応する参照画像の第２の格子点が画像の相関関係に基づいたポテンシャルエネルギーを考え、第１の格子点と第２の格子点の位置に基づいたポテンシャルエネルギーの勾配によって第２の格子点が受ける画像エネルギー力と、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力と、第２の格子点に生じる摩擦力とによる第２の格子点に関する第１の運動方程式を画像マッチング状態遷移モデルとして、この運動方程式に対して解法処理を行い、第２の格子点の平衡状態を求め、局所的な動きを求める。
【００１７】
更に本実施の形態にかかる画像マッチング装置１００では、複数個のパラメータ（動きパラメータ）によって写像関係を規定するパラメトリック動きモデルを導入している。すなわち、参照画像の画面内のすべての格子点の動き（第１の格子点との写像関係）が有限個のパラメトリック動きモデルによって決定されると仮定し、参照画像の画面を複数個の領域にクラスタリング（分割）し、動きモデルを規定するパラメータである動きパラメータを推定して、領域内の第２の格子点に最適の動きパラメータを割り当てる。そして、第１の運動方程式の解法処理結果である局所的な動きと、上記動きモデルによる動きとのずれを補正する力を補正項として第１の運動方程式に導入した第２の運動方程式を解法処理することによって、局所的に動きを更新する。このような、動きモデルの動きパラメータの推定処理、動きパラメータの割り当て処理、局所的な動きの更新（第２の運動方程式の解法処理）を繰り返すことによって、画像の正確な動きを推定して、第１の格子点と第２の格子点との写像関係を求めて画像マッチングを行っている。
【００１８】
第１解法処理部１０１は、対象画像上の複数の第１の格子点の各々と参照画像上で第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーの勾配により第２の格子点が受ける画像エネルギー力と、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力と、第２の格子点に生じる摩擦力とによる第２の格子点に関する運動方程式を、離散変数法としてのオイラー法による数値解析処理によって解法する処理部である。
【００１９】
画像エネルギー力計算部１０３は、第１の格子点と第２の格子点の位置に基づいた上記ポテンシャルエネルギーの勾配によって第２の格子点が受ける画像エネルギー力を計算する処理部である。
【００２０】
弾性エネルギー力計算部１０４は、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける第１の弾性エネルギー力を計算する処理部である。摩擦力計算部１０５は、第２の格子点に生じる摩擦力を計算する処理部である。また、弾性エネルギー力計算部１０４は、後述するラベルが同一の領域における第２の格子点と最適な動きパラメータによって定められる第２の格子点の動き後の第３の格子点に関しては第２の弾性エネルギー力を計算し、ラベルが異なる領域における第２の格子点と第３の格子点に関しては、ラベルが同一の領域における第２の弾性エネルギー力よりも小さくなるように第２の弾性エネルギー力を計算する。本実施の形態では、ラベルが異なる領域における第２の格子点と第３の格子点に関して、摩擦係数を０として第２の弾性エネルギー力を０とする処理を行っている。
【００２１】
動きクラスタリング処理部１０９は、第１解法処理部１０１による解法処理によって求めた局所的な動きから、クラスタリング手法により参照画像を複数の領域にクラスタリング（分割）して、分割した領域ごとに領域を識別するラベルを割り当てる処理部である。
【００２２】
動きパラメータ推定部１１０は、ラベルが割り当てられた領域における第２の格子点と第１の格子点の写像関係を示す動きモデルを規定するパラメータである動きパラメータを数値解析により推定する。本実施の形態では、最小二乗法による数値解析を行って、動きパラメータ（第１の動きパラメータ）を推定している。
【００２３】
動きパラメータ修正部１１３は、最急降下法により、推定した動きパラメータによって定まる対象画像上の点および参照画像上の点の間の画像の相関関係に基づくポテンシャルエネルギーが最小になるように動きパラメータを修正するものである。
【００２４】
動きパラメータ割り当て部１１１は、推定された動きパラメータから、ラベルに対応する領域内の各第２の格子点に対して最適なラベルを決定して最適な動きパラメータ（第２の動きパラメータ）を求め、この最適な動きパラメータを、ラベルに対応する領域内の各第２の格子点に対して割り当てる処理部である。
【００２５】
補正項計算部１１２は、割り当てられた最適な動きパラメータによる動きベクトルと第１の解法処理部による解法処理の結果である動きベクトルとの差分から第１の運動方程式を補正する補正項を算出する処理部である。
【００２６】
第１解法処理部１０８は、画像エネルギー力と、第２の格子点と最適な動きパラメータによって定められる第２の格子点の動き後の第３の格子点との間の弾性エネルギーから受ける第２の弾性エネルギー力と、摩擦力と、補正項計算部１１２で求めた補正項とによって第２の格子点の動きを規定した第２の運動方程式を数値解析によって解法する処理部である。
【００２７】
マッピング処理部１０７は、第２の解法処理部１０８による解法処理の結果から第２の格子点の平衡状態を求めることにより対象画像上の第１の格子点と参照画像上の第２の格子点との対応関係である写像関係を求める処理部である。
【００２８】
フレームメモリ１０６は、入力された対象画像と参照画像と保存するメモリである。なお、各部の詳細な処理については後述する。
【００２９】
次に、本実施の形態にかかる画像マッチング装置１００で導入する画像マッチング状態遷移モデルについて説明する。
【００３０】
２次元ユークリッド空間上の点ｘ∈Ｅ2における対象画像、参照画像の画素値をそれぞれ（１）式で示す。
【００３１】
【数１】
2007257518.tif 000003
（１）式は実数空間における連続画像モデルであり、本実施の形態ではデジタル画像を扱うため、２次元格子空間上の点ｎ∈Ｎ2について（２）式に示す離散画像モデルを考える。
【００３２】
【数２】
2007257518.tif 000004
従って、連続画像モデルの画素値はデータ補間関数ｇを用いて、（３）式に示すように定義される。ここで、データ補間関数ｇとしては、例えば双曲形補間関数を使用することができる。
【００３３】
【数３】
2007257518.tif 000005
本実施の形態では、対象画像から参照画像への画像マッチング問題を、対象画像上の曲面Ｘ⊂Ｅ2から参照画像上の曲面Ｙ⊂Ｅ2への写像として求めることにする。図２は、対象画像上の曲面から参照画像上の曲面への写像を示す説明図である。図３は、対象画像上の曲面Ｘの点ｘから参照画像上の曲面Ｙの点ｙへの写像を示す説明図である。すなわち、図２および３に示すように、画像マッチングは、対象画像上の曲面Ｘの点が参照画像上の曲面Ｙのどの点に対応しているのかを表す写像を求める問題として扱う。
【００３４】
このような写像に関する条件を導出すると以下のようになる。パラメータ（ｕ，ｖ）を用いて、曲面Ｘ上の点をｘ（ｕ，ｖ）∈Ｘ、曲面Ｙ上の点をｙ（ｕ，ｖ）∈Ｙと表し、パラメータ（ｕ，ｖ）と点ｘが同じスケールを有する場合、曲面Ｘ上の点ｘ（ｕ，ｖ）は、パラメータ（ｕ，ｖ）によって、（４）式で表すことができる。
【００３５】
【数４】
2007257518.tif 000006
ここで、曲面Ｘから曲面Ｙへの求めるべき写像を（５）式で定義する。写像の対応関係は曲面Ｘ上の点ｘが一意に曲面Ｙに対応づけることができればよいので、写像ｇはすべての点が一意に決定する全単射とする。よって同一パラメータ（ｕ，ｖ）を有する点ｘ（ｕ，ｖ）、ｙ（ｕ，ｖ）は，（６）式を満たすことになる。
【００３６】
【数５】
2007257518.tif 000007
このことは、すなわち、写像ｇについて考えることは、曲面Ｙそのものを考えることと同じである。一つの対象物が対象画像と参照画像に射影されて映っていると考えると、射影された対象画像上の像と参照画像上の像は、画素値が同じであると考えることができる。従って、対応づけるべき同一点の画素値が等しいと仮定すると、写像ｇに関して（７）式の条件式が導出される。かかる（７）式をパラメータｕ＝（ｕ，ｖ）を用いて表示すると、（８）式のようになる。
【００３７】
【数６】
2007257518.tif 000008
従って、画像マッチング問題は、（８）式の条件式を満たす写像ｇを求めるという問題に帰着することになる。
【００３８】
次に、上記画像マッチング問題の解法処理について説明する。（８）式は、（９）式のように変形することができる。
【００３９】
【数７】
2007257518.tif 000009
しかしながら、（９）式は、画像に含まれるノイズなどにより厳密には成立しないため、かかる式から直接写像ｇを求めることは困難である。このため、本実施の形態では、最小二乗法を利用している。最小二乗法によれば、（１０）式に示すエネルギー関数を考え、かかるエネルギー関数の最小化を行うことによりノイズが含まれていても条件を満足する写像ｇを求めることができる。
【００４０】
【数８】
2007257518.tif 000010
大数の法則が成立する程度の標本数があれば、最小二乗法により解が安定に求まると考えられるが、本実施の形態では、標本数が一つであるため、やはりノイズの影響を受けてしまう。解を安定して求めるためには、ベイズの定理により、近傍では滑らかに変化するという事前知識により事後分布の精度を高めるアプローチ、例えば、写像ｇに対して構造を規定して力学的なアプローチを考える。
【００４１】
写像ｇは全単射であるので、写像ｇについて考えることは、曲面Ｙそのものを考えることと同義である。力学的構造としては、例えば、剛体のような構造が考えられるが、画像の場合には変形等が考えられるため、剛体の構造は適切でない。そこで、本実施の形態では、曲面Ｙに対して、画像の変形等に対して柔軟に対応可能な弾性体の構造を導入する。
【００４２】
弾性体のエネルギーは、（１１）式で示される。
【００４３】
【数９】
2007257518.tif 000011
また、（１０）式で示される画像エネルギーを曲面Ｙ上の点を用いて表すと、（１２）式で表される。
【００４４】
【数１０】
2007257518.tif 000012
ポテンシャルエネルギーは、（１２）式の画像エネルギーと（１１）式の弾性エネルギーの線形結合により（１３）式のように定義される。
【００４５】
【数１１】
2007257518.tif 000013
このポテンシャルエネルギーを最小にすることは、画像マッチングの条件式である（７）式をできるだけ満たすように、かつ弾性体としての安定形態が取れるような曲面Ｙを探索することを意味する。このことをエネルギー最小化問題として定式化すると（１４）式で表される。
【００４６】
【数１２】
2007257518.tif 000014
従って、条件式（７）式を満たす写像ｇを求める画像マッチング問題は、エネルギー最小化問題の（１４）式を解法することに帰着する。
【００４７】
また（１４）式より最適性の必要条件は、（１５）式で求めることができる。
【００４８】
【数１３】
2007257518.tif 000015
従来の手法では、この最適性必要条件の（１５）式を満たすような曲面Ｙ上の点ｙを求めることにより解法している。これらは静的なエネルギーの最小化であり、静的な釣り合いの位置を求める静的な最適化である。図４は、ポテンシャルエネルギーの最小化の点を探索する場合の静的な最適点の説明をするための模式図である。図４に示すように、静的な最適化では初期値に最も近い位置で質点のポテンシャルエネルギーの局所最適点（釣り合い位置）が求まる。従って、この局所最適点が大域的な最適点になるかどうかは初期値に依存することになる。言い換えれば、初期値の選定が好ましくない場合には、大域的に最適な点が得られない場合がある。
【００４９】
このため、大域的な最適点を求めて（１４）式を解法するため、本実施の形態では、対象画像および参照画像上の曲面を格子状に分割して分割した各格子点を質点と考え、参照画像上の曲面Ｙの格子点（質点）における運動エネルギーを導入し、大域的な最適化を行っている。すなわち、曲面Ｙに対して時間軸を導入し、曲面Ｙが時間に応じて変形できるような構造とし、曲面に対して運動の構造を適用する。図５は、ポテンシャルエネルギーの最小化の点を探索する場合の大域的な最適点の説明をするための模式図である。図５に示すように、運動エネルギーの導入により、最初のポテンシャルエネルギーが運動エネルギーに変換されるため、格子点（質点）が局所最適点からさらに広い範囲で（１５）式を満たす大域的な最適点を探索することができる。これにより、静的な最適化よりもより広い範囲を探索できるので、初期値への依存性が低くなり、ノイズ等の影響に対してよりロバストになる可能性がある。
【００５０】
次に、（１５）式の大域的最適化のための運動方程式の導入について説明する。まず、運動方程式に導入する時間τを用いて曲面Ｙ上の点ｙを（１６）式のように拡張する。（１６）式で表した曲面Ｙ上の点ｙの全微分は（１７）式となるが、パラメータ（ｕ，ｖ）Tは時間τに対して独立であると仮定すると、（１８）式が成立するので、結局（１７）式は（１９）式で表現される。また、（１９）式から曲面Ｙ上の点ｙの２階微分は（２０）式で表される。（１９）式で示される曲面Ｙ上の点ｙの全微分と（２０）式で示される２階微分をそれぞれ（２１）式で表現する。
【００５１】
【数１４】
2007257518.tif 000016
（１９）式は、曲面Ｙ上の点ｙの速度を表し、（２０）式は、曲面Ｙ上の点ｙの加速度を表している。従って、曲面Ｙ上の点ｙの運動エネルギーは、（２２）式で定義することができる。
【００５２】
【数１５】
2007257518.tif 000017
このような運動エネルギーを導入したことは、曲面Ｙ上の点ｙを質点と考えて、曲面上の運動を定義したことになり、曲面Ｙ上の点ｙの運動という力学的な構造を規定したことを意味する。かかる曲面Ｙ上の点ｙの運動は運動方程式によって記述される。
【００５３】
次に、ラグランジュの方法に従って次のように曲面Ｙ上の点ｙの運動方程式を導出する。曲面Ｙは、エネルギー最小化の（１４）式に従うので、ポテンシャルエネルギーＵを最小にする方向に動く必要がある。ラグランジアンは（２３）式で定義される。
【００５４】
（２３）式は、位置のポテンシャルエネルギーである重力場における物体の自由落下などと同様に、物体は重力場を最小化するように落下することを意味している。ラグランジュ方程式は（２４）式で示されるので、（２３）および（２４）式から、曲面Ｙ上の点ｙの運動方程式は、（２５）式で表される。
【００５５】
【数１６】
2007257518.tif 000018
【００５６】
【数１７】
2007257518.tif 000019
（２５）式の運動方程式に従った運動は保存系であり、エネルギー保存則により解は収束しない。このため、系を非保存系にするために（２６）式で示される摩擦エネルギーを導入する。
【００５７】
【数１８】
2007257518.tif 000020
（２６）式で示される摩擦エネルギーを導入した場合のラグランジュ方程式は、（２７）式で表される。
【００５８】
【数１９】
2007257518.tif 000021
この（２７）式を解法することによって、摩擦エネルギーを導入した運動方程式が（２８）式で得られる。
【００５９】
【数２０】
2007257518.tif 000022
次に、動的な探索、（２８）式の運動方程式に従った解がエネルギー最小化の（１４）式に対する最適性必要条件である（１５）式を満足することについて説明する。
【００６０】
まず、探索が完了した場合を考える。この場合、点ｙは完全に停止するので、点ｙの速度（ｙの全微分）および加速度（ｙの２階微分）はともに０となるので、（２８）式の運動方程式により（２９）式となり、これにより（３０）式が成立する。
【００６１】
【数２１】
2007257518.tif 000023
（３０）式が成立することは、最適性必要条件の（１５）式を満足することを意味している。従って、探索が完了した場合、そのときの点ｙは最適点であることがわかる。
【００６２】
次に、曲面Ｙ上の点ｙが静止した場合を考える。この場合、点ｙの速度（ｙの全微分）は０となるが加速度（ｙの２階微分）は０にならないので、（２８）式の運動方程式により（３１）式が得られる。
【００６３】
【数２２】
2007257518.tif 000024
ここで、点ｙが最適点であると仮定すると、最適性条件の（１５）式により、点ｙの加速度（ｙの２階微分）は０となり、最適点の探索が完了する。一方、点ｙが最適点でない場合には、点ｙの加速度（ｙの２階微分）が０とならないため、最適点の探索は完了していないことになる。すなわち、点ｙが静止した位置が最適点の位置である場合に（３１）式により最適性必要条件の（１５）式が満足することを意味している。
【００６４】
次に、点ｙに摩擦エネルギーが作用している場合を考える。摩擦エネルギーの性質から（３２）式が導かれ、これは点ｙの運動が停止して時間の経過により最適点の探索が完了することを意味している。
【００６５】
【数２３】
2007257518.tif 000025
このように、最適点の探索が完了した場合、点ｙが静止した場合、点ｙに摩擦エネルギーが作用している場合の上記説明により、（２８）式の運動方程式に従った点ｙの運動が、（１５）式の最適性必要条件を満足することがわかる。
【００６６】
一方、（１５）式は１階微分の形式であるが、２階微分の最適性必要条件については、図４または５に示すポテンシャルエネルギーの山の頂点で探索が完了する場合があるため、満足されない場合がある。この場合には、ポテンシャルエネルギーの山の頂点では、最適点の探索が完了しないための条件として（３３）式を仮定することにより、２階の最適性必要条件を満足することができる。この（３３）式の仮定は、具体的には（３４）式の条件を探索アルゴリズムで実現することにより達成することができる。
【００６７】
以上により、（２８）式の運動方程式に従った最適点の探索は（１４）式で示される最適化問題に対する最適解を与えることがわかる。このような探索は、静的な探索に比べてロバストなものとなる。
【００６８】
ここで、（２８）式に示す運動方程式において、弾性エネルギーＥkの中に点ｙおよびパラメータ（ｕ，ｖ）に関する偏微分が含まれているため、有限差分法によって偏微分方程式を差分方程式によって置換する。弾性エネルギーＥkの右辺を有限差分法によって離散化すると、（３５）式のようになる。
【００６９】
【数２４】
2007257518.tif 000026
ここで、ｎは離散化されたパラメータ（ｕ，ｖ）の格子点、ｙnはｎに対応する離散化されたＹ上の点、ｙn+(-1,0)，ｙn+(1,0)，ｙn+(0,-1)，ｙn+(0,1)は点ｙnの偏微分に関連する上下左右の格子点を示している。
【００７０】
この（３５）式をｙnで偏微分すると（３６）式が得られ、かかる式は、弾性エネルギーによって発生する力である弾性エネルギー力を示している。
【００７１】
【数２５】
2007257518.tif 000027
この弾性エネルギー力を（３７）式で示し、離散パラメータΔｕ，Δｖをそれぞれ１とすると、弾性エネルギー力は（３８）式で表されることになる。
【００７２】
【数２６】
2007257518.tif 000028
ここで、画像エネルギーＥiは、パラメータ（ｕ，ｖ）に関して独立と考え、（２８）式の運動方程式もそのまま離散可能となり、（３９）式で表される。
【００７３】
【数２７】
2007257518.tif 000029
この（３９）式の右辺第１項は、点ｙnにおける画像エネルギーの勾配の逆向きの力である画像エネルギー力であり、（４０）式のようにおく。
【００７４】
【数２８】
2007257518.tif 000030
次に、この画像エネルギー力Ｆi（ｎ）、すなわち画像エネルギーの勾配の計算について説明する。本実施の形態では、画像エネルギーの勾配を、局所最適化の手法によって求めている。図６および図７は、画像エネルギー力の概念を示す説明図である。ここで、画像モデルＳcは連続画像モデルであり、実際にはサンプリングされた画像モデルＳpを利用しているため、局所最適化においてもサンプリングされた画像モデルに基づいておこなっている。
【００７５】
ｙn（τ）にもっとも近いサンプリング点を局所空間中心ｙcとして、ｙcを（４１）式で求める。
【００７６】
【数２９】
2007257518.tif 000031
そして、隣接空間Ｌを（４２）式で定義する。このとき、局所探索集合Ωは、（４３）式のように定義することができる。
【００７７】
【数３０】
2007257518.tif 000032
そして、局所最適化の後、その方向へのベクトルを求め、求めたベクトルを正規化し勾配の大きさを乗算すると、画像エネルギー力として（４４）〜（４６）式が得られる。本実施の形態では、画像エネルギー力として（４４）〜（４６）式を使用している。
【００７８】
【数３１】
2007257518.tif 000033
なお、実装上の扱いやすさ等を考慮してエネルギー関数を（４７）式のように定義して、（４７）式により導かれる画像エネルギー力を用いることもできる。
【００７９】
【数３２】
2007257518.tif 000034
以上説明した画像エネルギーＦi（ｎ）、弾性エネルギー力Ｆk（ｎ）、摩擦力を使用して、格子点ｙn（τ）に対する運動方程式を画像マッチング状態遷移モデルとして生成すると、運動方程式は、（４８）式で表される。また、ｙn（０）が参照画像に等しいとし、運動エネルギーは０、すなわちｙn（０）の速度が０として考えると、運動方程式の初期条件は、（４９）式となる。
【００８０】
【数３３】
2007257518.tif 000035
本実施の形態では、画像エネルギー力計算部１０３によって、（４４）〜（４６）式で示される画像エネルギー力Ｆiを計算し、弾性エネルギー力計算部１０４によって（３８）式で示される弾性エネルギー力Ｆkを計算し、摩擦力計算部１０５によって（４８）式の右辺第３項で示される摩擦力を計算する。そして、第１解法処理部１０１によって（４８）、（４９）式で示される運動方程式を数値解析処理によって解法している。
【００８１】
次に、この第１解法処理部１０１による運動方程式の解法について説明する。（４８）、（４９）式で示される運動方程式、すなわち常微分方程式は、一般的に解析的に解法することができないため、この運動方程式の系が収束するのに十分大きな時間Ｔを考え、数値解析によってｔ＝（０，Ｔ）区間を計算することによって運動方程式の収束状態を推定する。
【００８２】
本実施の形態では、常微分方程式は初期値が定まれば離散変数法によって一意に解が求まることを利用して、オイラー法による数値解析処理を行っている。
【００８３】
尚、本実施の形態ではオイラー法による数値解析処理を行っているが、これに限定されるものではない。離散変数法には、オイラー法以外にも、ルンゲクッタ法、ブリルシュ・ストア法、予測子・修正子法、隠的ルンゲクッタ法など種々の手法があり、このため、数値解析処理として、これらのいずれの手法を用いてもよい。以下はオイラー法を用いた数値解析処理を例にあげて説明する。
【００８４】
オイラー法は、一階の常微分方程式に対する数値解法であるため、（５０）式の変数変換を（４８）、（４９）式を運動方程式に施すことにより、（４８）、（４９）式を一階の常微分方程式に変換する。これにより、（５１）、（５２）式が得られる。
【００８５】
【数３４】
2007257518.tif 000036
【数３５】
2007257518.tif 000037
一方、常微分方程式（５３）式をオイラー法で数値解析処理を行って解法する場合には、（５４）式のようにｔ(m)からｔ(m+1)=ｔ(m)+ｈへ解を進展させる。
【００８６】
【数３６】
2007257518.tif 000038
【数３７】
2007257518.tif 000039
ここで、ｘ(m)はｍステップであることを示し、ｈはステップ幅である。オイラー法のスキームを（５１）、（５２）式に適用すると（５５）、（５６）式の更新式が得られる。
【００８７】
【数３８】
2007257518.tif 000040
本実施の形態では、第１解法処理部１０１によって、（５１）、（５２）式の更新式を繰り返して解法することによって第１の運動方程式（４８）、（４９）式の系の収束状態を求めている。そして、マッピング処理部１０７では、この収束状態から、曲面Ｙをｙn＝ｙn(T)のように求め、さらに、（６）式よりｙn＝ｇ（ｎ）の関係から、写像ｇを求めている。
【００８８】
（５５）、（５６）式の更新式は、各格子点の動きの状態が局所的な関係によって決定されるものであり、局所的な状態が平衡状態になりながら動きを推定していくため、より大域的な構造を導入し精度の向上を図る。
【００８９】
本実施の形態にかかる画像マッチング装置１００では、複数個のパラメータ（動きパラメータ）によって写像関係を規定するパラメトリック動きモデルを導入している。以下、パラメトリック動きモデルを利用した動き推定について概要を説明する。
【００９０】
参照画像の画面内のすべての格子点の動き（第１の格子点との写像関係）が有限個のパラメトリック動きモデルによって決定されると仮定している。動きパラメータが既知であると仮定すると、画面内の各格子点の動きを求める問題は、各格子点がどのパラメトリック動きモデルに属するかを決定すること、すなわち、最適な動きパラメータを選択する問題となる。このことは、複数個の動きパラメータの最適化を行うことであり、画面全体で取り得る動きに関して最適化問題を解法することに比べてきわめて容易な手法となる。
【００９１】
しかしながら、動きパラメータは既知ではないため、局所的な動的マッチング処理の（５５）、（５６）式の更新式に従った動き推定の結果を用い、かかる動き推定の結果をクラスタリング手法を利用して、複数個のクラスタ（グループ）に分割する。本実施の形態では、動きクラスタリング処理部１０９によって、k-means法によって数値解析を行うことにより、（５５）、（５６）式の更新式に従った動き推定の結果から第２の格子点を複数個の領域にクラスタリングしている。なお、クラスタリングの手法としては、k-means法の手法以外の手法を用いてもよい。
【００９２】
そして、動きパラメータ推定部１１０によって、各クラスタ内（領域内）の局所的な動きから最小二乗法を用いて、パラメトリック動きモデルの動きパラメータ（第１の動きパラメータ）を推定している。なお、動きパラメータを推定する手法としては、最小二乗法以外の手法を用いることも可能である。
【００９３】
動きパラメータの推定を行った後、動きパラメータ修正部１１３によって、推定された動きパラメータを修正する。その後、動きパラメータ割り当て部１１１は、各格子点に、各格子点が属する動きモデルの動きパラメータを割り当てる。具体的には、動きパラメータ割り当て部１１１は、各格子点がどの動きパラメータの動きモデルに属するかを示すラベルを導入し、ラベルを割り当てることにより、最適な動きパラメータ（第２の動きパラメータ）の割り当てを行う。このため、動きパラメータ割り当て部１１１は、動きモデルに対する評価関数を設定し、その評価関数が最小になる動きパラメータが示すラベルを各格子点に設定する。
【００９４】
各格子点に対して、第１解法処理部１０１によって（５５）、（５６）式の更新式により、局所的な動きが求められるが、同時に各格子点には、パラメトリック動きモデルからも動きが求められる。このため、各格子点の局所的な動きと動きモデルによる動きのズレを補正する力を示す補正項を導入し、かかる補正項を補正項計算部１１２によって求め、第２解法処理部１０８によって、求めた補正項を第１の運動方程式に加えた第２の運動方程式を数値解析により解法することによって局所的な更新をおこなっている。図８は、動きモデルを利用した動き推定の画像上における状態を示す説明図である。
【００９５】
このように本実施の形態にかかる画像マッチング装置１００では、上述した動きパラメータ処理部１１０による動きモデルの動きパラメータの推定処理、動きパラメータ割り当て部１１１による動きパラメータの割り当て処理、補正項計算部１１２および第２解法処理部１０８による局所的な動きの更新（第２の運動方程式の解法処理）を繰り返すという動き分割を用いた処理によって、画像の正確な動きを推定して、第１の格子点と第２の格子点との写像関係を求めて画像マッチングを行っている。図１０は、実施の形態１にかかるパラメトリック動きモデルを導入して動き分割を用いた画像マッチングの概念的な流れを示す説明図である。
【００９６】
例えば、図９に示すように、画像上のオブジェクトの一部が誤った動きで画像マッチングされている状態を考える。図９に示すような誤った動きがあると、境界の歪み等として視認されてしまう。このような場合には、誤った動きの領域に対して正確な動きの情報を割り当てることによりオブジェクト全体に正しい画像マッチングを行うことができる。本実施の形態における最適な動きパラメータの割り当て、すなわち格子点が属するパラメトリック動きモデルを決定することは、このような誤った動きの領域に対して正確な動きの情報を割り当てることに相当し、これにより、画像上で文字の一部が崩れるといった不具合を防止することができる。
【００９７】
次に、このようなパラメトリック動きモデルを導入した動き推定の詳細について説明する。まず、パラメトリック動きモデルを規定する。各パラメトリック動きモデルは、任意の数のラベルα∈Ｌ⊂Ｚによってラベル付けされているものとする。格子点ｎのラベルをｚ（ｎ）∈Ｌとし、ラベルαを有する領域の集合（各動き領域）をＮα＝｛ｎ｜ｚ（ｎ）＝α｝とする。ラベルαの領域ごとのパラメトリック動きモデルを（５７）式で定義する。
【００９８】
【数３９】
2007257518.tif 000041
動きパラメータａα(m)は種々の形式を採用することができる。例えば、動きパラメータａα(m)として、アファインモデルを用いた場合には（５８）式、２次形式モデルを用いた場合には（５９）式で示される。ここで、ｍは、イテレーションの変数である。
【００９９】
【数４０】
2007257518.tif 000042
【数４１】
2007257518.tif 000043
この他、動きパラメータとして、回転変換モデル、ユークリッド変換モデル、相似変換モデル、射影変換モデル、Ｌｉｅ変換モデルなどを用いてもよい。
【０１００】
次に、動きパラメータ推定について説明する。動きパラメータ推定の処理は、動きパラメータ推定部１１０によって行われる。パラメトリック動きモデルによって定まる動きベクトルは、（６０）式で示されるので、これと（５７）式を用いると、パラメトリック動きモデルによって定まる動きベクトルは（６１）式で示される。
【０１０１】
【数４２】
2007257518.tif 000044
一方、対応点により求まる動きベクトルは、（６２）式で示される。
【０１０２】
【数４３】
2007257518.tif 000045
パラメトリック動きモデルによって定まる動きベクトルと対応点により求まる動きベクトルとの誤差を（６３）式で定義すると、誤差エネルギーは（６４）式で表される。従って、各動きベクトルの定義式（６１）、（６２）式により、誤差エネルギーは（６５）式のようになり、結局（６６）式で示されることになる。
【０１０３】
【数４４】
2007257518.tif 000046
【数４５】
2007257518.tif 000047
一方、ラベルα領域内の誤差は、（６７）式で示され、このラベルα領域内の誤差が（６８）式に示すように最小になるように動きパラメータが推定される。
【０１０４】
【数４６】
2007257518.tif 000048
【数４７】
2007257518.tif 000049
本実施の形態では、動きパラメータと得られた動きベクトルの誤差がガウス分布に従っていると仮定して、最適な動きパラメータは最小二乗法を使用した数値解析によって求めている。すなわち、最小二乗法における正規方程式は、（６９−１）、（６９−２）、（６９−３）式に示されるため、動きパラメータは（７０−１）式で示される。
【０１０５】
【数４８】
2007257518.tif 000050
【数４９】
2007257518.tif 000051
かかる（７０−１）式は連立一次方程式であり、このため、動きパラメータは、（７０）式により特異値分解やＬＵ分解などの手法で求めることができる。
【０１０６】
次に、パラメータ修正部１１３によるパラメータ修正処理について説明する。上記動きパラメータ推定部１１０による動きパラメータの推定は、局所動きベクトルから動きパラメータを推定するが、推定された動きパラメータによって規定される動きモデルは画像に適合するかどうか分からない。すなわち、動きパラメータの画像に対する尤度は考慮されていないことになる。
【０１０７】
このことは、次に示すような問題として顕在化する。局所動き推定には推定誤差が入り込み、動きパラメータ推定にも推定誤差が入り込む。このため、得られた動きモデルは必ずしも画像に適合しているとは限らないため、ラベル割り当ての精度が低下する場合があり得る。その結果、正しいラベル割り当てが行われず、その情報を用いて次の局所動き推定をおこなうため、さらに誤差が拡大してしまう。これはフィードバックによる誤差の発散効果と見ることができるので、いずれかの段階で誤差のフィードバックをキャンセルする必要がある。上記の動きパラメータの画像に対する尤度を最適化することにより、動きパラメータの画像に対する誤差を減少することができ、十分上記フィードバックのキャンセルを行うことができる。
【０１０８】
動きパラメータの画像に対する尤度は、（７０−２）、（７０−３）式で示される。
【数５０】
2007257518.tif 000052
【０１０９】
（７０−２）、（７０−３）式は、ラベリングの尤度と等価であるが、ラベリングにおいてはラベルによって最適化がおこなわれるため、動きパラメータに関しては最適化されない。動きパラメータの最適化問題としては、（７０−４）式を解法すればよい。
【数５１】
2007257518.tif 000053
【０１１０】
（７０−４）式において、ノイズの標準偏差項は最尤推定の結果には影響を与えないため省略している。この問題は、Ｌｕｃａｓ−Ｋａｎａｄｅアルゴリズムによって安定して解くことができることが知られている。ここでは、各ラベルαに対して個別に最尤推定をおこないたいため、ラベル毎にＬｕｃａｓ−Ｋａｎａｄｅアルゴリズムを適用する。
【０１１１】
ラベルを有する領域の集合、すなわち各動き領域を（７０−５）式で示すと、動きパラメータの最適化問題は各領域毎に定義され、（７０−６）式で示される。
【数５２】
2007257518.tif 000054
【０１１２】
Ｌｕｃａｓ−Ｋａｎａｄｅアルゴリズムは、真の解（７０−７）式の十分近い近傍点（７０−８）式は既知として、（７０−９）式のように解を反復的に漸化させることにより近似解を求める手法であり、（７０−１０）式による最小化を行えばよい。
【数５３】
2007257518.tif 000055
【０１１３】
【数５４】
2007257518.tif 000056
【０１１４】
（７０−１０）式を、ａαの周りでテーラー展開すると、（７０−１１）式のようになる。かかる式をΔａαで微分すると（７０−１２）式が得られる。
【数５５】
2007257518.tif 000057
【０１１５】
かかる式に対し、ニュートン法を適用すると解の更新は（７０−１３）式のようになる。ここで、Ｈは（７０−１４）式で示されるHessianである。
【数５６】
2007257518.tif 000058
【０１１６】
従って、解の更新式は（７０−１５）式で示される。
【数５７】
2007257518.tif 000059
【０１１７】
ここで、例えば十分小さい定数ε（＞０）に対して、（７０−１６）式を満たしたら終了等のように終了条件を定めておく。なお、所定のイテレーション回数を繰り返した後に終了するように構成してもよい。
【０１１８】
（７０−１０）式は、各画素を等価に扱っているために、ラベル内に外れ値が入っていると、その影響をうけて解（平均）がずれてしまう。そこで、（７０−１０）式に重み付けを行い、（７０−１７）式で示す問題を考える。
【数５８】
2007257518.tif 000060
【０１１９】
この重みｗ（ｎ）を外れ値の場合には低くなるように設定することより、外れ値の影響に対してロバストになる。解の更新（７０−１３）式の導出と同様に更新式を導出すると、（７０−１８）式のようになる。なお、重みｗ（ｎ）の決定はＭ推定の枠組みなどを使うことできる。
【数５９】
2007257518.tif 000061
【０１２０】
従って、本実施の形態では、動きパラメータ修正部１１３によって動きパラメータ推定部１１０で推定した動きパラメータを、（７０−１８）、（７０−１４）、（７０−１５）式によって修正すればよい。
【０１２１】
次に、弾性エネルギー力計算部１０４による不連続対応弾性力の計算について説明する。互いに異なるパラメトリック動きモデルは、一般に互いに異なる動きパラメータを有している。このため、両者の境界部分の弾性定数は０に設定した方が歪みが少ない。その境界部分は互いに異なるラベルを有する部分となる。すなわち、格子点ｎ1，ｎ2間の弾性定数は、（７１）式で示されることになる。不連続に対応した弾性力（以下、「不連続対応弾性力」という）は、（７１）式と（３８）式により（７２）式で定義することができる。
【０１２２】
【数６０】
2007257518.tif 000062
【数６１】
2007257518.tif 000063
このため、弾性エネルギー力計算部１０４では、（７１）式に基づき弾性係数を決定し、後述する第２の運動方程式の項となる不連続対応弾性力を（７２）式により計算している。
【０１２３】
次に、補正項計算部１１２による補正項の計算について説明する。パラメトリック動きモデルによって定まる動きと各格子点の動きが異なる場合には、各格子点をパラメトリック動きモデルに近似させる補正をおこなう必要がある。このため、後述する第２の運動方程式にかかる近似を行うための補正項を加え、補正項計算部１１２によって補正項の計算を行っている。かかる近似の補正が線形な弾性エネルギーに従う場合とすれば、補正項は、（７３）式のように、パラメトリック動きモデルによって定まる動きベクトルと各格子点の動きベクトルの差を入力とした線形関数の出力で示される補正力として定義される。ここで、ｋ2は、補正の強度を示している。
【０１２４】
【数６２】
2007257518.tif 000064
また、パラメトリック動きモデルによって定まる動きと各格子点の動きとの差が所定値以上となる場合には、（７３）式に示す補正項では補正力が必要以上に大きくなってしまう。このような場合には、パラメトリック動きモデルによって定まる動きベクトルと各格子点の動きベクトルの差を入力とした非線形関数の出力で示される補正項、例えば、（７４）式のような対数をとった補正項を採用することにより、補正項の値が膨大になることを防止することができる。
【０１２５】
【数６３】
2007257518.tif 000065
次に、動きクラスタリング処理部１０９による動きクラスタリング処理について説明する。更新式（５５）、（５６）式による局所的な動的マッチング処理が完了した時点では、ラベルの割り当て処理が行われていないため、動きパラメータの推定を行うことができない。このため、更新式（５５）、（５６）式による局所的な動的マッチング処理が完了して最初に動きパラメータの推定処理を行う場合には、動きパラメータの推定処理の実行前に、動きクラスタリング処理部１０９によって、更新式（５５）、（５６）式によって求められた動き推定をクラスタリングすることによりラベルの初期値として使用している。ここで、更新式（５５）、（５６）式による動き推定のクラスタリング処理としては、輝度等の画素値に基づいて画像をｋ−Ｍｅａｎｓ法によって複数個にクラスタリングを行って画像を分割し、そのクラスタリング結果によって複数個の動きパラメータを求め、求めた動きパラメータを初期値としてｋ−Ｍｅａｎｓ法による数値解析処理により動きパラメータのラベルを求める。
【０１２６】
次に、動きパラメータ割り当て部１１１による動きパラメータ割り当て処理、すなわちラベル割り当て処理について説明する。動きパラメータ割り当て部１１１による動きパラメータ割り当て処理では、格子点ｎに対する動きパラメータの割り当て問題として、（７５）式に示す最適化問題を解放する。
【０１２７】
【数６４】
2007257518.tif 000066
ここで、Ｇはコスト関数であり、本実施の形態では、非特許文献「A. M. Tekalp, “digital video processing”, Prentice Hall, 1995，pp. 103」に示されている、（７６）式のようなＭＰＣ（Maximum Pel Count）を用いている。
【０１２８】
【数６５】
2007257518.tif 000067
ここで、Ｂ⊂Λ2は、格子点ｎを中心とした局所領域(ブロック)を表す集合であり、Ｎ（Ｂ）は集合の要素数である。また、非特許文献「A. M. Tekalp, “digital video processing”, Prentice Hall, 1995，pp. 102-103」に示されているＭＳＥ(Mean Squared Error), ＭＡＤ(Mean Absolute Error)などを用いても良い。
【０１２９】
ここで、動きパラメータを正確に割り当てるためには、動きパラメータの精度が高い必要があることが推察される。図１１に示すように、マッチング対象領域に対して画像マッチングが正しい領域（図１１中の白い領域）が広いほど、その領域に対する動きパラメータの推定精度も高くなる。これは最小二乗法の特性である（中心極限定理によりサンプル数が多いほど推定精度は向上する）。このため、画像マッチングが行われていない領域（図中の灰色領域）に対して、この動きパラメータを割り当てると、正確な画像マッチングとなりやすい。逆に、マッチング領域が狭い場合には、動きパラメータの推定精度は低下する。この場合、低精度の動きパラメータを拡張した場合でも、真の動きパラメータとのズレが大きくなり、正確な画像マッチングを行うことは困難である。
【０１３０】
ここで、画面全体での最適化を行うことも可能である。画像処理の有する局在性に着目し、画像の局在性をマルコフ性と等価なギブス分布によって定式化し、それを事前分布としてベイズの定理により、求めたい最適化問題を事後分布最大化問題として定式化すればよい。またその事後分布最大化問題は基本的にＮＰ困難であり、大域的最適解を得るのは難しく、通常はシミュレイティッドアニーリングなどにより求解する。さらに効率よく解法するため、非特許文献「Y. Boykov, O. Veksler and R. Zabih, "Fast Approximate Energy Minimization via Graph Cuts", IEEE Trans. on Pattern Analysis and Machine Int., Vol. 23, No. 11, pp. 1222-1239, 2001」に詳述されているグラフ理論のグラフカットを用いて定式化することも可能である。
【０１３１】
次に、第２解法処理部１０８によって解法処理を行う第２の運動方程式について説明する。第２の運動方程式は、（７７）式に示されるように、（５５）、（５６）式の第１の運動方程式に、（７３）式または（７４）式の補正項が加わり、さらに弾性エネルギー力として、（７１）、（７２）式に示す不連続対応弾性力の項を有している。第２解法処理部１０８では、かかる第２の運動方程式を、第１の運動方程式と同様に、オイラー法等の離散変数法を用いた数値解析処理により解法している。
【０１３２】
【数６６】
2007257518.tif 000068
次に、以上のように構成された実施の形態１にかかる画像マッチング装置１００による画像マッチング処理について説明する。図１２は、実施の形態１にかかる画像マッチング装置１００による画像マッチングの全体処理の手順を示すフローチャートである。まず、カウント回数Ｌを１に初期化し（ステップＳ１２０１）、第１解法処理部１０１、画像エネルギー力計算部１０３、弾性エネルギー力計算部１０４、摩擦力計算部１０５による（５５）、（５６）式による局所的な動的マッチング処理を行う（ステップＳ１２０２）。
【０１３３】
動的マッチング処理が完了したら、動きクラスタリング処理１０９によって動きクラスタリング処理を行い（ステップＳ１２０３）、参照画像の画面を複数個に分割して動きパラメータのラベルの初期値を求める。次に、動きパラメータ推定部１１０による動きパラメータ推定処理を行い（ステップＳ１２０４）、各ラベルの領域に対する動きパラメータを推定する。そして、動きパラメータ修正部１１３によって、動きパラメータ推定部１１０で推定した動きパラメータを、（７０−１８）、（７０−１４）、（７０−１５）式によって修正する（ステップＳ１２０５）。そして、動きパラメータ割り当て部１１１によって、修正した動きパラメータから最適な動きパラメータを決定して、決定された最適な動きパラメータを、ラベルに対応する領域内の各第２の格子点に対して割り当てる動きパラメータ割り当て処理を行う（ステップＳ１２０６）。
【０１３４】
次に、補正項を付加した動的マッチング処理、すなわち補正項計算部１１２および第２解法処理部１０８による（７７）式の第２の運動方程式の解法処理を行ってｙn(LT2)）を求め（ステップＳ１２０７）、フレームメモリ１０６に保存する。そして、カウント回数Ｌを１増加して（ステップＳ１２０８）、カウント回数Ｌが予め定められた所定回数を超えたか否かを調べる（ステップＳ１２０９）。
【０１３５】
カウント回数Ｌが予め定められた所定回数を超えていない場合には（ステップＳ１２０９：Ｎｏ）、ステップＳ１２０４からＳ１２０８までの処理を繰り返す。すなわち、動きパラメータ処理部１１０による動きモデルの動きパラメータの推定処理、動きパラメータ修正部１１３による動きパラメータの修正処理、動きパラメータ割り当て部１１１による動きパラメータの割り当て処理、補正項計算部１１２および第２解法処理部１０８による局所的な動きの更新（第２の運動方程式の解法処理）を繰り返す。
【０１３６】
一方、ステップＳ１２０９において、カウント回数Ｌが予め定められた所定回数を超えた場合には（ステップＳ１２０９：Ｙｅｓ）、第２解法処理部１０８によって、すべての格子点ｎに対して、ｙnにｙn(T)を設定する（ステップＳ１２１０）。そして、マッピング処理部１０７により、対象画像と参照画像との対応関係、すなわち写像を求める（ステップＳ１２１１）。このような動き分割を用いた処理によって、画像の正確な動きを推定して、第１の格子点と第２の格子点との写像関係を求めて画像マッチングが行われる。
【０１３７】
次に、ステップＳ１２０２における局所的な動的マッチング処理について説明する。図１３は、実施の形態１にかかる局所的な動的マッチング処理の手順を示すフローチャートである。具体的には、先に示した（５５）、（５６）式を数値解析により解法することにより画像マッチング処理を実現している。
【０１３８】
まず、第１解法処理部１０１は、時間τ(0)＝０に設定し（ステップＳ１３０１）、初期値ｙn(0)＝ｎ、ｖn(0)＝０を設定する（ステップＳ１３０２）。これにより、（５６）式が実行される。
【０１３９】
次に、画像エネルギー力計算部１０３によって、すべての格子点ｎに対してｍステップにおける画像相関ポテンシャル力Ｆi(m)(n)を計算する（ステップＳ１３０３）。かかる画像エネルギー力Ｆi(m)(n)の計算処理については後述する。
【０１４０】
次いで、弾性エネルギー力計算部１０４によって、すべての格子点ｎに対してｍステップにおける弾性エネルギー力Ｆk(m)(n)を（３８）式により算出する（ステップＳ１３０４）。そして、摩擦力計算部１０５によって、すべての格子点ｎに対してｍステップにおける摩擦力［−μｖn(m)］を計算する（ステップＳ１３０５）。
【０１４１】
次に、第１解法処理部１０１によって、（５５）式の更新式を、ステップＳ１３０３〜Ｓ１３０５で求めた画像エネルギー力Ｆi(m)(n)、弾性エネルギー力Ｆk(m)(n)、摩擦力［−μｖn(m)］で更新する（ステップＳ１３０６）。
【０１４２】
次に、第１解法処理部１０１によって、ｙn(m)の値をフレームメモリ１０６に保存する（ステップＳ１３０７）。そして、第１解法処理部１０１によって、τ(m+1)=τ(m)+ｈと更新し（ステップＳ１３０８）、ｍを１だけ増加する（ステップＳ１３０９）。そして、τ(m+1)があらかじめ定められた時間Ｔを越えたか否かを判断し（ステップＳ１３１０）、越えていない場合には、上記ステップＳ１３０３からＳ１３０９を繰り返し実行する。
【０１４３】
一方、τ(m+1)がＴを越えた場合には、第１解法処理部１０１によって、すべての格子点ｎに対して、ｙnにｙn(T)を設定する（ステップＳ１３１１）。
【０１４４】
次に、ステップＳ１３０３における画像エネルギー力の計算処理について説明する。図１４は、画像エネルギー力計算部１０３による画像エネルギー力の計算処理の手順を示すフローチャートである。
【０１４５】
まず、画像エネルギー力計算部１０３は、（４１）式により、ｙn（τ）にもっとも近いサンプリング点を局所空間中心ｙcとして計算する（ステップＳ１４０１）。そして次に、画像エネルギー力計算部１０３は、（４２）式で定義される隣接空間Ｌを設定し（ステップＳ１４０２）、局所探索集合Ωを（４３）式によって計算する（ステップＳ１４０３）。
【０１４６】
次に、画像エネルギー力計算部１０３は、（４６）式により局所最適化計算としてｙminを計算し（ステップＳ１４０４）、正規化のために、ステップＳ１４０４で求めた局所最適化のｙminを用いて（４５）式によりｄ＝ｙmin−ｙn(m)を計算する（ステップＳ１４０５）。そして、画像エネルギー力計算部１０３は、ｙmin、ｄ等を用いて、（４４）式により、画像エネルギー力Ｆi(m)（ｎ）を計算する（ステップＳ１４０６）。
【０１４７】
次に、ステップＳ１２０３における動きクラスタリング処理部１０９により動きクラスタリング処理について説明する。図１５は、実施の形態１にかかる動きクラスタリング処理部１０９により動きクラスタリング処理の手順を示すフローチャートである。
【０１４８】
まず、動きクラスタリング処理部１０９は、局所的な動的マッチング処理の結果としての輝度等の画素値に基づいて、参照画像をｋ−Ｍｅａｎｓ法によって複数個にクラスタリングを行って画像を分割しラベルを求める（ステップＳ１５０１）。
【０１４９】
そして、そのクラスタリング結果によるラベルに基づき、初期動きパラメータを推定し、複数個の初期動きパラメータを求める（ステップＳ１５０２）。次に、求めた初期動きパラメータを初期値としてｋ−Ｍｅａｎｓ法による数値解析処理により動きパラメータのラベルαを求める（ステップＳ１５０３）。ここで、ラベルαは、複数個の初期動きパラメータに対応して複数個求められる。かかる動きクラスタリング処理は、図１２に示すように、局所的な動的マッチング処理が完了した直後に１回だけ実行される。
【０１５０】
次に、ステップＳ１２０４における動きパラメータ推定部１１０による動きパラメータ推定処理について説明する。図１６は、実施の形態１にかかる動きパラメータ推定部１１０による動きパラメータ推定処理の手順を示すフローチャートである。
【０１５１】
まず、動きパラメータ推定部１１０は、動きクラスタリング処理部１０９によって求めたラベルαに対して、全ての格子点ｎに対して、（６９−２）式によりＡを、（６９−３）式によりｂを計算する（ステップＳ１６０１）。そして、（６９−１）式を特異値分解法（もしくはＬＵ法）による数値解析により解法処理を行う（ステップＳ１６０２）。そして、全てのラベルαに対して、上記ステップＳ１６０１およびＳ１６０２の処理を繰り返して実行する（ステップＳ１６０３）。このようにして推定された動きパラメータは、ステップＳ１２０５における動きパラメータ修正部１１３によって、（７０−１８）、（７０−１４）、（７０−１５）式により修正される。
【０１５２】
次に、ステップＳ１２０６における動きパラメータ割り当て部１１１による動きパラメータ割り当て処理について説明する。図１７は、実施の形態１にかかる動きパラメータ割り当て部１１１による動きパラメータ割り当て処理の手順を示すフローチャートである。
【０１５３】
まず、動きパラメータ割り当て部１１１は、全てのラベルαに対して、（７６）式に示すコスト関数Ｇ(n)が最小になるα0を求める（ステップＳ１７０１）。かかる処理は、すなわち、（７５）式を実行することである。そして、動きパラメータ割り当て部１１１は、α0を格子点ｎに対するラベルとして設定する（ステップＳ１７０２）。かかる処理により、格子点ｎに最適な動きパラメータが割り当てられたことになる。そして、動きパラメータ割り当て部１１１は、上記ステップＳ１７０１およびＳ１７０２の処理を全ての格子点ｎに対して繰り返し実行する（ステップＳ１７０３）。これにより、全ての格子点ｎに対して最適な動きパラメータが割り当てられる。
【０１５４】
次に、ステップＳ１２０７における補正項を付加した動的マッチング処理について説明する。図１８は、実施の形態１にかかる補正項を付加した動的マッチング処理の手順を示すフローチャートである。具体的には、先に示した（７７）式を数値解析により解法することにより画像マッチング処理を実現している。
【０１５５】
まず、第２解法処理部１０８は、時間τ(0)＝０に設定し（ステップＳ１８０１）、初期値ｙn(0)＝ｎ、ｖn(0)＝０を設定する（ステップＳ１８０２）。これにより、（７７）式の第２式が実行される。
【０１５６】
次に、画像エネルギー力計算部１０３によって、すべての格子点ｎに対してｍステップにおける画像相関ポテンシャル力Ｆi(m)(n)を計算する（ステップＳ１８０３）。かかる画像エネルギー力Ｆi(m)(n)の計算処理については、図１４で説明した処理と同様に行われる。
【０１５７】
次いで、弾性エネルギー力計算部１０４によって、すべての格子点ｎに対してｍステップにおける不連続対応弾性力Ｆdk(m)(n)を（７１）、（７２）式により算出する（ステップＳ１８０４）。不連続対応弾性力Ｆdk(m)(n)の計算処理については後述する。
【０１５８】
次に、補正項計算部１１２によって、すべての格子点ｎに対してｍステップにおける補正項Ｆd(m)(n)を（７３）式により算出する（ステップＳ１８０５）。そして、摩擦力計算部１０５によって、すべての格子点ｎに対してｍステップにおける摩擦力［−μ(τ(m))ｖn(m)］を計算する（ステップＳ１８０６）。
【０１５９】
次に、第２解法処理部１０８によって、（７７）式の更新式を、ステップＳ１８０３〜Ｓ１８０６で求めた画像エネルギー力Ｆi(m)(n)、不連続対応弾性力Ｆdk(m)(n)、補正項Ｆd(m)(n)、摩擦力［−μ(τ(m))ｖn(m)］で更新する（ステップＳ１８０７）。
【０１６０】
次に、第２解法処理部１０８によって、ｙn(m)の値をフレームメモリ１０６に保存する（ステップＳ１８０８）。そして、第２解法処理部１０８によって、τ(m+1)=τ(m)+ｈと更新し（ステップＳ１８０９）、ｍを１だけ増加する（ステップＳ１８１０）。そして、τ(m+1)があらかじめ定められた時間Ｔ1を越えたか否かを判断し（ステップＳ１８１１）、越えていない場合には、上記ステップＳ１８０３からＳ１８１０を繰り返し実行する。
【０１６１】
次に、ステップＳ１８０４における弾性エネルギー力計算部１０４による不連続対応弾性力Ｆdk(m)(n)の計算処理について説明する。図１９は、実施の形態１にかかる弾性エネルギー力計算部１０４による不連続対応弾性力Ｆdk(m)(n)の計算処理の手順を示すフローチャートである。
【０１６２】
まず、弾性エネルギー力計算部１０４は、ｍステップにおける格子点ｎ1のラベルｚ(m)（ｎ1）と格子点ｎ2のラベルｚ(m)（ｎ2）が等しいか否かを調べる（ステップＳ１９０１）。これは、格子点ｎ1と格子点ｎ2は同一の動きパラメータの動きモデルに属するか否かを調べている。
【０１６３】
そして、格子点ｎ1のラベルｚ(m)（ｎ1）と格子点ｎ2のラベルｚ(m)（ｎ2）が等しい場合（ステップＳ１９０１：Ｙｅｓ）、すなわち、格子点ｎ1と格子点ｎ2は同一の動きパラメータの動きモデルに属する場合には、格子点ｎ1と格子点ｎ2の間の弾性係数ｋ(m)（ｎ1,ｎ2）をk1,に設定する（ステップＳ１９０２）。そして、不連続対応弾性力Ｆdk(m)(n)を（７２）式により計算する（ステップＳ１９０４）。
【０１６４】
一方、ステップＳ１９０１において、格子点ｎ1のラベルｚ(m)（ｎ1）と格子点ｎ2のラベルｚ(m)（ｎ2）が等しくない場合（ステップＳ１９０１：Ｎｏ）、すなわち、格子点ｎ1と格子点ｎ2が異なる動きパラメータの動きモデルに属する場合には、弾性係数ｋ(m)（ｎ1,ｎ2）を０に設定する（ステップＳ１９０３）。これにより、（７２）式によって不連続対応弾性力Ｆdk(m)(n)は０となる。なお、ステップＳ１９０１からＳ１９０３までの処理は、（７１）式を実行する処理である。
【０１６５】
このようにして計算された不連続対応弾性力は、ステップＳ１８０７で、更新式（７７）式の項に入力されることになる。
【０１６６】
このように実施の形態１にかかる画像マッチング装置１００では、更新式（５５）、（５６）式による局所的な動的マッチング処理が終了した後、参照画像を複数の領域に分割して各領域ごとに、第２の格子点の動きを規定した動きモデルを分類するための動きパラメータの最適な動きパラメータを求めて動きモデルを示す第２の運動方程式（７７）式の数値解析による解法処理を行って、対象画像と前記参照画像との対応関係を求めているので、領域境界の不連続性と領域内の一様性をロバストに表現できるようになり、より高精度な画像マッチングを行うことができ、その結果、より高画質な補間フレームを生成することができる。
【０１６７】
（実施の形態２）
実施の形態２にかかる画像マッチング装置は、動きパラメータの推定処理において、ロバスト統計の手法であるＭ推定を利用したものである。
【０１６８】
実施の形態１では、動きパラメータ推定処理において、動きパラメータと得られた動きベクトルの誤差がガウス分布に従っていると仮定して最小二乗推定法を用いた動きパラメータの推定を行っていた。しかしながら、動きパラメータと動きベクトルの誤差が単純にガウス分布に従っているといえない場合には、ロバスト統計の手法であるＭ推定を利用することが好ましい。このような場合としては、例えば、画像の一部に大きな外れ値が含まれているような場合であり、動的マッチングにおいて、動きと動きの境界は歪みによってこのような状況になりやすい。
【０１６９】
実施の形態２にかかる画像マッチング装置の構成は実施の形態１と同様である。また、実施の形態２における画像マッチングの全体処理、局所的な動的マッチング処理、動きクラスタリング処理、動きパラメータ割り当て処理、補正項を付した動的マッチング処理については、実施の形態１と同様に行われる。
【０１７０】
Ｍ推定における標準正規変数ｚ(ｎ)は、（７８）式で示される。TukeyのBiweight法によると重み関数ｗ(ｎ)は、（７９）式で示される
【数６７】
2007257518.tif 000069
ここで、ｃ（＞０）は定数であり、正規分布の場合にはｃ＝６．０となる。動きパラメータ推定部１１０は、重み関数ｗ(ｎ)を用いて標本平均ｚ(i)、標本偏差ｓ(i)を（８０）式のように計算する。
【０１７１】
【数６８】
2007257518.tif 000070
ここで上付の(i)は、イテレーションを示している。また、動きパラメータ推定部１１０では、母分散は、標本偏差標本偏差ｓ(i)を（８１−１）、（８１−２）式で示すようにプラグインして推定することとする。
【０１７２】
【数６９】
2007257518.tif 000071
そして、動きパラメータ推定部１１０によって、母標準偏差の計算と重み関数の計算を交互にＩ回繰り返すことによってＭ推定による動きパラメータの推定が行われる。
【０１７３】
図２０は、実施の形態２にかかる動きパラメータ推定部１１０による動きパラメータ推定処理の手順を示すフローチャートである。まず、動きパラメータ推定部１１０は、重みｗ（ｎ）を１に初期化する（ステップＳ２００１）。そして、動きクラスタリング処理部１０９によって求めたラベルαに対して、全ての格子点ｎに対して、（８１−２）式によりＡを、（８１−３）式によりｂを計算する（ステップＳ２００２）。そして、（８１−１）式を特異値分解法（もしくはＬＵ法）による数値解析により解法処理を行う（ステップＳ２００３）。
【０１７４】
次に、動きパラメータ推定部１１０は、標準正規変数ｚ（ｎ）を（７８）式に従って計算する（ステップＳ２００４）。そして、重み関数ｗ(ｎ)を用いて標本平均ｚ(i)、標本偏差ｓ(i)を（８０）式に従って計算する（ステップＳ２００５）。次いで、動きパラメータ推定部１１０は、重みｗ（ｎ）を更新し（ステップＳ２００６）、イテレーションｉを１増加する（ステップＳ２００７）。そして、ｉが所定回数Ｉを越えたか否かを調べ（ステップＳ２００８）、越えていない場合には、ステップＳ２００２からＳ２００７までの処理を繰り返し実行する。
【０１７５】
一方、ステップＳ２００８において、ｉが所定回数Ｉを越えた場合には、全てのラベルαに対して、上記ステップＳ１６０１およびＳ１６０２の処理を繰り返して実行する（ステップＳ２００９）。
【０１７６】
このように実施の形態２にかかる画像マッチング装置では、動きパラメータの推定処理において、ロバスト統計の手法であるＭ推定を利用しているので、動きパラメータと動きベクトルの誤差がガウス分布に従っていない場合にも、領域境界の不連続性と領域内の一様性をロバストに表現できるようになり、より高精度な画像マッチングを行うことができ、その結果、より高画質な補間フレームを生成することができる。
【０１７７】
（実施の形態３）
実施の形態３にかかる画像マッチング装置は、第２の格子点が対象画像上には存在するが参照画像上には存在しない領域であるオクルージョン領域に含まれる場合には、当該格子点において動きパラメータの割り当てをおこなわず、オクルージョン領域に含まれない格子点に対してのみ動きパラメータの割り当てを行うものである。
【０１７８】
対象画像上には存在するが参照画像上には存在しない領域は、coveredなオクルージョン領域といい、本来的に画像マッチングを行うことができない領域である。実施の形態１にかかる画像マッチング装置１００では、coveredオクルージョン領域にも動きパラメータの割り当てがおこなわれるため、動きパラメータ推定の結果が歪んでしまう場合がある。
【０１７９】
実施の形態３にかかる画像マッチング装置では、coveredオクルージョン領域をあらかじめ検出しておき、そのオクルージョン領域の格子点には動きパラメータの割り当てをおこなわないようにしている。
【０１８０】
図２１は、実施の形態３にかかる画像マッチング装置２１００の構成を示すブロック図である。実施の形態３にかかる画像マッチング装置２１００では、オクルージョン検出部２１０１を備えている点が実施の形態１と異なり、他の構成については実施の形態１にかかる画像マッチング装置の構成と同様である。また、実施の形態３における画像マッチングの全体処理、局所的な動的マッチング処理、動きクラスタリング処理、動きパラメータ推定処理、補正項を付した動的マッチング処理については、実施の形態１と同様に行われる。
【０１８１】
オクルージョン検出部２１０１は、coveredオクルージョン領域を検出する処理部である。図２２および図２３は、coveredオクルージョン領域の状態を示す説明図である。coveredオクルージョン領域は、図２２に示すように、典型的に移動物体の前方に現れる。図２２に示す状態では、バネが縮んでいる領域がcoveredオクルージョン領域となる。
【０１８２】
このようなcoveredオクルージョン領域は、次のように定式化することができる。図２３に示すように、格子点の周囲の８格子を考える。そして、かかる８格子から形成される三角形の面積の和により、元の正方格子からの変形の度合いを調べる。かかる面積Ｓ（ｎ）をベクトルの外積により、（８２）式のように定義する。面積比ｒ（ｎ）は、正方格子の面積を４として、（８３）式で表される。なお、（８３）式において、×はベクトルの外積計算をあらわす。
【０１８３】
【数７０】
2007257518.tif 000072
そして、この面積比が閾値Ｔよりも小さければ、図２２に示すような縮んでいる領域である、すなわちcoveredオクルージョン領域の可能性が高い。本実施の形態にかかるオクルージョン検出部２１０１では、（８３）式に示す面積比を閾値Ｔと比較することにより、オクルージョン領域か否かを判定している。そして、動きパラメータ割り当て部１１１では、オクルージョン検出部２１０１によりcoveredオクルージョン領域であると判定された格子点に関しては、動きパラメータの割り当てをおこなわないように構成している。
【０１８４】
図２４は、実施の形態２にかかる画像マッチング装置２１００による動きパラメータ割り当て処理の手順を示すフローチャートである。
【０１８５】
まず、オクルージョン検出部２１０１は、格子点の周囲の８格子から形成される三角形の面積の総和Ｓ（ｎ）を（８２）式により計算する（ステップＳ２４０１）。次に、オクルージョン検出部２１０１は、Ｓ（ｎ）を用いて面積比ｒ（ｎ）を（８３）式により計算する（ステップＳ２４０２）。
【０１８６】
次に、オクルージョン検出部２１０１は、面積比ｒ（ｎ）が所定の閾値Ｔより小さいか否か、すなわちオクルージョン領域か否かを判定する（ステップＳ２４０３）。そして、面積比ｒ（ｎ）が所定の閾値Ｔより小さくない場合（ステップＳ２４０３：Ｎｏ）、すなわち、オクルージョン領域でない場合には、動きパラメータ割り当て部１１１は、全てのラベルαに対して、（７６）式に示すコスト関数Ｇ(n)が最小になるα0を求める（ステップＳ２４０４）。そして、動きパラメータ割り当て部１１１は、α0を格子点ｎに対するラベルとして設定する（ステップＳ２４０５）。かかる処理により、格子点ｎに最適な動きパラメータが割り当てられたことになる。
【０１８７】
一方、ステップＳ２４０３において、面積比ｒ（ｎ）が所定の閾値Ｔより小さい場合（ステップＳ２４０３：Ｙｅｓ）、すなわち、オクルージョン領域である場合には、ステップＳ２４０４およびＳ２４０５の最適な動きパラメータが割り当て処理は行わない。
【０１８８】
そして、動きパラメータ割り当て部１１１は、上記ステップＳ２４０１からＳ２４０５の処理を全ての格子点ｎに対して繰り返し実行する（ステップＳ２４０６）。これにより、オクルージョン領域でない全ての格子点ｎに対して最適な動きパラメータが割り当てられることになる。
【０１８９】
このように実施の形態３にかかる画像マッチング装置では、coveredオクルージョン領域をあらかじめ検出しておき、そのオクルージョン領域の格子点には動きパラメータの割り当てを行わないので、coveredオクルージョンによる歪みを最小限に抑えることができ、より高精度な画像マッチングを行うことができ、その結果、より高画質な補間フレームを生成することができる。
【０１９０】
（実施の形態４）
実施の形態４にかかる画像マッチング装置では、信頼性の低い領域に対しては動きパラメータの割り当てを行わないものである。
【０１９１】
実施の形態１にかかる画像マッチング装置では、動きパラメータ割り当て処理の最適化を全て信頼していた。しかしながら、実際には動きパラメータ割り当て処理の最適化の信頼度が低い領域も存在する。例えば、動きパラメータの最適値に近似する値が二つ以上存在するような場合には、信頼性が低いと言える。すなわち、異なる動きパラメータに対して、最適値に近似する値が二つ以上存在するような場合は、コスト関数面が多峰性になっていると考えられ、動きパラメータの最適値の信頼度は低いと考えられる。本実施の形態では、このような低信頼領域に対しては動きパラメータの割り当て処理を行わないように構成している。
【０１９２】
実施の形態４かかる画像マッチング装置の構成は、実施の形態１と同様である。また、実施の形態４における画像マッチングの全体処理、局所的な動的マッチング処理、動きクラスタリング処理、動きパラメータ推定処理、補正項を付した動的マッチング処理については、実施の形態１と同様に行われる。
【０１９３】
本実施の形態にかかる動きパラメータ割り当て部１１１では、格子点が低信頼領域に含まれるか否かを判定して低信頼性領域に含まれる場合には、動きパラメータの割り当てをおこなわない処理を行っている。
【０１９４】
低信頼領域の判定は、（８４）式により定式化することができる。
【０１９５】
【数７１】
2007257518.tif 000073
ここで、α0は、求めた最適なラベルであり、α1は次に最適なラベルである。また、ε（＞０）は定数である。ただし、動きパラメータの値が近似する場合には、（８４）式に該当する場合があるため、（８５）式を満たすほど動きパラメータが近似する場合には、（８４）式は適用しないこととする。
【０１９６】
【数７２】
2007257518.tif 000074
図２５は、実施の形態４にかかる動きパラメータ割り当て部１１１による動きパラメータ割り当て処理の手順を示すフローチャートである。まず、動きパラメータ割り当て部１１１は、全てのラベルαに対して、（７６）式に示すコスト関数Ｇ(n)が最小になるα0を求める（ステップＳ２５０１）。次に、動きパラメータ割り当て部１１１は、α0より一つ大きいα1、すなわち次に最適なラベルを求める（ステップＳ２５０２）。
【０１９７】
次に、動きパラメータが近似するか否か、すなわち（８５）式が成立するか否かを調べる（ステップＳ２５０３）。そして、（８５）式が成立しない場合には、（８６）式が成立するか否かにより低信頼性領域か否かを調べる（ステップＳ２５０４）。
【０１９８】
そして、（８６）式が成立せず低信頼性領域でないと判断された場合には（ステップＳ２５０４：Ｎｏ）、ステップＳ２５０１で求めたα0を格子点ｎに対する最適なラベルと設定する（ステップＳ２５０５）。一方、ステップＳ２５０４において、（８６）式が成立して低信頼性領域であると判断された場合には（ステップＳ２５０４：Ｙｅｓ）、ステップＳ２５０５の処理は行わず、最適なラベルの割り当て、すなわち最適な動きパラメータの割り当ては行われない。
【０１９９】
また、ステップＳ２５０３において、動きパラメータが近似する場合、すなわち（８５）式が成立する場合には（ステップＳ２５０３：Ｙｅｓ）、ステップＳ２４０４の低信頼性領域の判定処理を行わず、ステップＳ２５０１で求めたα0を格子点ｎに対する最適なラベルと設定する（ステップＳ２５０５）。
【０２００】
かかる処理により、格子点ｎに最適な動きパラメータが割り当てられたことになる。そして、動きパラメータ割り当て部１１１は、上記ステップＳ２５０１からＳ２５０５までの処理を全ての格子点ｎに対して繰り返し実行する（ステップＳ２５０６）。
【０２０１】
このように実施の形態４にかかる画像マッチング装置では、信頼性の低い領域に対しては動きパラメータの割り当てを行わないので、領域境界の不連続性と領域内の一様性をロバストに表現できるようになり、より高精度な画像マッチングを行うことができ、その結果、より高画質な補間フレームを生成することができる。
【０２０２】
実施の形態１〜４の画像マッチング装置は、ＣＰＵなどの制御装置と、ＲＯＭ（ReadOnlyMemory）やＲＡＭなどの記憶装置と、ＨＤＤ、ＣＤドライブ装置などの外部記憶装置と、ディスプレイ装置などの表示装置と、キーボードやマウスなどの入力装置を備えており、通常のコンピュータを利用したハードウェア構成となっている。
【０２０３】
実施の形態１〜４の画像マッチング装置で実行される画像マッチングプログラムは、インストール可能な形式又は実行可能な形式のファイルでＣＤ−ＲＯＭ、フレキシブルディスク（ＦＤ）、ＣＤ−Ｒ、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）等のコンピュータで読み取り可能な記録媒体に記録されて提供される。
【０２０４】
また、実施の形態１〜４の画像マッチング装置で実行される画像マッチングプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成しても良い。また、本実施形態の画像マッチング装置で実行される画像マッチングプログラムをインターネット等のネットワーク経由で提供または配布するように構成しても良い。
【０２０５】
また、本実施の形態の画像マッチングプログラムを、ＲＯＭ等に予め組み込んで提供するように構成してもよい。
【０２０６】
実施の形態１〜４の画像マッチング装置で実行される画像マッチングプログラムは、上述した各部（第１解法処理部１０１と、画像エネルギー力計算部１０３と、弾性エネルギー力計算部１０４と、摩擦力計算部１０５と、第２解法処理部１０８と、動きクラスタリング処理部１０９と、動きパラメータ推定部１１０と、動きパラメー割り当て部１１１と、補正項計算部１１２と、マッピング処理部１０７、オクルージョン検出部２１０１）を含むモジュール構成となっており、実際のハードウェアとしてはＣＰＵ（プロセッサ）が上記記憶媒体から画像マッチングプログラムを読み出して実行することにより上記各部が主記憶装置上にロードされ、上記各部（第１解法処理部１０１と、画像エネルギー力計算部１０３と、弾性エネルギー力計算部１０４と、摩擦力計算部１０５と、第２解法処理部１０８と、動きクラスタリング処理部１０９と、動きパラメータ推定部１１０と、動きパラメー割り当て部１１１と、補正項計算部１１２と、マッピング処理部１０７、オクルージョン検出部２１０１）が主記憶装置上に生成されるようになっている。
【０２０７】
なお、本発明は、上記実施の形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化することができる。また、上記実施の形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成することができる。例えば、実施の形態に示される全構成要素からいくつかの構成要素を削除してもよい。さらに、異なる実施の形態にわたる構成要素を適宜組み合わせても良い。
【図面の簡単な説明】
【０２０８】
【図１】実施の形態１にかかる画像マッチング装置の構成を示すブロック図である。
【図２】対象画像上の曲面から参照画像上の曲面への写像を示す説明図である。
【図３】対象画像上の曲面Ｘの点ｘから参照画像上の曲面Ｙの点ｙへの写像を示す説明図である。
【図４】ポテンシャルエネルギーの最小化の点を探索する場合の静的な最適点の説明をするための模式図である。
【図５】ポテンシャルエネルギーの最小化の点を探索する場合の大域的な最適点の説明をするための模式図である。
【図６】画像エネルギー力の概念を示す説明図である。
【図７】画像エネルギー力の概念を示す説明図である。
【図８】動きモデルを利用した動き推定の画像上における状態を示す説明図である。
【図９】画像上のオブジェクトの一部が誤った動きで画像マッチングされている状態を示す説明図である。
【図１０】実施の形態１にかかるパラメトリック動きモデルを導入して動き分割を用いた画像マッチングの概念的な流れを示す説明図である。
【図１１】マッチング対象領域と動きパラメータの推定精度の関係を説明するための模式図である。
【図１２】実施の実施の形態１にかかる画像マッチング装置１００による画像マッチングの全体処理の手順を示すフローチャートである。
【図１３】実施の形態１にかかる局所的な動的マッチング処理の手順を示すフローチャートである。
【図１４】画像エネルギー力計算部１０３による画像エネルギー力の計算処理の手順を示すフローチャートである。
【図１５】実施の形態１にかかる動きクラスタリング処理部１０９により動きクラスタリング処理の手順を示すフローチャートである。
【図１６】実施の形態１にかかる動きパラメータ推定部１１０による動きパラメータ推定処理の手順を示すフローチャートである。
【図１７】実施の形態１にかかる動きパラメータ割り当て部１１１による動きパラメータ割り当て処理の手順を示すフローチャートである。
【図１８】実施の形態１にかかる補正項を付加した動的マッチング処理の手順を示すフローチャートである。
【図１９】実施の形態１にかかる弾性エネルギー力計算部１０４による不連続対応弾性力Ｆdk(m)(n)の計算処理の手順を示すフローチャートである。
【図２０】実施の形態２にかかる動きパラメータ推定部１１０による動きパラメータ推定処理の手順を示すフローチャートである。
【図２１】実施の形態３にかかる画像マッチング装置２１００の構成を示すブロック図である。
【図２２】coveredオクルージョン領域を説明するための模式図である。
【図２３】coveredオクルージョン領域を説明するための模式図である。
【図２４】実施の形態２にかかる画像マッチング装置２１００による動きパラメータ割り当て処理の手順を示すフローチャートである。
【図２５】実施の形態４にかかる動きパラメータ割り当て部１１１による動きパラメータ割り当て処理の手順を示すフローチャートである。
【符号の説明】
【０２０９】
１００，２１００  画像マッチング装置
１０１  第１解法処理部
１０３  画像エネルギー力計算部
１０４  弾性エネルギー力計算部
１０５  摩擦力計算部
１０６  フレームメモリ
１０７  マッピング処理部
１０８  第２解法処理部
１０９  動きクラスタリング処理部
１１０  動きパラメータ推定部
１１１  動きパラメー割り当て部
１１２  補正項計算部
２１０１  オクルージョン検出部
【図１】
2007257518.tif 000075
【図２】
2007257518.tif 000076
【図３】
2007257518.tif 000077
【図４】
2007257518.tif 000078
【図５】
2007257518.tif 000079
【図６】
2007257518.tif 000080
【図７】
2007257518.tif 000081
【図８】
2007257518.tif 000082
【図９】
2007257518.tif 000083
【図１０】
2007257518.tif 000084
【図１１】
2007257518.tif 000085
【図１２】
2007257518.tif 000086
【図１３】
2007257518.tif 000087
【図１４】
2007257518.tif 000088
【図１５】
2007257518.tif 000089
【図１６】
2007257518.tif 000090
【図１７】
2007257518.tif 000091
【図１８】
2007257518.tif 000092
【図１９】
2007257518.tif 000093
【図２０】
2007257518.tif 000094
【図２１】
2007257518.tif 000095
【図２２】
2007257518.tif 000096
【図２３】
2007257518.tif 000097
【図２４】
2007257518.tif 000098
【図２５】
2007257518.tif 000099