(57)【特許請求の範囲】
【請求項１】
移動対象を含む画像を入力する画像入力手段と、
入力された前記画像を時系列に蓄積する画像蓄積手段と、
蓄積された前記画像を読み出して、前記画像毎に異なる空間解像度を有する階層画像を生成して階層画像蓄積手段に蓄積する階層画像生成手段と、
画像上の位置の輝度値を自己回帰モデルまたは多項式回帰モデルに適用し、時系列画像間の輝度変動量（＝（∂Ｅ（ｘ）／∂ｔ）＋∇Ｅ（ｘ）・Δｔ・ｖ（ｘ，ｔ））（但し、ｖは位置ｘにおける時刻ｔでの速度ベクトル、Ｅは位置ｘの輝度値、Δｔは離散化した時間幅）が当該自己回帰モデルまたは当該多項式回帰モデルに従うとした輝度変動方程式Ｆ１を導出し、当該式Ｆ１に対して前記画像上の位置の輝度の発散度および渦度を拘束条件に用いた式Ｆ２を定義し、前記式Ｆ１と前記式Ｆ２との和を最小化することを目的とした目的関数を導出して目的関数格納手段に格納する目的関数導出手段と、
格納された前記目的関数を読み出して、連続する２枚の前記画像に対する前記階層画像のうち空間解像度の低い２枚の前記階層画像の輝度値を前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の速度ベクトルを仮推定して速度ベクトル格納手段に格納する速度ベクトル仮推定手段と、
格納された前記目的関数と前記速度ベクトルとを読み出して、当該速度ベクトルと次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の前記速度ベクトルを推定して速度ベクトル格納手段に格納する速度ベクトル推定手段と、
を有することを特徴とする移動ベクトル検出装置。
【請求項２】
前記速度ベクトル推定手段は、当該速度ベクトル推定手段により格納された前記速度ベクトルと前記速度ベクトル推定手段で用いた前記２枚の階層画像の次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入して前記速度ベクトルを推定することを繰り返し、当該速度ベクトルを前記速度ベクトル格納手段に格納することを特徴とする請求項１に記載の移動ベクトル検出装置。
【請求項３】
前記速度ベクトル仮推定手段および前記速度ベクトル推定手段は、前記目的関数を演算する際に用いる非線形ロバスト関数として、ローレンツ関数を用いることを特徴とする請求項１又は２に記載の移動ベクトル検出装置。
【請求項４】
移動対象を含む画像を画像入力手段により入力する第１ステップと、
入力された前記画像を画像蓄積手段により時系列に蓄積する第２ステップと、
蓄積された前記画像を読み出して、前記画像毎に異なる空間解像度を有する階層画像を生成して階層画像生成手段により階層画像蓄積手段に蓄積する第３ステップと、
画像上の位置の輝度値を自己回帰モデルまたは多項式回帰モデルに適用し、時系列画像間の輝度変動量（＝（∂Ｅ（ｘ）／∂ｔ）＋∇Ｅ（ｘ）・Δｔ・ｖ（ｘ，ｔ））（但し、ｖは位置ｘにおける時刻ｔでの速度ベクトル、Ｅは位置ｘの輝度値、Δｔは離散化した時間幅）が当該自己回帰モデルまたは当該多項式回帰モデルに従うとした輝度変動方程式Ｆ１を導出し、当該式Ｆ１に対して前記画像上の位置の輝度の発散度および渦度を拘束条件に用いた式Ｆ２を定義し、前記式Ｆ１と前記式Ｆ２との和を最小化することを目的とした目的関数を導出して目的関数導出手段により目的関数格納手段に格納する第４ステップと、
格納された前記目的関数を読み出して、連続する２枚の前記画像に対する前記階層画像のうち空間解像度の低い２枚の前記階層画像の輝度値を前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の速度ベクトルを仮推定して速度ベクトル仮推定手段により速度ベクトル格納手段に格納する第５ステップと、
格納された前記目的関数と前記速度ベクトルとを読み出して、当該速度ベクトルと次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の前記速度ベクトルを推定して速度ベクトル推定手段により速度ベクトル格納手段に格納する第６ステップと、
を有することを特徴とする移動ベクトル検出方法。
【請求項５】
前記第６ステップは、前記速度ベクトル推定手段により格納された前記速度ベクトルと前記第６ステップで用いた前記２枚の階層画像の次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入して前記速度ベクトルを推定することを繰り返し、当該速度ベクトルを前記速度ベクトル推定手段により前記速度ベクトル格納手段に格納することを特徴とする請求項４に記載の移動ベクトル検出方法。
【請求項６】
前記第５ステップおよび前記第６ステップは、前記目的関数を演算する際に用いる非線形ロバスト関数として、ローレンツ関数を用いることを特徴とする請求項４又は５に記載の移動ベクトル検出方法。
【請求項７】
請求項４乃至６のいずれか１項に記載の移動ベクトル検出方法における各ステップをコンピュータによって実行させることを特徴とする移動ベクトル検出プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、輝度変動が著しい連続する時系列画像から移動対象の移動ベクトルを安定かつロバストに推定する技術に関する。
【背景技術】
【０００２】
現在、時系列画像から物体の動き（オプティカルフロー）を推定する研究は、コンピュータビジョン分野において重要な位置を占めており、日々盛んに研究が進められている。
【０００３】
動きを検出する最も基本的な方法は、連続する２枚の画像間における高い類似性を利用する相互相関法（ＣＣ法）であり、ＭＰＥＧに代表される符号化の分野において良く知られている。具体的には、２枚の時系列画像において輝度変化がないことを前提とし、画像間の類似性を相互相関関数などにより評価して類似する領域を画像間の移動ベクトルとみなす方法である。推定精度を向上させるために、画像を複数のブロックに分割してサブブロック単位で移動ベクトルを求める方法も用いられている。
【０００４】
一方、コンピュータビジョン分野では、オプティカルフロー法を利用して動きを検出する方法が知られている。このオプティカルフロー法は、動画像から各点の動きを求める手法の一つで、画像輝度の空間的変化と時間的変化を調べ、各点がどちらへどんな速さで移動したかを求める方法である。具体的には、連続する２枚以上の画像から正則化と呼ばれる枠組みで対象の動き速度を推定する方法であり、２枚の時系列画像に基づいて目的関数として定義された式から変分法（もしくは、オイラー・ラグランジュ法）に従って非線形連立方程式を得る。この非線形連立方程式を計算する動作を緩和法により反復誤差が小さくなるまで繰り返す。このようにして収束した解が、２枚の画像間における移動対象の速度ベクトルとなる。なお、画像の輝度値，強度値または濃淡値の意味するところは全て等価であり、計算する際に使用される画像データは画像の濃淡値である。
【０００５】
従来のオプティカルフロー法では、画像の一定領域内における流量の出入が保存されること、つまり、輝度変動は一定であることを前提に式（１）を導き、この式（１）を計算することにより速度ベクトルが求められていた。
【数１】
0004538426.tif 000002
【０００６】
式（１）中、Ｅは２次元画像（ｘ，ｙ）の時間ｔにおける輝度（濃淡）を表すものである。ｖは速度ベクトルである。左辺第一項は、微小時間における輝度変化の関数、第二項は、輝度の発散量の関数を示すものであり、右辺をゼロとすることで輝度変動を一定とすることを表している。
【０００７】
次に、この式（１）を微分演算の規則により展開することで、式（２）を得ることができる。
【数２】
0004538426.tif 000003
【０００８】
さらに、画像中の流動体の移動変位ベクトルｄは、離散時間幅Δｔと速度ベクトルｖにより、式（３）で定義することができる。
【数３】
0004538426.tif 000004
【０００９】
式（２）および式（３）を用いて、時間に関する積分計算を行うことにより、式（４）および式（５）を容易に得ることができる。
【数４】
0004538426.tif 000005
【数５】
0004538426.tif 000006
【００１０】
式（４）は、時間に関する一次微分方程式であり、この式（４）を解析的に解くことにより一般解である式（５）を得ることができる。なお、式（４）の右側に記述された式（∀ｕ∈（ｔ，ｔ＋Δｔ））は、離散時間Δｔの間にｘ軸方向の速度ｕが全て含まれていることを表している。
【００１１】
この式（５）に対し、テイラー展開を通じて輝度に関する輝度変動方程式が式（６）のように導かれる。この式（６）を最小化問題として収束計算することにより、速度ベクトルを推定する。
【数６】
0004538426.tif 000007
【００１２】
しかしながら、式（６）のみによる計算では解を収束させることができないので、式（６）に対して一定の条件を導入して計算を行う。ここで、流体力学的な拘束条件として、発散度および渦度を考慮した拘束条件を式（７）で定義する。
【数７】
0004538426.tif 000008
【００１３】
右辺の第一項および第二項が発散度方程式であり、第一項のξおよび第二項は収束計算を行うための発散度に関する誤差を考慮したものである。右辺の第三項と第四項が渦度方程式であり、第三項のζおよび第四項は収束計算上を行うための渦度に関する誤差を考慮したものである。
【００１４】
式（６）および式（７）から目的関数を式（８）を導出する。αは、正則化パラメータである。なお、ｄ＝（ｕ，ｖ）の式における左辺のｄは移動速度ベクトルであり、右辺の（ｕ，ｖ）は単位時間あたりのｘ軸およびｙ軸の速度ベクトルである。
【数８】
0004538426.tif 000009
【００１５】
従い、速度ベクトル（ｕ，ｖ）は、式（８）の発散度の誤差ξ，および渦度の誤差ζが小さくなるまで収束計算を行うことで計算することができる。
【００１６】
具体的には、式（８）は汎関数であるので、変分法の枠組みによりオイラー・ラグランジュ方程式を介して３つの変数（Δｄ，ξ，ζ）についての非線形連立方程式を容易に得ることができる。ここで、離散化された計算格子として、２次元画像の画素の位置を（ｉ，ｊ）で定義した場合、発散度および渦度は、有限差分法により離散化された式（９）で記述することができる。
【数９】
0004538426.tif 000010
【００１７】
また、この非線形連立方程式を計算する際に用いる非線形ロバスト関数として、Ｌｅｃｌｅｒｃ関数を用いる。これは、急峻な輝度変化による不安定な計算を抑制するためである。σは分散を表している。
【数１０】
0004538426.tif 000011
【００１８】
故に、速度ベクトル（ｕ，ｖ）は、上述した非線形連立方程式，式（９），式（１０），および２次元画像の位置（ｉ，ｊ）における画像の濃淡値を利用することにより、誤差ξ，ζが小さくなるまで収束計算を行うことで計算することができる。計算効率上、３つの変数（Δｄ，ξ，ζ）の初期値をいずれもゼロと設定し、これら３つのうち１つを変数、残り２つを固定値とみなして繰り返し収束計算をすることで、速度ベクトル（ｕ，ｖ）を推定する。
【００１９】
なお、ここで用いている定式化は、非特許文献１乃至３に従っている。
【非特許文献１】J.L. Barron、他２名、「Performance of Optical Flow Techniques」、IJCV、1994、vol.12、no.1、p.43-77
【非特許文献２】Richard P. Wildes、他３名、「Recovering Estimates of Fluid Flow from Image Sequence Data」、CVIU、2000、vol.80、p.246-266
【非特許文献３】Thomas Corpetti、他２名、「Dense Estimation of Fluid Flows」、IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE、2002年３月、vol.24、no.3、p.365-381
【発明の開示】
【発明が解決しようとする課題】
【００２０】
オプティカルフロー法を用いて動きを検出する場合、上記にて説明したように、一般的には輝度変動を一定とする前提条件を考慮した場合が多い。つまり、突発的に発生する輝度変化を考慮することなく、画像の動きは常に滑らかであると仮定し、速度の１次微分や２次微分が近似的にゼロになるように設定して動きを推定している。
【００２１】
しかしながら、我々の周辺に存在する実際の環境では、オクルージョンや陰影による局所的な輝度変動や、気象レーダパターンおよび衛星画像パターンから観測される雲や雨などの急峻な発達，衰退，分裂，融合，消失など複雑な輝度変動が常に発生している。
【００２２】
つまり、従来のオプティカルフロー法では、渦度や発散度と呼ばれる数値流体力学の分野で定式化されている方程式を拘束条件として適用しているが、渦の発達や衰退といった流体本来の特徴を有する輝度変動を考慮していないので、気象レーダパターンから推測される雲などの動きは不自然な乱れとして検出され、動きを推定する精度は低いという問題がある。
【００２３】
また、衛星画像のような流体状パターンでは、気象レーダパターンと比べて流体のエッジ構造が明瞭ではなく、画像間の輝度変動は比較的小さいので、渦などが観測された場合であっても一定の動きに対しては一定の精度を得ることができるが、集中豪雨や雷雲などのように対流現象の周期が極めて短い場合には大きな輝度変化が発生するので、輝度変動が常に発生しないと仮定して流体状の動きを推定すると、上記と同様に不自然なものとなり、結果として推定精度が低いという問題がある。
【００２４】
本発明は、上記を鑑みてなされたものであり、輝度変動を生ずる映像から移動対象の動き（速度）ベクトルを推定する精度の安定化と高精度化を行うことを課題とする。
【課題を解決するための手段】
【００２５】
第１の本発明に係る移動ベクトル検出装置は、移動対象を含む画像を入力する画像入力手段と、入力された前記画像を時系列に蓄積する画像蓄積手段と、蓄積された前記画像を読み出して、前記画像毎に異なる空間解像度を有する階層画像を生成して階層画像蓄積手段に蓄積する階層画像生成手段と、画像上の位置の輝度値を自己回帰モデルまたは多項式回帰モデルに適用し、時系列画像間の輝度変動量（＝（∂Ｅ（ｘ）／∂ｔ）＋∇Ｅ（ｘ）・Δｔ・ｖ（ｘ，ｔ））（但し、ｖは位置ｘにおける時刻ｔでの速度ベクトル、Ｅは位置ｘの輝度値、Δｔは離散化した時間幅）が当該自己回帰モデルまたは当該多項式回帰モデルに従うとした輝度変動方程式Ｆ１を導出し、当該式Ｆ１に対して前記画像上の位置の輝度の発散度および渦度を拘束条件に用いた式Ｆ２を定義し、前記式Ｆ１と前記式Ｆ２との和を最小化することを目的とした目的関数を導出して目的関数格納手段に格納する目的関数導出手段と、格納された前記目的関数を読み出して、連続する２枚の前記画像に対する前記階層画像のうち空間解像度の低い２枚の前記階層画像の輝度値を前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の速度ベクトルを仮推定して速度ベクトル格納手段に格納する速度ベクトル仮推定手段と、格納された前記目的関数と前記速度ベクトルとを読み出して、当該速度ベクトルと次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の前記速度ベクトルを推定して速度ベクトル格納手段に格納する速度ベクトル推定手段と、を有することを特徴とする。
【００２６】
本発明にあっては、輝度変動量を自己回帰モデルまたは多項式回帰モデルとする方程式から輝度変動方程式を導いて速度ベクトルを推定するので、実環境を映像でセンシングした時に観測される渦などの強い非線形な動きと大きな輝度変動を伴っている場合であっても、移動対象の動きを確実に追随することができ、その動き（速度）ベクトルの安定化かつ高精度化を実現することができる。
【００２７】
また、異なる空間解像度を有する複数の階層画像を連続する２枚の画像に対してそれぞれ生成し、生成された複数の階層画像のうち空間解像度の低い２枚の階層画像を用いて速度ベクトルを仮推定し、この速度ベクトルと次に空間解像度の高い２枚の階層画像とを用いて速度ベクトルを推定する。これにより、階層画像の面積および空間解像度を低くして計算量を減らすことができるので、速度ベクトルの収束計算の速度を向上させることができ、確実に収束させることができる。
【００２８】
上記移動ベクトル検出装置において、前記速度ベクトル推定手段は、当該速度ベクトル推定手段により格納された前記速度ベクトルと前記速度ベクトル推定手段で用いた前記２枚の階層画像の次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入して前記速度ベクトルを推定することを繰り返し、当該速度ベクトルを前記速度ベクトル格納手段に格納することを特徴とする。
【００２９】
本発明にあっては、速度ベクトル推定手段により推定された速度ベクトルと速度ベクトル推定手段で用いた２枚の階層画像の次に空間解像度の高い２枚の階層画像の輝度値とを用いて速度ベクトルを推定することを繰り返すことにより、画像間の移動対象の速度ベクトルを推定することができるので、上記と同様に、階層画像の面積および空間解像度を低くして計算量を減らすことができ、速度ベクトルの収束計算の速度を向上させることができる。
【００３０】
上記移動ベクトル検出装置において、前記速度ベクトル仮推定手段および前記速度ベクトル推定手段は、前記目的関数を演算する際に用いる非線形ロバスト関数として、ローレンツ関数を用いることを特徴とする。
【００３１】
本発明にあっては、ローレンツ関数を非線形ロバスト関数として適用することにより、速度ベクトルの収束計算の速度をより高速に行うことができる。
【００３２】
第２の本発明に係る移動ベクトル検出方法は、移動対象を含む画像を画像入力手段により入力する第１ステップと、入力された前記画像を画像蓄積手段により時系列に蓄積する第２ステップと、蓄積された前記画像を読み出して、前記画像毎に異なる空間解像度を有する階層画像を生成して階層画像生成手段により階層画像蓄積手段に蓄積する第３ステップと、画像上の位置の輝度値を自己回帰モデルまたは多項式回帰モデルに適用し、時系列画像間の輝度変動量（＝（∂Ｅ（ｘ）／∂ｔ）＋∇Ｅ（ｘ）・Δｔ・ｖ（ｘ，ｔ））（但し、ｖは位置ｘにおける時刻ｔでの速度ベクトル、Ｅは位置ｘの輝度値、Δｔは離散化した時間幅）が当該自己回帰モデルまたは当該多項式回帰モデルに従うとした輝度変動方程式Ｆ１を導出し、当該式Ｆ１に対して前記画像上の位置の輝度の発散度および渦度を拘束条件に用いた式Ｆ２を定義し、前記式Ｆ１と前記式Ｆ２との和を最小化することを目的とした目的関数を導出して目的関数導出手段により目的関数格納手段に格納する第４ステップと、格納された前記目的関数を読み出して、連続する２枚の前記画像に対する前記階層画像のうち空間解像度の低い２枚の前記階層画像の輝度値を前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の速度ベクトルを仮推定して速度ベクトル仮推定手段により速度ベクトル格納手段に格納する第５ステップと、格納された前記目的関数と前記速度ベクトルとを読み出して、当該速度ベクトルと次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入し、当該目的関数の前記発散度の誤差および前記渦度の誤差が収束するまで繰り返し演算することで前記移動対象の前記速度ベクトルを推定して速度ベクトル推定手段により速度ベクトル格納手段に格納する第６ステップと、を有することを特徴とする。
【００３３】
上記移動ベクトル検出方法において、前記第６ステップは、前記速度ベクトル推定手段により格納された前記速度ベクトルと前記第６ステップで用いた前記２枚の階層画像の次に空間解像度の高い２枚の前記階層画像の輝度値とを前記目的関数に代入して前記速度ベクトルを推定することを繰り返し、当該速度ベクトルを前記速度ベクトル推定手段により前記速度ベクトル格納手段に格納することを特徴とする。
【００３４】
上記移動ベクトル検出方法において、前記第５ステップおよび前記第６ステップは、前記目的関数を演算する際に用いる非線形ロバスト関数として、ローレンツ関数を用いることを特徴とする。
【００３５】
第３の本発明に係る移動ベクトル検出プログラムは、上記移動ベクトル検出方法に記載の移動ベクトル検出方法における各ステップをコンピュータによって実行させることを特徴とする。
【発明の効果】
【００３６】
本発明によれば、輝度変動を生ずる映像から移動対象の動き（速度）ベクトルを推定する精度の安定化と高精度化を行うことができる。
【発明を実施するための最良の形態】
【００３７】
以下、本発明の実施の形態について図面を用いて説明する。
【００３８】
図１は、本発明の実施の形態に係る移動ベクトル検出装置を示した構成図である。図１に示すように、移動ベクトル検出装置１０は、画像を入力する画像入力部１００、画像を時系列画像として蓄積する画像蓄積部１１０、画像毎に異なる空間解像度を有する複数の階層画像を生成する階層画像生成部１２０、目的関数を導出する目的関数導出部１３０、速度ベクトルを推定する速度ベクトル推定部１４０、および移動ベクトル表示部１５０を備えた構成である。
【００３９】
上記移動ベクトル検出装置１０は、コンピュータ本体を構成するＣＰＵ，メモリ，ハードディスクなどにより実現される。また、次に説明する各ステップ（Ｓ１〜Ｓ６）の動作処理はプログラムで実行され、画像蓄積部１１０，階層画像蓄積部，目的関数格納部，および速度ベクトル格納部は、メモリまたはハードディスクにより構成される。
【００４０】
上記にて説明した移動ベクトル検出装置１０の動作処理について説明する。
【００４１】
最初に、画像入力部１００には、移動対象が記録された複数の画像が入力される（Ｓ１）。この移動対象とは、例えば、気象レーダなどにより撮影された雲などの流動体を示すものである。
【００４２】
次に、画像蓄積部１１０には、画像入力部１００により入力された複数の画像が時系列に蓄積される（Ｓ２）。時系列とは、時間とともに変動する現象の観測値を時間とともに記録した並びを意味するものである。
【００４３】
続いて、階層画像生成部１２０は、画像蓄積部１１０により蓄積された複数の時系列画像のうち、連続する２枚の蓄積画像を画像蓄積部１１０から読み出して、その２枚の蓄積画像それぞれに対して、蓄積画像よりも空間解像度の低い（粗い）複数の階層画像を生成し、階層画像蓄積部に蓄積する（Ｓ３）。
【００４４】
この階層画像生成について、図２を用いて説明する。図２は、ある時間ｔにおける蓄積画像２００と時間ｔに連続する時間ｔ＋１における蓄積画像２１０とにおける階層画像を示した模式図である。階層画像生成部１２０は、蓄積画像２００を用いて、この蓄積画像２００よりも空間解像度の低い解像度を有する縮小された階層画像２０１を生成し、さらに、階層画像２０１よりも空間解像度の低い解像度を有する縮小された階層画像２０２を生成する。蓄積画像２１０についても同様に、階層画像生成部１２０は、この蓄積画像２１０よりも空間解像度の低い解像度を有する縮小された階層画像２１１を生成し、さらに、階層画像２１１よりも空間解像度の低い解像度を有する縮小された階層画像２１２を生成する。なお、１つの蓄積画像から３つ以上の階層画像を生成することも可能である。
【００４５】
そして、目的関数導出部１３０は、輝度変動量を自己回帰モデルまたは多項式回帰モデルとする方程式から輝度変動方程式を導き、この輝度変動方程式と発散度に関する発散度方程式と渦度に関する渦度方程式とから速度ベクトルと発散度の誤差と渦度の誤差に関する目的関数を導出し、目的関数格納部に格納する（Ｓ４）。目的関数の具体的な導出方法については後述する。なお、ここで説明したＳ４は、Ｓ３の後の動作に限られるものではなく、上記Ｓ１〜Ｓ３の間もしくはＳ１より前であってもよい。
【００４６】
その後、速度ベクトル推定部１４０は、目的関数格納部から目的関数を読み出して、階層画像生成部１２０により生成された複数の階層画像をその目的関数に代入することで速度ベクトルを推定し、速度ベクトル格納部に格納する（Ｓ５）。
【００４７】
具体的には、初期値をゼロとする速度ベクトルと複数の階層画像のうち最も空間解像度の低い階層画像２０２および階層画像２１２とを用いて速度ベクトルを仮推定する。これにより仮推定されたその速度ベクトルと次に空間解像度の高い階層画像２０１および階層画像２１１とを用いて第１の速度ベクトルを推定する。そして、この第１の速度ベクトルと階層画像を生成する元の画像である蓄積画像２００および蓄積画像２１０を用いて、第２の速度ベクトルを推定する。この第２の本速度ベクトルが、連続する画像間における移動対象の速度ベクトルとなる。
【００４８】
なお、１つの蓄積画像から３つ以上の階層画像が生成される場合であっても、最も空間解像度の低い２枚の階層画像から順番に、上記と同様に、階層画像を生成した元の２枚の画像（上記の場合、蓄積画像２００および蓄積画像２１０）を用いて速度ベクトルを推定するまで同様の計算を繰り返し行うことにより、連続する２枚の画像間における移動対象の速度ベクトルを推定することができる。なお、具体的な速度ベクトルの推定方法については後述する。
【００４９】
最後に、移動ベクトル表示部１５０は、速度ベクトル推定部１４０により推定された速度ベクトルを用いて、移動対象の移動ベクトルを矢印記号として画面に表示する（Ｓ６）。
【００５０】
次に、本実施の形態におけるオプティカルフロー法を用いた速度ベクトルの推定方法について説明する。
【００５１】
本来、背景技術にて説明した式（１）の右辺はゼロではなく、単位流量Ｑとして定義されている。つまり、従来の推定方法の場合、対流セルに見られるような発達と衰退を考慮しているが、一定領域内での流量の出入が保存されることを前提とする不自然なものであり、急峻に発達と衰退のそれぞれの量が異なる場合でのモデル化とはなっていないものであった。
【００５２】
ここで、流動体の急峻な発達と衰退は、時系列的に連続する画像上の移動対象の輝度変化に相当すること仮定することができる。つまり、従来の推定方法は、式（１）の右辺をゼロとする輝度変化一定モデルであったが、本実施の形態は、式（１）に代えて、自己回帰モデルを適用する式（１１）、または、多項式モデルを適用する式（１２）を利用して明瞭に画像間の輝度変動を表現することにより速度ベクトルを推定する。
【数１１】
0004538426.tif 000012
【数１２】
0004538426.tif 000013
【００５３】
αは自己回帰定数、βは多項式回帰定数、εは残差である。式（１１）に対し、テイラー展開を通じて輝度に関する輝度変動方程式が式（１３）のように導かれる。この式（１３）を最小化問題として解析することにより、速度ベクトルを推定する。
【数１３】
0004538426.tif 000014
【００５４】
式（１３）は、従来の式（６）に相当するものであり、この式（１３）と式（７）の拘束条件式から式（８）の目的関数を導出する。速度ベクトル（ｕ，ｖ）は、式（８）の発散度の誤差ξ，および渦度の誤差ζが小さくなるまで収束計算を行うことで計算することができる。式（８）を用いた具体的な計算方法については、背景技術にて説明した方法と同様であるので、ここでは説明を省略する。但し、従来において用いられていた非線形ロバスト関数としてのＬｅｃｌｅｒｃ関数の式（１０）を、ローレンツ関数の式（１４）に代えて収束計算を行う。
【数１４】
0004538426.tif 000015
【００５５】
次に、連続した２枚の時系列画像における画像パターンの変化について説明する。図３は、気象レーダにより撮影された国内上空の反射強度の時系列画像を示した図である。ある時間ｔにおける画像３００の画像パターンおよび時間ｔに連続する時間ｔ＋１における画像３１０の画像パターンは、ともに流体状の模様を有しており、時間ｔから時間ｔ＋１への時間経過にともない、輪郭および模様が不定形に著しく変化している。なお、画像３１０に示された矢印は、本実施の形態における移動ベクトル検出装置１０により、画像３００から画像３１０へ変化する際の移動ベクトルを推定したものである。
【００５６】
図４は、従来の速度ベクトルの推定方法を用いた場合（ａ）と本実施の形態の速度ベクトルの推定方法を用いた場合（ｂ）の流動体の移動ベクトルを示した比較図である。本来、両画像の中心付近に存在する流動体は、衰退するとともに右回転しているが、図４（ａ）に示す従来の推定方法を用いた場合には渦成分の移動ベクトルが検出されず、流動体の中心から外側に衰退する成分を移動ベクトルとして推定する不自然なものとなっている。一方、図４（ｂ）に示す本実施の形態の推定方法を用いた場合は、渦による右回転の回転ベクトルについても正しく検出されており、中心付近の流動体は右回転をしながら外側に衰退する移動ベクトルをより正確に表示している。
【００５７】
図５は、従来の速度ベクトルの推定方法を用いた場合（破線）と本実施の形態の速度ベクトルの推定方法を用いた場合（実線）の移動ベクトルの推定精度の評価結果を示した図である。図５（ａ）は、横軸を輝度変動量とし、縦軸を流動体の動き誤差とする測定結果を示したものであり、図５（ｂ）は、横軸を渦度および発散度とし、縦軸を流動体の動き誤差とする測定結果を示したものである。図５は、気象レーダにより観測された時系列に連続する人工衛星画像を用いて、従来の速度ベクトルの推定方法を用いた場合と本実施の形態の速度ベクトルの推定方法を用いた場合との輝度変動量もしくは渦度および発散度に対する移動ベクトルを測定し、その測定結果と人工衛星画像から得られる実際の移動ベクトルとの差異を動き誤差として表示したものである。図５（ａ）および図５（ｂ）から、従来の場合と比較して本実施の形態の場合の動き誤差は小さいと認識することができる。従い、この測定結果により、本実施の形態における速度ベクトルの推定方法における推定精度は従来の方法に比べてより高精度であり、本実施の形態における移動ベクトル検出装置１０は有効なものであることが示されている。
【００５８】
図６は、本実施の形態を用いた場合の血管内部に流れる血流の流れの移動ベクトルを示した図である。本実施の形態による移動ベクトル検出装置１０によれば、血管の内部に小さい血栓が存在する場合、その血栓の周辺には血流の乱れや渦度が発生していることを検出することができる。
【００５９】
本実施形態によれば、輝度変動量を自己回帰モデルまたは多項式回帰モデルとする方程式から輝度変動方程式を導いて速度ベクトルを推定するので、実環境を映像でセンシングした時に観測される渦などの強い非線形な動きと大きな輝度変動を伴っている場合であっても、移動対象の動きを確実に追随することができ、その動き（速度）ベクトルの安定化かつ高精度化を実現することができる。
【００６０】
これにより、例えば、血流の流れを移動ベクトルとして可視化し、渦度や発散度といった特異点を自動的に解析することで、血栓などによる血流の不規則な流れを容易に把握することができ、診断時における病巣の見逃しを防ぐことができる。医学の分野においては、超音波画像による血流の流れの可視化が行われているが、超音波画像から血栓周辺の不規則な血流の流れを瞬時に判断することは難しいため経験者による診察が中心であるが、実務経験が比較的浅い医者などであっても容易に病巣を把握することができる。その他、流体実験における粒子を流し込むことによる流れの可視化を通じて、動き解析を必要とする流体力学，医学，航空力学分野においても適用することができる。
【００６１】
また、異なる空間解像度を有する複数の階層画像を連続する２枚の画像に対してそれぞれ生成し、生成された複数の階層画像のうち空間解像度の低い２枚の階層画像を用いて速度ベクトルを仮推定し、この速度ベクトルを次に空間解像度の高い階層画像を用いて速度ベクトルを推定する。この速度ベクトルの推定を、階層画像を生成した元の画像である２枚の蓄積画像を用いるまで繰り返すことにより、連続する２枚の画像間における移動対象の速度ベクトルを推定する。これにより、階層画像の面積および空間解像度を低くして計算量を減らすことができるので、速度ベクトルの収束計算の速度を向上させることができ、確実に収束させることができる。例えば、非線形連立方程式を計算する際の演算量は、画像の面積にほぼ比例するので、階層画像の１辺が蓄積画像の半分である場合には、約４（２×２）倍で高速演算することが可能となる。
【００６２】
さらに、目的関数を演算する際に用いる非線形ロバスト関数として、ローレンツ関数を用いることにより、速度ベクトルの収束計算の速度をより高速に行うことができる。
【図面の簡単な説明】
【００６３】
【図１】本発明の実施の形態に係る移動ベクトル検出装置を示した構成図である。
【図２】ある時間ｔにおける蓄積画像と時間ｔに連続する時間ｔ＋１における蓄積画像とにおける階層画像を示した模式図である。
【図３】気象レーダにより撮影された国内上空の反射強度の時系列画像を示した図である。
【図４】従来の速度ベクトルの推定方法を用いた場合と本実施の形態の速度ベクトルの推定方法を用いた場合の流動体の移動ベクトルを示した比較図である。
【図５】従来の速度ベクトルの推定方法を用いた場合と本実施の形態の速度ベクトルの推定方法を用いた場合の移動ベクトルの推定精度の評価結果を示した図である。
【図６】本実施の形態を用いた場合の血管内部に流れる血流の流れの移動ベクトルを示した図である。
【符号の説明】
【００６４】
１０…移動ベクトル検出装置
１００…画像入力部
１１０…画像蓄積部
１２０…階層画像生成部
１３０…目的関数導出部
１４０…速度ベクトル推定部
１５０…表示ベクトル推定部
２００，２１０…蓄積画像
２０１，２０２，２１１，２１２…階層画像
３００，３１０…画像
Ｓ１〜Ｓ６…ステップ
【図１】
0004538426.tif 000016
【図２】
0004538426.tif 000017
【図４】
0004538426.tif 000018
【図５】
0004538426.tif 000019
【図６】
0004538426.tif 000020
【図３】
0004538426.tif 000021