(57)【要約】【課題】大域的な最適点の探索と収束性を両立させて、より高精度な画像マッチングを行うこと。【解決手段】画像マッチング装置において、参照画像の各第２の格子点の位置と第２の格子点に対応する対象画像の各第１の格子点の位置とに基づいたポテンシャルエネルギーの勾配により第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算部１０３と、各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する弾性エネルギー力計算部１０４と、摩擦係数を可変制御する摩擦係数制御部１０２と、各第２の格子点に作用する摩擦力を摩擦係数を利用して計算する摩擦力計算部１０５と、画像エネルギー力と弾性エネルギー力と摩擦力とによる各第２の格子点に関する運動方程式の解法処理を行う解法処理部１０１と、対象画像と参照画像との対応関係を求めるマッピング処理部１０７とを備えた。【選択図】        図１
2007188329.tif 000002
【特許請求の範囲】
【請求項１】
対象画像と参照画像との間の対応関係を求める画像マッチング装置であって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算部と、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する弾性エネルギー力計算部と、
前記各第２の格子点に作用する摩擦力についての摩擦係数を可変制御する摩擦係数制御部と、
前記各第２の格子点に作用する摩擦力を、前記摩擦係数を利用して計算する摩擦力計算部と、
前記画像エネルギー力と前記弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する運動方程式の解法処理を行う解法処理部と、
前記解法処理部による解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理部と、
を備えたことを特徴とする画像マッチング装置。
【請求項２】
前記摩擦係数制御部は、前記運動方程式における時間の経過に基づいて増加するように前記摩擦係数を計算することを特徴とする請求項１に記載の画像マッチング装置。
【請求項３】
前記摩擦係数制御部は、前記運動方程式における時間についての複数の線形関数に従って増加するように前記摩擦係数を計算することを特徴とする請求項２に記載の画像マッチング装置。
【請求項４】
前記摩擦係数制御部は、前記時間の経過に従って段階的に増加するように前記摩擦係数を計算することを特徴とする請求項２に記載の画像マッチング装置。
【請求項５】
前記摩擦係数制御部は、前記時間についての指数関数に従って増加するように前記摩擦係数を計算することを特徴とする請求項２に記載の画像マッチング装置。
【請求項６】
前記摩擦係数制御部は、前記運動方程式における時間の開始時において前記摩擦係数を負値とし、前記時間の経過に従って増加するように前記摩擦係数を計算することを特徴とする請求項２に記載の画像マッチング装置。
【請求項７】
前記摩擦係数制御部は、前記運動方程式における時間の開始時において前記摩擦係数を負値とし、前記時間についての線形関数に従って増加するように前記摩擦係数を計算することを特徴とする請求項６に記載の画像マッチング装置。
【請求項８】
前記摩擦係数制御部は、前記運動方程式における時間の開始時において前記摩擦係数を負値とし、前記時間の経過に従って段階的に増加するように前記摩擦係数を計算することを特徴とする請求項６に記載の画像マッチング装置。
【請求項９】
前記ポテンシャルエネルギーが平衡状態に達しているか否かを判定する平衡状態判定部を更に備え、
前記摩擦係数制御部は、前記平衡状態判定部による判定結果に基づいて前記摩擦係数を計算することを特徴とする請求項１に記載の画像マッチング装置。
【請求項１０】
前記摩擦係数制御部は、前記ポテンシャルエネルギーが平衡状態に達している場合には、平衡状態に達していない場合の前記摩擦係数の値よりも大きい値となるように前記摩擦係数を計算することを特徴とする請求項９に記載の画像マッチング装置。
【請求項１１】
前記平衡状態判定部は、前記第２の格子点の座標の分布範囲の分散が予め定められた範囲内にある場合には、前記ポテンシャルエネルギーが平衡状態に達していると判断することを特徴とする請求項１０に記載の画像マッチング装置。
【請求項１２】
対象画像と参照画像との間の対応関係を求める画像マッチング方法であって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算ステップと、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する弾性エネルギー力計算ステップと、
前記各第２の格子点に作用する摩擦力についての摩擦係数を可変制御する摩擦係数制御ステップと、
前記各第２の格子点に作用する摩擦力を、前記摩擦係数を利用して計算する摩擦力計算ステップと、
前記画像エネルギー力と前記弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する運動方程式の解法処理を行う解法処理ステップと、
前記解法処理ステップによる解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理ステップと、
を含むことを特徴とする画像マッチング方法。
【請求項１３】
対象画像と参照画像との間の対応関係を求める画像マッチング方法であって、
前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算ステップと、
前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する弾性エネルギー力計算ステップと、
前記各第２の格子点に作用する摩擦力についての摩擦係数を可変制御する摩擦係数制御ステップと、
前記各第２の格子点に作用する摩擦力を、前記摩擦係数を利用して計算する摩擦力計算ステップと、
前記画像エネルギー力と前記弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する運動方程式の解法処理を行う解法処理ステップと、
前記解法処理ステップによる解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理ステップと、
をコンピュータに実行させることを特徴とする画像マッチングプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象画像と参照画像の２つの画像からそれぞれの対応点を検出して対応付ける画像マッチング装置、画像マッチング方法および画像マッチングプログラムに関する。
【背景技術】
【０００２】
一つの対象画像から他方の参照画像の各画素それぞれの対応点を検出して対応関係を求める画像マッチング技術は、例えば、動画像における動き検出、ステレオマッチング、画像モーフィング、画像認識、動画像符号化等の画像処理技術の分野において利用されている。
【０００３】
このような画像マッチング技術としては、オプティカルフロー手法、ブロックベース手法、勾配法、ベイジアンメソッドという主として４種類に分類することができる（例えば、非特許文献１参照）。
【０００４】
オプティカルフロー手法は、「輝度の変化は一定である」というオプティカルフロー式を導出しそのオプティカルフロー式を拘束条件としてフローを求めるものである。ブロックベース手法は、画像を所定のブロックに分割し、ブロック毎のテンプレートマッチングによって動きを求める手法である。勾配法は、画像の輝度勾配が減少する方向にマッチングをおこなう手法である。ベイジアンメソッドは、確率的に尤もらしいマッチングを求める手法である。
【０００５】
また、従来の画像マッチング技術として、複数の多重解像度フィルタを用い、複数の多重解像度画像ピラミッドを生成し、生成された画像ピラミッドを上位の階層から下位の階層に順にマッチング処理を行うことによって、大きな動きから小さな動きまで画像の対応付けを行えるロバスト性の高い画像マッチング技術がある（例えば、特許文献１参照）。
【０００６】
【非特許文献１】A. Murat Tekalp, “Digital Video Processing”, Prentice Hall, 1995
【特許文献１】特許第２９２７３５０号公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、このような従来の画像マッチング技術には次のような問題がある。オプティカルフロー手法では、ノイズに敏感で速い動きに対応することが本質的に困難であるという問題があり、ブロックベースの手法では、画像の中のブロック毎に画像マッチングを行っているので、画像中のオブジェクトの平行移動などの動きには高い信頼性を有するが、画像中のオブジェクトが変形したり、回転する等の動きに対して本質的に対応することが難しいという問題がある。勾配法では、画像の輝度勾配が減少する方向にマッチングを行うため、安定してオブジェクトの動きを探索することが困難であるという問題がある。ベイジアンメソッドでは、大域的最適点の求め方が困難であるという問題がある。
【０００８】
一方、特許文献１に開示された技術では、原理的に複数の多重解像度フィルタを用い、多重解像度画像ピラミッドの最上位階層から最下位階層までのマッチング処理を行っているが、基本的には局所的な最適化処理であり大域的に最適になるとは限らない。
【０００９】
本発明は、上記に鑑みてなされたものであって、大域的な最適点の探索と収束性を両立させて、より高精度な画像マッチングを行うことができる画像マッチング装置、画像マッチング方法および画像マッチングプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
上述した課題を解決し、目的を達成するために、本発明は、対象画像と参照画像との間の対応関係を求める画像マッチング装置であって、前記対象画像上の複数の第１の格子点の各々と前記参照画像上で当該第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーを計算するとともに、前記各第２の格子点の位置と当該第２の格子点に対応する前記各第１の格子点の位置とに基づいて前記ポテンシャルエネルギーの勾配により前記第２の格子点が受ける画像エネルギー力を計算する画像エネルギー力計算部と、前記各第２の格子点と当該第２の格子点に隣接する他の第２の格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する弾性エネルギー力計算部と、前記各第２の格子点に作用する摩擦力についての摩擦係数を可変制御して計算する摩擦係数制御部と、前記各第２の格子点に作用する摩擦力を、前記摩擦係数を利用して計算する摩擦力計算部と、前記画像エネルギー力と前記弾性エネルギー力と前記摩擦力とによる前記各第２の格子点に関する運動方程式の解法処理を行う解法処理部と、前記解法処理部による解法処理の結果から前記対象画像と前記参照画像との対応関係を求めるマッピング処理部と、を備えたことを特徴とする。
【００１１】
また、本発明は、上記画像マッチング装置で実行することができる画像マッチング方法および画像マッチングプログラムである。
【発明の効果】
【００１２】
本発明によれば、摩擦係数を可変制御して計算して、参照画像上における各第２の格子点に作用する摩擦力を計算することで、大域的な最適点の探索と収束性を両立させて、よりより高精度な画像マッチングを行うことができるという効果を奏する。
【発明を実施するための最良の形態】
【００１３】
以下に添付図面を参照して、この発明にかかる画像マッチング装置、画像マッチング方法および画像マッチングプログラムの最良な実施の形態を詳細に説明する。
【００１４】
なお、以下の実施の形態にかかる画像マッチング装置は、動画像の動き補償を行って補償画像を生成する動き補償装置に適用する他、ステレオマッチング、画像モーフィング、画像認識、動画像符号化等の画像処理技術に適用することもできる。
【００１５】
（実施の形態１）
図１は、実施の形態１にかかる画像マッチング装置の構成を示すブロック図である。本実施の形態にかかる画像マッチング装置１００は、図１に示すように、解法処理部１０１と、画像エネルギー力計算部１０３と、弾性エネルギー力計算部１０４と、摩擦力計算部１０５と、摩擦係数制御部１０２と、マッピング処理部１０７と、フレームメモリ１０６とを主に備えている。
【００１６】
実施の形態１にかかる画像マッチング装置１００では、対象画像および参照画像のそれぞれを所定の格子状に分割した各格子点を力学上の質点として扱い、対象画像の第１の格子点と第１の格子点に１対１に対応する参照画像の第２の格子点が画像の相関関係に基づいたポテンシャルエネルギーを考え、第１の格子点と第２の格子点の位置に基づいたポテンシャルエネルギーの勾配によって第２の格子点が受ける画像エネルギー力と、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける弾性エネルギー力と、第２の格子点に生じる摩擦力とによる第２の格子点に関する運動方程式を画像マッチング状態遷移モデルとして、この運動方程式に対して解法処理を行い、第２の格子点の平衡状態を求めることにより第１の格子点と第２の格子点との写像関係を求めて画像マッチングを行っている。
【００１７】
そして、実施の形態１にかかる画像マッチング装置１００では、摩擦力を計算する際に、摩擦力の摩擦係数を可変制御することによって、大域的にも最適化処理を行うことを可能としている。
【００１８】
解法処理部１０１は、対象画像上の複数の第１の格子点の各々と参照画像上で第１の格子点に一対一に対応する複数の第２の格子点の各々との間で画像の相関関係に基づくポテンシャルエネルギーの勾配により第２の格子点が受ける画像エネルギー力と、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける弾性エネルギー力と、第２の格子点に生じる摩擦力とによる第２の格子点に関する運動方程式を、離散変数法としてのオイラー法による数値解析処理によって解法する処理部である。
【００１９】
マッピング処理部１０７は、解法処理部１０１による解法処理の結果から第２の格子点の平衡状態を求めることにより対象画像上の第１の格子点と参照画像上の第２の格子点との対応関係である写像関係を求める処理部である。
【００２０】
画像エネルギー力計算部１０３は、第１の格子点と第２の格子点の位置に基づいた上記ポテンシャルエネルギーの勾配によって第２の格子点が受ける画像エネルギー力を計算する処理部である。
【００２１】
弾性エネルギー力算出部１０４は、第２の格子点に隣接する格子点との間の弾性エネルギーから受ける弾性エネルギー力を計算する処理部である。
【００２２】
摩擦係数制御部１０２は、各第２の格子点に作用する摩擦力の計算に使用する摩擦係数を可変制御する処理部である。具体的には、摩擦係数制御部１０２は、運動方程式における時間の経過に従って、当該時間の線形関数に従って増加するように摩擦係数を計算することで摩擦係数を可変制御する。
【００２３】
摩擦力算出部１０５は、第２の格子点に生じる摩擦力を摩擦係数を用いて計算する処理部である。
【００２４】
フレームメモリ１０６は、入力された対象画像と参照画像と保存するメモリである。なお、各部の詳細な処理については後述する。
【００２５】
次に、本実施の形態にかかる画像マッチング装置１００で導入する画像マッチング状態遷移モデルについて説明する。
【００２６】
２次元ユークリッド空間上の点ｘ∈Ｅ2における対象画像、参照画像の画素値をそれぞれ（１）式で示す。
【００２７】
【数１】
2007188329.tif 000003
（１）式は実数空間における連続画像モデルであり、本実施の形態ではデジタル画像を扱うため、２次元格子空間上の点ｎ∈Ｎ2について（２）式に示す離散画像モデルを考える。
【００２８】
【数２】
2007188329.tif 000004
従って、連続画像モデルの画素値はデータ補間関数ｇを用いて、（３）式に示すように定義される。ここで、データ補間関数ｇとしては、例えば双曲形補間関数を使用することができる。
【００２９】
【数３】
2007188329.tif 000005
本実施の形態では、対象画像から参照画像への画像マッチング問題を、対象画像上の曲面Ｘ⊂Ｅ2から参照画像上の曲面Ｙ⊂Ｅ2への写像として求めることにする。図２は、対象画像上の曲面から参照画像上の曲面への写像を示す説明図である。図３は、対象画像上の曲面Ｘの点ｘから参照画像上の曲面Ｙの点ｙへの写像を示す説明図である。すなわち、図２および３に示すように、画像マッチングは、対象画像上の曲面Ｘの点が参照画像上の曲面Ｙのどの点に対応しているのかを表す写像を求める問題として扱う。
【００３０】
このような写像に関する条件を導出すると以下のようになる。パラメータ（ｕ，ｖ）を用いて、曲面Ｘ上の点をｘ（ｕ，ｖ）∈Ｘ、曲面Ｙ上の点をｙ（ｕ，ｖ）∈Ｙと表し、パラメータ（ｕ，ｖ）と点ｘが同じスケールを有する場合、曲面Ｘ上の点ｘ（ｕ，ｖ）は、パラメータ（ｕ，ｖ）によって、（４）式で表すことができる。
【００３１】
【数４】
2007188329.tif 000006
ここで、曲面Ｘから曲面Ｙへの求めるべき写像を（５）式で定義する。写像の対応関係は曲面Ｘ上の点ｘが一意に曲面Ｙに対応づけることができればよいので、写像ｇはすべての点が一意に決定する全単射とする。よって同一パラメータ（ｕ，ｖ）を有する点ｘ（ｕ，ｖ）、ｙ（ｕ，ｖ）は，（６）式を満たすことになる。
【００３２】
【数５】
2007188329.tif 000007
このことは、すなわち、写像ｇについて考えることは、曲面Ｙそのものを考えることと同じである。一つの対象物が対象画像と参照画像に射影されて映っていると考えると、射影された対象画像上の像と参照画像上の像は、画素値が同じであると考えることができる。従って、対応づけるべき同一点の画素値が等しいと仮定すると、写像ｇに関して（７）式の条件式が導出される。かかる（７）式をパラメータｕ＝（ｕ，ｖ）を用いて表示すると、（８）式のようになる。
【００３３】
【数６】
2007188329.tif 000008
従って、画像マッチング問題は、（８）式の条件式を満たす写像ｇを求めるという問題に帰着することになる。
【００３４】
次に、上記画像マッチング問題の解法処理について説明する。（８）式は、（９）式のように変形することができる。
【００３５】
【数７】
2007188329.tif 000009
しかしながら、（９）式は、画像に含まれるノイズなどにより厳密には成立しないため、かかる式から直接写像ｇを求めることは困難である。このため、本実施の形態では、最小二乗法を利用している。最小二乗法によれば、（１０）式に示すエネルギー関数を考え、かかるエネルギー関数の最小化を行うことによりノイズが含まれていても条件を満足する写像ｇを求めることができる。
【００３６】
【数８】
2007188329.tif 000010
大数の法則が成立する程度の標本数があれば、最小二乗法により解が安定に求まると考えられるが、本実施の形態では、標本数が一つであるため、やはりノイズの影響を受けてしまう。解を安定して求めるためには、ベイズの定理により、近傍では滑らかに変化するという事前知識により事後分布の精度を高めるアプローチ、例えば、写像ｇに対して構造を規定して力学的なアプローチを考える。
【００３７】
写像ｇは全単射であるので、写像ｇについて考えることは、曲面Ｙそのものを考えることと同義である。力学的構造としては、例えば、剛体のような構造が考えられるが、画像の場合には変形等が考えられるため、剛体の構造は適切でない。そこで、本実施の形態では、曲面Ｙに対して、画像の変形等に対して柔軟に対応可能な弾性体の構造を導入する。
【００３８】
弾性体のエネルギーは、（１１）式で示される。
【００３９】
【数９】
2007188329.tif 000011
また、（１０）式で示される画像エネルギーを曲面Ｙ上の点を用いて表すと、（１２）式で表される。
【００４０】
【数１０】
2007188329.tif 000012
ポテンシャルエネルギーは、（１２）式の画像エネルギーと（１１）式の弾性エネルギーの線形結合により（１３）式のように定義される。
【００４１】
【数１１】
2007188329.tif 000013
このポテンシャルエネルギーを最小にすることは、画像マッチングの条件式である（７）式をできるだけ満たすように、かつ弾性体としての安定形態が取れるような曲面Ｙを探索することを意味する。このことをエネルギー最小化問題として定式化すると（１４）式で表される。
【００４２】
【数１２】
2007188329.tif 000014
従って、条件式（７）式を満たす写像ｇを求める画像マッチング問題は、エネルギー最小化問題の（１４）式を解法することに帰着する。
【００４３】
また（１４）式より最適性の必要条件は、（１５）式で求めることができる。
【００４４】
【数１３】
2007188329.tif 000015
従来の手法では、この最適性必要条件の（１５）式を満たすような曲面Ｙ上の点ｙを求めることにより解法している。これらは静的なエネルギーの最小化であり、静的な釣り合いの位置を求める静的な最適化である。図４は、ポテンシャルエネルギーの最小化の点を探索する場合の静的な最適点の説明をするための模式図である。図４に示すように、静的な最適化では初期値に最も近い位置で質点のポテンシャルエネルギーの局所最適点（釣り合い位置）が求まる。従って、この局所最適点が大域的な最適点になるかどうかは初期値に依存することになる。言い換えれば、初期値の選定が好ましくない場合には、大域的に最適な点が得られない場合がある。
【００４５】
このため、大域的な最適点を求めて（１４）式を解法するため、本実施の形態では、対象画像および参照画像上の曲面を格子状に分割して分割した各格子点を質点と考え、参照画像上の曲面Ｙの格子点（質点）における運動エネルギーを導入し、大域的な最適化を行っている。すなわち、曲面Ｙに対して時間軸を導入し、曲面Ｙが時間に応じて変形できるような構造とし、曲面に対して運動の構造を適用する。図５は、ポテンシャルエネルギーの最小化の点を探索する場合の大域的な最適点の説明をするための模式図である。図５に示すように、運動エネルギーの導入により、最初のポテンシャルエネルギーが運動エネルギーに変換されるため、格子点（質点）が局所最適点からさらに広い範囲で（１５）式を満たす大域的な最適点を探索することができる。これにより、静的な最適化よりもより広い範囲を探索できるので、初期値への依存性が低くなり、ノイズ等の影響に対してよりロバストになる可能性がある。
【００４６】
次に、（１５）式の大域的最適化のための運動方程式の導入について説明する。まず、運動方程式に導入する時間τを用いて曲面Ｙ上の点ｙを（１６）式のように拡張する。（１６）式で表した曲面Ｙ上の点ｙの全微分は（１７）式となるが、パラメータ（ｕ，ｖ）Tは時間τに対して独立であると仮定すると、（１８）式が成立するので、結局（１７）式は（１９）式で表現される。また、（１９）式から曲面Ｙ上の点ｙの２階微分は（２０）式で表される。（１９）式で示される曲面Ｙ上の点ｙの全微分と（２０）式で示される２階微分をそれぞれ（２１）式で表現する。
【００４７】
【数１４】
2007188329.tif 000016
（１９）式は、曲面Ｙ上の点ｙの速度を表し、（２０）式は、曲面Ｙ上の点ｙの加速度を表している。従って、曲面Ｙ上の点ｙの運動エネルギーは、（２２）式で定義することができる。
【００４８】
【数１５】
2007188329.tif 000017
このような運動エネルギーを導入したことは、曲面Ｙ上の点ｙを質点と考えて、曲面上の運動を定義したことになり、曲面Ｙ上の点ｙの運動という力学的な構造を規定したことを意味する。かかる曲面Ｙ上の点ｙの運動は運動方程式によって記述される。
【００４９】
次に、ラグランジュの方法に従って次のように曲面Ｙ上の点ｙの運動方程式を導出する。曲面Ｙは、エネルギー最小化の（１４）式に従うので、ポテンシャルエネルギーＵを最小にする方向に動く必要がある。ラグランジアンは（２３）式で定義される。
【００５０】
（２３）式は、位置のポテンシャルエネルギーである重力場における物体の自由落下などと同様に、物体は重力場を最小化するように落下することを意味している。ラグランジュ方程式は（２４）式で示されるので、（２３）および（２４）式から、曲面Ｙ上の点ｙの運動方程式は、（２５）式で表される。
【００５１】
【数１６】
2007188329.tif 000018
【数１７】
2007188329.tif 000019
（２５）式の運動方程式に従った運動は保存系であり、エネルギー保存則により解は収束しない。このため、系を非保存系にするために（２６）式で示される摩擦エネルギーを導入する。
【００５２】
【数１８】
2007188329.tif 000020
（２６）式で示される摩擦エネルギーを導入した場合のラグランジュ方程式は、（２７）式で表される。
【００５３】
【数１９】
2007188329.tif 000021
この（２７）式を解法することによって、摩擦エネルギーを導入した運動方程式が（２８）式で得られる。
【００５４】
【数２０】
2007188329.tif 000022
次に、動的な探索、（２８）式の運動方程式に従った解がエネルギー最小化の（１４）式に対する最適性必要条件である（１５）式を満足することについて説明する。
【００５５】
まず、探索が完了した場合を考える。この場合、点ｙは完全に停止するので、点ｙの速度（ｙの全微分）および加速度（ｙの２階微分）はともに０となるので、（２８）式の運動方程式により（２９）式となり、これにより（３０）式が成立する。
【００５６】
【数２１】
2007188329.tif 000023
（３０）式が成立することは、最適性必要条件の（１５）式を満足することを意味している。従って、探索が完了した場合、そのときの点ｙは最適点であることがわかる。
【００５７】
次に、曲面Ｙ上の点ｙが静止した場合を考える。この場合、点ｙの速度（ｙの全微分）は０となるが加速度（ｙの２階微分）は０にならないので、（２８）式の運動方程式により（３１）式が得られる。
【００５８】
【数２２】
2007188329.tif 000024
ここで、点ｙが最適点であると仮定すると、最適性条件の（１５）式により、点ｙの加速度（ｙの２階微分）は０となり、最適点の探索が完了する。一方、点ｙが最適点でない場合には、点ｙの加速度（ｙの２階微分）が０とならないため、最適点の探索は完了していないことになる。すなわち、点ｙが静止した位置が最適点の位置である場合に（３１）式により最適性必要条件の（１５）式が満足することを意味している。
【００５９】
次に、点ｙに摩擦エネルギーが作用している場合を考える。摩擦エネルギーの性質から（３２）式が導かれ、これは点ｙの運動が停止して時間の経過により最適点の探索が完了することを意味している。
【００６０】
【数２３】
2007188329.tif 000025
このように、最適点の探索が完了した場合、点ｙが静止した場合、点ｙに摩擦エネルギーが作用している場合の上記説明により、（２８）式の運動方程式に従った点ｙの運動が、（１５）式の最適性必要条件を満足することがわかる。
【００６１】
一方、（１５）式は１階微分の形式であるが、２階微分の最適性必要条件については、図４または５に示すポテンシャルエネルギーの山の頂点で探索が完了する場合があるため、満足されない場合がある。この場合には、ポテンシャルエネルギーの山の頂点では、最適点の探索が完了しないための条件として（３３）式を仮定することにより、２階の最適性必要条件を満足することができる。この（３３）式の仮定は、具体的には（３４）式の条件を探索アルゴリズムで実現することにより達成することができる。
【００６２】
以上により、（２８）式の運動方程式に従った最適点の探索は（１４）式で示される最適化問題に対する最適解を与えることがわかる。このような探索は、静的な探索に比べてロバストなものとなる。
【００６３】
ここで、（２８）式に示す運動方程式において、弾性エネルギーＥkの中に点ｙおよびパラメータ（ｕ，ｖ）に関する偏微分が含まれているため、有限差分法によって偏微分方程式を差分方程式によって置換する。弾性エネルギーＥkの右辺を有限差分法によって離散化すると、（３５）式のようになる。
【００６４】
【数２４】
2007188329.tif 000026
ここで、ｎは離散化されたパラメータ（ｕ，ｖ）の格子点、ｙnはｎに対応する離散化されたＹ上の点、ｙn+(-1,0)，ｙn+(1,0)，ｙn+(0,-1)，ｙn+(0,1)は点ｙnの偏微分に関連する上下左右の格子点を示している。
【００６５】
この（３５）式をｙnで偏微分すると（３６）式が得られ、かかる式は、弾性エネルギーによって発生する力である弾性エネルギー力を示している。
【００６６】
【数２５】
2007188329.tif 000027
この弾性エネルギー力を（３７）式で示し、離散パラメータΔｕ，Δｖをそれぞれ１とすると、弾性エネルギー力は（３８）式で表されることになる。
【００６７】
【数２６】
2007188329.tif 000028
ここで、画像エネルギーＥiは、パラメータ（ｕ，ｖ）に関して独立と考え、（２８）式の運動方程式もそのまま離散可能となり、（３９）式で表される。
【００６８】
【数２７】
2007188329.tif 000029
この（３９）式の右辺第１項は、点ｙnにおける画像エネルギーの勾配の逆向きの力である画像エネルギー力であり、（４０）式のようにおく。
【００６９】
【数２８】
2007188329.tif 000030
次に、この画像エネルギー力Ｆi（ｎ）、すなわち画像エネルギーの勾配の計算について説明する。本実施の形態では、画像エネルギーの勾配を、局所最適化の手法によって求めている。図６および図７は、画像エネルギー力の概念を示す説明図である。ここで、画像モデルＳcは連続画像モデルであり、実際にはサンプリングされた画像モデルＳpを利用しているため、局所最適化においてもサンプリングされた画像モデルに基づいておこなっている。
【００７０】
ｙn（τ）にもっとも近いサンプリング点を局所空間中心ｙcとして、ｙcを（４１）式で求める。
【００７１】
【数２９】
2007188329.tif 000031
そして、隣接空間Ｌを（４２）式で定義する。このとき、局所探索集合Ωは、（４３）式のように定義することができる。
【００７２】
【数３０】
2007188329.tif 000032
そして、局所最適化の後、その方向へのベクトルを求め、求めたベクトルを正規化し勾配の大きさを乗算すると、画像エネルギー力として（４４）〜（４６）式が得られる。本実施の形態では、画像エネルギー力として（４４）〜（４６）式を使用している。
【００７３】
【数３１】
2007188329.tif 000033
なお、実装上の扱いやすさ等を考慮してエネルギー関数を（４７）式のように定義して、（４７）式により導かれる画像エネルギー力を用いることもできる。
【００７４】
【数３２】
2007188329.tif 000034
以上説明した画像エネルギーＦi（ｎ）、弾性エネルギー力Ｆk（ｎ）、摩擦力を使用して、格子点ｙn（τ）に対する運動方程式を画像マッチング状態遷移モデルとして生成すると、運動方程式は、（４８）式で表される。また、ｙn（０）が参照画像に等しいとし、運動エネルギーは０、すなわちｙn（０）の速度が０として考えると、運動方程式の初期条件は、（４９）式となる。
【００７５】
【数３３】
2007188329.tif 000035
本実施の形態では、画像エネルギー力計算部１０３によって、（４４）〜（４６）式で示される画像エネルギー力Ｆiを計算し、弾性エネルギー力計算部１０４によって（３８）式で示される弾性エネルギー力Ｆkを計算し、摩擦力計算部１０５によって（４８）式の右辺第３項で示される摩擦力を計算する。そして、解法処理部１０１によって（４８）、（４９）式で示される運動方程式を数値解析処理によって解法している。
【００７６】
次に、この解法処理部１０１による運動方程式の解法について説明する。（４８）、（４９）式で示される運動方程式、すなわち常微分方程式は、一般的に解析的に解法することができないため、この運動方程式の系が収束するのに十分大きな時間Ｔを考え、数値解析によってｔ＝（０，Ｔ）区間を計算することによって運動方程式の収束状態を推定する。
【００７７】
本実施の形態では、常微分方程式は初期値が定まれば離散変数法によって一意に解が求まることを利用して、オイラー法による数値解析処理を行っている。
【００７８】
尚、本実施の形態ではオイラー法による数値解析処理を行っているが、これに限定されるものではない。離散変数法には、オイラー法以外にも、ルンゲクッタ法、ブリルシュ・ストア法、予測子・修正子法、隠的ルンゲクッタ法など種々の手法があり、このため、数値解析処理として、これらのいずれの手法を用いてもよい。以下はオイラー法を用いた数値解析処理を例にあげて説明する。
【００７９】
オイラー法は、一階の常微分方程式に対する数値解法であるため、（５０）式の変数変換を（４８）、（４９）式を運動方程式に施すことにより、（４８）、（４９）式を一階の常微分方程式に変換する。これにより、（５１）、（５２）式が得られる。
【００８０】
【数３４】
2007188329.tif 000036
【数３５】
2007188329.tif 000037
一方、常微分方程式（５３）式をオイラー法で数値解析処理を行って解法する場合には、（５４）式のようにｔ(m)からｔ(m +1)=ｔ(m)+ｈへ解を進展させる。
【００８１】
【数３６】
2007188329.tif 000038
【数３７】
2007188329.tif 000039
ここで、ｘ(m)はｍステップであることを示し、ｈはステップ幅である。オイラー法のスキームを（５１）、（５２）式に適用すると（５５）、（５６）式の更新式が得られる。
【００８２】
【数３８】
2007188329.tif 000040
本実施の形態では、解法処理部１０１によって、（５１）、（５２）式の更新式を繰り返して解法することによって運動方程式（４８）、（４９）式の系の収束状態を求めている。そして、マッピング処理部１０７では、この収束状態から、曲面Ｙをｙn＝Ψ1,n(T)のように求め、さらに、（６）式よりｙn＝ｇ（ｎ）の関係から、写像ｇを求めている。
【００８３】
次に、摩擦係数制御部１０２の詳細について説明する。摩擦係数制御部１０５は、摩擦力計算部１０５によって（４８）式の右辺第３項の摩擦力を計算する際に、摩擦係数μを可変制御して計算する。
【００８４】
摩擦係数μは探索の仕方を決定づけるパラメータであり、例えばμ→∞のときに、運動エネルギーはすべて摩擦エネルギーとして吸収されてしまうため、最適点の探索は、動的な探索とはならず、図４で示したような静的な探索にとどまる。一方、μ＝０のときには、図５に示すような動的な探索となり最適点の広範囲な探索を行うことができるが、エネルギー保存則により解が収束しないという問題がある。
【００８５】
例えば、図８に示すような対象画像と参照画像のテスト画像を使用して摩擦係数μの値によるポテンシャルエネルギーの収束を実験により求めた。図８は、対象画像の例とこの対象画像から画面全体が左から右に８ピクセル移動した参照画像の例を示す説明図である。この図８に示す対象画像と参照画像を用いて、摩擦係数μを一定として上記で説明した（５１）、（５２）式の更新式を繰り返して解法することによって運動方程式（４８）、（４９）式の系の収束状態を求め、この収束状態から写像ｇを求めることを試みた。
【００８６】
図９−１は、摩擦係数μ＝５における対象画像と参照画像の写像関係を示す説明図である。図９−２は、摩擦係数μ＝０の場合における対象画像と参照画像の写像関係を示す説明図である。図９−１および９−２において灰色線が写像関係を示し、白線が方向を示している。
【００８７】
図８では、参照画像が対象画像から画面全体が左から右に８ピクセル移動したものであるため、図９−１、９−２の画像中に現れる写像関係の線が直線状である場合に、画面全体の写像が正確と判断する。
【００８８】
図９−１に示すように、摩擦係数μが大きい場合には、画面全体の写像が正しく得られないことがわかる。一方、図９−２に示すように、摩擦係数μが小さい場合には、画面全体の方向性はほぼ正しいが、写像関係が正確に表れていないことがわかる。
【００８９】
これは次のように考えることができる。図１０は、摩擦係数の相違による画像エネルギー力と時間との関係を示す説明図である。図１１は、摩擦係数の相違による弾性エネルギー力と時間との関係を示す説明図である。図１２は、摩擦係数の相違によるポテンシャルと時間との関係を示す説明図である。
【００９０】
図１２に示すポテンシャルエネルギーの時間変化を見ると、μ＝０の場合はエネルギー保存則により、ポテンシャルエネルギーが収束していないことがわかる。一方、μ＝５の場合は、ポテンシャルエネルギーの変化が非常に遅いことがわかる。図１０に示す画像エネルギーの変化と図１１に示す弾性エネルギーの変換を参酌すると、このようなポテンシャルエネルギーの変化の原因は、画像エネルギーは使用されず、弾性エネルギーのみで最適点の探索が行われていると推測される。
【００９１】
従って、最適点の探索自体は、摩擦係数が小さい方がよいが、このままでは収束しないので、本実施の形態の摩擦係数制御部１０２では、摩擦係数を可変制御し、運動方程式における時間が早い時期は摩擦係数を小さくし、後半の収束させる時期に摩擦係数を大きくするようにしている。これにより、ロバストに低エネルギーを求めることができる。
【００９２】
具体的には、摩擦係数制御部１０２は、（５７）式に示すように、摩擦係数μ（τ）を、運動方程式における時間τの２段階の線形関数により、時間τが経過するに従って増加するように計算している。
【００９３】
【数３９】
2007188329.tif 000041
図１３は、摩擦係数μの時間τによる変化の状態を示すグラフである。図１４は、（５７）式により摩擦係数μを計算して、画像マッチング処理を行った場合におけるポテンシャルエネルギーの時間的変化の状態を示す説明図である。μ制御の線が、摩擦係数を（５７）式により計算した場合である。なお、図１４において、μ1＝０、μmax＝０．５、Ｔ＝１００としている。図１４に示すように、ポテンシャルエネルギーが低い状態で収束していることがわかる。
【００９４】
次に、以上のように構成された本実施の形態にかかる画像マッチング装置１００による画像マッチング処理について説明する。図１５は、画像マッチング装置１００による画像マッチング処理の手順を示すフローチャートである。具体的には、先に示した（５１）、（５２）式を数値解析により解法することにより画像マッチング処理を実現している。
【００９５】
まず、解法処理部１０１は、時間τ(0)＝０に設定し（ステップＳ１５０１）、初期値Ψ1,n（０）＝ｎ、Ψ2,n（０）＝０を設定する（ステップＳ１５０２）。これにより、（５２）式が実行される。
【００９６】
次に、画像エネルギー力計算部１０３によって、すべてのｎに対してｍステップにおける画像相関ポテンシャル力Ｆi(m)(n)を計算する（ステップＳ１５０３）。かかる画像エネルギー力Ｆi(m)(n)の計算処理については後述する。
【００９７】
次いで、弾性エネルギー力計算部１０４によって、すべてのｎに対してｍステップにおける弾性エネルギー力Ｆk(m) (n)を（３８）式により算出する（ステップＳ１５０４）。そして、摩擦力計算部１０５によって、すべてのｎに対してｍステップにおける摩擦力［−μ(τ(m))Ψ2,n(m)］を計算する（ステップＳ１５０５）。
【００９８】
次に、解法処理部１０１によって、（５５）式の更新式を、ステップＳ１５０３〜Ｓ１５０５で求めた画像エネルギー力Ｆi(m) (n)、弾性エネルギー力Ｆk(m) (n)、摩擦力［−μ(τ(m))Ψ2,n(m)］で更新する（ステップＳ１５０６）。
【００９９】
次に、解法処理部１０１によって、Ψ1,n(m)の値をフレームメモリ１０６に保存する（ステップＳ１５０７）。そして、解法処理部１０１によって、τ(m+1)=τ(m)+ｈと更新し（ステップＳ１５０８）、ｍを１だけ増加する（ステップＳ１５０９）。そして、τ(m+1)があらかじめ定められた時間Ｔを越えたか否かを判断し（ステップＳ１５１０）、越えていない場合には、上記ステップＳ１５０３からＳ１５０９を繰り返し実行する。一方、τ(m+1)がＴを越えた場合には、解法処理部１０１によって、すべてのｎに対して、ｙnにΨ1,n(T)を設定する（ステップＳ１５１１）。そして、マッピング処理部１０７により、対象画像と参照画像との対応関係、すなわち写像を求める（ステップＳ１５１２）。
【０１００】
次に、ステップＳ１５０３における画像エネルギー力の計算処理について説明する。図１６は、画像エネルギー力計算部１０３による画像エネルギー力の計算処理の手順を示すフローチャートである。
【０１０１】
まず、画像エネルギー力計算部１０３は、（４１）式により、ｙn（τ）にもっとも近いサンプリング点を局所空間中心ｙcとして計算する（ステップＳ１６０１）。そして次に、画像エネルギー力計算部１０３は、（４２）式で定義される隣接空間Ｌを設定し（ステップＳ１６０２）、局所探索集合Ωを（４３）式によって計算する（ステップＳ１６０３）。
【０１０２】
次に、画像エネルギー力計算部１０３は、（４６）式により局所最適化計算としてｙminを計算し（ステップＳ１６０４）、正規化のために、ステップＳ１６０４で求めた局所最適化のｙminを用いて（４５）式によりｄ＝ｙmin−Ψ1,n(m)を計算する（ステップＳ１６０５）。そして、画像エネルギー力計算部１０３は、ｙmin、ｄ等を用いて、（４４）式により、画像エネルギー力Ｆi(m)（ｎ）を計算する（ステップＳ１６０６）。
【０１０３】
次に、ステップＳ１５０５における摩擦力［−μ(τ(m))Ψ2,n(m)］の計算処理について説明する。図１７は、摩擦力の計算処理の手順を示すフローチャートである。まず、現在の時間τ(m)における摩擦係数μ（τ(m)）を（５７）式に従って求めるため、摩擦係数制御部１０２は、現在の時間τ(m)が０以上Ｔ／２以下否かを判断する（ステップＳ１７０１）。
【０１０４】
そして、現在の時間τ(m)が０以上Ｔ／２以下である場合には（ステップＳ１７０１：Ｙｅｓ）、摩擦係数制御部１０２は、（５７）式に従い摩擦係数μ（τ(m)）を（２μ1／Ｔ）τ(m)により計算する（ステップＳ１７０２）。
【０１０５】
一方、現在の時間τ(m)がＴ／２を経過している場合には（ステップＳ１７０１：Ｎｏ）、摩擦係数制御部１０２は、（５７）式に従い摩擦係数μ（τ(m)）を[２(μmax−μ1)／Ｔ］τ(m)＋２μ1−μmaxにより計算する（ステップＳ１７０３）。そして、摩擦力計算部１０５は、ステップＳ１７０２またはＳ１７０３で計算した摩擦係数を用いて摩擦力［−μ(τ(m))Ψ2,n(m)］を計算する（ステップＳ１７０４）。
【０１０６】
このようにして求めた画像エネルギー力Ｆi(m) (n)、弾性エネルギー力Ｆk(m) (n)、摩擦力［−μ(τ(m))Ψ2,n(m)］により、上述したように、ステップＳ１５０６で（５５）式の更新式が更新されることになる。
【０１０７】
図１８は、（５７）式により摩擦係数を可変に制御して場合における図８で示した対象画像と参照画像の写像関係を示す説明図である。図１８に示すように、画面全体に渡って正確に写像が求まっている。この場合、図１４と同様にポテンシャルエネルギーＵが低い状態で収束することができる。
【０１０８】
このように実施の形態１にかかる画像マッチング装置１００では、摩擦係数制御部１０２によって摩擦力計算の際の摩擦係数を、運動方程式における時間に応じて線形関数に従って増加するように可変制御しているので、大域的な最適点の探索と収束性を両立させて、より高精度な画像マッチングを行うことができる。
【０１０９】
なお、実施の形態１では、摩擦係数を、運動方程式における時間に応じて線形関数に従って増加するように可変制御しているが、これに限定されるものではなく、摩擦係数を、運動方程式における時間に応じて段階的に増加するように可変制御してもよい。例えば、摩擦係数μを（５８）式に従って段階的に増加するように制御することができる。図１９は、摩擦係数を運動方程式における時間に応じて（５８）式に従い段階的に増加する場合における摩擦係数μの時間的変化を示す説明図である。
【０１１０】
【数４０】
2007188329.tif 000042
また、摩擦係数制御部１０２によって、摩擦係数を、運動方程式における時間に応じて指数関数に従って増加するように可変制御してもよい。例えば、摩擦係数μを（５９）式で示す指数関数に従って増加するように制御することができる。図２０は、摩擦係数を、運動方程式における時間に応じて（５９）式の指数関数に従って増加する場合における摩擦係数μの時間的変化を示す説明図である。
【０１１１】
【数４１】
2007188329.tif 000043
【０１１２】
（実施の形態２）
実施の形態２にかかる画像マッチング装置１００では、摩擦係数μの初期値として負値を設定可能で、摩擦係数をかかる初期値から運動方程式の時間に応じて増加するように可変制御するものである。
【０１１３】
実施の形態２にかかる画像マッチング装置１００の構成は実施の形態１と同様である。また、画像マッチングの全体処理、画像エネルギー力計算処理は、弾性エネルギー力計算処理は、実施の形態１と同様である。実施の形態２では、摩擦力計算処理が実施の形態１と異なっている。
【０１１４】
通常、摩擦係数μは正値をとるが、実施の形態２にかかる摩擦係数制御部１０２では、摩擦係数μの初期値として負値をとり、かかる初期値から運動方程式の時間に応じて増加するように可変制御している。これにより最適点のさらなる広範囲な探索が可能となる。ただし、この場合には、摩擦係数が負値に固定されていると、加速して発散してしまうため、可変制御の条件が重要となる。具体的には、摩擦係数μを（６０）式に従って、負値から時間τに応じて線形関数に従って増加させている。摩擦係数μが負値の場合には、加速エネルギーととらえることができる。
【０１１５】
【数４２】
2007188329.tif 000044
図２１は、実施の形態２における摩擦係数μの時間τによる変化の状態を示すグラフである。ここで、μ0は負値となっており、このため、時間τ＝０における摩擦係数μは負値の初期値をとることになる。
【０１１６】
図２２は、実施の形態２における摩擦力計算処理の手順を示すフローチャートである。まず、現在の時間τ(m)における摩擦係数μ（τ(m)）を（５７）式に従って求めるため、摩擦係数制御部１０２は、現在の時間τ(m)が０以上Ｔ／２以下否かを判断する（ステップＳ２２０１）。
【０１１７】
そして、現在の時間τ(m)が０以上Ｔ／２以下である場合には（ステップＳ２２０１：Ｙｅｓ）、摩擦係数制御部１０２は、（６０）式に従い摩擦係数μ（τ(m)）を（２μ1／Ｔ）τ(m)＋μ0により計算する（ステップＳ２２０２）。
【０１１８】
一方、現在の時間τ(m)がＴ／２を経過している場合には（ステップＳ２２０１：Ｎｏ）、摩擦係数制御部１０２は、（６０）式に従い摩擦係数μ（τ(m)）を[２(μmax−μ1)／Ｔ］τ(m)＋２μ1−μmaxにより計算する（ステップＳ２２０３）。そして、摩擦力計算部１０５は、ステップＳ２２０２またはＳ２２０３で計算した摩擦係数を用いて摩擦力［−μ(τ(m))Ψ2,n(m)］を計算する（ステップＳ２２０４）。
【０１１９】
図２３は、対象画像の例とこの対象画像から画面全体が左から右に１２ピクセル移動した参照画像の例を示す説明図である。この図２３に示す対象画像と参照画像を用いて、摩擦係数μに負値の初期値を用いない場合と負値を用いた場合に（５１）、（５２）式の更新式を繰り返して解法することによって運動方程式（４８）、（４９）式の系の収束状態を求め、この収束状態から写像ｇを求めることを試みた。
【０１２０】
図２４−１は、摩擦係数の初期値に負値でない値（０以上の値）を用いて時間に応じて増加させる制御を行った場合における対象画像と参照画像の写像関係を示す説明図である。図２４−２は、摩擦係数の初期値に負値を用いて（６０）式に従って時間に応じて増加させる制御を行った場合における対象画像と参照画像の写像関係を示す説明図である。図２４−１および２４−２においては、図９−１，９−２と同様に、灰色線が写像関係を示し、白線が方向を示している。
【０１２１】
図２４−１に示すように、摩擦係数の初期値に負値でない値を用いて時間に応じて増加させる制御を行った場合には、画面半分の写像が正確に求めることができないことがわかる。一方、図２４−１に示すように、摩擦係数の初期値に負値を用いて時間に応じて増加させる制御を行った場合には、画面全体の写像が正確に求められていることがわかる。
【０１２２】
図２５は、（５７）式により摩擦係数μを計算して、画像マッチング処理を行った場合におけるポテンシャルエネルギーの時間的変化の状態を示す説明図である。
【０１２３】
図２５は、摩擦係数の初期値に負値でない値を用いて時間に応じて増加させる制御を行った場合と摩擦係数の初期値に負値を用いて時間に応じて増加させる制御を行った場合におけるポテンシャルエネルギーの時間的変化の状態を示す説明図である。図２５に示すように、初期値に負値を用いて摩擦係数の可変制御を行った場合の方が、より大域的な最適点に達しエネルギーが低いことがわかる。
【０１２４】
このように実施の形態２にかかる画像マッチング装置では、摩擦係数μの初期値として負値を設定可能で、摩擦係数をかかる初期値から運動方程式の時間に応じて増加するように可変制御しているので、より大域的な最適点の探索と収束性を両立させて、より高精度な画像マッチングを行うことができる。
【０１２５】
なお、実施の形態２では、摩擦係数の初期値として負値を用い、線形関数に従って増加させているが、これに限定されるものではなく、時間に応じて段階的に負値の初期値から増加させるように可変制御することもできる。例えば、摩擦係数μを（６１）式に従って負値の初期値から段階的に増加するように制御することができる。図２６は、摩擦係数を運動方程式における時間に応じて（６１）式に従い負値の初期値から段階的に増加する場合における摩擦係数μの時間的変化を示す説明図である。
【０１２６】
【数４３】
2007188329.tif 000045
【０１２７】
（実施の形態３）
実施の形態１および２にかかる画像マッチング装置では、摩擦係数を運動方程式における時間に応じて増加させるように可変制御していたが、この実施の形態３の画像マッチング装置は、ポテンシャルエネルギーの平衡状態に応じて摩擦係数の可変制御を行うものである。
【０１２８】
図２７は、実施の形態３にかかる画像マッチング装置の構成を示すブロック図である。本実施の形態にかかる画像マッチング装置２７００は、図２７に示すように、解法処理部１０１と、画像エネルギー力計算部１０３と、弾性エネルギー力計算部１０４と、摩擦力計算部１０５と、摩擦係数制御部２７０２と、マッピング処理部１０７と、平衡状態判定部２７０３と、フレームメモリ１０６とを主に備えている。
【０１２９】
ここで、解法処理部１０１、画像エネルギー力計算部１０３、弾性エネルギー力計算部１０４、摩擦力計算部１０５、マッピング処理部１０７については実施の形態１と同様である。
【０１３０】
平衡状態判定部２７０３は、ポテンシャルエネルギーが平衡状態にあるか否かを判定する処理部である。図１２より、摩擦係数μが０の場合、エネルギー保存則によりポテンシャルエネルギーの振動が止まらないが、ある一定の状態の周りを振動していることがわかる。これは、ポテンシャルエネルギーが平衡点周りを振動している運動ととらえることができる。そして、ポテンシャルエネルギーがかかる状態に移行している場合には、既に平衡状態に到達しているため、最適点の大域的な探索は完了していると判断することができる。このため、本実施の形態では、平衡状態判定部２７０３によって、ポテンシャルエネルギーが平衡状態であるか否かを判断し、平衡状態にある場合には、平衡状態への移行前よりも摩擦係数を大きく設定している。
【０１３１】
平衡状態に到達しているか否かは、例えばポテンシャルエネルギーの分布状態により判断することができる。図２８は、ポテンシャルエネルギーが平衡状態に達している場合の、ポテンシャルエネルギーの分布状態を示す説明図である。図２９は、ポテンシャルエネルギーが平衡状態に達していない場合の、ポテンシャルエネルギーの分布状態を示す説明図である。
【０１３２】
平衡状態に達しておらず状態が分散しているときには、図２９のような分布を示している。一方、平衡状態に達している場合には、状態がある程度まとまりを示すので、図２８のような分布を示す。これらの違いは分布の分散が大きいか、小さいかである。そこで、平衡状態判定部２７０３は、ポテンシャルエネルギーの分布の分散に大きな変化がなくなったら平衡状態に達したと判断することとしている。
【０１３３】
また、平衡状態に達している場合には、ある点ｙが新しい移動領域に移動していくことがない状態であり、同じ領域内を移動している状態である。このため、平衡状態判定部２７０３は、点ｙの移動領域の変化率を求め、変化率が小さくなった場合には平衡状態と判断してもよい。具体的には、点ｙの座標の分布範囲の分散が一定範囲内にある場合には、点ｙの移動領域の変化率が小さくなったとして、平衡状態に達したと判断すればよい。
【０１３４】
図３０は、実施の形態３における摩擦力計算処理の手順を示すフローチャートである。まず、平衡状態判定部２７０３は、すべてのｎに対してポテンシャルエネルギーＵ（ｎ）の標本分散を求め（ステップＳ３００１）、求めた標本分散から母分散を求める（ステップＳ３００２）。ここで、ポテンシャルエネルギーＵ（ｎ）の標本分散は、標本平均を（６２）式により求め、標本平均を使用する標本分散を（６３）式により計算することにより求める。
【０１３５】
【数４４】
2007188329.tif 000046
次に、平衡状態判定部２７０３は、ｍステップについて推定値がｍ−１ステップについての推定値より所定値以上変化しているか否かを調べる（ステップＳ３００３）。そして、所定値以上変化している場合には、まだ、平衡状態に達していないと判断して、摩擦係数μ（τ(m)）にμ0を設定する（ステップＳ３００４）。一方、ステップｍについて推定値がステップｍ−１についての推定値より所定値以上変化していない場合には、平衡状態に達したと判断して、摩擦係数μ（τ(m)）にμ1を設定する（ステップＳ３００５）。ここで、μ0＜μ1である。このため、平衡状態に達した場合には摩擦係数が大きく設定されることになる。
【０１３６】
そして、摩擦力計算部１０５は、ステップＳ３００４またはＳ３００５で設定した摩擦係数を用いて摩擦力［−μ(τ(m))Ψ2,n(m)］を計算する（ステップＳ３００６）。
【０１３７】
このように実施の形態３にかかる画像マッチング装置では、ポテンシャルエネルギーの平衡状態に応じて摩擦係数の可変制御を行っているので、画像により適応した画像マッチング処理を実現することができる。
【０１３８】
実施の形態１〜３の画像マッチング装置は、ＣＰＵなどの制御装置と、ＲＯＭ（Read Only Memory）やＲＡＭなどの記憶装置と、ＨＤＤ、ＣＤドライブ装置などの外部記憶装置と、ディスプレイ装置などの表示装置と、キーボードやマウスなどの入力装置を備えており、通常のコンピュータを利用したハードウェア構成となっている。
【０１３９】
実施の形態１〜３の画像マッチング装置で実行される画像マッチングプログラムは、インストール可能な形式又は実行可能な形式のファイルでＣＤ−ＲＯＭ、フレキシブルディスク（ＦＤ）、ＣＤ−Ｒ、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）等のコンピュータで読み取り可能な記録媒体に記録されて提供される。
【０１４０】
また、実施の形態１〜３の画像マッチング装置で実行される画像マッチングプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成しても良い。また、本実施形態の画像マッチング装置で実行される画像マッチングプログラムをインターネット等のネットワーク経由で提供または配布するように構成しても良い。
【０１４１】
また、本実施の形態の画像マッチングプログラムを、ＲＯＭ等に予め組み込んで提供するように構成してもよい。
【０１４２】
実施の形態１〜３の画像マッチング装置で実行される画像マッチングプログラムは、上述した各部（解法処理部、画像エネルギー力計算部、弾性エネルギー力計算部、摩擦力計算部、摩擦係数制御部、平衡状態判定部）を含むモジュール構成となっており、実際のハードウェアとしてはＣＰＵ（プロセッサ）が上記記憶媒体から画像マッチングプログラムを読み出して実行することにより上記各部が主記憶装置上にロードされ、上記各部（解法処理部、画像エネルギー力計算部、弾性エネルギー力計算部、摩擦力計算部、摩擦係数制御部、平衡状態判定部）が主記憶装置上に生成されるようになっている。
【０１４３】
なお、本発明は、上記実施の形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化することができる。また、上記実施の形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成することができる。例えば、実施の形態に示される全構成要素からいくつかの構成要素を削除してもよい。さらに、異なる実施の形態にわたる構成要素を適宜組み合わせても良い。
【図面の簡単な説明】
【０１４４】
【図１】実施の形態１にかかる画像マッチング装置の構成を示すブロック図である。
【図２】対象画像上の曲面から参照画像上の曲面への写像を示す説明図である。
【図３】対象画像上の曲面Ｘの点ｘから参照画像上の曲面Ｙの点ｙへの写像を示す説明図である。
【図４】ポテンシャルエネルギーの最小化の点を探索する場合の静的な最適点の説明をするための模式図である。
【図５】ポテンシャルエネルギーの最小化の点を探索する場合の大域的な最適点の説明をするための模式図である。
【図６】画像エネルギー力の概念を示す説明図である。
【図７】画像エネルギー力の概念を示す説明図である。
【図８】対象画像の例とこの対象画像から画面全体が左から右に８ピクセル移動した参照画像の例を示す説明図である。
【図９−１】摩擦係数μ＝５における対象画像と参照画像の写像関係を示す説明図である。
【図９−２】摩擦係数μ＝０の場合における対象画像と参照画像の写像関係を示す説明図である。
【図１０】摩擦係数の相違による画像エネルギー力と時間との関係を示す説明図である。
【図１１】摩擦係数の相違による弾性エネルギー力と時間との関係を示す説明図である。
【図１２】摩擦係数の相違によるポテンシャルと時間との関係を示す説明図である。
【図１３】摩擦係数μの時間τによる変化の状態を示すグラフである。
【図１４】（５７）式により摩擦係数μを計算して、画像マッチング処理を行った場合におけるポテンシャルエネルギーの時間的変化の状態を示す説明図である。
【図１５】画像マッチング装置１００による画像マッチング処理の手順を示すフローチャートである。
【図１６】画像エネルギー力計算部１０３による画像エネルギー力の計算処理の手順を示すフローチャートである。
【図１７】摩擦力の計算処理の手順を示すフローチャートである。
【図１８】（５７）式により摩擦係数を可変に制御して場合における図８で示した対象画像と参照画像の写像関係を示す説明図である。
【図１９】摩擦係数を運動方程式における時間に応じて（５８）式に従い段階的に増加する場合における摩擦係数μの時間的変化を示す説明図である。
【図２０】摩擦係数を、運動方程式における時間に応じて（５９）式の指数関数に従って増加する場合における摩擦係数μの時間的変化を示す説明図である。
【図２１】実施の形態２における摩擦係数μの時間τによる変化の状態を示すグラフである。
【図２２】実施の形態２における摩擦力計算処理の手順を示すフローチャートである。
【図２３】対象画像の例とこの対象画像から画面全体が左から右に１２ピクセル移動した参照画像の例を示す説明図である。
【図２４−１】摩擦係数の初期値に負値でない値を用いて時間に応じて増加させる制御を行った場合における対象画像と参照画像の写像関係を示す説明図である。
【図２４−２】摩擦係数の初期値に負値を用いて（６０）式に従って時間に応じて増加させる制御を行った場合における対象画像と参照画像の写像関係を示す説明図である。
【図２５】（５７）式により摩擦係数μを計算して、画像マッチング処理を行った場合におけるポテンシャルエネルギーの時間的変化の状態を示す説明図である。
【図２６】摩擦係数を運動方程式における時間に応じて（６１）式に従い負値の初期値から段階的に増加する場合における摩擦係数μの時間的変化を示す説明図である。
【図２７】実施の形態３にかかる画像マッチング装置の構成を示すブロック図である。
【図２８】ポテンシャルエネルギーが平衡状態に達している場合の、ポテンシャルエネルギーの分布状態を示す説明図である。
【図２９】ポテンシャルエネルギーが平衡状態に達していない場合の、ポテンシャルエネルギーの分布状態を示す説明図である。
【図３０】実施の形態３における摩擦力計算処理の手順を示すフローチャートである。
【符号の説明】
【０１４５】
１００，２７００  画像マッチング装置
１０１  解法処理部
１０２，２７０２  摩擦係数制御部
１０３  画像エネルギー力計算部
１０４  弾性エネルギー力計算部
１０５  摩擦力計算部
１０６  フレームメモリ
１０７  マッピング処理部
２７０３  平衡状態判定部
【図１】
2007188329.tif 000047
【図２】
2007188329.tif 000048
【図３】
2007188329.tif 000049
【図４】
2007188329.tif 000050
【図５】
2007188329.tif 000051
【図６】
2007188329.tif 000052
【図７】
2007188329.tif 000053
【図８】
2007188329.tif 000054
【図９−１】
2007188329.tif 000055
【図９−２】
2007188329.tif 000056
【図１０】
2007188329.tif 000057
【図１１】
2007188329.tif 000058
【図１２】
2007188329.tif 000059
【図１３】
2007188329.tif 000060
【図１４】
2007188329.tif 000061
【図１５】
2007188329.tif 000062
【図１６】
2007188329.tif 000063
【図１７】
2007188329.tif 000064
【図１８】
2007188329.tif 000065
【図１９】
2007188329.tif 000066
【図２０】
2007188329.tif 000067
【図２１】
2007188329.tif 000068
【図２２】
2007188329.tif 000069
【図２３】
2007188329.tif 000070
【図２４−１】
2007188329.tif 000071
【図２４−２】
2007188329.tif 000072
【図２５】
2007188329.tif 000073
【図２６】
2007188329.tif 000074
【図２７】
2007188329.tif 000075
【図２８】
2007188329.tif 000076
【図２９】
2007188329.tif 000077
【図３０】
2007188329.tif 000078