(57)【要約】【課題】従来の動き推定法では得られなかった分岐のある流れ等の不連続的な動きを伴うような流れのある時系列画像における動きを推定する。【解決手段】本発明は、移動する物体の連続及び不連続な動きを同時に推定する方法であって、移動する物体の連続及び不連続な波現象の特性を表現する数理モデルの目的関数によりフローモデルのパラメータを２つの段階により推定する。推定は、画像データのノイズの影響を緩和するためフローベクトル等のモデルパラメータを粗く推定する第一の最適化と、第一の最適化で推定される周波数などの統計量に基づいてモデルパラメータを更に最適化する第二の最適化とを含む。その後、推定されたモデルパラメータと統計モデルとの類似性を判定し、類似性が高いと判断した場合、前記モデルパラメータが所定の精度に達したものと判定し、類似性が低いと判定した場合、モデルパラメータを変更した上で、再度第一の最適化を行う。【選択図】図１
2016066184.tif 000002
【特許請求の範囲】
【請求項１】
移動する物体の連続及び不連続な動きを同時に推定する方法であって、
データ入力部により移動する物体の時系列画像を取り込むステップと、
前記時系列画像の画像データをデータ蓄積部に保存するステップと、
前記移動する物体のモデルパラメータを推定するステップであって、前記モデルパラメータは、画像の輝度の時系列的な変化を表す波生成方程式と、画像処理における動き推定モデルであるオプティカルフローモデルとを統合した、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルの目的関数により推定され、前記目的関数は、波の分散関係式による拘束条件を含み、前記推定するステップは、
第一の最適化部において、前記画像データのノイズの影響を緩和するためのロバスト関数を前記目的関数に適用して、前記目的関数の最小化問題の枠組みで前記モデルパラメータを粗く推定する第一の最適化ステップと、
第二の最適化部において、前記波生成方程式の項の数及び目的関数を適用する画像領域の大きさであるモデルパラメータを第一の最適計算で推定されるパラメータの統計量に基づいて最適化する第二の最適化ステップと
を含む、ステップと、
収束判定部において、前記推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求め、事前に設定された統計モデルと照らし合わせを行うことにより前記モデルパラメータと前記統計モデルとの類似性を判定するステップと、
前記モデルパラメータと前記統計モデルとの類似性が高いと判断した場合、前記モデルパラメータが所定の精度に達したものと判定して、前記モデルパラメータの最終結果を表示部において表示するステップと、
前記モデルパラメータと前記統計モデルとの類似性が低いと判定した場合、前記波生成方程式の項の数及び前記目的関数を適用する画像領域の大きさであるモデルパラメータを変更した上で、前記第一の最適化ステップに戻るステップと
を含むことを特徴とする方法。
【請求項２】
前記モデルパラメータは、オプティカルフローの縦及び横の成分と、波数と、周波数と、位相とが含まれることを特徴とする請求項１に記載の方法。
【請求項３】
前記拘束条件の式の項数は、前記目的関数の項数に連動して変化することを特徴とする請求項１又は２に記載の方法。
【請求項４】
前記収束判定部において、移動する対象ごとに特化した動きに関する特徴を統計分布モデルとして用いることを特徴とする請求項１乃至３のいずれか１項に記載の方法。
【請求項５】
移動する物体の連続及び不連続な動きを同時に推定する装置であって、
移動する物体の時系列画像を取り込むデータ入力部と、
前記時系列画像の画像データを保存するデータ蓄積部と、
前記移動する物体のモデルパラメータを推定する最適化部であって、前記モデルパラメータは、画像の輝度の時系列的な変化を表す波生成方程式と、画像処理における動き推定モデルであるオプティカルフローモデルとを統合した、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルの目的関数により推定され、前記目的関数は、波の分散関係式による拘束条件を含み、前記最適化部は、
画像データのノイズの影響を緩和するためのロバスト関数を前記目的関数に適用して、前記目的関数の最小化問題の枠組みで前記モデルパラメータを粗く推定する第一の最適化部と、
前記波生成方程式の項の数及び目的関数を適用する画像領域の大きさであるモデルパラメータを第一の最適計算で推定されるパラメータの統計量に基づいて最適化する第二の最適化部と
を含む、最適化部と、
前記推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求め、事前に設定された統計モデルと照らし合わせを行うことにより前記モデルパラメータと前記統計モデルとの類似性を判定する収束判定部であって、前記モデルパラメータと前記統計モデルとの類似性が低いと判定した場合、前記波生成方程式の項の数及び前記目的関数を適用する画像領域の大きさであるモデルパラメータを変更した上で、前記第一の最適化部に前記移動する物体のモデルパラメータを推定させる、収束判定部と、
収束判定部において、前記モデルパラメータと前記統計モデルとの類似性が高いと判断した場合、前記モデルパラメータが所定の精度に達したものと判定して最終結果を表示する表示部と、
を含むことを特徴とする装置。
【請求項６】
前記モデルパラメータは、オプティカルフローの縦及び横の成分と、波数と、周波数と、位相とが含まれることを特徴とする請求項５に記載の装置。
【請求項７】
前記拘束条件の式の項数は、前記目的関数の項数に連動して変化することを特徴とする請求項５又は６に記載の装置。
【請求項８】
前記収束判定部において、移動する対象ごとに特化した動きに関する特徴を統計分布モデルとして用いることを特徴とする請求項５乃至７のいずれか１項に記載の装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、移動する物体の動きを推定する方法及び装置に関し、より詳細には、人、車、川及び海などの監視映像において、移動する物体の動きの特徴と状況変化を画像処理により推定する方法および装置に関する。
【背景技術】
【０００２】
実環境において、環境の状況を理解するために、遠隔で移動する物体の移動及び変形の時間変化を、移動対象物体を撮影したカメラ画像などを通じて、解析することは重要な手段の一つとなっている。一方で、移動対象物体の中には、センサを装着できないもの、不定特定多数の人、車及び雲、波といった自然現象、体内の細胞、赤血球など、多くの対象が存在する。特に、剛体系が重心の追従が可能であるのに対して、流体系は重心をもたないため、表面の形状、テクスチャの変化を精度よく検出することが課題となっている。剛体系についても単体から多数の移動体が視線上、折り重なって見える場合は個々の動きよりも全体としての流れに関心が移るため、ダイナミック・テクスチャと呼ばれるように、流体系と同様に扱われることが多い。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】J.L. Barron, D.J. Fleet, and S.S. Beauchemin, “Performance of optical flow techniques”, IJCV, vol. 12, no. 1, pp. 43-77, 1994.
【非特許文献２】堀川清司，[新編]“海岸工学”，東京大学出版会，1981.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
推定される画像上の見かけの動きベクトルはオプティカルフローあるいはフロー（ｆｌｏｗ）とも言われている。オプティカルフローの多くのモデルは滑らかな流れにおいてフローを推定することは可能であるが、分岐、すれ違い、及び衝突などにおいて生じる不連続的な流れについては、フローの平滑化処理のために、分岐などがわかるようなオプティカルフローモデルが得られにくいという課題があった。
【０００５】
従来より流体モデルに基づいたオプティカルフローの推定方法が数多く示されてきている。具体的には、数値シミュレーションで一般的に用いられているナビエ・ストークス（ＮＳ）方程式や連続方程式などがある。しかし、ＮＳ式は連続的な流れ、不連続的な流れの速度を計算の解として導くことができる一方で、フローを推定するときの目的関数の最小化問題においては、解として不連続的なフローが得られない。これは、最小化問題において、反復計算において解が収束していくための条件に、ＮＳ式の解がなるべく滑らかになることが含まれているためである。ＮＳ式は、特に、流体特性を示すレイノイズ数（Ｒｅ）が数千程度までを弱い乱流までを許容するモデルである。そのため、自然現象に見られるようなＲｅが５０００以上の強い乱流を伴った動き（オプティカルフロー）を安定かつ精度よく推定することが課題の一つであった。また、広く用いられている画像の輝度の時系列的な変化が一定という拘束条件は、適用範囲が限定的であるという問題があった。
【０００６】
移動する物体の動き推定は、コンピュータビジョン（ＣＶ）においても重要な位置を占めている。別の基本的な動き検出方法として、ＭＰＥＧに代表される符号化の分野で用いられている相互相関（ＣＣ）法（パターンマッチング法）が知られている（例えば非特許文献１）。ＣＣ法は、連続する２枚の画像間の類似性を相互相関関数などにより評価し、類似する領域を画像間の移動距離（ベクトル）とみなすものである。精度をあげるために、画像を複数のブロックに分割して、それぞれのブロックの移動距離を求める。ここで、ＣＣ法は、類似する画像間で輝度変動がないことが前提条件となっている。しかし実際には、屋外において移動対象の表面が見かけの輝度が変形とともに変化するため、画像間において単純な対応づけができなくなることが知られている。
【０００７】
移動する物体の動き（オプティカルフロー）を推定する際の問題において、数多くの拘束条件が提案され、用いられている（非特許文献１参照）。その多くは、フローの滑らかさに関するものであり、滑らかさの程度は、適用する対象の違い、画像領域、テクスチャに関わらず、同一のもの、固定化されたものである。具体的には、拘束条件として、オプティカルフローモデルの空間一次微分、二次微分を行うこと等が代表的である。そのため、従来において、画像内における動きの推定精度を向上させるには、限界があった。特に、不連続的な動き推定に難があった。
【課題を解決するための手段】
【０００８】
このような問題及び課題を解決するために、本発明の第１の態様は、移動する物体の連続及び不連続な動きを同時に推定する方法であって、データ入力部により移動する物体の時系列画像を取り込むステップと、前記時系列画像の画像データをデータ蓄積部に保存するステップと、前記移動する物体のモデルパラメータを推定するステップであって、前記モデルパラメータは、画像の輝度の時系列的な変化を表す波生成方程式と、画像処理における動き推定モデルであるオプティカルフローモデルとを統合した、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルの目的関数により推定され、前記目的関数は、波の分散関係式による拘束条件を含み、前記推定するステップは、第一の最適化部において、前記画像データのノイズの影響を緩和するためのロバスト関数を前記目的関数に適用して、前記目的関数の最小化問題の枠組みで前記モデルパラメータを粗く推定する第一の最適化ステップと、第二の最適化部において、前記波生成方程式の項の数及び目的関数を適用する画像領域の大きさであるモデルパラメータを第一の最適計算で推定されるパラメータの統計量に基づいて最適化する第二の最適化ステップとを含む、ステップと、収束判定部において、前記推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求め、事前に設定された統計モデルと照らし合わせを行うことにより前記モデルパラメータと前記統計モデルとの類似性を判定するステップと、前記モデルパラメータと前記統計モデルとの類似性が高いと判断した場合、前記モデルパラメータが所定の精度に達したものと判定して、前記モデルパラメータの最終結果を表示部において表示するステップと、前記モデルパラメータと前記統計モデルとの類似性が低いと判定した場合、前記波生成方程式の項の数及び前記目的関数を適用する画像領域の大きさであるモデルパラメータを変更した上で、前記第一の最適化ステップに戻るステップとを含むことを特徴とする。
【０００９】
また、本発明の第２の態様は、第１の態様の方法であって、前記モデルパラメータは、オプティカルフローの縦及び横の成分と、波数と、周波数と、位相とが含まれることを特徴とする。
【００１０】
また、本発明の第３の態様は、第１又は第２の態様であって、前記拘束条件の式の項数は、前記目的関数の項数に連動して変化することを特徴とする。
【００１１】
また、本発明の第４の態様は、第１乃至第３のいずれか１つの態様であって、前記収束判定部において、移動する対象ごとに特化した動きに関する特徴を統計分布モデルとして用いることを特徴とする。
【００１２】
また、本発明の第５の態様は、移動する物体の連続及び不連続な動きを同時に推定する装置であって、移動する物体の時系列画像を取り込むデータ入力部と、前記時系列画像の画像データを保存するデータ蓄積部と、前記移動する物体のモデルパラメータを推定する最適化部であって、前記モデルパラメータは、画像の輝度の時系列的な変化を表す波生成方程式と、画像処理における動き推定モデルであるオプティカルフローモデルとを統合した、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルの目的関数により推定され、前記目的関数は、波の分散関係式による拘束条件を含み、前記最適化部は、画像データのノイズの影響を緩和するためのロバスト関数を前記目的関数に適用して、前記目的関数の最小化問題の枠組みで前記モデルパラメータを粗く推定する第一の最適化部と、前記波生成方程式の項の数及び目的関数を適用する画像領域の大きさであるモデルパラメータを第一の最適計算で推定されるパラメータの統計量に基づいて最適化する第二の最適化部とを含む、最適化部と、前記推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求め、事前に設定された統計モデルと照らし合わせを行うことにより前記モデルパラメータと前記統計モデルとの類似性を判定する収束判定部であって、前記モデルパラメータと前記統計モデルとの類似性が低いと判定した場合、前記波生成方程式の項の数及び前記目的関数を適用する画像領域の大きさであるモデルパラメータを変更した上で、前記第一の最適化部に前記移動する物体のモデルパラメータを推定させる、収束判定部と、収束判定部において、前記モデルパラメータと前記統計モデルとの類似性が高いと判断した場合、前記モデルパラメータが所定の精度に達したものと判定して最終結果を表示する表示部と、を含むことを特徴とする。
【００１３】
また、本発明の第６の態様は、第５の態様であって、前記モデルパラメータは、オプティカルフローの縦及び横の成分と、波数と、周波数と、位相とが含まれることを特徴とする。
【００１４】
また、本発明の第７の態様は、第５又は第６の態様であって、前記拘束条件の式の項数は、前記目的関数の項数に連動して変化することを特徴とする。
【００１５】
また、本発明の第８の態様は、第５乃至第７のいずれかひとつの態様であって、前記収束判定部において、移動する対象ごとに特化した動きに関する特徴を統計分布モデルとして用いることを特徴とする。
【発明の効果】
【００１６】
本発明によれば、従来の動き推定法では得られなかった分岐のある流れ、衝突時の反射のある流れ、及び渦の回転のある流れなど、連続的な動きに加えて、不連続的な動きを伴うような流れのある時系列画像における動きを推定することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明の１実施形態に係る動き推定装置の構成を示すブロック図である。
【図２】本発明の１実施形態に係る動き推定方法の処理の流れを示すフローチャートである。
【図３】移動対象を写した連続した２枚の画像（（ａ）及び（ｂ））から波物理モデルに基づいた対応づけを行う方法に関する概念を示す図である。（ａ）及び（ｂ）は移動対象の連続する２つの画像であり、（ｃ）は（ａ）の画像におけるフローベクトルの基底を示す図であり、（ｄ）は（ｂ）の画像におけるフローベクトルを示す画像である。
【図４】本発明における移動する物体の画像において分割された画像を示す図である。（ａ）物体の画像において分割された画像を示し、（ｂ）は分割部分のフローベクトルを示す画像である。
【図５】本実施形態における収束判定において使用されるヒストグラム分布を示す図である。（ａ）はブレッドシュナイダーのエネルギースペクトラム分布の実験データ（ａ）を示す図である。（ｂ）はフロー推定においてＭを変えたときのヒストグラム分布を示す図である。
【図６】左方より右方へ数千以上の高いレイノイズ数を乱流を画像で撮影し、流れが途中の障害物に衝突したときの様子を示す図である。（ａ）は従来法による動き推定、（ｂ）は本発明による動き推定を示している。
【図７】渦状の時系列画像における動き推定を示す図である。（ａ）は従来法における動き推定、（ｂ）は本発明における動き推定を示している。
【図８】従来法と本発明法とによるフローの推定誤差評価の結果を示す図である。
【発明を実施するための形態】
【００１８】
以下、図面を参照して本発明の実施形態について説明する。
【００１９】
図１は、本発明の一実施形態に係る動き推定装置１００の構成を示すブロック図である。動き推定装置１００は、画像データ入力部１１０と、画像データ入力部１００からの画像データを蓄積する画像データ蓄積部１２０とを備える。また、動き推定装置１００は、画像データ蓄積部１２０に蓄積された画像データと、移動する物体の動きを推定するための流体モデルとを使用して、動き推定を行うために流体モデルのモデルパラメータを粗く推定する第一の最適化部１３０と、第一の最適化部１３０によリ推定されたモデルパラメータを基にモデルパラメータの更なる最適化を行う第二の最適化部１４０とを備える。また、動き推定装置１００は、第二の最適化部１３０において最適化されたモデルパラメータが所定の制度に達したか否かを判定する収束判定部１５０と、所定の精度に達したモデルパラメータの値を表示する表示部１６０と、第二の最適化部１３０において最適化されたモデルパラメータが所定の制度に達していない場合に、再度動き推定を行うためにモデルパラメータの値を変更するパラメータ変更部１７０とを備える。なお、画像データ入力部１１０は、カメラ１１５に接続されている。カメラ１１５は、動き推定装置１００内に含まれていても、外部に接続されていてもよい。
【００２０】
図２は本発明の１実施形態に係る動き推定方法の処理の流れを示すフローチャートである。本実施形態は、まず、ステップ２１０において、カメラやセンサから取得して得られる対象物の時系列画像、又は動画等の画像データを、データ入力部１１０から取り込む。ステップ２２０において、取り込んだ画像データをデータ蓄積部１２０に保存する。
【００２１】
次に、ステップ２３０において、移動する物体の動きを推定するための流体モデルを定義する。ここで、本実施形態における流体モデルは、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルであり、画像の輝度の時系列的な変化が複数の三角関数の項の和からなる波生成方程式と、画像処理における動き推定モデルである、オプティカルフローモデルとを統合して目的関数として定義した。流体モデルについては、まず、画像の輝度の時系列的な変化が複数の三角関数の項の和からなる波生成方程式を定義する。次に、定義した波生成方程式をオプティカルフローモデルに適用し、解の安定性、フローベクトル推定の精度の向上を図るための波の分散関係式による拘束条件と、画像データのノイズや外れ値の影響を緩和するロバスト関数を適用して、目的関数を定義する。
【００２２】
次に、ステップ２４０において、第一段階の最適化を行う。第一段階の最適化においては、第一の最適化部１３０において、ステップ２２０において取り込んだ画像データに目的関数を適用して、目的関数の最小化問題の枠組みにより、画像データのノイズや外れ値の影響を緩和して、モデルパラメータを粗く推定する。具体的には、オプティカルフローモデルのフローベクトルの縦横の成分、更には一定画像領域の波の波数、周波数、位相、及び見かけの厚さを推定する。
【００２３】
ステップ２４０においてモデルパラメータを粗く推定した後、ステップ２５０において、第二段階の最適化を行う。第二段階の最適化は、第二の最適化部１４０において波生成方程式の項の数や目的関数を適用する画像領域の大きさなどといったモデルパラメータを第一段階の最適化において推定されるモデルパラメータ、特に周波数などの統計量に基づいて更に最適化する。
【００２４】
ステップ２６０において、推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求める。
【００２５】
その後、ステップ２７０において、収束判定を行う。収束判定部１５０において、ステップ２６０において求められたパラメータに関する頻度分布と、事前に設定された統計モデルと照らし合わせを行って、推定されたモデルパラメータと事前に設定された統計モデルとの類似性を判定する。
【００２６】
ステップ２７０において、推定されたモデルパラメータと統計モデルとの類似性が高いと判定された場合、ステップ２８０において、モデルパラメータが所定の精度に達したものと判定してモデルパラメータの最終結果を表示部１６０に表示する。
【００２７】
一方で、ステップ２７０において推定されたモデルパラメータと統計モデルとの類似性が低いと判定された場合は、ステップ２９０において画像領域の大きさや項数といったモデルパラメータを変更した上で、再度、ステップ２４０に戻り、第一の最適化及び第二の最適化を行い、モデルパラメータを最適化する。繰り返し最適化計算を行うことで、連続的及び不連続的な波を含むフローを同時に推定することができる。
【００２８】
上記において、モデルパラメータには、オプティカルフローの縦横の成分に加えて、波数、周波数及び位相が含まれている。また、拘束条件の項数が目的関数の項数に連動して大小に変化すること、さらに、収束判定部において、移動する対象ごとに特化した動きに関する特徴を統計分布モデルとして用いることができる。
【００２９】
図３は、移動対象を写した連続した２枚の画像から波物理モデルに基づいた対応づけを行う方法に関する概念を示す図である。図３において、移動対象は水面における波を例としている。図２のステップ２１０において、カメラやセンサから取得して得られる対象物の時系列画像、又は動画等の画像データを、データ入力部から取り込む。取り込んだ画像のうち、連続する２枚の画像が、図３の（ａ）（時刻ｔ−１）及び（ｂ）（時刻ｔ）である。図２のステップ２２０において、取り込んだ連続する２枚の画像である図３の（ａ）及び（ｂ）の画像データをデータ蓄積部に保存する。
【００３０】
次に、図２のステップ２３０において、取り込んだ画像内の移動する物体の動きベクトル（フローベクトル）を推定するための流体モデルを定義する。本実施形態の流体モデルは、連続及び不連続な動きを含む波現象の特性を表現する波に関する数理モデルであり、波物理現象（非特許文献２参照）に関する数理モデル（波生成方程式）と、画像処理又はコンピュータビジョン（ＣＶ)の枠組みにおける動き推定モデルである、オプティカルフローモデルとを統合させた目的関数であることを特徴としている。カメラ視線からみた移動対象については、時間の変化と共に変形を伴うことも少なくない。特に、自然現象において、雲、波及び煙などは、変形を伴う移動対象の典型的な例である。これらの移動対象は、見かけの表面のテクスチャ、輪郭線などが時間ごとに変化するために、移動対象を画像間の対応づけ、フローベクトルの推定を行うためには、モデルが必要である。
【００３１】
そこで、画像内における移動、変形をする対象について、剛体、流体に関わらず、２次元位置（ｘ，ｙ）、時間tにおける見かけの移動対象の高さ（画像の輝度）を、Ｈ（ｘ，ｙ，ｔ）とする。また、ある画像領域における波は、異なった波数、周波数、位相をさまざまな波が線形に重なり合ったものと近似する。すると、ある時間ｔにおける波の画像の輝度は、
【００３２】
【数１】
2016066184.tif 000003
【００３３】
という波生成方程式により表現される。ここで、Ｍは線形和の項数である。また、ａｍは第ｍ番目の波の振幅、（ｋｘ，ｍ，ｋｙ，ｍ）は第ｍ番目の波の波数、ｆｍは第ｍ番目の波の周波数、θｍは第ｍ番目の波の位相である。項数Ｍの決定方法は後述する。
【００３４】
図３において、連続した２枚の波画像（時刻ｔ−１及びｔ）を入力としたとき、波画像内の局所的なフローベクトルを推定するために、２枚の画像間において対応づけを行う。まず、時刻ｔ−１において、取得した波画像（図３（ａ））の一部の領域内に複数の基底を与える（多重基底：図３（ｃ））。次に時刻ｔにおいて取得した波画像（図３（ｂ））において、図３（ｃ）の基底に対応する位置における物体は、図３（ａ）に対して移動、変形する。ここで、図３（ａ）の基底と、図３（ｂ）の基底に対応する位置との間のベクトルを推定されるフローベクトルとみなす（図３（ｄ））。図３の波のような流体現象の場合は、移動対象の移動と変形のため、対応付けが乱れる問題があったが、式（１）に従って輝度が変化することを考慮した目的関数を加味することで、画像間の対応付けの乱れは緩和される。局所的な画像領域においても、式（１）に示すように複数の三角関数基底から構成されている数理モデルにより表現される。
【００３５】
次に、式（１）の波生成方程式の特性をオプティカルフローモデルに取り入れるために、ＣＶにおけるオプティカルフロー法の枠組みについて述べる。入力する画像１枚の輝度をＩ（ｘ，ｙ，ｔ）で定義し、単位時間当たりの輝度の変化が一定ではないとする。すなわち、
【００３６】
【数２】
2016066184.tif 000004
【００３７】
である。続いて、輝度の変化量については、式（１）の数理モデルを用いると、
【００３８】
【数３】
2016066184.tif 000005
【００３９】
と表される。ここで、推定されるフローベクトルを、
【００４０】
【数４】
2016066184.tif 000006
【００４１】
とおく。“ｔｒ”はベクトルの転置を示す。式（３）の左辺において、テイラー展開をすると、
【００４２】
【数５】
2016066184.tif 000007
【００４３】
となる。ここで、Ｉｘ、Ｉｙ、Ｉｔはそれぞれｘ、ｙ、ｔに関する一次微分である。
【００４４】
以上より、式（３）を整理すると、
【００４５】
【数６】
2016066184.tif 000008
【００４６】
を得る。ここで、
【００４７】
【数７】
2016066184.tif 000009
【００４８】
である。
【００４９】
数理モデル（波生成方程式）を下記の目的関数に組み入れる準備として、式（５）を２次形式の誤差関数として定義する。
【００５０】
【数８】
2016066184.tif 000010
【００５１】
フローベクトルの推定において、フローベクトルに関する拘束条件を付加することにより、解の安定性やフローベクトル推定の精度の向上がなされることが知られている。そのため、本発明においてもフローベクトルに関する拘束条件を加味して目的関数を定義する。なお、式（１）の数理モデルにおいては、波の連続性と不連続性の両方の表現が含まれていることに注意する。そこで、速度に関する波物理特性として分散関係式を目的関数に導入する。この分散関係式は、波数が増加すると角周波数が増加することを示すものである。分散関係式は、
【００５２】
【数９】
2016066184.tif 000011
【００５３】
で定義される。ここで、ωは角周波数、ｋは波数、ｈは見かけの厚さ、ｇは重力である。
【００５４】
なお、直接的に見かけの深さを見いだすのは困難であるが、後述する第一段階の最適化において、式（７）のｋ、ｈ及びｇの３つのパラメータがバランスをもつようにフローベクトルを推定していく。
【００５５】
他の基本的な変数について、１／ｆは周期、Ｌ＝２π／ｋは波長である。また、スカラー速度ｃはω＝ｃｋと表される。本発明では、式（７）においてｔａｎｈ（ｘ）を扱いやすくするために、テイラー展開
【００５６】
【数１０】
2016066184.tif 000012
【００５７】
により近似した。ただし、
【００５８】
【数１１】
2016066184.tif 000013
【００５９】
である。
【００６０】
これにより、式（７）は、
【００６１】
【数１２】
2016066184.tif 000014
【００６２】
と定義できる。ここで、ｃはフローベクトルにより、
【００６３】
【数１３】
2016066184.tif 000015
【００６４】
と関係づけられる。以上、式（８）、（９）より、新しい速度に関する拘束条件として次を得た。
【００６５】
【数１４】
2016066184.tif 000016
【００６６】
ただし、
【００６７】
【数１５】
2016066184.tif 000017
【００６８】
ここで、
【００６９】
【数１６】
2016066184.tif 000018
【００７０】
は、ユークリッドノルムである。
【００７１】
第一段階の最適化（図２のステップ２４０）において、モデルパラメータを推定するが、データにはノイズが含まれている。したがって、この外れ値の影響を少なくするために、ロバスト推定法を適用して目的関数を定義する。本実施形態においては、ロバスト関数ρとしてローレンツ関数を選択した。σは形状パラメータと呼ばれており、５．０とおいた。また、式（６）と式（１０）の影響具合について重み係数をηとおく。本実施形態では、このηを最適化する。さらに、２次元画像領域をΩとおく。以上より、本発明の最終的な目的関数は、
【００７２】
【数１７】
2016066184.tif 000019
【００７３】
となる。ただし、
【００７４】
【数１８】
2016066184.tif 000020
【００７５】
である。
【００７６】
第一段階の最適化（図２のステップ２４０）において、画像データのノイズや外れ値の影響を緩和するためにロバスト関数を適用した目的関数（式（１１））の最小化問題の枠組においてモデルパラメータを粗く推定する。具体的には、式（６）と式（１０）より定義した目的関数（式（１１））を使用して、８個の未知数を推定する。ここで、第一の最適化において推定する８つのパラメータは、ｕ、ｖ、ｋｘ、ｋｙ、ｆ、θ、ｈである。連続した２枚の画像について、２枚の異なる時間の画像の対応する一定画像領域より複数の極所的なパラメータを算出し、それぞれ目的関数に代入していくことにより目的関数を最小化する。目的関数の最小化により、８つのパラメータがバランスよく推定されることになる。式（１１）については、Ｍ、Ωの初期値を与えることで、最小化問題として解くことができ、画像データからフローベクトルをΩごとに推定することができる。
【００７７】
式（１１）の第一項のρｉｍｇ、第二項のρｖｅｌはそれぞれロバスト関数である。式（１１）は、η＝０でも解くことはできる。ｅ０及びｅ１，ｍはＭ項からなる複数の基底をもつ。ｅ０のＦＷは一定画像領域内において、異なる波数
【００７８】
【数１９】
2016066184.tif 000021
【００７９】
と見かけの厚さｈｍをＭ項含むため、ｅ１，ｍの項数もＭである。一方で、項数Ｍが推定された場合は、
【００８０】
【数２０】
2016066184.tif 000022
【００８１】
の項数Ｍも変わる。画像領域の大きさは｜Ω｜は、次のように最適化した。第一の最適化において、式（１１）は、初期値｜Ω｜とＭが与えられると、共役勾配法（Ｃｏｎｊｕｇａｔｅ  Ｇｒａｄｉｅｎｔ（ＣＧ））より最小化させることができる。収束性を高めるため、ロバスト関数中のσは反復回数とともに、一定値から小さくした。式（１１）において、９つの未知数に関してその一次微分値を用いた。反復計算は、反復誤差が０．００１未満になるまで続けられた。第一段階の最適化においては、反復計算に際して、初期値として、｜Ω｜＝１０×１０ｐｉｘｅｌと、Ｍ＝３とを与えた。また、９つの未知数の初期値は０とおいた。
【００８２】
計算を進めるためには、画像Ｉ（ｘ，ｙ，ｔ）を計算空間内で離散化した。ここで、
【００８３】
【数２１】
2016066184.tif 000023
【００８４】
とおいた。ただし、ｎ＝１，２，・・・と整数である。
【００８５】
【数２２】
2016066184.tif 000024
【００８６】
Δｘ＝Δｙ＝Δｔ＝１．０
Ｉ（ｘ，ｙ，ｔ）のｘ，ｙ，ｔの偏微分値は有限差分法により、
【００８７】
【数２３】
2016066184.tif 000025
【００８８】
とおいた。Ｍを除いて、すべての未知数は実数である。大きさ｜Ω｜ｐｉｘｅｌｓは１枚の画像よりも小さく、計算を進めるときは、Ωを重複させながら行った。
【００８９】
第二段階の最適化（図２のステップ２５０）は、波生成方程式の項の数や目的関数を適用する画像領域の大きさなどといったモデルパラメータを、第一段階の最適化において推定される周波数などの統計量に基づいて最適化する。第二段階の最適化においては、Ｍ及び｜Ω｜の最適化のために、第１の最適化により推定される周波数の画像内での生起頻度数を与える。
【００９０】
図４は、本発明における移動する物体の画像において分割された画像を示す図である。図４において、渦を巻いているハリケーン画像を画像分割の例としている。
【００９１】
Ｍが大きい場合は、式（１）の多項式の項数が増える。多項式の項数の増加により、ｃｏｓ-ｓｉｎｅ基底の線形の重ね合わせを行うことができ、異なる画像の表面やテクスチャについて幅広く近似することが可能である。一方で、２次元画像領域Ωの大きさ｜Ω｜が大きくなる場合は、フローベクトルはより平滑化されたものになる。したがって、局所的な動きを得るためには、｜Ω｜を小さくする必要がある。それゆえに、Ｍ及び｜Ω｜のそれぞれの大きさを最適化することによりフローベクトルの推定精度は向上する。図４（ａ）のように渦を巻いているハリケーン画像においては、小さい画像｜Ω｜を複数重ね合わせている。また、図４（ｂ）には、推定されるフローの一部を例として示している。
【００９２】
式（１１）の目的関数において、波物理特性が２つの拘束条件として含められたが、Ｍ、｜Ω｜については、式（１１）の目的関数（第一段階の最適化）だけでは最適化することが困難である。そこで本実施形態では、推定されたモデルパラメータについて、収束判定（図２のステップ２７０）を行う。収束判定は、波物理現象の枠組みにおいて、実環境のさまざまな方向性や大きき、スケールをもつ波現象を観測された統計モデルを利用して、推定されたモデルパラメータに基づいてパラメータに関する頻度分布を求めて（図２のステップ２６０）、事前に設定された統計モデルと比較することにより行う。このような統計モデルの一つに、ブレッドシュナイダー（Ｂｒｅｄｔｓｃｈｎｅｉｄｅｒ）のエネルギースペクトラム分布があり、
【００９３】
【数２４】
2016066184.tif 000026
【００９４】
と記述される。ここで、Ｈ＊及びＴ＊は定数であり、Ｈ＊は有義波高と、Ｔ＊は有義周期と呼ばれている。
【００９５】
図５は、本実施形態における収束判定において使用されるヒストグラム分布を示す図である。図５（ａ）はブレッドシュナイダーのエネルギースペクトラム分布の実験データ（ａ）を示す図である。また、図５（ｂ）はフロー推定においてＭを変えたときのヒストグラム分布を示す図である。式（１１）より推定された図５（ａ）が示すように、ブレッドシュナイダーのエネルギースペクトラム分布は単峰性であり、縦軸は無次元化周波数を、横軸はエネルギースペクトラムを示す。分布の帯域は広く、中心ピークはやや低周波数側にあるのが特徴である。時系列画像と式（１１）の目的関数からは、複数のパラメータが推定されるが、周波数ごとにエネルギーが変わるものと仮定した。これについては、画像中、周波数ごとに得られる画素数より、ヒストグラムをつくる。これを｜Ω｜の局所的な画像領域ごとに得る。続いて、大局的なヒストグラムを得るために、すべての｜Ω｜に関して平均値を求めた。このとき、ヒストグラムの平均値の面積は１になるように正規化した。
【００９６】
これにより、式（１２）のブレッドシュナイダー分布（図５（ａ））と、式（１１）より推定されたヒストグラム分布（図５（ｂ））との類似性を比較することができる。図５（ｂ）には、Ｍ＝１〜６まで変えた場合のヒストグラム分布をそれぞれ重ねて表示してある。図５（ｂ）において、Ｍが大きくなるほど、図５（ａ）のグラフとのマッチングが良くなる。このとき、｜Ω｜は１０×１０ｐｉｘｅｌｓの例である。図２のステップ２７０において、モデルパラメータと統計モデルとの類似性が高いと判定された場合、つまり、モデルパラメータと統計モデルとのズレが所定の範囲内に収まっている場合、所定の精度に達したものと判定する。一方で、モデルパラメータと統計モデルとの類似性が低いと判定された場合、つまり、モデルパラメータと統計モデルとのズレが所定の範囲内に収まっている場合、は、画像領域の大きさや項数といったモデルパラメータを変更した上で（図２のステップ２９０）、再度、第一段階の最適化（図２のステップ２４０）に戻り、目的関数を最適化する。このようにして、Ｍと｜Ω｜を一定範囲において変化させて、ヒストグラム分布の類似性をすべての組み合わせについて求める。
【００９７】
Ｅｓｔ（ｆ）Ｍ，｜Ω｜を、式（１１）で推定されたヒストグラム分布とすると、ヒストグラム分布は、正規化相関係数ＣＲ（・,・）により評価できる。即ち、
【００９８】
【数２５】
2016066184.tif 000027
【００９９】
ただし、
【０１００】
【数２６】
2016066184.tif 000028
【０１０１】
ＣＲ（・，・）が１．０に近づくほど、Ｂ（ｆ）とＥｓｔ（ｆ）Ｍ，｜Ω｜は類似する値になる。一組の時系列画像ごとに、ＣＲ値が最大となるときの、Ｍと｜Ω｜が選択され、局所的にも大局的にも最適化され、フローの推定精度は向上する。
【０１０２】
なお、収束判定においては、他の統計モデルを適用することもできる。特定の対象については、事前情報（知識）から統計モデルを作成して適用することもできる。
【０１０３】
図６及び図７は、従来法と本発明法による時系列画像から動き（フロー）推定した結果を示す。このうち図６は、左方より右方へ数千以上の高いレイノイズ数を乱流を画像で撮影し、流れが途中の障害物に衝突したときの様子を示す図であり、（ａ）は従来法による動き推定、（ｂ）は本発明による動き推定を示している。図６（ａ）に示す従来法の動き推定においては、やや滑らかなフローが推定されているが、図６（ｂ）に示す本発明の動き推定においては、障害物の影響を受けたやや不連続なフローが得られた。
【０１０４】
次の例として、渦状の時系列画像における動き推定を図７に示す。図７の画像において、（ａ）は従来法における動き推定、（ｂ）は本発明における動き推定を示している。図７（ａ）に示す従来法の動き推定よりも図７（ｂ）における本発明における動き推定により推定されたフローベクトルの方が渦の流れ構造をよりはっきりと得るものとなった。渦の中心部ほど、局所的に不連続な動きを示すため、図７（ａ）に示す従来法の滑らかな動き拘束条件では十分にその特徴を細かく捉えることができなかったと考えられる。この渦のテクスチャのように、それぞれの場所において、表面の周波数が異なったパターンを示しているが、本発明では、三角関数の項数による近似が最適化において変化するために、周辺から渦中心までフローが得られたと考えられる。
【０１０５】
図８は、従来法と本発明法とによるフローの推定誤差評価の結果を示す。渦や不連続性の動きにおいて時系列画像を数値シミュレーションにより生成し、従来法の平均誤差（正解フローと推定フローの一画素当りの平均誤差）と、本発明による比較評価をした結果、定量的にフロー推定誤差が改善された。なお、推定誤差を求めるときの正解のオプティカルフローは、数値シミュレーションで事前に動き情報が１００％分かっているもの、又は実際の画像において目視で特定の狭い領域内の動きを拾ったものを用いている。
【産業上の利用可能性】
【０１０６】
本発明は、画像工学、医用画像工学、生物学、防災学、流体力学、土木工学、海洋学、河川学、及び災害監視等の分野に関して、移動する物体の動きの推定を行うことができる。
【符号の説明】
【０１０７】
１１０  画像データ入力部
１２０  画像データ蓄積部
１３０  第一の最適化部
１４０  第二の最適化部
１５０  収束判定部
１６０  表示部
１７０  モデルパラメータ変更部
【図１】
2016066184.tif 000029
【図２】
2016066184.tif 000030
【図５】
2016066184.tif 000031
【図６】
2016066184.tif 000032
【図８】
2016066184.tif 000033
【図３】
2016066184.tif 000034
【図４】
2016066184.tif 000035
【図７】
2016066184.tif 000036