JP5874996B2 - 物体運動推定装置、物体運動推定方法及びプログラム - Google Patents

Description

本発明は、物体運動推定装置、物体運動推定方法及びプログラムに関し、特に、画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する物体運動推定装置等に関する。

オプティカルフローは、画像における物体の動作を表すベクトルである。オプティカルフロー方程式は、式(1)の保存則により表される（非特許文献１参照）。ここで、Ｉ(x,y,t)は輝度であり、(x,y)は画像上の座標を示し、ｔは時間である。そして、[ｕ ｖ]^Tは、オプティカルフローを表す。

変数ｕ及びｖを決定するには、拘束条件として、式(1)のみでは十分でない。そのため、従来、オプティカルフローを計算するために、付加的な拘束条件を導入することがなされてきた。

付加的な拘束条件を表す一つのアプローチは、すべてのオプティカルフローは同じと仮定して、(x,y)の近傍に空間的に分布する複数の点を考えることである（非特許文献２参照）。オプティカルフローは、各点に対する保存則を構成する連立代数方程式を解くことによって計算される。このアプローチの基本的な推定原理は、式(2)により与えられる線形方程式によって表すことができる（非特許文献３及び４参照）。ここで、添字は、偏微分を表す。すなわち、例えば、I_tx＝∂²Ｉ／∂t∂xである。輝度のヘッセ行列（式(2)の行列の下部分）は、画像における物体の角に深く関係している（非特許文献５参照）。このヘッセ行列の行列式（線形方程式の解法の数値的正確性に関係する）も、推定されたオプティカルフローの信頼性の度合いとして使用される（非特許文献６参照）。

拘束条件を表す他のアプローチは、勾配法（大域法）と呼ばれるものである。これは、付加的な拘束条件として、オプティカルフローの滑らかさを導入するものである（非特許文献１、５、７及び８参照）。勾配法では、ヘッセ行列の正則性は失われ、オプティカルフロー推定の信頼性が低下する画像中の領域が存在する。このような画像中の領域に対しては、信頼できる領域のオプティカルフローを使って補間によって求める。

勾配法の関連研究には、輝度の一定性に代えて、勾配（すなわち、式(2)の第２及び第３行）の一定性を使うものもある（非特許文献８参照）。この修正により、結果として、推定法において、照明変動に対して、より頑健になる。勾配の一定性は、同様に、推定された領域において、並進的な動きを仮定する。

他方、カルマンフィルタは、観測出力の時系列データから動的システムの状態を逐次的に推定する。さらに、カルマンフィルタを使用することにより、結果として、時間解像度を犠牲にするものの、時間の導関数の明示的な計算を必要としない。それゆえ、カルマンフィルタを用いることで、より頑健でかつ正確なオプティカルフロー推定を構築することができる。さらに、近年、時間的な逐次処理は、動画像（ムービー）のような時系列データを扱う場合に、より重要なものとなってきた。カルマンフィルタは、潜在的に、時系列データの処理に対して、強力な方法になりうる。

オプティカルフロー推定にカルマンフィルタを使う試みは、様々になされてきた（非特許文献９−１２参照）。これらの研究は、すべて、共通のフレームワークがある。すなわち、
(i)輝度の時空間勾配（すなわち、Ｉ_x、Ｉ_y、Ｉ_t）を計算する。
(ii)線形方程式又は最小化問題を解くことにより、オプティカルフローを推定する。
(iii)スムーザーとして、カルマンフィルタを適用する。

第３のステップ(iii)において、オプティカルフローの動的モデルは、式(3)として与えられる。ここで、ｕ_(k)及びｖ_(k)は、時刻ｋにおけるオプティカルフローベクトルを示し、Ｎ(0,s)は、標準偏差ｓのガウス性白色ノイズを示す。

発明者らは、オプティカルフロー方程式を、オプティカルフローを推定するためにカルマンフィルタを適用する動的システムとして考察した（特許文献１、非特許文献１３参照）。これにより、他の従来技術と比較して、画像の輝度とオプティカルフローの振る舞いをシミュレートする動的システムモデルを提案した。オプティカルフローを、式(2)のような線形方程式を解くことなく、カルマンフィルタ処理中に推定することができる。これは、Lucas-Kanade法（非特許文献２参照）やHorn-Schunck法（非特許文献１参照）とは異なり、オプティカルフローの並進性の仮定も、オプティカルフローの滑らかさに対する人工的なペナルティパラメータも必要としない。

世界公開第２０１０／０９８１６０号パンフレット

B.K.P.Horn、外１名，"Determining Optical Flow,"Artificial Intelligence，vol.17，pp.185-203，1981. B.D.Lucas、外１名，"An iterative image registration technique with an application to stereo vision,"in IJCAI(P.J.Hayes，ed.)，pp.674-679，William Kaufmann，1981. S.Uras、外３名,"A computational approach to motion perception,"Biological Cybernetics，vol.60，pp.79-87，1988. A.Verri、外２名,"Differential techniques for optical flow,"Journal of optical Society of America A，vol.7，no.5，pp.912-922，1990. H.H.Nagel、外１名,"An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences,"IEEE transactions on Pattern Analysis and Machine Intelligence，vol.8，pp.565-593，1986. N.Ohta,"Image movement detection with reliability indices,"IEICE Transactions，vol.E74，no.10，pp.3379-3388，1991. T.Brox、外３名,"High accuracy optical flow estimation based on a theory for warping,"in Computer Vision−ECCV 2004(T.Pajdla and J.Matas,eds.)，vol.3024 of Lecture Notes in Computer Science，pp.205-216，Berlin：Springer，2004. J.Weickert、外３名,"A survey on variational optic flow methods for small displacement,"Tech.Rep.152，Universit¨at des Saarlandes，2005. A.Singh,"Incremental estimation of image flow using a Kalman filter,"Journal of Visual Communication and Image Representation，vol.3，no.1，pp.39-57，1992. T．Chin、外２名,"Probabilistic and sequential computation of optical flow using temporal coherence,"Image Processing，IEEE Transactions on，vol.3，no.6，pp.773-788，1994. J.Cooper、外１名,"Uncertainty in the optical flow using the Kalman filter,"in Proceedings.Fourth international conference on Control，Automation，Robotics and Vision,(Westin Stamford，Singapore)，pp.1681-1685，1996. M.Elad、外１名,"Recursive optical flow estimation-adaptive filtering approach,"Journal of Visual Communication and Image Representation，vol.9，no.2，pp.119-138，1998. N.Sebe、外２名,"Dynamical system approach to optical flow estimation,"39th SICE Symposium on Control Theory，pp.357-362，2009. M.Otte、外１名,"Estimation of optical flow based on higher-order spatiotemporal derivatives in interlaced and non-interlaced image sequences,"Artificial Intelligence，vol.78，pp.5-43，1995. S．Baker、外５名,"A database and evaluation methodology for optical flow,"International Journal of Computer Vision，vol.92，no.1，pp.1-31，2011.

付加的な拘束条件を表すアプローチ（非特許文献１−８参照）は、ほとんどが、輝度の時間的かつ空間的な導関数の明示的計算が必要であり、画像データから導関数を得るために、差分の計算が不可避である。しかしながら、このような差分の計算は、外乱に対して弱い。この外乱は、推定の信頼性に重大な影響を与える。

また、従来のオプティカルフロー推定にカルマンフィルタを使う試み（非特許文献９−１２参照）は、式(3)の動的モデルを採用し、オプティカルフロー推定後のローパスフィルタとしてカルマンフィルタを適用した。しかしながら、これらの研究は、オプティカルフロー方程式の動特性を扱っていない。オプティカルフロー方程式の動特性は、輝度の時間変化と空間勾配を結びつける、本質的な特性である。

発明者らの提案（特許文献１、非特許文献１３参照）は、回転及び／又は拡大を伴う画像の系列に対するオプティカルフロー推定に適用可能ではある。しかしながら、輝度とオプティカルフローのより高次の勾配推定が必要である。より高次の勾配に対する共分散行列は、正確にモデル化し、設定することは難しい。さらに、このモデルの観測変数は、基本的に、サイズ３×３ピクセルの領域の輝度である。

ゆえに、本発明は、第１の目的として、従来のオプティカルフロー推定法と比較して、推定精度を向上させ、かつ、例えば照明変動などに対して頑健な物体運動推定装置等を提供する。

さらに、物体運動推定の立場からは、推定されたオプティカルフローから、複数の異なる運動をする物体の分離（つまり境界の検出）も重要である。本発明は、第２の目的として、異なる運動物体の境界等も同時に検出可能な物体運動推定装置等を提供する。

本願発明の第１の観点は、画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する物体運動推定装置であって、前記動的システムモデルの状態は、前記画像データの推定領域における各ピクセルの輝度、及び、前記各ピクセルの全部又は一部でのオプティカルフローであり、前記動的システムモデルの前記状態を記憶する状態記憶手段と、前記状態記憶手段に記憶された前記状態における前記輝度及び前記オプティカルフローから、新たな時刻の前記状態における前記輝度を演算して、前記状態記憶手段に記憶する輝度時間発展手段と、少なくとも前記状態記憶手段に記憶された前記状態における前記オプティカルフローから、新たな時刻の前記状態における前記オプティカルフローを演算して、前記状態記憶手段に記憶するオプティカルフロー時間発展手段を備えるものである。

本願発明の第２の観点は、第１の観点において、前記動的システムモデルは、ノイズを含み、前記輝度時間発展手段及び／又は前記オプティカルフロー時間発展手段は、カルマンフィルタを用いて演算するものである。

本願発明の第３の観点は、第１又は第２の観点において、前記各ピクセルでの前記オプティカルフローには、前記状態に含まれるものと前記状態に含まれないものが存在し、前記輝度情報発展手段は、前記状態記憶手段に記憶された前記状態に含まれるオプティカルフローを補間して前記状態に含まれないオプティカルフローを求めて、前記状態記憶手段に記憶された前記状態における前記輝度及び前記オプティカルフロー、並びに、前記補間して得られた前記状態に含まれないオプティカルフローから、新たな時刻の前記状態における前記輝度を演算するものであり、前記オプティカルフロー時間発展手段は、少なくとも前記状態記憶手段に記憶された前記状態に含まれる前記オプティカルフローから、新たな時刻の前記状態に含まれる前記オプティカルフローを演算するものである。

本願発明の第４の観点は、第１から第３のいずれかの観点において、前記動的システムモデルの観測変数と予測推定値との観測残差から、前記オプティカルフローの不連続性を判断する不連続検出手段を備えるものである。

本願発明の第５の観点は、物体運動推定装置において、画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する物体運動推定方法であって、前記物体運動推定装置は、状態記憶手段と、時間発展手段を備え、前記状態記憶手段は、前記動的システムモデルの状態を記憶するものであり、前記動的システムモデルの状態は、前記画像データの推定領域における各ピクセルの輝度、及び、前記各ピクセルの全部又は一部での前記オプティカルフローであり、前記時間発展手段が、前記状態記憶手段に記憶された前記状態における前記輝度及び前記オプティカルフローから、新たな時刻の前記状態における前記輝度を演算して前記状態記憶手段に記憶し、少なくとも前記状態記憶手段に記憶された前記状態における前記オプティカルフローから、新たな時刻の前記状態における前記オプティカルフローを演算して前記状態記憶手段に記憶する時間発展ステップを含むものである。

本願発明の第６の観点は、コンピュータにおいて第５の観点の物体運動推定方法を実現するためのプログラムである。

なお、本願発明を、第６の観点のプログラムを（定常的に）記録するコンピュータ読み取り可能な記録媒体として捉えてもよい。

また、本願発明において、輝度時間発展手段が、状態における輝度の空間差分を近似する空間差分オペレータを用いて、新たな時刻の状態における輝度を演算するものであってもよい。

本願発明によれば、第１の目的について、輝度の観測量（画像データにおける輝度）を直接使用するのではなく、推定された輝度を使用することにより、頑健なオプティカルフロー推定が行うことが可能になる。すなわち、本願発明によれば、画像データの系列とは別に、推定領域における多くの画素の輝度及びオプティカルフローの情報を利用して、新たな時刻における各画素の輝度を計算する。よって、動的システムモデルの状態として、推定領域における多くの画素の輝度とオプティカルフローを利用することにより、推定精度を向上させることができる。さらに、例えば照明変動などは、画像データ全体に現れる。そのため、推定領域における多くの画素の情報を利用することにより、頑健な推定を実現することができる。

すなわち、発明者らが提案した手法（特許文献１、非特許文献１３参照）と同様に、本願発明は、各ピクセルにおけるオプティカルフロー方程式を動的システムとして考える。本願発明は、並進性の過程を要求しない。そして、背景技術と比較して、本願発明は、任意の領域における各ピクセルでのオプティカルフロー推定に適用可能である。そして、特に照明変動に対して、観測された領域を拡張することによって、頑健な推定を成し遂げることができる。

さらに、本願発明は、輝度とオプティカルフローのより高次の勾配推定を要しなくすることができる。よって、本願発明によれば、推定精度を決定する可調整パラメータである共分散行列を調整するための労力は、軽減される。

さらに、他の技術分野における研究結果を利用することにより、潜在的に、推定性能を改善することができうる。例えば、流体解析での移流方程式の差分近似、システム制御での故障検出などである。

さらに、本願発明の第２の観点にあるように、カルマンフィルタを利用して予測推定誤差共分散を得ることによって、推定したオプティカルフローの信頼性を、通常は知りえない真値が得られなくても評価することが可能になる。

さらに、本願発明の第３の観点にあるように、オプティカルフローがゼロに近い場合でも、疑似稠密推定モデルを使用することにより、原理的に推定を行うことが可能になる。

さらに、本願発明の第４の観点にあるように、第２の目的について、観測残差を利用して不連続点を検出することにより、例えば異なる運動物体の分離に役立つ、物体の境界などを検出することが可能になる。

本願発明の実施例である物体運動推定装置の概略ブロック図である。 第１の例における入力画像を示す図である。 第１の例において推定された(120,200)でのオプティカルフローを示す図である。 補間オペレータＭ_uを示す図である。 第２の例のオリジナル画像を示す図である。 第２の例のGround truth flowを示す図である。 第２の例の推定されたオプティカルフロー[^q_u(k|k) ^q_v(k|k)]^T（ｋ＝５）を示す図である。 第２の例の共分散行列の対角要素Ｐ_u(k|k)＋Ｐ_v(k|k)（ｋ＝５）を示す図である。 第２の例の測定の二乗推定誤差|ｚ(k)−Ｈ^q_(k|k)|²（ｋ＝５）を示す図である。 フローのカラーコーディングを示す図である。

以下、図面を参照して、本願発明の実施例について説明する。ただし、本願発明は、記載した実施例に限定されるものではない。

図１は、本願発明の実施例である物体運動推定装置１の概略ブロック図である。物体運動推定装置１は、画像データ記憶部３と、状態記憶部５（本願発明の「状態記憶手段」の一例）と、時間発展部７（本願発明の「時間発展手段」の一例）と、不連続検出部９（本願発明の「不連続検出手段」の一例）を備える。時間発展部７は、輝度時間発展部１１（本願発明の「輝度時間発展手段」の一例）と、オプティカルフロー時間発展部１３（本願発明の「オプティカルフロー時間発展手段」の一例）を備える。

物体運動推定装置１は、画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する。画像データ記憶部３は、画像データの系列を記憶する。状態記憶部５は、動的システムモデルの状態を記憶する。

物体運動推定装置１の動的システムモデルの一例について、稠密推定の場合を説明する。簡単のため、推定領域（画像データの一部又は全部の領域で、複数のピクセルを含む。本願発明の「推定領域」の一例である。）は、幅と高さがそれぞれＭとＮ（ピクセル）の長方形の領域であると仮定する。

動的システムモデルの時刻ｋにおける状態ｑ_(k)は、式(4)で定義される。ここで、ｑ_I(k)、ｑ_u(k)及びｑ_v(k)は、それぞれ、式(5a)、(5b)及び(5c)である。この例では、状態は、推定領域の各ピクセルの輝度及びオプティカルフローである。状態記憶部５は、少なくとも、時刻ｋにおける状態ｑ_(k)を記憶している。

この例での動的システムモデルは、式(6)で表される（以下、「稠密推定モデル」ともいう。）。ここで、η_I(k+1)、η_u(k+1)及びη_v(k+1)はシステムノイズ（共分散行列Ｑ）であり、ζ_(k)は、観測ノイズ（共分散行列Ｒ）である。

式(7a)及び(7b)にあるように、ｑ_Iの空間微分ｑ_Ix及びｑ_Iyを近似する空間差分オペレータを、それぞれ、Ｍ_x及びＭ_yとする。

輝度時間発展部１１は、式(8a)に従い、時刻ｋにおける状態の輝度ｑ_I(k)及びオプティカルフローｑ_u(k)及びｑ_v(k)から、時刻ｋ＋１における状態の輝度ｑ_I(k+1)を演算する。ただし、η_I(k)=0とする。演算結果は、状態記憶部５に記憶する。

オプティカルフロー時間発展部１３は、式(8b)及び(8c)に従い、時刻ｋにおける状態のオプティカルフローｑ_u(k)及びｑ_v(k)から、時刻ｋ＋１における状態のオプティカルフローｑ_u(k+1)及びｑ_v(k+1)を演算する。ただし、η_u(k)，η_v(k)は0とする。演算結果は、状態記憶部５に記憶する。

不連続検出部９は、予測推定値^ｑ_(k+1|k)が得られることから、故障検出技術を利用して、観測残差ｚ_(k)−^ｑ_I(k|k-1)を用いてオプティカルフローの不連続点を検出し、例えば物体の境界等を検出する。観測残差が大きいことは、Ｉ_t＝−（Ｉ_xｕ＋Ｉ_yｖ）が間違っていることを意味する。これは、例えばある物体が別の物体に隠れた場合（隠れ）などに起因する。そのため、オプティカルフローの不連続点を検出することにより、物体の境界等を検出することが可能になる。

さらに、輝度時間発展部１１及びオプティカルフロー時間発展部１３は、例えば、カルマンフィルタを使用することにより、濾過推定値^ｑ_u(k|k)、^ｑ_v(k|k)だけでなく、濾過推定誤差共分散^Ｐ_(k|k)も得られる。そのため、オプティカルフローについて、真値がなくても、推定精度の評価が可能になる。

物体運動推定装置１の動的システムモデルの他の一例について、疑似稠密推定の場合を説明する。簡単のため、推定領域は、同様に、幅と高さがそれぞれＭとＮ（ピクセル）の長方形の領域であると仮定する。

疑似稠密推定では、推定領域内の点のうち、一部のもののみ状態とする。間引いた推定点のオプティカルフローの状態を、~ｑ_u(k)及び~ｑ_v(k)と表す。

疑似稠密推定において、動的システムモデルの時刻ｋにおける状態ｑ_(k)は、式(9)で定義される。この例では、状態は、推定領域の各ピクセルの輝度、及び、推定領域の一部のピクセルのオプティカルフローである。状態記憶部５は、少なくとも、時刻ｋにおける状態ｑ_(k)を記憶している。

この例での動的システムモデルは、式(10)で表される（以下、「疑似稠密推定モデル」ともいう。）。ここで、η_I(k+1)、η_u(k+1)及びη_v(k+1)はシステムノイズ（共分散行列Ｑ）であり、ζ_(k)は、観測ノイズ（共分散行列Ｒ）である。

オプティカルフロー時間発展部１３は、式(11a)及び(11b)にあるように、補間オペレータＭ_u及びＭ_vを使って、すべてのピクセルでのオプティカルフローｑ_u(k)及びｑ_v(k)を求める。

輝度時間発展部１１は、式(12a)に従い、時刻ｋにおける状態の輝度ｑ_I(k)及びオプティカルフローｑ_u(k)及びｑ_v(k)から、時刻ｋ＋１における状態の輝度ｑ_I(k+1)を演算する。ただし、η_I(k)=0とする。演算結果は、状態記憶部５に記憶する。

オプティカルフロー時間発展部１３は、式(12b)及び(12c)に従い、時刻ｋにおける状態のオプティカルフロー~ｑ_u(k)及び~ｑ_v(k)から、時刻ｋ＋１における状態のオプティカルフロー~ｑ_u(k+1)及び~ｑ_v(k+1)を演算する。ただし、η_u(k)，η_v(k)は0とする。演算結果は、状態記憶部５に記憶する。

不連続検出部９は、同様に、故障検出技術を利用して、観測残差を用いてオプティカルフローの不連続点を検出し、例えば境界等を検出する。

さらに、輝度時間発展部１１及びオプティカルフロー時間発展部１３は、例えばカルマンフィルタを使用することにより、真値を知らなくても推定精度の評価が可能になる。

発明者らは、動的システムモデルに基づくオプティカルフローの点推定を提案した（特許文献１及び非特許文献１３参照）（以下、「点推定モデル」という。）。点推定モデルでは、各ピクセルに関してオプティカルフローに使用されたとき、推定処理は、他のピクセルから独立である。点推定モデルの推定の変数は、３×３領域の輝度である。よって、これらは、近傍のピクセルとオーバーラップする。しかしながら、近傍のピクセルの推定処理を併合するフレームワークがない。

本実施例では、稠密推定モデルによるオプティカルフロー推定について説明する。この稠密推定モデルは、推定領域の各ピクセルの輝度とオプティカルフローを、一つの動的システムとして考える。輝度の時間更新は、オプティカルフローによって定義される。このモデルの観測変数は、推定領域の各ピクセルの輝度である。

輝度の偏微分（例えば、Ｉ_x、Ｉ_y)に差分近似を適用したとき、このモデルは、輝度とオプティカルフローの勾配推定を必要としない。それゆえ、輝度とオプティカルフローの空間勾配は、点推定モデルでは状態を構成するものではあったが、稠密推定モデルの状態には含まれない。

稠密推定モデルの詳細について、説明する。簡単のため、推定領域は、幅と高さがそれぞれＭとＮ（ピクセル）の長方形の領域であると仮定する。整数Ｌは、オプティカルフローが推定されるピクセルの数を示す。稠密推定モデルでは、オプティカルフローは、各ピクセルで推定される。すなわち、Ｌ＝ＭＮである。

式(13a)及び(13b)で与えられる非線形離散時間システムを考える。推定領域における各ピクセルの輝度Ｉとオプティカルフロー[ｕ ｖ]^Tに関して、状態ｑ_(k)は、式(14)で定義される。ここで、ｑ_I(k)，ｑ_u(k)及びｑ_v(k)は、それぞれ、式(15a)、(15b)及び(15c)である。

各ピクセルの輝度の空間勾配は、式(16a)、(16b)及び(16c)で定義される。ｑ_Ix及びｑ_Iyは、式(17a)及び(17b)で近似される。ここで、Ｍ_xとＭ_yは、離散偏微分オペレータ（すなわち、偏微分オペレータを離散近似する。）を示す。これらのオペレータの具体例については、後に説明する。

状態空間方程式(1)を前進差分近似で離散化すると、差分方程式は式(18)で与えられる。式(18)の右辺は、式(13)のf(q_(k))の具体的な表現の一例である。式(18)は、移流方程式の近似モデルである。よって、他にも、時間微分近似の手法が多数存在する。

オプティカルフローの時間更新が、式(19a)及び(19b)として与えられるとする。式(19)は、暗に、オプティカルフローが時不変であることを仮定する。カルマンフィルタは、この仮定にもかかわらず、時変オプティカルフローを推定することができる。

観測値z_(k)は、式(20)のように、推定領域の各ピクセルの輝度とする。よって、式(13b)の行列Ｈは、式(21)で与えられる。ここで、Ｅ_m×nは、m×nの単位行列を表し、O_m×nは、m×nの零行列を表す。

続いて、空間勾配の近似について説明する。輝度の空間勾配Ｉ_x及びＩ_yは、式(17)によって近似される。

数値流体運動の分野では、風上差分スキームがよく使用されている。しかし、私たちは、未知のフローを推定することを目的とする。よって、式(22)により与えられる２次の中心差分スキームを適用する。変数Ｉ_y(x,y)も、同様に定義される。行列Ｍ_x及びＭ_yは、これらの近似を使うことで容易に計算することができる。

偏差分オペレータ(22)の主要部が２次の中心差分近似であるとしても、他の差分スキーム（例えば、風上差分スキーム）を、中心差分近似の代わりに適用することができる。同様に、画像処理の分野でよく用いられる、ソーベルフィルタや他の空間の偏差分フィルタを適用することができる。行列Ｍ_x及びＭ_yは、このような空間的差分フィルタを使っても定義される。

続いて、稠密推定モデルの局所的可観測性について説明し、オプティカルフロー推定の可能性を検討する。

オプティカルフロー推定の可能性は、推定モデルの可観測性を解析することで、簡単に決定することができる（特許文献１及び非特許文献１３参照）。Ｍ_x及びＭ_yとして、通常の偏微分オペレータを選択すれば、稠密推定モデル(13)は、可観測指数は３であり、局所可観測である。よって、一般的な状態に対して、観測された輝度データから、オプティカルフローを推定することができる。

他方、非線形システムが局所可観測だとしても、オプティカルフローは、すべての状態（すなわち、画像の輝度と物体の動きのすべての組み合わせ）に対しては推定できるとは限らない。特別な状態では、システムが不可観測になり、このような場合、オプティカルフローは、観測値を使って推定することはできない。可観測性解析は、このような特別な状態を同定できる。

瞬時推定に対する正確な線形方程式は、推定領域のサイズと選択された偏微分オペレータに依存する。しかし、離散化する前の連続時間モデルにおける本質的な推定方程式は、式(23)で与えられる。各ピクセルにおいて、推定方程式(23)は、オプティカルフローの瞬時推定に使用される。オプティカルフローを一意に決定するには、輝度の空間勾配ベクトル[Ｉ_x Ｉ_y]^Tとその時間微分[Ｉ_tx Ｉ_ty]^Tの線形独立性が必要である。一般には、これらのベクトルは、線形独立であり、オプティカルフローは一意に決定することができる。しかし、静的な（動きのない）状態（すなわち、[Ｉ_tx Ｉ_ty]^T＝Ｏ_1×2）は、輝度の空間勾配ベクトル[Ｉ_x Ｉ_y]^Tとその時間微分[Ｉ_tx Ｉ_ty]^Tが線形従属になる特別な場合である。よって、稠密推定モデル(13)は、静的状態を決定するためには使うことができない。

静的状態において、稠密推定モデルの局所的可観測性が失われる本質的な理由は、以下のものである。
(i)オプティカルフロー方程式(1)は、輝度変数とオプティカルフロー変数に関して双線形である。
(ii)観測出力は、輝度データのみで構成される。

オプティカルフロー方程式の双線形構造は、オプティカルフロー方程式のより高次の時間微分でも維持される。(ii)より、推定方程式の形式から、左辺は輝度の時間微分を表し、右辺は輝度及びオプティカルフローの時間微分の積を表す。推定方程式の右辺の係数（２より大きい可観測指数に対応する。）は、輝度の時間微分からなり、静的状態ではゼロになる。よって、輝度のより高次の微分を使って、静的な状態を決定することができない。静的な状態を決定するために、可観測指数が２の推定モデルが必要になる。

本実施例では、疑似稠密推定モデルによるオプティカルフロー推定について説明する。

観測式(13b)と(21)の特別な構造が示すように、可観測指数が２であるモデルを得るためには、推定される変数の数は、観測変数の数よりも小さくなければならない。この要求を満たすため、オプティカルフローの空間解像度を犠牲にする。オプティカルフローを推定するピクセルの数を減らすため、すべてのピクセルのオプティカルフローを推定することはしない。オプティカルフローが推定されないピクセルのオプティカルフローは、他のピクセルの値を補間して推定する。これらのピクセルの状態空間方程式は、補間されたオプティカルフローを使ったオプティカルフロー方程式によって定義される。

推定領域は、同様に、幅と高さがそれぞれＭとＮ（ピクセル）の長方形の領域を仮定する。式(13)と同様のシステムを考える。状態q_(k)は、式(24)及び(25)で定義される。ベクトル~ｑ_uと~ｑ_vは、間引かれた推定点におけるオプティカルフローを示す。式(15b)と(15c)で定義されるオプティカルフローベクトルｑ_uとｑ_vは、式(26a)及び(26b)で与えられるように~ｑ_uと~ｑ_vから補間によって生成される。ここで、Ｍ_uとＭ_vは、補間オペレータを示す。間引きと補間オペレータの明示的な例は、後述する。推定エリアにおけるオプティカルフローを推定するピクセルの数を、Ｌで定義する。可観測指数を２に減らすために、式(27)の条件を満たす必要がある。式(18)と(21)の場合にあるように、このモデルの状態空間と観測方程式は、式(28a)及び(28b)で与えられる。ここで、ｆ（ｑ_(k)）は式(29)であり、Ｈは式(21)で定義される。

続いて、オプティカルフローの間引きと補間について説明する。間引きと補間の方法には、多くの組み合わせがある。ここでは、間引きと補間オペレータの一例について説明する。

ＭとＮが偶数とする。ＸとＹの座標が共に偶数であるピクセルのオプティカルフローは、稠密推定モデルと同じ方法で推定する。よって、オプティカルフローの状態は、式(30)で与えられる。

他のピクセル（例えば、(1,2,k)）のオプティカルフローは、式(31)にあるように、~ｑ_uと~ｑ_vから補間によって生成する。

変数ｖも、同様に補間する。行列Ｍ_uとＭ_vは、この近似を使って容易に計算することができる。式(31)で定義される補間オペレータを使って、推定される変数の数は、稠密推定モデルにおける推定される変数の約１／４になる。そして、疑似稠密推定モデルの可観測指数は２である。疑似稠密推定モデルは、ゼロベクトルのオプティカルフロー（つまり静止状態）を推定することができる。

さらに、可観測指数を減少させると、推定が速くなる。加えて、オプティカル推定点をピクセルの中心にそろえる必要がない。例えば、推定点を２×２ピクセル領域の中心とすることも可能である。

本実施例では、カルマンフィルタと疑似稠密モデルを使った推定について説明する。稠密モデルについても、同様に、カルマンフィルタを使って推定できる。

まず、ノイズを加えた疑似稠密推定モデル(32)及び(33)を考える。ここで、ｆ(q_(k))とＨは、それぞれ、式(29)と(21)で与えられる。システムノイズη_(k)と観測ノイズζ_(k)は、それぞれ、共分散行列Q_(k)とR_(k)のガウス性白色ノイズであると仮定する。

続いて、カルマンフィルタを使った実装について説明する。

時刻ｋにおける観測値z_(k)を用いて推定された時刻ｍの状態ベクトルを^ｑ_(m|k)とする。そして、P_(m|k)は、誤差共分散行列を示す。疑似稠密推定モデルに対して、式(32)のｆ（ｑ_(k)）のヤコビ行列は、式(34)で与えられる。ここで、Ｎ(ｑ_(k))、Ｊ_u（ｑ_(k)）及びＪ_v（ｑ_(k)）は、それぞれ、式(35a)、(35b)及び(35c)である。

式(32)と(33)の離散時間システムへ適用された拡張カルマンフィルタの更新手続きは、時間更新が式(36)、(37)及び(38)で与えられ、観測更新が式(39)、(40)及び(41)で与えられる。この手続を実行することで、オプティカルフローは、観測した輝度だけを使って推定することができる。

点推定モデルと比較すると、このモデルの状態空間方程式は、差分計算を含む。この差分計算は、推定された状態に適用される。これは、カルマンフィルタによって平滑化されたものである。そのため、本実施例では、ノイズの影響は、観測された輝度を直接用いる従来法に比べて、より抑制される。

続いて、カルマンフィルタのパラメータの設定について説明する。カルマンフィルタにおいて、推定性能（例えば、推定精度と推定値への収束の速さ）は、共分散パラメータに依存する。

ノイズの共分散Ｑ_(k)とＲ_(k)が、時不変であり、共分散行列が式(42)の構造を持つと仮定する。ここで、Ｑ_I、Ｑ_u及びＱ_vは、それぞれ、η_I、η_u及びη_vの共分散行列である。

システムノイズη_u及びη_vは、各時間ステップでオプティカルフローを更新する仮想的なノイズである。そのため、その共分散行列Ｑ_u及びＱ_vは、時間領域の性能を決定する可調整パラメータである。もしＱ_uとＱ_vを小さく設定すれば、推定結果は、時間方向に滑らかで安定である。しかし、結果は、大きな時間遅れを伴う可能性がある。このように、もしオプティカルフローの時間的な変動が大きいと想定されるならば、Ｑ_uとＱ_vは、大きな値に設定するべきである。

オプティカルフロー方程式にシステムノイズη_Iが含まれることから、Ｑ_Iが大きい場合は、オプティカルフロー方程式が正確には成立しないことを表す。一般的に、Ｑ_Iは、小さな値に設定できる。

観測誤差ζは、撮像環境を反映する。その共分散Ｒは、撮像環境に従って設定されるべきである。もし、（例えば量子化ノイズのように、）各ピクセルに、平均がゼロで分散λ²である白色ノイズが現れるのであれば、Ｒは、λ²Ｅとなる。平滑化プレフィルタとして標準偏差σであるガウスフィルタを適用すると、Ｒの(i,j)要素は、式(43)である。ここで、ｄは、インデックスｉとｊに対応するピクセルの間の距離である。推定精度の観点からは、より小さなＲが好ましい。それゆえ、（オプティカルフローの不連続性が存在する場合を除き、）プレフィルタリングによる平滑化を行うことがとても重要である。

撮像環境に起因する一般的な問題の多くは、照明変動である。照明変動が存在すると、勾配法の基本的な前提条件である輝度一定の仮定が成立しないが、Ｒを適当に設定することにより、照明変動の影響を減らしうる。照明変動は、すべてのピクセルに同様に影響する。すなわち、各ピクセルに存在する観測ノイズ間に強い相関が存在する。よって、もし照明変動が、平均がゼロで分散がξ²のランダムな変数としてモデル化されるならば、式(44)が成り立つ。上述の各ピクセルに個別に加わるノイズと共通に加わるノイズの両方が存在する場合には、式(43)と(44)により定義される２つのタイプのＲを加算すればよい。

ここで、共分散行列の設定について、さらに説明する。

ほとんどの場合、初期推定誤差の共分散行列は、未知である。それゆえ、共分散行列は、10¹〜10³を乗じ、初期共分散行列として使うことが好ましい。

付加的な照明変動に対して、観測ノイズ共分散Ｒの設定は、理論的には、式(44)で与えられる。しかしながら、各ピクセルの照明変動の実際の影響は、撮像環境と撮像デバイスの特性に依存する。それゆえ、照明変動の影響を軽減するためには、式(45)のガウスタイプの共分散行列を使うことが好ましい。ここで、σ_ilは、平滑化プレフィルタのものよりも十分大きな値であり、ｄは、インデックスｉとｊに対応するピクセル間の距離である。

本実施例では、疑似稠密推定モデルの効果について、具体例を通じて説明する。

まず、照明変動を伴う回転する画像の系列のアニメーションに疑似稠密推定法を適用する。

従来の勾配法（非特許文献８参照）では、輝度一定に代えて勾配一定を使用する。推定領域において、勾配一定は、並進の流れの存在を仮定するから、このアプローチは、回転及び／又は拡大を伴う画像系列に対して適切ではない。それに対して、本願発明は、並進のフローを仮定せずに、照明変動に対応した推定手法を与える。

輝度Ｉ(r,θ,k)が式(46)で与えられる入力画像を考える。ここで、(r,θ)は、極座標を表す。アニメーションは、図２に示される入力画像を、画像の中心に位置する点(0,0)を回転の中心として、一定の角速度で回転することによって生成される。式(46)の右辺の第２項は、照明変動を示す。パラメータの値は、I₀=128、k_r=16、k_θ=20、k_t=1/10、β_il=3、及び、ω_il=π/10である。そして、時間、空間及び輝度の量子化がなされる。特に、整数時間のｋに対して、(x,y)は、整数の点である。そして、各画像は、標準偏差σ=2のガウスフィルタを使ってプレフィルタ処理がなされる。推定領域は、長方形である。ここで、110≦ｘ≦130であり、190≦ｙ≦210である（幅も高さも21ピクセルである)。

偏差分オペレータＭ_x及びＭ_yは、式(22)で定義される。補間オペレータＭ_u及びＭ_vは、式(31)で定義される。カルマンフィルタ処理で使用される共分散行列は、時不変である。具体的な設定は、式(47)及び(48)である。共分散行列Ｒ_(k)は、式(49)で与えられる。ここで、Ｒ_nsとＲ_ilは、それぞれ、式(43)と(44)で定義される個別と共通のノイズの共分散行列を示す。これらのノイズの源は、主として、量子化ノイズと照明変動である。パラメータ値は、λ²=1×10^-1及びξ²=100である。ｋ＝１での初期推定値は、式(50)、(51)及び(52)で与えられる。

図３は、疑似稠密推定モデルをつかって得られる推定されたオプティカルフローベクトルを示す。線Ｌ₁及びＬ₂は、それぞれ、ｕ及びｖの理論値（theoretical value）である。線Ｌ₃及びＬ₄は、それぞれ、ｕ及びｖの推定値（estimated value）である。この疑似稠密推定モデルは、照明変動の影響を取り除き、フローに関する事前知識を必要とせず、頑健な推定を達成する。

続いて、実世界の画像系列について説明する。ここでは、“RubberWhale”データセット（インターネット＜URL：http://vision.middlebury.edu/flow/＞）で、疑似稠密推定法を適用する。このデータセットは、オプティカルフローの真値（実測）を持つ。そのため、真値（実測）を使って、提案した推定法の性能を評価することができる。このデータセットの詳細は、非特許文献１５と上記のウェブサイトに紹介されている。

オプティカルフロー推定には、ＲＧＢカラーデータセットの最初の５フレームを使う。最初に、カラー画像を、グレースケール画像に変換する。続いて、各画像を、標準偏差σ=2のガウスフィルタでプレフィルタ処理する。与えられた画像は、12×12ピクセルの領域に分割される。ここで、各推定領域は、隣接した推定領域と４ピクセルずつオーバーラップする。そして、各領域に、推定法を適用する。推定後、各領域の中心（８×８ピクセル）での結果を統合し、画像全体の推定値とする。

なお、一般的に、偏微分方程式を解くための適当な境界条件を得ることはできない。そのため、推定領域において、境界の推定を適切に扱うことは、とても重要である。実際の実装では、推定領域を１又は２ピクセル拡張し、推定処理を適用する。拡張した領域は、推定後に、取り除く。

本発明によれば、推定誤差の共分散行列も得られる。しかしながら、この特徴は、計算量を増加させる。そのため、例えば、大きな画像系列を１辺が5〜30ピクセルの矩形の領域に分割し、オプティカルフローを並列によって推定することが好ましい。

偏差分オペレータＭ_x及びＭ_yは、式(22)で定義される。補間オペレータＭ_u及びＭ_vは、式(31)で定義されるものとは異なり、式(53)で決定される。ここで、~ｕ_jは、2×2ピクセルの中心でのフローである。ｕ_*は、ピクセルの中心での補間フロー（interpolated flow）である。図４に、これらの変数を並べて示す。

カルマンフィルタ処理で使用される共分散行列は、時不変である。具体的な設定は、式(54)で与えられる。適当な境界条件を設定することは難しい。よって、推定領域の境界に対応するＱ_I(k)の対角要素は、式(55)及び(56)にあるように、大きい値に設定する。共分散行列Ｒ_(k)は、式(57)で与えられる。ここで、Ｒ_nsとＲ_ilは、それぞれ、式(43)及び(45)で定義される個別と共通のノイズの共分散行列を示す。これらのノイズの源は、主として、量子化ノイズと照明変動である。パラメータ値は、λ²=1×10^-1、ξ²=1、及び、σ_il=20である。ｋ＝１での初期推定値は、式(58)、(59)及び(60)で与えられる。

推定誤差Ｐ_(k|k)の共分散行列の対角要素が初期値の10⁴倍よりも大きくなったならば、対応するピクセルの推定は、失敗したとする。このような場合、^ｑ_I(k|k)を、ｋ番目のフレームの輝度にリセットする。そして、失敗したピクセルに対応する^ｑ_u(k|k)、^ｑ_v(k|k)並びにＰ_(k|k)の行及び列を、式(59)及び(60)で与えられる初期値にリセットする。

なお、オプティカルフロー方程式から数式(8a)，(12a)を導出する際には、簡単のため、時間軸に沿った離散化に、オイラー法を適用した。オイラー法に代えて、推定の正確性を改善するために、ルンゲクッタ法のようなより正確な数値積分法を使用することもできる。

移流方程式に対して解法を数値的に安定にするためは、式(63)の不等式を満たすことが望ましい。ここで、ｃは、フローの大きさ（速度）を示す。δtとδxは、それぞれ、時間と空間の離散化ステップを示す。式(63)の左辺の不等式の値は、クーラン数と呼ばれる。観測更新処理において、離散化により引き起こされる数値的な誤差が小さい方が望ましいことから、条件(63)を満たすようにδｔを小さくする方がよい。よって、は、時間更新処理において、フレーム間を分割し、離散時間ステップを小さくすべきである。

フレーム間の時間ステップは、８つのステップに分割される。そして、各時間ステップで、４次のルンゲクッタ法を適用する。

拡張カルマンフィルタに代えて、共分散行列の時間更新の正確性を改善するため、無香料カルマンフィルタを使用することができる。しかしながら、無香料カルマンフィルタの使用は、計算量を増加させる。

データベースで与えられる真値（実測）を使って推定結果を評価する。非特許文献１５を参照して、統計的な評価を得ることができる。表１及び２は、フローの統計的な評価を示す。エンドポイントエラー（EE）とアンギュラーエラー（AE）は、それぞれ、式(61)及び(62)で定義される。ここで、[ｕ_true ｖ_true]^Tは、オプティカルフローの真値（実測）を示す。AVとSDは、それぞれ、エラー観測の平均と標準偏差を示す。RXは、頑健性の尺度（すなわち、誤差値がＸ以下となるピクセルの個数のパーセンテージ）を示す。AXは、正確性の尺度（すなわち、エラーを低い方から高い方へ並べた後の百分位数でＸ番目のピクセルで誤差値（非特許文献１５参照））を示す。本願発明のパフォーマンスは、現時点では、最先端の手法に比べれば劣っているものの、多重解像度解析と可調整パラメータの適応的な設定を統合することにより推定性能を改善する可能性を有する。

図５〜９に、推定結果を示す。図５は、オリジナル画像の５番目のフレームを示す。図６及び７は、それぞれ、真値（実測）と推定されたフローを示す。そして、図１０は、フローのカラーコーディングを示す。図６において、黒のピクセルは、隠れなどのため、信頼できる真値（実測）を決定できない領域である。図８は、ｋ＝５での共分散行列Ｐ_u(k|k)+Ｐ_v(k|k)の対角成分を示す。これらの値は、各ピクセルのフローの推定値のエンドポイントエラーの意味で信頼性の尺度として使うことができる。すなわち、共分散が大きくなれば、推定結果の信頼性は小さくなる。この信頼性の尺度は、真のオプティカルフローを知ることなく、得られる。これは、本発明の特徴の一つである。図８は、テクスチャが乏しい、黄色の箱（図５のＰ₁）と、オレンジのリング状のもの（図５のＰ₂）と、木の格子（図５のＰ₃）の領域で、推定されたオプティカルフローの信頼性が、低いことを示している。

図９は、ｋ＝５での観測残差|ｚ_(k)-H^ｑ_(k|k-1)|²を示す。この観測残差は、実際のシステムの振る舞いが、動的システムモデルの振る舞いとは異なる場合に増加する。この値は、動的システムの故障を検出するためにも使用される。オプティカルフロー推定の場合、この観測残差を利用して、オプティカルフローの不連続性を検出することができる。図９より、オプティカルフローの不連続を正確に検出できていることがわかる。

１ 物体運動推定装置，３ 画像データ記憶部，５ 状態記憶部，７ 時間発展部，９ 不連続検出部，１１ 輝度時間発展部，１３ オプティカルフロー時間発展部

Claims (6)

1. 画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する物体運動推定装置であって、
   前記動的システムモデルの状態は、前記画像データの推定領域における複数のピクセルの輝度、及び、前記各ピクセルの全部又は一部でのオプティカルフローであって、輝度とオプティカルフローの空間勾配を含まず、
   前記動的システムモデルの前記状態を記憶する状態記憶手段と、
   前記状態記憶手段に記憶された前記状態における前記輝度及び前記オプティカルフローから、新たな時刻の前記状態における前記輝度を演算して、前記状態記憶手段に記憶する輝度時間発展手段と、
   少なくとも前記状態記憶手段に記憶された前記状態における前記オプティカルフローから、新たな時刻の前記状態における前記オプティカルフローを演算して、前記状態記憶手段に記憶するオプティカルフロー時間発展手段を備える物体運動推定装置。
2. 前記動的システムモデルは、ノイズを含み、
   前記輝度時間発展手段及び／又は前記オプティカルフロー時間発展手段は、カルマンフィルタを用いて前記状態に適用される差分計算を含む状態空間方程式により前記ノイズの影響を抑制するとともに、濾過推定値に加えて予測推定誤差共分散を得ることにより、信頼性の尺度として用いて推定されたオプティカルフローの信頼性の評価を可能にする、請求項１記載の物体運動推定装置。
3. 前記状態に含まれるオプティカルフローの数は、前記推定領域におけるピクセル数よりも少なく、前記推定領域における各ピクセルでのオプティカルフローには、前記状態に含まれないものが存在し、
   前記輝度時間発展手段は、
   前記状態記憶手段に記憶された前記状態に含まれるオプティカルフローを補間して前記状態に含まれないオプティカルフローを求めて、
   前記状態記憶手段に記憶された前記状態における前記輝度、並びに、前記状態に含まれるオプティカルフローの全部若しくは一部と共に又は前記状態に含まれるオプティカルフローに代えて前記補間して得られたオプティカルフローから、新たな時刻の前記状態における前記輝度を演算するものであり、
   前記オプティカルフロー時間発展手段は、少なくとも前記状態記憶手段に記憶された前記状態に含まれる前記オプティカルフローから、新たな時刻の前記状態に含まれる前記オプティカルフローを演算する、請求項１又は２記載の物体運動推定装置。
4. 前記動的システムモデルの観測変数と予測推定値との観測残差から、前記オプティカルフローの不連続性を判断する不連続検出手段を備える請求項１から３のいずれかに記載の物体運動推定装置。
5. 物体運動推定装置において、画像データの系列における物体の運動を表すオプティカルフローを、動的システムモデルに基づき推定する物体運動推定方法であって、
   前記物体運動推定装置は、状態記憶手段と、時間発展手段を備え、
   前記状態記憶手段は、前記動的システムモデルの状態を記憶するものであり、
   前記動的システムモデルの状態は、前記画像データの推定領域における複数のピクセルの輝度、及び、前記各ピクセルの全部又は一部での前記オプティカルフローであって、輝度とオプティカルフローの空間勾配を含まず、
   前記時間発展手段が、
   前記状態記憶手段に記憶された前記状態における前記輝度及び前記オプティカルフローから、新たな時刻の前記状態における前記輝度を演算して前記状態記憶手段に記憶し、
   少なくとも前記状態記憶手段に記憶された前記状態における前記オプティカルフローから、新たな時刻の前記状態における前記オプティカルフローを演算して前記状態記憶手段に記憶する時間発展ステップを含む物体運動推定方法。
6. コンピュータにおいて請求項５記載の物体運動推定方法を実現するためのプログラム。