(57)【要約】      （修正有）【課題】原画像情報上の指定された点の微小な移動量を計測することができる動体の移動量計測装置、方法及びプログラムを提供する。【解決手段】動体の移動量計測装置１０は、入力インタフェース部２０が、観測画像情報を時系列に複数フレームを取得する。システム状態量算出部が、取得されたフレームにおける観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出する。観測量算出部が、取得されたフレームにおける観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出する。推定部が、処理領域ｎのシステム状態量の最適値を、処理領域ｎのエッジ抽出画像情報として推定する。追跡波形処理部７０が、推定部により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、入力された着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測する。【選択図】図１
2015152975.tif 000002
【特許請求の範囲】
【請求項１】
モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、
前記観測画像情報を時系列に複数フレームを取得する画像取得手段と、
前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出するシステム状態量算出手段と、
前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出する観測量算出手段と、
前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記処理領域ｎのシステム状態量の最適値を、前記処理領域ｎのエッジ抽出画像情報として推定する推定手段と、
ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力する入力手段と、
前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測する追跡手段と、
を含む動体の移動量計測装置。
【請求項２】
前記処理領域ｎは、前記処理領域（ｎ−１）と一部重複するように設定される請求項１記載の動体の移動量計測装置。
【請求項３】
モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、
前記観測画像情報を時系列に複数フレームを取得する画像取得手段と、
前記画像取得手段により取得された連続するフレームの画像間の位置合わせを行い、前記連続するフレーム間における画像内の物体の位置変化量を推定する位置変化量推定手段と、
時刻ｎでの前記観測画像情報におけるシステム状態量を、前記位置変化量推定手段により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移を用いて算出するシステム状態量算出手段と、
時刻ｎでの前記観測画像情報の観測量を算出する観測量算出手段と、
時刻ｎでの前記観測画像情報における、前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記時刻ｎでのシステム状態量の最適値をエッジ抽出画像情報として推定する推定手段と、
ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力する入力手段と、
前記推定手段によって推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測する追跡手段と、
を含む動体の移動量計測装置。
【請求項４】
前記状態遷移を、時刻（ｎ−１）におけるシステム状態量と時刻ｎにおけるシステム状態量との相関を表す忘却係数を用いて表した請求項３記載の動体の移動量計測装置。
【請求項５】
前記観測対象を臓器とした請求項１〜請求項４の何れか１項記載の動体の移動量計測装置。
【請求項６】
前記観測対象を腸とした請求項５記載の動体の移動量計測装置。
【請求項７】
モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、かつ画像取得手段、システム状態量算出手段、観測量算出手段、推定手段、入力手段、及び追跡手段を含む動体の移動量計測装置における動体の移動量計測方法であって、
前記画像取得手段が、前記観測画像情報を時系列に複数フレームを取得すると、
システム状態量算出手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出するステップと、
前記観測量算出手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出するステップと、
前記推定手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記処理領域ｎのシステム状態量の最適値を、前記処理領域ｎのエッジ抽出画像情報として推定するステップと、
前記入力手段が、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力するステップと、
前記追跡手段が、前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測するステップと、
を含む動体の移動量計測方法。
【請求項８】
モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、かつ、画像取得手段、位置変化量推定手段、システム状態量算出手段、観測量算出手段、推定手段、入力手段、及び追跡手段を含む動体の移動量計測装置における動体の移動量計測方法であって、
前記画像取得手段が、前記観測画像情報を時系列に複数フレームを取得するステップと、
前記位置変化量推定手段が、前記画像取得手段により取得された連続するフレームの画像間の位置合わせを行い、前記連続するフレーム間における画像内の物体の位置変化量を推定するステップと、
前記システム状態量算出手段が、時刻ｎでの前記観測画像情報におけるシステム状態量を、前記位置変化量推定手段により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移を用いて算出するステップと、
前記観測量算出手段が、時刻ｎでの前記観測画像情報の観測量を算出するステップと、
前記推定手段が、時刻ｎでの前記観測画像情報における、前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記時刻ｎでのシステム状態量の最適値をエッジ抽出画像情報として推定するステップと、
前記入力手段が、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力するステップと、
前記追跡手段が、前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測するステップと、
を含む動体の移動量計測方法。
【請求項９】
コンピュータを、請求項１〜請求項６のいずれか１項記載の動体の移動量計測装置の各手段として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、動体の移動量計測装置、動体の移動量計測方法及びそのプログラムに関する。
【背景技術】
【０００２】
従来、動体の動きを検出する技術としては、動体の動き速度のｎ倍の高速度の分解能を有する撮影装置で撮影することにより、動きに対するぼけを防止する技術が提案されている。しかし、動体の動きが２次元空間であれば対応が可能であるが、撮影対象の動体の動き速度が必ずしも定まっておらず、動体の想定される動き速度に対して全て対応することは難しかった。更に、動体の動きが３次元空間で移動する場合には、撮影装置の撮影深度問題から生ずるぼけを防止することが出来なかった。また、撮影装置の固有の雑音による撮影した画像に含まれることになる影響は避けることは出来なかった。
【０００３】
近年、画像工学の分野において、動画像のオプティカルフローに関する技術の研究開発が多くなされている。すなわち、動画像に含まれる対象物の動きを追跡する技術が提案されている。
【０００４】
例えば、Lucas-Kanade法では、画像の小領域（例えば５×５ピクセル）のフローベクトルについては変化量が一定であると仮定し、最小二乗法に基づいてオプティカルフローを推定している（例えば、非特許文献１参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】B. D. Lucas and T. Kanade，“An Iterative Image Registration Technique with an Application to Stereo Vision”，1981，Proceedings of Imaging Understanding Workshop，p.121-130
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、非特許文献１に記載の技術では、時刻間の大きな変動への追従は困難である。また、画素の明度変化が小さい領域では推定精度が低下する。
【０００７】
上記非特許文献１に記載の技術では、小領域１つに対して１つのフローベクトルを設定するため、ピクセル単位での移動量（変化量）の導出は困難である、という問題がある。
【０００８】
また、例えば、動画像からピクセル単位の移動量（変化量）の導出を必要とする課題の一例として、新生児の手の動きや足の動きの検知がある。新生児は、脳が発達すると、手で顔を触ったり、体を触ったり、柱を触ったり、また足を開いたり閉じたりする。このような動作を行うことで、新生児は空間把握や機能を確認している。逆にいえば、発達障害又は脳の発達に問題がある新生児の手は、空中をさまようだけになる。また発達障害又は脳の発達に問題がある新生児の足は全く動かない状態になる。
【０００９】
現在は看護師が時間を定めて、新生児の手が顔を触った回数や足の動きの回数を数えており、看護師が忙しいときは、新生児をビデオに録画して、後からビデオを見ながら回数を記録している。このように新生児など発言しないものについては、人間が目視で確認することで異常を早期発見しているのが現場の現状である。
【００１０】
また、例えば、不整脈は、心臓の心拍を医師が聞くことで判明するが、心臓のどの部分に異常があるか、又はないかの判定は、目視で確認しているのが現状である。その他、診断を行う際には、患者の体の一部を切って患部等の動きを見るのは簡単だが、患者の体の一部を切ると患者に負担がかかってしまう。
【００１１】
上記のような分野において、上記非特許文献１等に記載の技術を適用しようとすると、様々な問題がある。
【００１２】
本発明は、上記の事情に鑑みなされたものであり、エッジ抽出画像情報上の指定された箇所の微小な移動量を計測することができる動体の移動量計測装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１３】
上記目的を達成するために、第１の発明に係る動体の移動量計測装置は、モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、前記観測画像情報を時系列に複数フレームを取得する画像取得手段と、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出するシステム状態量算出手段と、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出する観測量算出手段と、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記処理領域ｎのシステム状態量の最適値を、前記処理領域ｎのエッジ抽出画像情報として推定する推定手段と、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力する入力手段と、前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測する追跡手段と、を含んで構成されている。
【００１４】
第２の発明に係る動体の移動量計測方法は、モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、かつ画像取得手段、システム状態量算出手段、観測量算出手段、推定手段、入力手段、及び追跡手段を含む動体の移動量計測装置における動体の移動量計測方法であって、前記画像取得手段が、前記観測画像情報を時系列に複数フレームを取得すると、システム状態量算出手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出するステップと、前記観測量算出手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出するステップと、前記推定手段が、前記取得されたフレームにおける前記観測画像情報内の時刻ｎでの処理領域ｎの前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記処理領域ｎのシステム状態量の最適値を、前記処理領域ｎのエッジ抽出画像情報として推定するステップと、前記入力手段が、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力するステップと、前記追跡手段が、前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測するステップと、を含んで構成されている。
【００１５】
第１の発明に係る動体の移動量計測装置、及び第２の発明に係る動体の移動量計測方法における前記処理領域ｎは、前記処理領域（ｎ−１）と一部重複するように設定されるようにすることができる。
【００１６】
第３の発明に係る動体の移動量計測装置は、モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、前記観測画像情報を時系列に複数フレームを取得する画像取得手段と、前記画像取得手段により取得された連続するフレームの画像間の位置合わせを行い、前記連続するフレーム間における画像内の物体の位置変化量を推定する位置変化量推定手段と、時刻ｎでの前記観測画像情報におけるシステム状態量を、前記位置変化量推定手段により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移を用いて算出するシステム状態量算出手段と、時刻ｎでの前記観測画像情報の観測量を算出する観測量算出手段と、時刻ｎでの前記観測画像情報における、前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記時刻ｎでのシステム状態量の最適値をエッジ抽出画像情報として推定する推定手段と、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力する入力手段と、前記推定手段によって推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測する追跡手段と、を含んで構成されている。
【００１７】
第４の発明に係る動体の移動量計測方法は、モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、前記原画像情報にぼけ及び雑音が混在した観測画像情報から前記原画像情報としてエッジ抽出画像情報を推定復元し、推定復元した前記エッジ抽出画像情報上の指定された箇所の移動量を計測する動体の移動量計測装置であって、かつ、画像取得手段、位置変化量推定手段、システム状態量算出手段、観測量算出手段、推定手段、入力手段、及び追跡手段を含む動体の移動量計測装置における動体の移動量計測方法であって、前記画像取得手段が、前記観測画像情報を時系列に複数フレームを取得するステップと、前記位置変化量推定手段が、前記画像取得手段により取得された連続するフレームの画像間の位置合わせを行い、前記連続するフレーム間における画像内の物体の位置変化量を推定するステップと、前記システム状態量算出手段が、時刻ｎでの前記観測画像情報におけるシステム状態量を、前記位置変化量推定手段により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移を用いて算出するステップと、前記観測量算出手段が、時刻ｎでの前記観測画像情報の観測量を算出するステップと、前記推定手段が、時刻ｎでの前記観測画像情報における、前記システム状態量算出手段により算出されたシステム状態量、及び前記観測量算出手段により算出された観測量に基づいて、前記時刻ｎでのシステム状態量の最適値をエッジ抽出画像情報として推定するステップと、前記入力手段が、ユーザによって指定された、前記推定手段により推定されたエッジ抽出画像情報上の、前記移動量を計測する着目箇所を入力するステップと、前記追跡手段が、前記推定手段により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、前記入力手段から入力された前記着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、取得されたフレーム間で逐次移動量を計測するステップと、を含んで構成されている。
【００１８】
第３の発明に係る動体の移動量計測装置、及び第４の発明に係る動体の移動量計測方法における前記状態遷移を、時刻（ｎ−１）におけるシステム状態量と時刻ｎにおけるシステム状態量との相関を表す忘却係数を用いて表したようにすることができる。
【００１９】
本発明は、前記観測対象を臓器とすることができる。
【００２０】
本発明は、前記観測対象を腸とすることができる。
【００２１】
本発明に係るプログラムは、コンピュータを、上記の動体の移動量計測装置の各手段として機能させるためのプログラムである。
【発明の効果】
【００２２】
以上説明したように、本発明の動体の移動量計測装置、方法、及びプログラムによれば、原画像情報にぼけ及び雑音が混在した観測画像情報から、算出されたシステム状態量、及び算出された観測量に基づいて、状態空間モデルに基づく予測法により、システム状態量の最適値をエッジ抽出画像情報として推定し、推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、ユーザによって指定された着目箇所に対応するエッジ上の箇所を、隣接するフレーム間で追跡し、フレーム間で逐次移動量を計測することにより、エッジ抽出画像情報上の指定された箇所の微小な移動量を計測することができる、という効果が得られる。
【図面の簡単な説明】
【００２３】
【図１】第１及び第２の実施の形態に係る動体の移動量計測装置の概略構成を示すブロック図である。
【図２】第１の実施の形態における画像処理部の概略構成を示すブロック図である。
【図３】第１及び第２の実施の形態における状態空間モデルを説明するための図である。
【図４】第１の実施の形態における画像処理の内容を示すフローチャートである。
【図５】状態方程式を説明するための図である。
【図６】第１の実施の形態におけるシステム状態量算出処理の内容を示すフローチャートである。
【図７】観測過程を説明するための図である。
【図８】観測方程式を説明するための図である。
【図９】状態空間に基づく予測処理の内容を示すフローチャートである。
【図１０】追跡波形処理の内容を示すフローチャートである。
【図１１】追跡波形を生成するための処理を説明するための図である。
【図１２】観測画像内の２点の着目点について追跡波形を生成した例を示す図である。
【図１３】ユーザに表示する追跡画像の内容を説明するための図である。
【図１４】処理結果の一例を示す図である。
【図１５】第２の実施の形態における画像処理部の概略構成を示すブロック図である。
【図１６】第２の実施の形態における画像処理の内容を示すフローチャートである。
【図１７】フレーム間のブロックマッチングを説明するための図である。
【図１８】各フレームの時系列を説明するための図である。
【図１９】第２の実施の形態におけるシステム状態量算出処理の内容を示すフローチャートである。
【図２０】第２の実施の形態における時刻（ｎ−１）の状態ベクトルの構成を説明するための図である。
【図２１】第２の実施の形態における時刻ｎの状態ベクトルの構成を説明するための図である。
【図２２】同一のシステム状態量を持つ領域を説明するための図である。
【図２３】同一のシステム状態量を持つ領域を説明するための図である。
【図２４】第２の実施の形態における忘却係数を説明するための図である。
【図２５】第２の実施の形態における状態方程式の構成を示すイメージ図である。
【図２６】第２の実施の形態における状態方程式の構成の他の例を示すイメージ図である。
【図２７】本発明の実施の形態を医療分野に適用させた場合の例を説明する説明図である。
【図２８】本発明の実施の形態を安全・安心分野に適用させた場合の例を説明する説明図である。
【図２９】本発明の実施の形態を植物の生産・管理分野に適用させた場合の例を説明する説明図である。
【発明を実施するための形態】
【００２４】
＜概要＞
まず、本発明の実施の形態の概要について説明する。
【００２５】
本発明の実施の形態に係る動体の移動量計測装置は、モデル化可能な観測対象のシステムにおける状態空間モデルであって、原画像情報を駆動源として含む状態空間モデルに基づく予測法を用いて、原画像情報にぼけ及び雑音が混在した撮影された観測画像情報を時系列に複数フレーム取得し、算出されたシステム状態量、及び算出された観測量に基づいて、状態空間モデルに基づく予測法により、システム状態量の最適推定値を算出し、撮影された観測画像から原画像情報（エッジ抽出画像情報）を推定復元する。
【００２６】
また、本発明の実施の形態に係る動体の移動量計測装置は、推定原画像情報上に、ユーザによって指定された箇所について、フレーム間で当該箇所を追跡し、フレーム間各々での当該点の移動量を計測することにより、指定された動物体の動きを追跡する。本実施の形態では、医療画像として腹腔のＭＲＩ撮像画像から腸の動きを検出するため、ユーザにより指定された腸の特定箇所の動きを計測し、指定された腸の腸壁の動きを追跡することにより、指定された腸の時間的経過に対応した移動量を計測する装置に本発明を適用した場合について説明する。
尚、実施形態として、腸の動きの移動量を計測する態様について以下に説明するが、この実施形態に限定されるものではなく、他の用途に応用してもよい。
【００２７】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２８】
＜第１の実施の形態＞
まず、第１の実施の形態について説明する。第１の実施の形態では、腸を撮影した観測画像情報を対象とし、複数のフレームの各々について、当該フレームの処理領域１，・・・，処理領域（ｎ−１）に基づいて、処理領域ｎのエッジ抽出画像を推定し、推定されたエッジ抽出画像から、ユーザにより指定された任意の着目箇所を、撮影した観測画像情報上で追跡する動体の移動量計測装置に本発明を適用した場合について説明する。
【００２９】
図１に示すように、第１の実施の形態に係る動体の移動量計測装置１０は、ＣＰＵ、ＲＯＭ、ＲＡＭ、及びＨＤＤを含むコンピュータで構成されており、入力インタフェース部２０、操作部３０、記憶部４０、画像処理部６０、追跡波形処理部７０、追跡波形解析部７５、及び出力インタフェース部８０を備えている。入力インタフェース部２０は、画像取得手段の一例である。
【００３０】
入力インタフェース部２０は、時系列に観測された複数フレームの観測画像である画像データを取得する。すなわち、入力インタフェース部２０は、画像入力装置から提供される画像データをコンピュータで処理可能なデータ形式に変換するなどの入力処理を行う。
【００３１】
画像入力装置は、処理の対象となる画像データ（観測画像）をデジタルデータとしてコンピュータに入力するための入力装置であり、動画像を構成する時系列に連続した各フレームの各画像データを動体の移動量計測装置１０に入力する。画像入力装置としては、例えば、カメラ、記録メディア、モデムなどを用いることができる。カメラは、撮像機能を有する全ての装置を意味し、例えば、デジタルビデオカメラの他に、カメラ機能を搭載した携帯電話や、防犯カメラ（監視カメラ）、画像診断を行うための医療機器（内視鏡、レントゲン、エコー、ＣＴ、ＭＲＩなど）などを含むことができる。記録メディアは、画像データを記録可能な記録媒体を広く意味し、例えば、磁気ディスク（ＨＤＤやＦＤなど）や、光ディスク（ＣＤやＤＶＤ、ＢＤなど）、光磁気ディスク（ＭＯ）、フラッシュメモリ（メモリカードやＵＳＢメモリなど）などである。モデムは、外部の通信ネットワーク（例えば、電話回線やＬＡＮ、インターネットなど）と接続するための装置である。なお、画像入力装置の種類は、動体の移動量計測装置１０の用途に応じて適用可能なものを選択すればよい。本実施の形態では、処理の対象となる画像データが腸の動画像であるため、画像入力装置はＭＲＩを想定する。
【００３２】
図示は省略するが、入力インタフェース部２０は、画像入力装置の種類に応じて別個独立に設けられる。例えば、記録メディアの入力インタフェース部２０は、ドライブとも呼ばれ、記録メディアの種類に応じて様々な種類のドライブが使用可能である。なお、ドライブは、記録メディアを読み書きする装置であり、記録メディアに関する限り、通常は、入力インタフェース部２０と出力インタフェース部８０とは一体化されている。また、通常、モデムは、画像入力装置としても画像出力装置としても機能しうるため、モデムに関しても、通常、入力インタフェース部２０と出力インタフェース部８０とは一体化されている。入力インタフェース部２０は、コンピュータ本体の内部に格納された内蔵カード（ボード）であってもよいし、内部インタフェース部（図示省略）を介して接続された外部設置型機器であってもよい。
【００３３】
なお、画像情報をアナログデータとして入力される画像入力装置の場合、対応する入力インタフェース部２０は、サンプリング部及びＡ／Ｄ変換部（共に図示せず）を有する。サンプリング部は、所定のサンプリング周波数で、入力されたアナログ信号をサンプリング処理し、Ａ／Ｄ変換部に出力する。サンプリング周波数は、復元処理対象（情報源）の種類に応じて変更可能である。Ａ／Ｄ変換部は、サンプリングされた信号の振幅値を所定の分解能でＡ／Ｄ変換処理する。
【００３４】
操作部３０は、例えば、キーボートやマウス、タッチパネル等であるが、音声認識装置などを用いてもよい。使用者は、操作部３０を用い、本装置を操作することができる。また、操作部３０は、パラメータ設定部３１、領域指定部３２、処理フレーム指定部３３、画像情報設定部３４、及び着目点指定部３５を有する。パラメータ設定部３１は、使用者の入力操作により、本実施の形態における画像処理に必要な各種パラメータの値を設定する。領域指定部３２は、使用者の入力操作により、入力された画像に対して画像処理の対象となる処理領域（画像の特定の範囲）を指定する。処理フレーム指定部３３は、使用者の入力操作により、後述する処理フレーム数を指定する。画像情報設定部３４は、入力された画像を白黒画像、又はＲＧＢ画像の何れかとするかを設定する。着目点指定部３５は、後述する追跡波形生成処理において、ユーザの入力操作により、着目箇所の位置を指定する。着目箇所とは、後述する推定部６４により推定されたエッジ抽出画像情報上の、移動量を計測する箇所である。着目点指定部３５は、入力手段の一例である。
【００３５】
記憶部４０は、コンピュータの本体を構成する一要素であって、主記憶部４１、補助記憶部４７、及び表示メモリ４８を有する。主記憶部４１は、入力インタフェース部２０を介して取り込んだ画像データが記憶される取込画像メモリ４２と、後述する追跡波形処理部７０によって生成された追跡波形が記憶される追跡波形保存メモリ４３と、後述する推定部６４によって推定された処理結果が記憶される被写体領域メモリ４４と、後述する画像処理を実行するための画像処理プログラムや各種パラメータが記憶された処理パラメータメモリ４５と、追跡波形の考察結果が記憶される波形考察結果メモリ４６とを有する。取込画像メモリ４２は、少なくとも後述する画像処理において必要となる複数フレーム分の観測画像を含む画像データを記憶可能な容量とする。
【００３６】
表示メモリ４８は、ディスプレイ等で構成された画像出力装置（図示省略）に表示するための画像データが記憶される領域であり、取り込んだ画像（観測画像）を表示するために処理前の画像データ（取り込んだ画像データ）が記憶される取込画像メモリ４９と、画像処理において一時的に必要な記憶領域である処理メモリ５０と、処理された画像を表示するために処理後の画像データが記憶される処理画像メモリ５１とを有する。
【００３７】
補助記憶部４７は、主記憶部４１の容量不足を補う。補助記憶部４７は、例えば、ハードディスクドライブ（ＨＤＤ）であってもよいし、ＣＤ−ＲＯＭ、ＤＶＤ、ＳＳＤ（Solid State Drive）、フラッシュメモリなどの可搬性のあるものであってもよいし、また、それらの組み合わせであってもよい。
【００３８】
また、本実施の形態では、処理パラメータメモリ４５に予め画像処理プログラムが記憶されている場合について説明するが、画像処理プログラムは、記録メディアから入力インタフェース部２０を介して記憶部４０にインストールしたり、モデム及び入力インタフェース部２０を介して外部から記憶部４０にダウンロードしたりしてもよい。
【００３９】
画像処理部６０は、後述する画像処理を実行するための処理部であり、機能的には、図２に示すように、領域設定部６１と、システム状態量算出部６２と、観測量算出部６３と、推定部６４と、エッジ抽出画像出力部６５とを含んだ構成で表すことができる。
【００４０】
領域設定部６１は、連続するフレームの画像の各々について、Ｎ個の処理領域を設定する。領域設定部６１は、処理領域ｎを、処理領域（ｎ−１）と一部重複するように設定する。
【００４１】
システム状態量算出部６２は、取込画像メモリ４２に記憶された画像データの特定フレームにおける観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量を算出する。
具体的には、システム状態量算出部６２は、取込画像メモリ４２に記憶された画像データに含まれる観測画像について、フレーム毎に、当該観測画像内の処理領域ｎの各々について、当該処理領域ｎにおけるシステム状態量を、処理領域（ｎ−１）のシステム状態量、処理領域（ｎ−１）から処理領域ｎへの状態遷移、及び処理領域（ｎ−１）のシステム状態量に含まれない処理領域ｎのシステム状態量の要素を含む有色性駆動源を用いて算出する。なお、第１の実施の形態では、処理領域ｎを時々刻々と変化させて推定を行うため、「時刻ｎでの処理領域ｎ」と表現する。
【００４２】
観測量算出部６３は、取込画像メモリ４２に記憶された画像データの特定フレームにおける観測画像情報内の時刻ｎでの処理領域ｎの観測量を算出する。
具体的には、観測量算出部６３は、取込画像メモリ４２に記憶された画像データに含まれる観測画像について、フレーム毎に、当該観測画像内の処理領域ｎの各々について、観測画像の処理領域ｎの情報、処理領域ｎのシステム状態量、ぼけ及び雑音を用いて、エッジ抽出画像にぼけ及び雑音が加わって観測画像となる過程を観測量として算出する。
【００４３】
推定部６４は、取込画像メモリ４２に記憶された画像データの特定フレームにおける観測画像情報内の時刻ｎでの処理領域ｎのシステム状態量算出部６２により算出されたシステム状態量、及び観測量算出部６３により算出された観測量に基づいて、処理領域ｎのシステム状態量の最適値を、当該処理領域ｎのエッジ抽出画像情報として推定する。
具体的には、推定部６４は、取込画像メモリ４２に記憶された画像データに含まれる観測画像について、フレーム毎に、当該観測画像内の処理領域ｎの各々について、システム状態量算出部６２により算出されたシステム状態量、及び観測量算出部６３により算出された観測量に基づいて、状態空間モデルに基づく予測法により、当該処理領域ｎのシステム状態量の最適値を、当該処理領域ｎのエッジ抽出画像として推定する。
【００４４】
エッジ抽出画像出力部６５は、推定部６４により得られたフレーム毎の各処理領域ｎのエッジ抽出画像を結合して、フレーム毎にエッジ抽出画像を生成する。
【００４５】
なお、画像処理部６０の各部の処理については、後述する画像処理において詳述する。
【００４６】
追跡波形処理部７０は、推定部６４により推定されたエッジ抽出画像情報に含まれるエッジ上の箇所であって、かつ、着目点指定部３５から入力された着目箇所に対応するエッジ上の箇所を隣接するフレーム間で追跡し、フレーム間で逐次移動量を計測する。
ここで、取込画像メモリ４２に取込まれた画像は出力インタフェース部８０を介し、図示しない画像出力部に表示される。画像出力部に表示される画像上でユーザにより指定された箇所は、エッジ抽出画像出力部６５によって生成されたエッジ抽出画像上に対応して反映して表示され、ユーザによって指定された箇所は、エッジ画像上のフレーム間で当該点を追跡し、フレーム間各々での当該箇所の移動量を計測する。そして、追跡波形処理部７０は、エッジ上の箇所の各々の追跡結果及び移動量に基づいて、追跡波形を生成し、出力インタフェース部８０を介して、追跡波形を画像出力装置に表示させることにより、可視化する。追跡波形処理部７０で処理され、指定された箇所の追跡結果及び移動量は、取込画像メモリ４２の取り込み画像情報に取り込まれ、画像出力部に表示出力される。なお、追跡波形処理部７０の処理については、後述する追跡波形生成処理において詳述する。追跡波形処理部７０は、追跡手段の一例である。
【００４７】
追跡波形解析部７５は、追跡波形処理部７０によって得られた追跡波形に対して、ウェーブレット解析等の処理を実行することで、解析結果を出力する。
【００４８】
出力インタフェース部８０は、画像処理部６０で復元処理され、処理画像メモリ５１に記憶された画像データ（エッジ抽出画像）を、ディスプレイなどで構成された画像出力装置に出力可能なデータ形式に変換して出力処理を行うと共に、追跡波形処理部７０で生成された追跡波形及び移動量の出力処理を行う。
【００４９】
なお、処理画像の出力形態は画像出力装置に表示する場合に限定されず、他の画像出力装置に出力するようにしてもよい。例えば、画像出力装置として、プリンタ、記録メディア、モデムなどを用いることができる。記録メディア及びモデムは、画像入力装置としての記録メディア及びモデムとそれぞれ共用してもよい。なお、画像出力装置の種類は、動体の移動量計測装置１０の用途に応じて適用可能なものを選択すればよい。
【００５０】
また、図示は省略するが、出力インタフェース部８０は、画像出力装置の種類に応じて別個独立に設けられる。上記のように、記録メディア及びモデムに関しては、通常、入力インタフェース部２０と出力インタフェース部８０とは一体化されている。出力インタフェース部８０も、入力インタフェース部２０と同様に、コンピュータ本体の内部に格納された内蔵カード（ボード）であってもよいし、内部インタフェース部を介して接続された外部設置型機器であってもよい。
【００５１】
次に、第１の実施の形態における画像処理の理論となる観測画像モデル及び状態空間モデルについて説明する。
【００５２】
観測画像は、一般的に、原画像にぼけ及び雑音を加えたモデルによって得られる。雑音は、画素に関係なく発生し、暗電流や熱雑音などに起因して不可避的に発生する。従って、第１の実施の形態におけるエッジ抽出画像の生成は、観測画像からぼけ及び雑音を取り除き、エッジを抽出することである。
【００５３】
第１の実施の形態では、エッジ抽出画像情報のみを用いてシステム状態量を算出し、かつ、観測画像情報、エッジ抽出画像情報、及びぼけ及び雑音を用いて観測量を算出した状態空間モデルであって、かつ、有色性駆動源を含む状態空間モデルに基づく予測法により、システム状態量の最適推定値を求めることにより、観測画像からエッジ抽出画像を生成する。
【００５４】
上記の状態空間モデルは、下記（１）式の状態方程式と、下記（２）式の観測方程式とで表される。ただし、状態方程式において、ベクトルｘｐ１（ｎ）は時刻ｎの状態ベクトル（エッジ抽出画像のシステム状態量）、行列Φｐ１は状態遷移行列、ベクトルδｐ１（ｎ）は駆動源ベクトルである。また、観測方程式において、ベクトルｙｐ１（ｎ）は観測ベクトル（観測画像情報）、行列Ｍは観測行列、ベクトルεｐ１（ｎ）は観測雑音ベクトルである。また、添え字ｐ１は第１の実施の形態に係るものであることを表す。
【００５５】
【数１】
2015152975.tif 000003
【００５６】
上記（１）式、及び（２）式をブロック線図で表したものを図３に示す。図３のうち、状態遷移過程が上記（１）式に対応し、観測過程が上記（２）式に対応する。
【００５７】
（１）式の状態方程式は、観測対象のシステムを状態空間モデルで記述したものであり、内部状態つまり状態変数（ここでは、状態ベクトルｘｐ１（ｎ−１）、ｘｐ１（ｎ））の時間に対する生成過程を表している。また、（２）式の観測方程式は、何らかの観測装置を通じて観測する過程を記述したものであり、観測結果（ここでは、観測ベクトルｙｐ１（ｎ））が、被観測量つまり入力（ここでは、状態ベクトルｘｐ１（ｎ））に依存して生成される様子を示している。
【００５８】
次に、図４を参照して、第１の実施の形態に係る動体の移動量計測装置１０の作用について説明する。なお、図４に示す画像処理ルーチンは、記憶部４０に制御プログラムとして格納されており、図示しないＣＰＵによって実行される。
【００５９】
ステップ１００で、画像入力装置から入力インタフェース部２０を介して、所定フレームの観測画像を取り込み、主記憶部４１の取込画像メモリ４２に記憶する。また、表示メモリ４８の取込画像メモリ４９にも記憶し、出力インタフェース部８０を介して、画像出力装置に取込画像を表示する。
【００６０】
次に、ステップ１０２で、上記ステップＳ１００で記憶した所定フレームの画像情報のうち、処理対象とする特定フレーム数を設定する。
【００６１】
次に、ステップ１０４で、領域設定部６１が、画像処理の処理単位となる処理領域のサイズを設定する。処理領域のサイズの設定は、ユーザにより操作部３０の領域指定部３２から指定された値を設定するようにしてもよい。
【００６２】
ここでは、ステップ１０４において、処理領域のサイズｐ×ｑ＝３×３画素と設定されたものとする。この場合、図５に示すように、左から右へ流れている、処理領域（ｎ−１）、処理領域ｎという画像情報が処理対象となる。
【００６３】
次に、ステップ１０５で、システム状態量算出部６２が、処理領域ｎと処理領域（ｎ−１）との位置関係に基づいて、状態遷移行列Φｐ１を構成する。
【００６４】
具体的には、まず、処理領域（ｎ−１）から処理領域ｎ間に、水平に左から右へ１画素ずつ処理領域を状態変化させると仮定する。この場合、各処理領域において、上記図５に示すように、同一のシステム状態量を持つ領域が存在する。同一のシステム状態量を持つ領域は、上記図５において網掛けされた画素領域である。上記図５に示すように、状態（ｎ−１）のシステム状態量ｘ１，２，ｘ２，２，ｘ３，２，ｘ１，３，ｘ２，３，ｘ３，３と、状態ｎのシステム状態量ｘ１，１，ｘ２，１，ｘ３，１，ｘ１，２，ｘ２，２，ｘ３，２とが同一のシステム状態量を持つ領域となる。
【００６５】
そして、上述した同一のシステム状態量を持つ領域が状態（ｎ−１）から状態ｎへの変化によって一致した状態遷移行列Φｐ１を構成する。
【００６６】
次に、ステップ１０６で、領域設定部６１が、上記ステップＳ１０４で設定されたサイズの処理領域（ｎ−１）、処理領域ｎを、上記ステップ１０２で設定されたフレームの画像情報に対して設定する。
【００６７】
次に、ステップ１０８で、システム状態量算出部６２が、システム状態量算出処理を実行する。ここで、図６を参照して、システム状態量算出処理ルーチンについて説明する。
【００６８】
ステップ１０８０で、システム状態量算出部６２が、処理領域（ｎ−１）内の各画素のシステム状態量を要素とする時刻（ｎ−１）の状態ベクトル（ｎ−１）と、処理領域ｎ内の各画素のシステム状態量を要素とする時刻ｎの状態ベクトルｎとを導出する。以下、各時刻の状態ベクトルの導出方法を具体的に説明する。
【００６９】
まず、時刻（ｎ−１）の状態ベクトル（ｎ−１）は、上記図５に示すように、ｐ×ｑサイズの処理領域（ｎ−１）内の各画素のシステム状態量を組み合わせて構成される。すなわち、時刻（ｎ−１）のシステム状態ベクトル（ｎ−１）は９次元のベクトルｘｐ１（ｎ−１）＝［ｘ１，１（ｎ−１）  ｘ２，１（ｎ−１）  ・・・  ｘ３，３（ｎ−１）］Ｔとなる。次に、時刻ｎの状態ベクトルｎは、上記図５に示すように、ｐ×ｑサイズの処理領域ｎ内の各画素のシステム状態量を組み合わせて構成される。すなわち、時刻ｎの状態ベクトルｎは９次元のベクトルｘｐ１（ｎ）＝［ｘ１，１（ｎ）  ｘ２，１（ｎ）  ・・・  ｘ３，３（ｎ）］Ｔとなる。ただし、Ｔはベクトルの転置を表す。また、状態ベクトルｘｐ１（ｎ−１）及びｘｐ１（ｎ）のサイズ（次元数）は設定された処理領域ｐ×ｑのサイズによって変化する。
【００７０】
次に、ステップ１０８２で、システム状態量算出部６２が、駆動源ベクトルを構成する。上記図５の時刻ｎの状態ベクトルｎのシステム状態量ｘ１，３，ｘ２，３，ｘ３，３は、時刻（ｎ−１）の状態ベクトル情報からは構成することができないため、本実施の形態では、駆動源と呼ばれる領域にシステム状態量ｘ１，３，ｘ２，３，ｘ３，３を用いることで、駆動源ベクトルδｐ１（ｎ）を構成する。これは、物体がフレーム内を位置変化した際に、時刻（ｎ−１）の状態ベクトル（ｎ−１）のシステム状態量では表せられない領域のシステム状態量を、時刻ｎの状態ベクトルｎのシステム状態量で構成することを表しており、エッジ抽出画像の情報が含まれることから、駆動源ベクトルδｐ１（ｎ）は有色性駆動源となる。
【００７１】
次に、ステップ１０８４で、システム状態量算出部６２が、上記ステップ１０８０で導出した時刻（ｎ−１）の状態ベクトル（ｎ−１）及び時刻ｎの状態ベクトルｎの状態ベクトルｘｐ１（ｎ−１）及びｘｐ１（ｎ）、上記ステップ１０５で構成した状態遷移行列Φｐ１、及び上記ステップ１０８２で構成した駆動源ベクトルδｐ１（ｎ）に基づいて、システム状態量を算出する。上記図５に示すように、本実施の形態における状態方程式は、クリアな画像からなる状態ベクトルｘｐ１（ｎ）と、０及び１からなる状態遷移行列Φｐ１と、有色信号であるクリアな画像からなる駆動源ベクトルδｐ１（ｎ）とから構成されるため、所望のシステム状態量（エッジ抽出画像の画素情報）、つまり、クリアな画像情報（エッジ抽出画像情報）のみで構成されることになる。
【００７２】
システム状態量算出処理ルーチンが終了すると、画像処理ルーチンにリターンして、ステップ１１０へ移行し、観測量算出部６３が、観測画像モデルに基づいて、図７に示すように、エッジ抽出画像情報にぼけ及び雑音が加わって観測画像情報となる過程を、観測量（観測方程式）として算出する。観測ベクトルｙｐ１（ｎ）は、観測画像情報を含む観測量として時刻ｎの状態ベクトルｎの処理領域内の各画素の画素情報を要素とする９次元のベクトルである。
【００７３】
また、上記図７に示すように、ｚｐ１（ｎ）を原画像としたとき、原画像ｚｐ１（ｎ）とエッジ抽出画像ベクトルｘｐ１（ｎ）との関係は、以下のように表すことができる。
【００７４】
【数２】
2015152975.tif 000004
【００７５】
そして、上記（２）式をエッジ抽出画像ベクトルｘｐ１（ｎ）の式に書き換えると、以下のように表すことができる。
【００７６】
【数３】
2015152975.tif 000005
【００７７】
すなわち、観測行列Ｍは、エッジ抽出行列Ｍ−１の一般化逆行列により導出することができる。エッジ抽出行列Ｍ−１は、従来既知のエッジ抽出フィルタに基づいて予め定められるため、観測行列Ｍの各要素は、既知であって、予めデータ化して適切に設定されている。また、観測雑音ベクトルεｐ１（ｎ）は、処理領域に含まれる画素に対応する観測雑音を要素とする９次元のベクトルである。
【００７８】
上記ステップ１１０によって算出された観測量（観測方程式）の対応関係を図８に示す。図８に示すように、観測ベクトルｙｐ１（ｎ）は、観測行列Ｍとエッジ抽出画像ベクトルｘｐ１（ｎ）との積と、観測雑音ベクトルεｐ１（ｎ）との和によって表せる。観測行列Ｍの各要素は、予め定められている。
【００７９】
次に、ステップ１１２で、推定部６４が、上記ステップ１０８で算出されたシステム状態量（状態方程式）及び上記ステップ１１０で算出された観測量（観測方程式）により、下記に示す有色性駆動源を含む状態空間モデルに基づく予測法を導出する。
【００８０】
【数４】
2015152975.tif 000006
【００８１】
上記のアルゴリズムは、初期設定の過程［Initialization］と反復の過程［Iteration］とに大別され、反復の過程では、１〜５の手順を逐次繰り返す。以下、図９を参照して、有色性駆動源を含む状態空間モデルに基づく予測法のアルゴリズムの詳細について説明する。
【００８２】
ステップ１１２０で、アルゴリズムの処理が初回か否かを判定する。初回の場合には、ステップ１１２２へ移行し、２回目以降の場合には、ステップ１１２４へ移行する。
【００８３】
ステップ１１２２では、初期設定を行う。具体的には、状態ベクトル、つまり、システム状態量である所望信号（エッジ抽出画像信号）ベクトルの最適推定値（以下「所望信号の最適推定値ベクトル」という）の初期値ｘ＾ｐ１（０｜０）、所望信号ベクトルの推定誤差（以下「所望信号の推定誤差ベクトル」という）の相関行列の初期値Ｐｐ１（０｜０）、駆動源ベクトルの共分散の値Ｒδｐ１（ｎ）［ｉ，ｊ］、及び観測雑音ベクトルの共分散の値Ｒεｐ１（ｎ）［ｉ，ｊ］を、上述の初期設定の過程［Initialization］に示した初期状態に設定する。
【００８４】
ただし、ベクトル０Ｊは、Ｊ次元の零ベクトルであり、行列ＩＪは、Ｊ次の単位行列である。また、Ｊは、√Ｊ×√Ｊの大きさの処理領域の画素数である。本実施の形態では、処理領域のサイズｐ×ｑ＝３×３＝９が処理領域の画素数である。また、Ｅ［ｘ＾ｐ１（ｎ），ｘ＾ｐ１Ｔ（ｎ）］は、駆動源ベクトルｘ＾ｐ１（ｎ）の自己相関行列の期待値である。Ｅ［εｐ１（ｎ）εｐ１Ｔ（ｎ）］［ｉ，ｊ］は、観測雑音ベクトルεｐ１の自己相関行列の期待値であり、ここではσ２ε［ｉ，ｊ］と等しく、既知と仮定している。ここでいう「既知」とは、別の任意の方法（アルゴリズム）で求められて与えられることを意味する。もし観測雑音ベクトルεｐ１（ｎ）が白色雑音であり零平均であれば、σ２εは所定のサンプル数の分散で与えられる。なお、ベクトル及び行列の要素を示す場合、ベクトルａ（ｎ）のｉ番目の要素をａ（ｎ）［ｉ］と表記し、行列Ａ（ｎ）のｉ行ｊ列の要素をＡ（ｎ）［ｉ，ｊ］と表記することとする。
【００８５】
次に、ステップ１１２４で、上記ステップ１０５で定義した状態空間モデルにおける状態遷移行列Φｐ１、上記ステップ１１２２で設定した所望信号の推定誤差ベクトルの相関行列の初期値Ｐｐ１（０｜０）（ｎ＝２の場合）、または１時刻前に後述するステップ１１３０で更新された相関行列Ｐｐ１（ｎ｜ｎ）（ｎ＞２の場合）、及び駆動源ベクトルの共分散の値Ｒδｐ１（ｎ）［ｉ，ｊ］を用いて、状態（ｎ−１）までの情報により時刻ｎの状態ベクトルを推定した場合の誤差である相関行列Ｐｐ１（ｎ｜ｎ−１）を計算する（上述の反復の過程［Ｉｔｅｒａｔｉｏｎ］の手順１）。
【００８６】
次に、ステップ１１２６で、上記ステップ１１２４で計算した相関行列Ｐｐ１（ｎ｜ｎ−１）、上記ステップ１１０で定義した状態空間モデルにおける観測行列Ｍ、及び観測雑音ベクトルの共分散の値Ｒεｐ１（ｎ）［ｉ，ｊ］を用いて、重み係数行列Ｋｐ１（ｎ）を計算する（同手順２）。
【００８７】
次に、ステップ１１２８で、状態遷移行列Φｐ１、及び上記ステップ１１２２で設定した所望信号の最適推定値ベクトルの初期値ｘ＾ｐ１（０｜０）（ｎ＝２の場合）、または１時刻前に本ステップで得られた所望信号の最適推定値ベクトルｘ＾ｐ１（ｎ｜ｎ）（ｎ＞２の場合）を用いて、状態（ｎ−１）までの情報による状態ｎでの所望信号の最適推定値ベクトルｘ＾ｐ１（ｎ｜ｎ−１）を計算する（同手順３）。そして、計算した所望信号の最適推定値ベクトルｘ＾ｐ１（ｎ｜ｎ−１）、上記ステップ１１２６で計算した重み係数行列Ｋｐ１（ｎ）、観測ベクトルｙｐ１（ｎ）、及び観測行列Ｍを用いて、状態ｎまでの情報によるその時刻での所望信号の最適推定値ベクトルｘ＾ｐ１（ｎ｜ｎ）を計算する（同手順４）。
【００８８】
次に、ステップ１１３０で、単位行列Ｉ、重み係数行列Ｋｐ１（ｎ）、観測行列Ｍ、及び上記ステップ１１２４で計算された相関行列Ｐｐ１（ｎ｜ｎ−１）を用いて、状態ｎまでの情報によるその状態での相関行列Ｐｐ１（ｎ｜ｎ）を更新する（同手順５）。次に、ステップ１１３２で、上記ステップ１１２８で計算された所望信号の最適推定値ベクトルｘ＾ｐ１（ｎ｜ｎ）を、現在設定されている処理領域の処理結果として一旦被写体領域メモリ４４に記憶する。
【００８９】
有色性駆動源を含む状態空間モデルに基づく予測法のアルゴリズムが終了すると、画像処理ルーチンへリターンして、ステップ１１４へ移行し、上記ステップ１０２で設定されたフレームの全処理領域について画像処理が終了したか否かを判定する。全処理領域について処理が終了していない場合には、ステップ１１５へ移行してｎを１インクリメントして、ステップ１０６へ戻って、次の処理領域を設定して処理を繰り返す。全処理領域について処理が終了した場合には、ステップ１１６へ移行して、上記ステップ１００で取り込んだ全フレームについて処理が終了したか否かを判定する。全フレームについて処理が終了していない場合には、ステップ１１８へ移行してフレームＬを１インクリメントして、ステップ１０２へ戻って、次のフレームについて処理を繰り返す。全フレームについて処理が終了した場合には、ステップ１２０へ移行する。
【００９０】
ステップ１２０では、エッジ抽出画像出力部６５が、上記ステップ１１３０で被写体領域メモリ４４に記憶された各フレームの処理領域の処理結果を結合して、各フレームのエッジ抽出画像を生成し、被写体領域メモリ４４に記憶すると共に、各フレームのエッジ抽出画像を追跡波形処理部７０へ出力する。
【００９１】
次に、ステップ１２２で、追跡波形処理部７０が、追跡波形生成処理を実行する。ここで、図１０を参照して、追跡波形生成処理ルーチンについて説明する。
【００９２】
ステップ１２２０で、追跡波形処理部７０が、上記ステップ１２０で出力された各フレームのエッジ抽出画像を受け付ける。
【００９３】
次に、ステップ１２２２で、追跡波形処理部７０が、出力インタフェース部８０を介して、画像出力装置にフレーム１のエッジ抽出画像を表示する。
【００９４】
次に、上記ステップ１２２２で表示されたフレーム１のエッジ抽出画像をユーザが確認し、ユーザが操作部３０の着目点指定部３５を操作して着目点の位置が入力されると、ステップ１２２４で、追跡波形処理部７０が、着目点の位置の入力を受け付ける。着目点の位置の入力は、例えば、表示されたフレーム１のエッジ抽出画像の任意の点を、ユーザが指定することにより行われる。
【００９５】
次に、ステップ１２２６で、追跡波形処理部７０が、上記ステップ１２２４で受け付けた着目点の位置を用いて、上記ステップ１２２０で受け付けたフレーム１のエッジ抽出画像に含まれるエッジ上の追跡点を設定する。エッジ上の追跡点を設定する方法の一例としては、例えば、着目点の位置に最も近いエッジ上の点を、エッジ上の追跡点と設定することができる。
【００９６】
次に、ステップ１２２８で、追跡波形処理部７０が、Ｌに２を代入する。
【００９７】
次に、ステップ１２３０で、追跡波形処理部７０が、上記ステップ１２２６で設定されたフレーム１のエッジ抽出画像上での追跡点の座標、又は後述するステップ１２３２で前回記憶された追跡点の座標に基づいて、フレームＬ−１とフレームＬとの間のエッジ抽出画像上での追跡点の移動量を算出する。
【００９８】
具体的には、追跡波形処理部７０が、図１１に示すように、フレームＬのエッジ抽出画像上において、フレーム（Ｌ−１）の追跡点の座標の上方向に存在するエッジ上の点までの移動量と、フレーム（Ｌ−１）の追跡点の座標の下方向に存在するエッジ上の点までの移動量とを算出する。
【００９９】
そして、追跡波形処理部７０が、上記図１１に示すように、上記ステップ１２３０で算出された、上方向に存在するエッジ上の点までの移動量と下方向に存在するエッジ上の点までの移動量とを比較して、移動量が少ない方向に存在するエッジ上の点を、フレームＬにおける追跡点と判定し、フレームＬ−１とフレームＬとの間のエッジ抽出画像上での追跡点の移動量を算出する。なお、本実施の形態では、上方向に存在するエッジ上の点までの移動量、及び下方向に存在するエッジ上の点までの移動量が、観測対象とする動体の最大移動量を予め設定することにより、例えば、成人の小腸の伸縮量の場合には４ピクセル以上である場合にはエラーであると判定し、移動量は０とする。なお、本実施の形態では解析対象を成人の小腸とするため、４ピクセル以上である場合をエラーと判定する場合を例に説明したが、これに限られるものではない。解析対象が異なる場合には、解析対象に応じて定められた移動量の範囲を用いて、エラーと判定することができる。
【０１００】
次に、ステップ１２３２で、追跡波形処理部７０が、上記ステップＳ１２３０で算出された追跡点の移動量を、追跡波形保存メモリ４３に記憶させる。また、追跡波形処理部７０が、上記ステップＳ１２３０で求められたフレームＬにおける追跡点の座標を、追跡波形保存メモリ４３に記憶させる。
【０１０１】
次に、ステップ１２３４で、追跡波形処理部７０が、上記ステップ１２２０で受け付けた全てのフレームのエッジ抽出画像について、上記ステップ１２３０〜ステップ１２３２の処理を実行したか否かを判定する。全てのフレームのエッジ抽出画像について、上記ステップ１２３０〜ステップ１２３２の処理を実行した場合には、ステップ１２３８へ移行する。一方、上記ステップ１２３０〜ステップ１２３２の処理を実行していないフレームのエッジ抽出画像が存在する場合には、ステップ１２３６へ移行してＬを１インクリメントして、ステップ１２３０へ戻って、次のフレームを設定して処理を繰り返す。
【０１０２】
次に、ステップ１２３８で、追跡波形処理部７０が、上記ステップ１２２４で受け付けた全ての着目点について、上記ステップ１２２５〜ステップ１２３４の処理を実行したか否かを判定する。全ての着目点について、上記ステップ１２２５〜ステップ１２３４の処理を実行した場合には、ステップ１２４０へ移行する。一方、上記ステップ１２２５〜ステップ１２３４の処理を実行していない着目点が存在する場合には、ステップ１２２５へ戻って、次の着目点を設定して処理を繰り返す。
【０１０３】
次に、ステップ１２４０で、追跡波形処理部７０が、上記ステップ１２３４で追跡波形保存メモリに着目点毎に記憶された、フレーム間の各々の追跡点の移動量に基づいて、上記ステップ１２２４で受け付けた着目点毎に、追跡波形を生成し、追跡波形保存メモリ４３に記憶させる。
【０１０４】
追跡波形生成処理ルーチンが終了すると、画像処理ルーチンにリターンして、ステップ１２４へ移行し、追跡波形解析部７５が、上記ステップ１２２で得られた着目点毎の追跡波形に対して解析を行い、腸内の異常箇所を検出し、検出結果を波形考察結果メモリ４６へ記憶させる。具体的には、追跡波形解析部７５が、着目点毎の追跡波形に対してウェーブレット解析等の処理を実行することで解析結果を出力し、例えば、腸内の異常を検出する。追跡波形解析部７５による解析結果は、波形考察結果メモリ４６に記憶される。
【０１０５】
そして、ステップ１２６で、画像処理部６０が、上記ステップ１２３４で追跡波形保存メモリ４３に記憶された着目点毎の追跡座標の各々を、上記ステップ１００で取込画像メモリ４２に記憶された所定フレームの取り込み画像に転写して追跡画像を各々生成する。
【０１０６】
そして、ステップ１２８で、画像処理部６０が、追跡波形保存メモリ４３に記憶された追跡波形と、波形考察結果メモリ４６に記憶された解析結果とを、出力インタフェース部８０を介して出力して、画像出力装置に表示すると共に、フレーム毎に生成された追跡画像を連続して画像出力装置に表示して、画像処理ルーチンを終了する。例えば、図１２に示すように、上記ステップ１２２４で受け付けた着目点毎に、追跡波形を画像出力装置に表示する。
【０１０７】
なお、上記ステップ１２８で画像出力装置に表示される追跡画像の各々は、図１３に示すように、エッジ抽出画像上の追跡点の各々を、上記ステップ１００で受け付けた観測画像の各々に転写したものである。すなわち、フレーム毎に生成された追跡画像を連続して画像出力装置に表示することにより、ユーザ（医師等）は普段見慣れている観測画像上で、腸壁の変動を確認することができる。図１４に、画像出力装置に表示される処理結果の一例を示す。なお、追跡画像の代わりに、着目点毎の追跡座標の各々をエッジ抽出画像に転写して生成した画像を画像出力装置に表示してもよい。
【０１０８】
以上説明したように、第１の実施の形態の動体の移動量計測装置１０によれば、観測画像を時系列に複数フレーム取得し、フレーム毎に、当該観測画像内の処理領域ｎの各々について、当該処理領域ｎにおけるシステム状態量を、処理領域（ｎ−１）のシステム状態量、処理領域（ｎ−１）から処理領域ｎへの状態遷移、及び処理領域（ｎ−１）のシステム状態量に含まれない処理領域ｎのシステム状態量の要素を含む有色性駆動源を用いて算出し、算出されたシステム状態量、及び算出された観測量に基づいて、状態空間モデルに基づく予測法により、当該処理領域ｎのシステム状態量の最適値を、当該処理領域ｎのエッジ抽出画像として推定し、推定されたエッジ抽出画像に含まれるエッジ上の点であって、かつ、ユーザによって指定された点について、フレーム間で当該点を追跡し、フレーム間各々での当該点の移動量を計測することにより、腸壁上の点の微小な移動量を計測することができ、可視化することもできる。すなわち、ピクセル単位で腸壁の変動を追跡することができる。
【０１０９】
また、腸壁の変動に対する客観的な評価をすることができる。これにより、医師の判断を支援することができる。
【０１１０】
また、処理領域毎に処理していくことで、処理する行列のサイズを小さくして演算量を抑えることができ、並列処理を適用することで、画像処理を高速に実行することができる。
【０１１１】
また、有色性駆動源を含む状態空間モデルに基づく予測法を用いることにより、演算量を抑えることができ、画像処理を高速に実行することができる。
【０１１２】
また、着目点毎の追跡波形に対してウェーブレット解析等の処理を実行することで、解説結果を出力することができ、例えば、腸内の異常を検出することができる。これにより、医師の判断を支援することができる。
【０１１３】
また、ユーザによって指定された動体上の観測対象の箇所のみについて、フレーム間で当該箇所を追跡し、フレーム間各々での当該箇所の移動量を計測することにより、演算量を抑えることができ、画像処理を高速に実行することができる。
【０１１４】
＜第２の実施の形態＞
次に、第２の実施の形態について説明する。第２の実施の形態では、腸の動画像を対象とし、複数のフレームの各々について、時刻（ｎ−１）の（Ｌ−１）フレームまでに得られる画像情報に基づいて、時刻（ｎ）の（Ｌ）フレームのエッジ抽出画像を推定し、推定されたエッジ抽出画像からエッジ上の箇所を追跡する動体の移動量計測装置に本発明を適用した場合について説明する。なお、第２の実施の形態における「時刻（ｎ）」は、上記第１の実施の形態における「時刻（ｎ）」には対応せず、第２の実施の形態における「時刻（ｎ）」は、複数フレームの各々におけるフレーム（Ｌ）の時刻に対応する。
【０１１５】
第２の実施の形態に係る動体の移動量計測装置２１０の構成について、第１の実施の形態に係る動体の移動量計測装置１０と異なる部分について説明し、他の部分については第１の実施の形態と同一符号を付して詳細な説明を省略する。
【０１１６】
図１５に示すように、第２の実施の形態における画像処理部２６０は、位置変化量推定部２６１と、システム状態量算出部２６２と、観測量算出部２６３と、推定部２６４と、エッジ抽出画像出力部６５とを含んだ構成で表すことができる。
【０１１７】
位置変化量推定部２６１は、連続するフレームの画像間の位置合わせを行い、連続するフレーム間における画像内の物体の位置変化量を推定する。
【０１１８】
システム状態量算出部２６２は、時刻ｎでの観測画像情報におけるシステム状態量を、位置変化量推定部２６１により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移を用いて算出する。
具体的には、システム状態量算出部２６２は、画像内の処理領域毎に、時刻ｎにおけるシステム状態量を、時刻（ｎ−１）までに取得されたフレームから得られる時刻（ｎ−１）のシステム状態量、位置変化量推定部２６１により推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移、及び時刻（ｎ−１）のシステム状態量に含まれない時刻ｎのシステム状態量の要素を含む有色性駆動源を用いて算出する。
【０１１９】
観測量算出部２６３は、時刻ｎでの観測画像情報の観測量を算出する。
具体的には、観測量算出部２６３は、画像内の処理領域毎に、時刻ｎのフレームの観測画像の情報、時刻ｎのシステム状態量、ぼけ及び雑音を示す情報を用いて、エッジ抽出画像にぼけと雑音が加わって観測画像となる過程を観測量として算出する。
【０１２０】
推定部２６４は、時刻ｎでの観測画像情報における、システム状態量算出部２６２により算出されたシステム状態量、及び観測量算出部２６３により算出された観測量に基づいて、時刻ｎでのシステム状態量の最適値をエッジ抽出画像情報として推定する。
具体的には、推定部２６４は、画像内の処理領域毎に、システム状態量算出部２６２により算出されたシステム状態量、及び観測量算出部２６３により算出された観測量に基づいて、状態空間モデルに基づく予測法により、時刻ｎのシステム状態量の最適値を、時刻ｎのエッジ抽出画像として推定する。
【０１２１】
次に、第２の実施の形態における画像処理の理論となる状態空間モデルについて説明する。
【０１２２】
第２の実施の形態における状態空間モデルは、下記（５）式の状態方程式と、下記（６）式の観測方程式とで表される。ただし、状態方程式において、ベクトルｘｐ２（ｎ）は時刻ｎの状態ベクトル（エッジ抽出画像のシステム状態量）、行列Φｐ２は状態遷移行列、ベクトルδｐ２（ｎ）は駆動源ベクトルである。また、観測方程式において、ベクトルｙｐ２（ｎ）は観測ベクトル（観測画像情報）、行列Ｍは観測行列、ベクトルεｐ２（ｎ）は観測雑音ベクトルである。また、添え字ｐ２は第２の実施の形態に係るものであることを表す。
【０１２３】
【数５】
2015152975.tif 000007
【０１２４】
（５）式の状態方程式は、観測対象のシステムを状態空間モデルで記述したものであり、内部状態つまり状態変数（ここでは、状態ベクトルｘｐ２（ｎ−１）、ｘｐ２（ｎ））の時間に対する生成過程を表している。また、（６）式の観測方程式は、何らかの観測装置を通じて観測する過程を記述したものであり、観測結果（ここでは、観測ベクトルｙｐ２（ｎ））が、被観測量つまり入力（ここでは、状態ベクトルｘｐ２（ｎ））に依存して生成される様子を示している。
【０１２５】
次に、図１６を参照して、第２の実施の形態に係る動体の移動量計測装置２１０の作用について説明する。なお、図１６に示す画像処理ルーチンは、記憶部４０に制御プログラムとして格納されており、図示しないＣＰＵによって実行される。
【０１２６】
ステップ２００で、画像入力装置から入力インタフェース部２０を介して、所定フレームの観測画像を取り込み、主記憶部４１の取込画像メモリ４２に記憶する。また、表示メモリ４８の取込画像メモリ４９にも記憶し、出力インタフェース部８０を介して、画像出力装置に取込画像を表示する。なお、以下では、時刻ｎのフレームの画像情報をＬフレームと表記する。また、ここでは、時刻ｎの開始値を「０」とし、ｎ≧（ｋ−１）となる時刻のフレームから、以下の処理を行う。ｋは後述する処理フレーム数であり、本実施の形態では「３」である。ここでは、ｎ＝２から以下の処理を行う。
【０１２７】
次に、ステップ２０２で、位置変化量推定部２６１が、取込画像メモリ４２に記憶されている時刻ｎ−２、ｎ−１、ｎの連続するフレームの画像情報を２つずつ読み出し、各フレームの画像間の位置合わせを行い、連続するフレーム間における画像内の物体の位置変化量を推定する。画像の位置合わせには、ブロックマッチングやＬｕｃａｓ−Ｋａｎａｄｅ法といった既存の手法を用いることができる。本実施の形態では、図１７に示すようなブロックマッチングを用いる場合について説明する。
【０１２８】
ブロックマッチングとは、時刻（ｎ−１）の（Ｌ−１）フレームで設定した探索範囲内において、時刻ｎのＬフレームで類似する領域を探す手法であり、時刻（ｎ−１）の（Ｌ−１）フレームの小領域の画素値と時刻ｎのＬフレームで取り出した小領域の画素値との誤差の総和が最小となる領域を求め、対応する位置を探す手法である。画素値の誤差の総和を取る手段として、下記（７）式に示す各画素値の２乗距離の総和を取るＳＳＤ（Sum of Squared Distance）や、下記（８）式に示す各画素値の絶対値距離の総和を取るＳＡＤ（Sum of Absolute Distance）が挙げられる（非特許文献２「藤原  慎矢、田口  亮、“ブロックマッチング動き推定に基づくフレームレート向上に関する一手法”、IEICE Technical Report SIS2006-15, Jun. 2006」参照）。なお、（７）式及び（８）式におけるｘｉ，ｊは画素（ｉ，ｊ）の画素値である。
【０１２９】
【数６】
2015152975.tif 000008
【０１３０】
例えば、画像内の物体が左から右へ水平方向に移動していた場合には、時刻ｎのＬフレームにおいて、時刻（ｎ−１）の（Ｌ−１）フレームに設定した小領域の位置より右方向にずれた位置から取り出した小領域が、時刻（ｎ−１）の（Ｌ−１）フレームの小領域とマッチングすることになる。このときのずれ量を時刻（ｎ−１）の（Ｌ−１）フレームと時刻ｎのＬフレーム間の位置変化量として推定する。以下、画像内の物体が左から右へ水平方向に移動していた場合を例に説明する。
【０１３１】
次に、ステップ２０４で、位置変化量推定部２６１が、画像処理の処理単位となる処理領域のサイズを設定する。処理領域のサイズは、上記ステップ２０２で推定された位置変化量に応じて決定される。処理領域のサイズの設定は、ユーザにより操作部３０の領域指定部３２から指定された値を設定するようにしてもよいし、位置変化量と局所領域のサイズとの関係を予めテーブル等に定めておいて、推定された位置変化量に基づいて自動的に設定するようにしてもよい。
【０１３２】
また、位置変化量推定部２６１が、後述する状態方程式における時刻（ｎ−１）のシステム状態量として用いるフレーム数を設定する。フレーム数の設定は、ユーザにより操作部３０の処理フレーム指定部３３から指定された値を設定するようにしてもよいし、予め定めた値を設定するようにしてもよい。
【０１３３】
ここでは、ステップ２０４において、処理領域のサイズｐ×ｑ＝３×３画素、フレーム数３と設定されたものとする。この場合、図１８に示すように、左から右へ流れている、時刻（ｎ−２）の（Ｌ−２）フレーム、時刻（ｎ−１）の（Ｌ−１）フレーム、時刻ｎのＬフレームという時系列の３フレーム分の画像情報が処理対象となる。
【０１３４】
次に、ステップ１０６で、位置変化量推定部２６１が、上記ステップ２０４で設定されたサイズの処理領域を各フレームに対して同一の領域に設定する。
【０１３５】
次に、ステップ２０８で、システム状態量算出部２６２が、システム状態量算出処理を実行する。ここで、図１９を参照して、システム状態量算出処理ルーチンについて説明する。
【０１３６】
ステップ２０８０で、システム状態量算出部２６２が、時刻（ｎ−２）の（Ｌ−２）フレーム及び時刻（ｎ−１）の（Ｌ−１）フレームに設定された局所領域内の各画素のシステム状態量を要素とする時刻（ｎ−１）の状態ベクトルと、時刻ｎのＬフレームに設定された局所領域内の各画素のシステム状態量を要素とする時刻ｎの状態ベクトルとを導出する。以下、各時刻の状態ベクトルの導出方法を具体的に説明する。なお、以下では、状態変化の説明を分かり易くするため、人物が左から右へ歩行している様子を例に挙げて説明する。
【０１３７】
まず、時刻（ｎ−１）の状態ベクトルは、図２０に示すように、時刻（ｎ−１）の（Ｌ―１）フレーム及び時刻（ｎ−２）の（Ｌ―２）フレームにおけるｐ×ｑサイズの処理領域内の各画素のシステム状態量を組み合わせて構成される。すなわち、時刻（ｎ−１）のシステム状態ベクトルは１８次元のベクトルｘｐ２（ｎ−１）＝［ｘ１，１Ｌ−１（ｎ−１）  ｘ２，１Ｌ−１（ｎ−１）  ・・・  ｘ３，３Ｌ−１（ｎ−１）  ｘ１，１Ｌ−２（ｎ−１）  ｘ２，１Ｌ−２（ｎ−１）  ・・・  ｘ３，３Ｌ−２（ｎ−１）］Ｔとなる。次に、時刻ｎの状態ベクトルは、図２１に示すように、時刻ｎのＬフレームにおけるｐ×ｑサイズの処理領域内の各画素のシステム状態量を組み合わせて構成される。すなわち、時刻ｎの状態ベクトルは９次元のベクトルｘｐ２（ｎ）＝［ｘ１，１Ｌ（ｎ）  ｘ２，１Ｌ（ｎ）  ・・・  ｘ３，３Ｌ（ｎ）］Ｔとなる。ただし、Ｔはベクトルの転置を表す。また、状態ベクトルｘｐ２（ｎ−１）及びｘｐ２（ｎ）のサイズ（次元数）は設定されたフレーム数及び処理領域ｐ×ｑのサイズによって変化する。
【０１３８】
次に、ステップ２０８２で、システム状態量算出部２６２が、上記ステップ２０２で推定した位置変化量、及び時刻ｎのシステム状態量と時刻（ｎ−１）のシステム状態量との相関に基づく忘却係数を用いて状態遷移行列Φｐ１を構成する。
【０１３９】
具体的には、まず、画像内の物体が、時刻（ｎ−２）から時刻（ｎ−１）間、及び時刻（ｎ−１）から時刻ｎ間に、水平に左から右へ１画素ずつ移動していたとする。この場合、各フレームの処理領域において、図２２に示すように、同一のシステム状態量を持つ領域が存在する。ここで、同一のシステム状態量を持つ領域は、上記ステップ２０２で推定された位置変化量から判断できる。同一のシステム状態量を持つ領域は、図２２において網掛けされた画素領域である。図２３に、図２２の各フレームの処理領域部分を抜き出した図を示す。図２３に示すように、時刻（ｎ−２）の（Ｌ−２）フレームのシステム状態量ｘ１，１Ｌ−２，ｘ２，１Ｌ−２，ｘ３，１Ｌ−２と、時刻（ｎ−１）の（Ｌ−１）フレームのシステム状態量ｘ１，２Ｌ−１，ｘ２，２Ｌ−１，ｘ３，２Ｌ−１と、時刻ｎのＬフレームのシステム状態量ｘ１，３Ｌ，ｘ２，３Ｌ，ｘ３，３Ｌとが同一のシステム状態量を持つ領域となる。また、時刻（ｎ−１）の（Ｌ−１）フレームのシステム状態量ｘ１，１Ｌ−１，ｘ２，１Ｌ−１，ｘ３，１Ｌ−１と、時刻ｎのＬフレームのシステム状態量ｘ１，２Ｌ，ｘ２，２Ｌ，ｘ３，２Ｌとも同一のシステム状態量を持つ領域となる。
【０１４０】
また、時刻（ｎ−１）のシステム状態量と時刻ｎのシステム状態量とは相関を持ち、図２４に示すように、時刻（ｎ−１）の状態ベクトルを構成する時刻（ｎ−２）の（Ｌ−２）フレーム及び時刻（ｎ−１）の（Ｌ−１）フレームにおいて、時刻ｎのＬフレームとの時刻が近いほど強い相関を示し、遠いほど弱い相関を示す。そこで、時刻ｎのＬフレームからフレームが遠ざかるにつれて相関が低くなるようなパラメータ（忘却係数ρｉ，ｊ）を設定する。一般的に忘却係数は時刻ｎのＬフレームに近いほど相関が強まるように１に近い値となり、時刻ｎのＬフレームから遠ざかるほど相関が弱まるように０に近い値となるように構成する。ここで、忘却係数は推定して導出してもよいし、既知として扱ってもよい。
【０１４１】
そして、上述した同一のシステム状態量を持つ領域が時刻（ｎ−１）から時刻ｎへの変化によって一致し、かつ忘却係数を導入した状態遷移行列Φｐ２を構成する。
【０１４２】
次に、ステップ２０８４で、システム状態量算出部２６２が、駆動源ベクトルを構成する。図２３の時刻ｎのＬフレームのシステム状態量ｘ１，１Ｌ，ｘ２，１Ｌ，ｘ３，１Ｌは、時刻（ｎ−１）の情報からは構成することができないため、本実施の形態では、駆動源と呼ばれる領域にシステム状態量ｘ１，１Ｌ，ｘ２，１Ｌ，ｘ３，１Ｌを用いることで、駆動源ベクトルδｐ１（ｎ）を構成する。これは、物体がフレーム内を位置変化した際に、時刻（ｎ−１）のシステム状態量（（時刻（ｎ−２）のＬ−２）フレーム及び時刻（ｎ−１）の（Ｌ−１）フレーム）では表せられない領域のシステム状態量を、時刻ｎのシステム状態量で構成することを表しており、原画像の情報が含まれることから、駆動源ベクトルδｐ１（ｎ）は有色性駆動源となる。
【０１４３】
次に、ステップ２０８６で、システム状態量算出部２６２が、上記ステップ２０８０で導出した時刻（ｎ−１）及び時刻ｎの状態ベクトルｘｐ２（ｎ−１）及びｘｐ２（ｎ）、上記ステップ２０８２で構成した状態遷移行列Φｐ２、及び上記ステップ２０８４で構成した駆動源ベクトルδｐ２（ｎ）に基づいて、システム状態量を算出する。算出されたシステム状態量（状態方程式）の一例と、各ステップで構成された各情報との関係を図２５に示す。このように、本実施の形態における状態方程式は、クリアな画像からなる状態ベクトルｘｐ２（ｎ）と、０及び忘却係数ρｉ，ｊからなる状態遷移行列Φｐ２と、有色信号であるクリアな画像からなる駆動源ベクトルδｐ２（ｎ）とから構成されるため、所望のシステム状態量（エッジ抽出画像の画素情報）、つまり、クリアな画像情報（エッジ抽出画像情報）のみで構成されることになる。
【０１４４】
システム状態量算出処理ルーチンが終了すると、画像処理ルーチンにリターンして、ステップ１１０へ移行し、観測量算出部２６３が、第１の実施の形態と同様に、上記図７に示した観測画像モデルに基づいて、上記図８に示すように、エッジ抽出画像情報にぼけと雑音が加わって観測画像情報となる過程を、観測量（観測方程式）として算出する。
【０１４５】
なお、観測ベクトルｙｐ２（ｎ）は、観測画像情報を含む観測量として時刻ｎのＬフレームの処理領域内の各画素の画素情報を要素とする９次元のベクトルである。また、観測行列Ｍは、第１の実施の形態と同様に、エッジ抽出行列Ｍ−１の一般化逆行列により導出することができる。エッジ抽出行列Ｍ−１は、従来既知のエッジ抽出フィルタに基づいて予め定められるため、観測行列Ｍの各要素は、既知であって、予めデータ化して適切に設定されている。また、観測雑音ベクトルεｐ２（ｎ）は、処理領域に含まれる画素に対応する観測雑音を要素とする９次元のベクトルである。
【０１４６】
次に、ステップ１１２で、推定部２６４が、上記ステップ２０８で算出されたシステム状態量（状態方程式）及び上記ステップ１１０で算出された観測量（観測方程式）により、第１の実施の形態で示した有色性駆動源を含む状態空間モデルに基づく予測法を導出する。ステップ１１２は、上記図９に示した状態空間に基づく予測処理ルーチンによって実現される。なお、第２の実施の形態における状態空間モデルに基づく予測法においては、「時刻（ｎ）」は、上記第１の実施の形態における「時刻（ｎ）」には対応しない。第２の実施の形態における「時刻（ｎ）」は、複数フレームの各々におけるフレーム（Ｌ）の時刻に対応する。
【０１４７】
有色性駆動源を含む状態空間モデルに基づく予測法のアルゴリズムが終了すると、画像処理ルーチンへリターンして、ステップ１１４へ移行し、上記ステップ２０４で設定されたフレームの全処理領域について画像処理が終了したか否かを判定する。全処理領域について処理が終了していない場合には、ステップ１０６へ戻って、次の処理領域を設定して処理を繰り返す。全処理領域について処理が終了した場合には、ステップ１１６へ移行して、上記ステップ２００で取り込んだ全フレームについて処理が終了したか否かを判定する。全フレームについて処理が終了していない場合には、ステップ２１８へ移行してｎを１インクリメントして、ステップ２０２へ戻って、次のフレームについて処理を繰り返す。全フレームについて処理が終了した場合には、ステップ１２０へ移行する。
【０１４８】
以下、ステップ１２０〜１２８を実行して、画像処理ルーチンを終了する。
【０１４９】
以上説明したように、第２の実施の形態の動体の移動量計測装置２１０によれば、時刻ｎにおけるシステム状態量を、時刻（ｎ−１）までに取得されたフレームから得られる時刻（ｎ−１）のシステム状態量、推定された位置変化量に基づく時刻（ｎ−１）から時刻ｎへの状態遷移、及び時刻（ｎ−１）のシステム状態量に含まれない時刻ｎのシステム状態量の要素を含む有色性駆動源を用いて算出し、算出されたシステム状態量、及び算出された観測量に基づいて、状態空間モデルに基づく予測法により、時刻ｎのシステム状態量の最適値を、時刻ｎのエッジ抽出画像として推定し、推定原画像情報上にユーザによって指定された箇所について、フレーム間で当該箇所を追跡し、フレーム間各々での当該箇所の移動量を計測することにより、腸壁上の点の微小な移動量を計測することができ、可視化することができる。すなわち、ピクセル単位で腸壁の変動を追跡することができる。
【０１５０】
また、状態遷移行列を構成する際に、時刻ｎからの時間間隔に応じた忘却係数を導入することで、状態遷移行列をより適切に表現することができ、より精度良くエッジ抽出画像を得ることができる。
【０１５１】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１５２】
例えば、第１及び第２の実施の形態では、追跡波形処理部７０が、フレームＬのエッジ抽出画像上において、フレーム（Ｌ−１）の追跡箇所の座標の上方向に存在するエッジ上の箇所までの移動量と、フレーム（Ｌ−１）の追跡箇所の座標の下方向に存在するエッジ上の箇所までの移動量とを算出し、移動量が少ない方向に存在するエッジ上の箇所を、フレームＬにおける追跡点と判定する場合を例に説明したが、これに限定されるものではない。腸壁の移動方向が左右方向である場合には、追跡波形処理部７０が、フレームＬのエッジ抽出画像上において、フレーム（Ｌ−１）の追跡箇所の座標の右方向に存在するエッジ上の箇所までの移動量と、フレーム（Ｌ−１）の追跡箇所の座標の左方向に存在するエッジ上の点までの移動量とを算出し、移動量が少ない方向に存在するエッジ上の箇所を、フレームＬにおける追跡箇所と判定する。
【０１５３】
また、第２の実施の形態では、画像内の物体が水平方向に移動する場合を例に説明したがこれに限定されない。特に、位置変化量推定部２６１で推定された位置の変化に応じたシステム状態量を算出すればよい。例えば、図２６に示すように、画像内の物体が上方へ移動する場合でも、システム状態量が重なる領域を考慮した状態遷移行列Φｐ１及び駆動源ベクトルδｐ１を構成するとよい。
【０１５４】
また、第２の実施の形態では、忘却係数として時刻変化における相関値の変化に基づく値を用いる場合について説明したが、時刻変化によらないフレーム間の相関値に基づく値を用いるようにしてもよい。
【０１５５】
なお、上記実施の形態では、ＣＰＵによりプログラムを実行することで各部の処理を実行する場合について説明したが、各部をハードウエアにより構成するようにしてもよい。
【０１５６】
また、上記実施の形態では、医療画像として撮像された腸の動画像を対象とし、腸壁の動きを追跡する動体の移動量計測装置に本発明を適用した場合を例に説明したが、これに限定されるものではなく、他の臓器の医療画像に対しても適用可能である。
【０１５７】
また、上記実施の形態は、３次元表示するための複数枚の画像の時系列データに対しても応用可能である。
【０１５８】
なお、上記実施の形態では、医療画像として撮像された腸の動画像を対象とし、腸壁の動きを追跡する動体の移動量計測装置に本発明を適用した場合を例に説明したが、これに限定されるものではない。例えば、本発明は以下の分野に適用させることができる。
【０１５９】
（１）医療分野
本発明は医療分野に適用させることもでき、例えば、入院患者の行動検知システムに本発明を適用させることができる。
例えば、図２７に示すように、ベッドに寝ていた入院患者が、急な行動（起き上がる行動等）をした状態、もしくは入院患者がベッドから落ちそうな状態又はベッドから落ちてしまった状態を検知し、医療従事者に知らせるような行動検知システムに本発明を適用させることができる。
この場合には、観測対象である入院患者の輪郭（エッジ）上の着目箇所を追跡箇所として設定し、予め定められた輪郭の移動方向に当該追跡箇所を追跡することにより、入院患者の行動を検知することができる。これにより、画像情報の全てを医療従事者に知らせるのではなく、入院患者の行動の異常を検知して医療従事者に知らせることができる。
【０１６０】
また、例えば、新生児の手足の動きを撮影計測することにより脳の発達障害状況を把握すること、あるいは、不整脈症状が発症している患者に対し、心臓を動きを撮影し心臓の各部位の移動量を計測することにより発症部位の特定把握すること、更には、小腸の異常発症に対し、小腸を撮影し小腸の移動量を計測することにより、発症している小腸部位を特定把握することができる。
【０１６１】
（２）安全・安心分野
本発明は安全・安心分野に適用させることもでき、例えば、防犯システムや災害システムに本発明を適用させることができる。
例えば、図２８に示すように、人が防犯カメラの映像を監視し続けるのには無理があるため、防犯カメラの映像内で異常があった場合にのみ知らせるような防犯システムに適用させることができる。
また、災害発生時の初期行動において、止まっている人や何度も往復する人、又は避難方向と逆方向に向かう人等を検知し、そのような人々に対し、正しい避難方向へ誘導するような災害システムに本発明を適用させることができる。
この場合には、観測対象者の輪郭（エッジ）上の着目箇所を追跡箇所として設定し、予め定められた輪郭の移動方向に当該追跡箇所を追跡することにより、異常等を検知することができる。
【０１６２】
（３）植物の生産・管理分野
本発明は植物の生産・管理分野に適用させることもでき、例えば、植物の成長を監視するシステムに適用させることもできる。
例えば、図２９に示すように、観測対象植物（野菜）の輪郭（エッジ）上の着目箇所を追跡箇所として設定し、予め定められた輪郭の移動方向に追跡箇所を追跡することにより、植物の成長度合や、異常の検知をすることができる。
【符号の説明】
【０１６３】
１０，２１０  動体の移動量計測装置
２０   入力インタフェース部
３０   操作部
３１   パラメータ設定部
３２   領域指定部
３３   処理フレーム指定部
３４   画像情報設定部
３５   着目点指定部
４０   記憶部
４１   主記憶部
４２   取込画像メモリ
４３   追跡波形保存メモリ
４４   被写体領域メモリ
４５   処理パラメータメモリ
４６   波形考察結果メモリ
４７   補助記憶部
４８   表示メモリ
４９   取込画像メモリ
５０   処理メモリ
５１   処理画像メモリ
６０，２６０  画像処理部
６１   領域設定部
６２，２６２  システム状態量算出部
６３，２６３  観測量算出部
６４，２６４  推定部
６５   エッジ抽出画像出力部
７０   追跡波形処理部
７５   追跡波形解析部
８０   出力インタフェース部
【図１】
2015152975.tif 000009
【図２】
2015152975.tif 000010
【図４】
2015152975.tif 000011
【図６】
2015152975.tif 000012
【図９】
2015152975.tif 000013
【図１０】
2015152975.tif 000014
【図１５】
2015152975.tif 000015
【図１６】
2015152975.tif 000016
【図１９】
2015152975.tif 000017
【図３】
2015152975.tif 000018
【図５】
2015152975.tif 000019
【図７】
2015152975.tif 000020
【図８】
2015152975.tif 000021
【図１１】
2015152975.tif 000022
【図１２】
2015152975.tif 000023
【図１３】
2015152975.tif 000024
【図１４】
2015152975.tif 000025
【図１７】
2015152975.tif 000026
【図１８】
2015152975.tif 000027
【図２０】
2015152975.tif 000028
【図２１】
2015152975.tif 000029
【図２２】
2015152975.tif 000030
【図２３】
2015152975.tif 000031
【図２４】
2015152975.tif 000032
【図２５】
2015152975.tif 000033
【図２６】
2015152975.tif 000034
【図２７】
2015152975.tif 000035
【図２８】
2015152975.tif 000036
【図２９】
2015152975.tif 000037