JP3557982B2

JP3557982B2 - オプティカルフロー推定方法

Info

Publication number: JP3557982B2
Application number: JP2000022967A
Authority: JP
Inventors: ロイセバスチャン
Original assignee: NEC Corp
Current assignee: NEC Corp
Priority date: 1999-04-20
Filing date: 2000-01-31
Publication date: 2004-08-25
Anticipated expiration: 2020-01-31
Also published as: EP1047019A3; EP1047019A2; CA2297233C; US6507661B1; JP2000306108A; CA2297233A1

Description

【０００１】
【発明の属する技術分野】
本発明は、機械に人間の視覚系を持たせることを目的とするマシンビジョン（機械視覚）の分野に関するものであり、特に複数の画像間でのオプティカルフローを有効に推定する方法に関する。
【０００２】
【従来の技術】
運動推定は、ロボット工学（ナビゲーションおよび障害物回避を含む）、自律走行自動車、医学画像解析（血管造影等の非剛直運動を含む）等の多くの種類のマシンビジョン処理の際に生じる重要な問題である。２個以上の時系列連続画像間の動きが小さい場合、２個の異なる像間の２次元の動きベクトル場として定義されるオプティカルフローによって説明される。オプティカルフローは、画像中の対象物が、どのように運動し、どこに向かって運動し、どの程度の速さであるかを示すものである。
【０００３】
輝度一定の仮定（ＣｏｎｓｔａｎｔＢｒｉｇｈｔｎｅｓｓＡｓｓｕｍｐｔｉｏｎ：以下ＣＢＡと称する）下では、画素の動きは１次元方向に制限することができる。しかしながら、１個の画素におけるフローには２成分（すなわち、方向（向きおよび角度）と絶対値（すなわち速度））が存在するため、オプティカルフロー推定は固有の困難さを有する問題である。従って、その問題に対処すべく、いくつかの試みが行われてきた。
【０００４】
ほとんどの先行技術が、フロー場を「規則化」することで、すなわちフロー場に対して何らかの形の平滑化を行うことで、その問題を克服するものである（Ｈｏｒｎｅｔａｌ．， ”ＤｅｔｅｒｍｉｎｉｎｇＯｐｔｉｃａｌＦｌｏｗ， ”ＡｒｔｉｆｉｃｉａｌＩｎｔｅｌｌｉｇｅｎｃｅ，Ｖｏｌ．１７，ｐｐ．１８５−２０３（１９８１）；Ｈ．Ｎａｇｅｌｅｔａｌ．， ”ＯｎＴｈｅＥｓｔｉｍａｔｉｏｎＯｆＯｐｔｉｃａｌＦｌｏｗ：ＲｅｌａｔｉｏｎｓＢｅｔｗｅｅｎＤｉｆｆｅｒｅｎｔＡｐｐｒｏａｃｈｅｓＡｎｄＳｏｍｅＮｅｗＲｅｓｕｌｔｓ，” ＡｒｔｉｆｉｃｉａｌＩｎｔｅｌｌｉｇｅｎｃｅ，Ｖｏｌ．３３，ｐｐ．２９９−３２４（１９８７）参照）。ＣＢＡはまた、２個の画像間の最小二乗差を最小化することでフロー場を推定するエネルギー最小化とすることもできる（Ｐ．Ａｎａｎｄａｎ， ”ＡＣｏｍｐｕｔａｔｉｏｎａｌＦｒａｍｅｗｏｒｋＡｎｄＡｎＡｌｇｏｒｉｔｈｍＦｏｒＴｈｅＭｅａｓｕｒｅｍｅｎｔＯｆＳｔｒｕｃｔｕｒｅＦｒｏｍＭｏｔｉｏｎ， ”Ｉｎｔ’ｌＪｏｕｒｎａｌｏｆＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｖｏｌ．２，ｐｐ．２８３−３１０（１９８９）；Ａ．Ｓｉｎｇｈ，ＯｐｔｉｃＦｌｏｗＣｏｍｐｕｔａｔｉｏｎ：ＡＵｎｉｆｉｅｄＰｅｒｓｐｅｃｔｉｖｅ，ＩＥＥＥＣｏｍｐｕｔｅｒＳｏｃｉｅｔｙＰｒｅｓｓ（１９９２）参照）。オプティカルフローはまた、小さい画像区画全体にわたって局所輝度を分割することで計算することもできる（Ｂ．Ｌｕｃａｓｅｔａｌ．， ”ＡｎＩｔｅｒａｔｉｖｅＩｍａｇｅＲｅｇｉｓｔｒａｔｉｏｎＴｅｃｈｎｉｑｕｅＷｉｔｈＡｎＡｐｐｌｉｃａｔｉｏｎＴｏＳｔｅｒｅｏＶｉｓｉｏｎ， ”ＤＡＲＰＡＩＵＷｏｒｋｓｈｏｐ，ｐｐ．１２１−１３０（１９８１）参照）。平滑化の問題には、パラメータ化された画像全体の運動モデルを適合化することで対処することもできる（Ｓ．Ｓｒｉｎｉｖａｓａｎｅｔａｌ．， ”ＯｐｔｉｃａｌＦｌｏｗＵｓｉｎｇＯｖｅｒｌａｐｐｅｄＢａｓｉｓＦｕｎｃｔｉｏｎｓＦｏｒＳｏｌｖｉｎｇＧｌｏｂａｌＭｏｔｉｏｎＰｒｏｂｌｅｍｓ，” ＰｒｏｃｅｅｄｉｎｇｓｏｆＥｕｒｏｐｅａｎＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｆｒｅｂｕｒｇ，Ｇｅｒｍａｎｙ，ｐｐ．２８８−３０４（１９８８）参照）。
【０００５】
多くの先行技術による推定方法によれば、コスト関数を最小化することで、輝度の制約と平滑化との間の均衡が得られる。それらの方法は、反復非線形法に基づくものであるため、広域最小値に収束するとは限らず、従って、局所最小値に収束する際に満足できる結果を与えない。
【０００６】
【発明が解決しようとする課題】
本発明の方法は、フロー推定の問題を、マルコフランダム場（ＭａｒｋｏｖＲａｎｄｏｍＦｉｅｌｄ：ＭＲＦ」）の枠組みでのラベリング問題として公式化することで、上記の制限を克服するものである。従って本発明は、フロー場における不連続性を保持しながら、高密度でノンパラメトリックなフローを解くものである。
【０００７】
ある種のＭＲＦでは、グラフ上での最大フロー計算によって、正確な帰納的最大（ＭａｘｉｍｕｍＡＰｏｓｔｅｒｉｏｒｉ：ＭＡＰ）推定値を効率良く得ることができる。最適であることが保証されていることから、この計算によって、局所最小解の問題が回避される。ＭＲＦ公式化およびグラフ理論解を用いる最近の一部の方法について、各種文献等にその例が記載されている（Ｓ．Ｒｏｙｅｔａｌ．， ”ＡＭａｘｉｍｕｍ−ＦｌｏｗＦｏｒｍｕｌａｔｉｏｎＯｆＴｈｅｎ−ＣａｍｅｒａＳｔｅｒｅｏＣｏｒｒｅｓｐｏｎｄｅｎｃｅＰｒｏｂｌｅｍ，” Ｉｎｔ’ｌＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｍｕｍｂａｉ，Ｉｎｄｉａ，ｐｐ．４９２−４９９（１９９８）；ロイ（Ｓ．Ｒｏｙ）によって１９９７年１１月２６日に出願された米国特許出願０８／９７８、８３４号（発明の名称「ＭａｘｉｍｕｍＦｌｏｗＭｅｔｈｏｄＦｏｒＳｔｅｒｅｏＣｏｒｒｅｓｐｏｎｄｅｎｃｅ」。）；Ｈ．Ｉｓｈｉｋａｗａｅｔａｌ．， ”Ｏｃｃｌｕｓｉｏｎｓ，Ｄｉｓｃｏｎｔｉｎｕｉｔｉｅｓ，ａｎｄＥｐｉｐｏｌａｒＬｉｎｅｓＩｎＳｔｅｒｅｏ，” ＰｒｏｃｅｅｄｉｎｇｓｏｆＥｕｒｏｐｅａｎＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｆｒｅｉｂｕｒｇ，Ｇｅｒｍａｎｙ，ｐｐ．２３２−２３７（１９９８）；Ｙ．Ｂｏｙｋｏｗｅｔａｌ．， ”ＭａｒｋｏｖＲａｎｄｏｍＦｉｅｌｄｓＷｉｔｈＥｆｆｉｃｉｅｎｔＡｐｐｒｏｘｉｍａｔｉｏｎｓ，” ＰｒｏｃｅｅｄｉｎｇｓｏｆＩＥＥＥＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎａｎｄＰａｔｔｅｒｎＲｅｃｏｇｎｉｔｉｏｎ，ｐｐ．６４８−６５５（１９９８）参照）。
【０００８】
フロー推定における別の重要な問題は画像導関数の計算である。画像は空間、時間および強度の次元で識別されることから、空間時間導関数の離散的計算の正確さには制限がある。この問題は、複雑な導関数フィルターによってある程度解決される。実際には導関数は、照明の変化、輝度の尺度および反射などの輝度一定の仮定からの逸脱によっても信頼性が低下する。従って、輝度の制約が、「真の」厳密な制約と考えるべきではない。この不確定性の考え方について説明するため、本発明では、輝度の制約を確率的枠組みに入れる。オプティカルフローの確率的解釈についての関連する例が、シモンセリらの論文（Ｅ．Ｓｉｍｏｎｃｅｌｌｉｅｔａｌ， ”ＰｒｏｂａｂｉｌｉｔｙＤｉｓｔｒｉｂｕｔｉｏｎｓｏｆＯｐｔｉｃａｌＦｌｏｗ，”，ＰｒｏｃｅｅｄｉｎｇｓｏｆＩＥＥＥＣｏｎｆｅｒｅｎｃｅｏｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎａｎｄＰａｔｔｅｒｎＲｅｃｏｇｎｉｔｉｏｎ，ｐｐ．３１０−３１５（１９９１））に記載されている。この方法では、非確率的アプローチにおける問題の一部が克服されているが、オプティカルフローの確率の非線形特性について考慮されておらず、また画像導関数における誤差を適切に考慮せずに過度に単純化されたオプティカルフローモデルを用いているため充分な効果を得られるものではない。
【０００９】
そこで、輝度の制約を行ないながら、画像導関数の測定における誤差を適切にモデル化し、しかもそのモデルの環境においてオプティカルフローに対する画像全体の最適解を効率良く得られるオプティカルフローの推定方法が必要とされている。
【００１０】
従って本発明の目的は、複数画像間のオプティカルフローを有効かつ正確に推定するオプティカルフローの推定方法を提供することにある。
【００１１】
本発明の別の目的は、輝度一定の仮定の制約を行ないながら、画像導関数の測定における誤差を適切にモデル化するオプティカルフローの推定方法を提供することにある。
【００１２】
本発明のさらに別の目的は、そのモデルの環境下で、オプティカルフローに対する画像の全体的な最適解を効率良く与えるオプティカルフローの推定方法を提供することにある。
【００１３】
【課題を解決するための手段】
本発明によれば、複数の画像間でのオプティカルフローを推定する方法が提供される。本発明のオプティカルフロー推定方法は、運動の方向成分と運動速度成分を得るステップとから構成されている。運動方向成分を求める方法は、複数の画像の空間時間導関数を用いて第１のグラフを作成するステップと、第１のグラフで第１の最大フローについて解を求めることで、それから第１の最小カットを得るステップと、第１の最小カットから運動方向成分を計算する段階を行うステップとから構成される。また、運動速度成分を求めるステップは、複数の画像の空間時間導関数および前記運動方向成分を用いて第２のグラフを作成するステップと、第２のグラフで第２の最大フローについて解を求めることで、それから第２の最小カットを得るステップと、第２の最小カットから運動速度成分を計算する段階を行うステップとから構成される。そして、運動方向成分および運動速度成分とを組み合わせて、複数画像間のオプティカルフローが推定される。本発明のオプティカルフロー推定方法は、輝度の制約を行ないながら、画像導関数の測定における誤差を適切にモデル化し、そのモデルの環境において、オプティカルフローに対して画像全体の最適解を効率良く提供するものである。
【００１４】
【発明の実施の形態】
次に、本発明の実施の形態について詳細に説明する。
【００１５】
本発明の実施形態について説明する前に、本発明についての理解を深めるため、輝度の制約を行ないながら画像導関数における誤差をモデル化する場合の問題について説明および公式化する。
Ａ．問題についての公式化
輝度の制約は、画素の画像輝度が一定であると仮定して得られる。その結果、空間時間座標に関する画素の強度の全体的導関数はゼロである。従って、以下の式のようになる。
【００１６】
【数６】

式中、Ｉ_ｘ、Ｉ_ｙおよびＩ_ｔは空間時間画像導関数であり、υ_ｘ、υ_ｙは、ｘ方向およびｙ方向でのフロー成分である。この制約は、直線の式を説明するものである（図１（ａ）および図１（ｂ）参照）。図１（ａ）において、斜線を施した「許容」領域は、法線ベクトルυ_ｎについての全ての可能な運動を表している。
図１（ｂ）において、法線ベクトルυ_ｎを中心とする斜線を施した半円における全ての方向が同等の確率を有する。前述のように、輝度の制約は、画像導関数における固有の不確定性によって緩和されるはずである。以下に説明するように、ベイズ（Ｂａｙｅｓｉａｎ）の枠組みにおいて単純かつ直観的な前提（すなわち、問題の解を得る前にわかっている知見を表す、アプリオリ確率分布）を用いることで、ＣＢＡの有用なモデルを得ることができる。
【００１７】
便宜上のため、本明細書では、以下の表記を用いる。空間導関数Ｉ_ｘ、Ｉ_ｙを∇Ｉと称し、空間時間導関数Ｉ_ｘ、Ｉ_ｙ、Ｉ_ｔはＩ_ｄと称する。これらの画像導関数はいかなる方法でも求めることができ、そのための多くのアルゴリズムが公知である。ある画素でのフローは、υ、すなわちυ＝［υ_ｘ、υ_ｙ、１］と表記される。
【００１８】
図１（ａ）および図１（ｂ）に示したように、Ｉ_ｄ ^０と表記される真の空間時間導関数は、動きベクトルυを、Ｉ_ｄ ^０・υ＝０で示される直線上に来るように制限する。確率の形では、Ｐ（υ｜Ｉ_ｄ）は、ノイズ画像導関数Ｉ_ｄによって決まるフローの確率と定義される。画像導関数についての誤差モデルは、以下のように定義される。
【００１９】
【数７】

式中、ｎは観察誤差であり、平均値ゼロおよび何らかの共分散Σにてガウス分布していると仮定されている。Ｐ（υ｜Ｉ_ｄ）を得るために、ベイズの法則を用いて下記式が得られる。
【００２０】
【数８】

加法的ノイズモデルを考慮すると、条件的確率Ｐ（Ｉ_ｄ ^０｜Ｉ_ｄ）はガウス分布であり、平均はＩ_ｄ、共分散がΣである。従って、真の画像導関数Ｐ（υ｜Ｉ_ｄ ^０）によって決まるフローのアプリオリ分布を考慮すると、所望の条件的確率Ｐ（υ｜Ｉ_ｄ）を示すことができる。
【００２１】
式（３）と同じ条件的確率について記載している同様の確率的方法がこれまで用いられているが（シモンセリら（同上）参照）、先行技術の方法は、２つの重要な点で本実施形態と異なっている。第１に、先行技術の方法のノイズモデルは、画像導関数ではなくフローベクトルに誤差の原因があり、それは最初から誤差のあることが知られている。第２に、先行技術の方法では、フローベクトルＰ（υ）でのアプリオリ分布を選択する必要がある。この前提条件は説明が非常に困難であり、運動の種類、シーンでの奥行き分布などによって変わる。さらに、解析を容易にするためには、Ｐ（υ）についてゼロ平均のガウス分布を選択する必要があり、それは実際には実現できる場合が少ない。
【００２２】
それとは対照的に本実施形態では、条件的分布Ｐ（υ｜Ｉ_ｄ ^０）、すなわち画素の真の画像導関数を考慮したフロー確率を選択する必要があるだけである。そこで、本実施形態で使用される前提の方が扱い易く、画面全体の運動パターンＰ（υ）についての知識を必要としない。この前提の選択およびそれが解に与える影響について、以下のセクションで説明する。
Ｂ．輝度制約についての確率モデル
図１（ａ）および図１（ｂ）からわかる通り、動きベクトルυの未知成分はＣＢＡ直線上にあり、角度θによってパラメータ化することができる。これは、可能なθ値の空間を、許容（斜線）領域と非許容領域に分けるものである。許容領域は、Ｉ_ｄ ^０に関連する法線ベクトルυ_ｎを中心とする半円である。そこで、必要な前提条件Ｐ（υ｜Ｉ_ｄ ^０）をθの条件的前提条件と記載することができる。
【００２３】
最も弱い形では、θに関する前提は単に、許容領域でのフロー方向が同等の確率を有するというものである（図１のＰ（θ｜Ｉ_ｄ ^０）参照）。従来技術における“核”は、以下の通りである。
【００２４】
【数９】

式中、θ_ｎは、法線ベクトルυ_ｎの方向である。所望に応じて、フローについての具体的知見を用いて、フロー方向の条件的分布を変えることができる。例として、フローの速度が厳密に規定された場合に、θ_ｎからの許容される角度逸脱の範囲を縮小することができる。
【００２５】
真のフローはθによって十分にわかることから、条件的前提Ｐ（θ｜Ｉ_ｄ ^０）を選択することで、条件的前提Ｐ（υ｜Ｉ_ｄ ^０）が自動的に決まる。それは、以下のように示すことができる。
【００２６】
【数１０】

式中、υ_ｎは‖Ｉ_ｔ‖／‖∇Ｉ‖に等しい速度を有する。式（３）、（４）および（５）を比較することで、Ｐ（υ｜Ｉ_ｄ）はＰ（Ｉ_ｄ ^０｜Ｉ_ｄ）の関数として表すことができる。しかしながら、この関数は簡単な解析型を有するものではない。実際、それは数値的に評価するのが好ましい。
【００２７】
各画素が画像導関数Ｉ_ｄを生じる。次に、分布Ｐ（Ｉ_ｄ ^０｜Ｉ_ｄ）から誘導されて、一連の実現的な値が得られる。各実現的な値について、従来の核が所望の分布Ｐ（υ｜Ｉ_ｄ）上に累積される。真のフローＰ（υ｜Ｉ_ｄ）の条件的分布は、異なる方向を示す核の加重平均であり、その場合加重は、条件的分布Ｐ（Ｉ_ｄ ^０｜Ｉ_ｄ）によって決定される。
【００２８】
上記の確率分布について説明するため、図２には、３種類の画像導関数［２０、２０、１０］、［１０、１０、５］および［４、４、２］についての法線フロー分布および条件フロー分布Ｐ（υ｜Ｉ_ｄ）を示してある。これらの導関数は、各種量の画像テクスチャを特徴づける領域で認められる同じ法線フローベクトル［−０．３５、−０．３５］に相当する。画像導関数における誤差は、空間時間次元の各次元での標準偏差が１のガウス分布によってモデル化される。高レベルのテクスチャの場合（Ｉ_ｄ＝［２０、２０、１０］）、輝度制限と法線フローベクトルは信頼性が高い。従って、得られる法線フロー分布は非常にコンパクトであり、フロー分布全体は、輝度制限線方向のみが不確定である。中程度のテクスチャの場合（法線フローベクトルの位置および全フローの両方における不確定性が高くなる。画像テクスチャの量が低い場合（Ｉ_ｄ＝［４、４、２］）、法線フローおよび全フローの両方の値における不確定性の程度が大幅に高くなる。これは、法線フローおよび輝度制限の信頼性が局所区画にある画像テクスチャの量によって決まるという直観的事実に相当するものである。低テクスチャ領域ではこのモデルは、輝度制限線からの大幅な逸脱をもたらすものではない。
【００２９】
図３には、図２と同じＩ_ｄ値を用いて、フローの方向および速度について得られる分布を示してある。図３において各縦軸には、記載されている画像導関数についてのフローの方向（上図）および速度（下図）の条件的分布を示してある。図３からわかる通り、フロー方向の分布は本質的に、利用可能なテクスチャの量による影響を受けない。しかしながら、テクスチャの量は、フローの速度の確率に大きく影響する。高テクスチャの場合、法線フローが信頼性が高く、従って全フローの速度は法線フロー（垂線で示してある）の速度より大きいはずである。テクスチャの量が減少するに連れて、法線フローの速度の信頼性が低くなり、法線フローより小さいフローの速度の確率が高くなる。これは、信頼性の低い法線フローは全フロー値の範囲をさほど制限するものではないという直観的事実を裏付けるものである。極端な場合は、識別可能な運動がない場合、すなわちＩ_ｄ≒［０、０、０］の場合であると考えられる。その場合、シミュレーションされる方向分布は、［−π、π］の範囲で均一である。結果的に、そのような画素の方向は、強制的な平滑化のために、完全に隣接画素の方向によって決まることになる。
Ｃ．オプティカルフローの解法
ほとんどの先行技術の方法では、フロー場は局所的には平滑であるとの仮定のもとに、輝度の制約に対する忠実度を左右するコスト関数を最小とすることによりオプティカルフローの推定を行っている。奥行きの不連続性のため、フロー場は各区分ごとに平滑であるのが普通である（すなわちそれには、大きい不連続部によって分離された平滑運動区画がある）。平滑化を行うことにより、フローの推定が、それらの境界部を通って平滑化され、結果的にフロー推定が不正確になる。
【００３０】
概して、得られるコスト関数は、反復非線形最適化法を用いて最小化され、広域最小値に収束する保証はない。フロー推定を、制限がある種類のＭＲＦモデルに関するラベル問題として公式化することで、反復法を回避することができ、画像全体の最適解が保証される。グラフ上で最大フロー問題への変換を行うことで、このラベル問題の正確な帰納的最大（ＭＡＰ）推定値を得ることができる。この広域最小値は、大きい不連続部を保存する傾向を有する。
【００３１】
ＭＲＦのＭＡＰ推定値に対して最大のフロー解を得るには、ラベルが１次元である必要がある。残念ながら、全ての画素のフローが、２次元ベクトルによって説明される。そのため、フローを２個の１次元空間にパラメータ化する必要がある。本実施形態においては、２次元フロー場［υ_ｘ、υ_ｙ］は、相当する角度−速度表示［θ、ｍ］へとパラメータ化される。このパラメータ化の好ましい選択について、以下でさらに詳細に説明する。
【００３２】
一般的なＭＲＦの考え方は当業界では知られており、それについての詳細な説明がリー（Ｓ．ＬＩ）らの著作にある（Ｓ．Ｌｉｅｔａｌ．，ＭａｒｋｏｖＲａｎｄｏｍＦｉｅｌｄＭｏｄｅｌｉｎｇＩｎＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇｐｕｂｌ．（１９９５））。しかしながら、本実施形態の方法の公式化に先だって、ＭＲＦの基礎となる考え方を以下に簡単に説明する。
【００３３】
Ｓ＝｛０、．．．、ｍ−１｝で表される箇所（画素）の集合が与えられた場合には、個々のラベル問題は、ラベル集合Ｌ＝｛０、．．．、ｍ−１｝から引き出される固有のラベル（方向または速度）を各箇所に割り当てるという問題となる。ラベルの各構成は、確率変数Ｆ＝｛Ｆ_０、．．．、Ｆ_ｍ−１｝群から引き出される。ＭＲＦのマルコフ特性は、ある場所が一定のラベルｆ_ｉを取る確率がそれに隣接するものによってのみ決まるように決定される。概して、その確率は決定が困難であるが、ハンマースレー−クリフォード（Ｈａｍｍｅｒｓｌｅｙ−Ｃｌｉｆｆｏｒｄ）の定理により、その確率をギブズ分布を用いて「クリーク電位」Ｖ_ｃ（ｆ）に関連させ得ることが明らかである。すなわち下記式の通りである。
【００３４】
【数１１】

上記式において、Ｕ（ｆ）＝Σ_ｃ _∈Ｖ_ｃ（ｆ）である。すなわち、全クリークにわたって合計されたクリーク電位である。クリークは局所的近傍Ｎ全体にわたって考慮され、この近傍としては例えば、画素の４個の隣接画素（各画素が隣接する画素を４個のみ有すると考える場合）その他の隣接関数があると考えられる。ベイズの式では、事後確率Ｐ（Ｆ＝ｆ｜Ｘ＝ｘ）（ｘは観察されたデータである）を最大とすることが望ましい。ベイズ則を用いると、下記式のようになる。
【００３５】
【数１２】

ノイズが“ｉｉｄ”（ｉｄｅｎｔｉｃａｌｌｙａｎｄｉｎｄｅｐｅｎｄｅｎｔｌｙｄｉｓｔｒｉｂｕｔｅｄ：独立し同様に分布）であると仮定すると、確度の項は以下のように定義される。
【００３６】
【数１３】

上記式において、積は全箇所にわたるものである（すなわち、全ピクセル）。要約すると、ＭＡＰの推定は、エネルギーが下記式で表されるエネルギー最小化問題に書き換えることができる。
【００３７】
【数１４】

上記式は、ラベル構成からの寄与と得られるクリークポテンシャルからの寄与とを含んでいる。代表的には、クリークポテンシャルは、問題の事前の知見を反映するものであり、オプティカルフローの場合には、推定されたフロー場に平滑化を課すのに使用される。
【００３８】
前述のように、本実施形態は、非反復的広域最小化法を用いて、ラベル問題の解を得るものである。これは、Ｅ（ｆ）をフローグラフとして表し、そこで最大フローの計算を行うことで得られる。平均計算量を実験的に測定したところ、Ｏ（ｎ^１．１５、ｄ^１．３１）である（ｎは画像のサイズであり、ｄはラベル数である）。この環境では、クリークポテンシャルＶ（）は線形である必要があり、下記の形の平滑化項が得られる。
【００３９】
【数１５】

式中、βは解において望まれる平滑化の量を制御する比例定数である。
１．オプティカルフローについての最大フロー解
前セクションで説明したように、最大フロー計算を用いて最小化するコスト関数は以下の通りである。
【００４０】
【数１６】

最大フロー式およびＭＲＦ解釈の詳細は公知であり、上述の米国特許出願０８／９７８、８３４号ならびにロイ（Ｒｏｙ）らやイシカワ（Ｉｓｈｉｋａｗａ）らの論文（前出）などに記載されている。ＭＡＰ推定に関連する最小コストカットの広域最適性も保証されることが知られており、ボイコフ（Ｂｏｙｋｏｖ）ら（前出）やイシカワら（前出）の文献に記載されている。
【００４１】
前述のように、本実施形態におけるフロー場のパラメータ化は、（θ、ｍ）表示である。フローについての解を得るには、フロー速度分布Ｐ（υ｜Ｉ_ｄ）をそれの角度成分Ｐ（θ｜Ｉ_ｄ）および速度成分Ｐ（ｍ｜Ｉ_ｄ）に簡単に因数分解することで、上記の段落Ｂに記載の方法に従って、条件的確率Ｐ（θ｜Ｉ_ｄ）を計算する。方向フロー場θ（全画素についての方向の構成を示す）についての解を得るため、式１１は以下の形となる。
【００４２】
【数１７】

ＭＲＦ法は有限数のラベルを用いることから、θ＝［−π、π］の値の範囲は有限数の段階に区分する必要があることが明らかであろう。本実施形態を用いた実験では、段階のサイズは１°〜４°を用いた。画素の運動を区分することで、非離散的表現の場合と比較して大きい誤差を生じるように思われるかも知れないが、この実験から、それは当てはまらないことが明らかになった。
【００４３】
フロー方向についての解が得られたら次に、各画素についての速度ｍについての解を得る必要がある。速度は、フロー方向の解を求めた方法と同様にして解を得ることができる。しかしながら実際には、速度の計算は、フローの方向の計算よりかなり難しい。好ましくは、計算された方向推定値によって得られる追加データを利用することで、条件的分布Ｐ（ｍ｜Ｉ_ｄ）を修正する。それによってＰ（ｍ｜θ_ｓ、Ｉ_ｄ）が得られる（θ_ｓは、画素の方向についての解である）。そこで、運動速度を計算するためのコスト関数は、以下のようになる。
【００４４】
【数１８】

留意すべき点として、式（１２）および（１３）におけるβを、それぞれβ_１およびβ_２と表すことで、βの特定の値が任意であり、運動方向と運動速度の両方で両式において同じであっても、あるいは運動方向および運動速度について２式で異なっていても良いことを示している。
【００４５】
上記の修正により、得られる効果は大幅に向上される。それは、方向推定がフロー全体を直線に制限することで、速度の分布の不確定性を低減することで説明される。すなわち、輝度制約線方向の曖昧さがなくなっていることから、この新たな条件的分布Ｐ（ｍ｜θ_ｓ、Ｉ_ｄ）は、真のフローの速度を代表する程度がかなり高くなっている。２つの推定値（すなわち、θおよびｍ）を合わせることで、２個の画像間のオプティカルフローが得られる。
【００４６】
図４について説明すると、本実施形態の方法全体を描いたフローチャートを示してある。時系列連続画像４００が、本方法に対する入力として提供される。時系列連続画像４００は代表的には、７個以上の画像の連続ビデオ画像であるが、画像導関数の計算ができるだけの時間的密度を有する複数画像であればいかなるものであっても良い。時系列連続画像４００を、本実施形態の方法の２つの段階に対する入力として用いる。第１段階では運動方向を推定し（ステップ４０２、４０４、４０６および４０８）、第２段階では、運動速度を推定する（ステップ４１０、４１２、４１４および４１６）。運動方向の結果も運動速度を得るための段階への入力として提供されることから、運動方向を得るための段階を最初に行うのが普通である。
【００４７】
運動方向を得るための段階では、第１のフローグラフＧ_１がステップ４０２において作成される。第１のフローグラフＧ_１は、画像の空間時間導関数（式（１））を用いて作成され、コスト関数が得られ（式（１２））、それの最小値が、運動の方向成分となる。このフローグラフＧ_１は上述の米国特許出願０８／９７８、８３４号と同様の構成となっているが、式（１２）のコスト関数を用いて、エッジ容量関数（ｏｃｃ（ｕ、ｖ）＝βおよびｒｅｇ（ｕ、ｖ）＝−ｌｎ（Ｐ（θ｜Ｉ_ｄｉ）））を誘導している。次に、本実施形態の方法では、ステップ４０４において、第１のグラフＧ_１で中の最大フローの解を求め、上述の米国特許出願０８／９７８、８３４号に記載の方法と同様にして、第１のグラフＧ_１から最小カットを抜き出す。ステップ４０６では、該最小カットから運動方向を計算する。方向θ_ｉ（全画素について、ｉ∈Ｓ）は、最小カットにおける「ラベル」エッジによって直接得られる。結果として、運動方向４０８が得られ、それは運動の方向であることから、オプティカルフローの１成分を表す。
【００４８】
運動速度段階では、ステップ４１０で第２のフローグラフＧ_２が作成される。第２のフローグラフＧ_２は、画像の空間時間導関数（式（１））と前段階で計算された画素の運動方向４０８とを用いて作成される。コスト関数が得られ（式（１３））、その最小値が、Ｐ（ｍ｜θ、Ｉ_ｄ）に当てはめた場合に、運動速度成分を与える。このフローグラフＧ_２は上述の米国特許出願０８／９７８、８３４号と同様の構成となっているが、このコスト関数を用いて、エッジ容量関数（ｏｃｃ（ｕ、ｖ）＝βおよびｒｅｇ（ｕ、ｖ）＝−ｌｎ（Ｐ（ｍ｜θ_ｓｉ、Ｉ_ｄｉ）））を誘導している。次に本実施形態の方法では、ステップ４１２において、第２のフローグラフＧ_２中での最大フローの解を求め、上述の米国特許出願０８／９７８、８３４号に記載の方法と同様にして、第２のグラフＧ_２から最小カットを抜き出す。ステップ４１４では、その最小カットから運動速度を計算する。速度ｍ_ｉ（全画素について、ｉ∈Ｓ）は、最小カットにおける「ラベル」エッジによって直接得られる。結果として、運動速度４１６が得られ、それは運動速度であることから、オプティカルフローの別の成分を表す。
【００４９】
オプティカルフロー４１８は、運動方向成分４０８と運動速度成分４１６を合わせたものであることから、オプティカルフロー場全体となる。
２．２次元フローのパラメータ化
前述のように、オプティカルフローはパラメータ化されて、２個の１次元表現になる。これら２つのパラメータはできるだけ互いに独立であることが望ましい（すなわち、Ｐ（υ｜Ｉ_ｄ）＝Ｐ（ａ（υ）｜Ｉ_ｄ）Ｐ（ｂ（υ）｜Ｉ_ｄ）であって、式中ａ（υ）およびｂ（υ）はフローを表す新たな１次元パラメータである）。そこで、角度−速度表現（θ、ｍ）および速度成分（υ_ｘ、υ_ｙ）という２つの選択肢を検討した。最良の表現を決定するため、相互相関係数を実験的に測定した。多数の代表的画像導関数の場合（５００の実験）、相当する条件的分布Ｐ（υ｜Ｉ_ｄ）を得て、２つの異なるパラメータ化について相互相関係数を計算した。相互相関係数ρは以下のように定義される。
【００５０】
【数１９】

式中、Ｅは期待値であり、μは平均を示し、（ａ、ｂ）は（θ、ｍ）または（υ_ｘ、υ_ｙ）のいずれかである。ρの平均値は、（θ、ｍ）表現の場合は０．０４であり、（υ_ｘ、υ_ｙ）表現の場合は０．４である。（θ、ｍ）表現はほとんど独立であるが、（υ_ｘ、υ_ｙ）表現はそうではないことが明らかである。従って、角度−速度のパラメータ化を選択するのが適切である。
Ｄ．結果
本セクションでは、バロン（Ｂａｒｒｏｎ）らによる評価についての論文（Ｂａｒｒｏｎｅｔａｌ．， ”ＰｅｒｆｏｒｍａｎｃｅＯｆＯｐｔｉｃａｌＦｌｏｗＴｅｃｈｎｉｑｕｅｓ，” Ｉｎｔ’ｌＪｏｕｒｎａｌｏｆＣｏｍｐｕｔｅｒＶｉｓｉｏｎ，Ｖｏｌ．２，Ｎｏ．１，ｐｐ．４３−７７（１９９４））からの合成データ集合および実データ集合について本実施形態の方法を用い、さらにはその論文に記載の各種方法の結果と本実施形態の結果とを比較することで、本実施形態の方法の効果を評価する。
【００５１】
本実施形態の方法について試験を行う際には、画像導関数の計算は、空間−時間ガウスフィルター（σ＝１．５）の適用と、次に４点差演算子（１／１２）［−１、８、−８、１］の適用から成るものである。バロンら（同上）における修正ホーン−シュンク（ＨｏｒｎａｎｄＳｈｕｎｋ）アルゴリズムは、同じ導関数計算を使用するものである。ほとんどの実験に要する実行時間は、小さい画像の場合で数秒の範囲であり、高速ワークステーションでの大きい画像の場合で１０分以内である。これらの実行時間は、解にほとんど影響を与えることなく、運動パラメータについての比較的粗い離散化を用いることで、容易に短縮することができる。本セクションに示した結果はいずれも、事後処理を行わずに、本実施形態の方法によって得られた生のフロー場である。
１．合成画像
本実施形態のオプティカルフロー推定方法を、正しい結果が得られているバロンらの５種類の合成画像列について行った。この５種類の合成画像列は、バロンらの論文において、様々なアルゴリズムを比較するために用いられている画像列の例であり、それぞれ「Ｓｉｎｕｓｏｉｄ１」、「ｓｑｕａｒｅ２」、「ＴｒａｎｓｌａｔｉｎｇＴｒｅｅ」、「ＤｉｖｅｒｇｉｎｇＴｒｅｅ」、「Ｙｏｓｅｍｉｔｅ」というタイトルがつけられている。
【００５２】
本実施形態のオプティカルフロー推定方法による結果を、１００％のフロー場密度を与えるバロンらにおける５種類のアルゴリズムの結果と比較した。本実施形態は特に、高密度フロー場の推定を行うためのものであって、密度の低い場を与えるよう修正することは容易ではないことから、低密度法を直接比較することはできない。誤差の測定は、バロンらにおいて用いられている方法と同じである。２つの動き［ｕ_０、υ_０］および［ｕ_１、υ_１］の場合、誤差の測定値は、２個のベクトル［ｕ_０、υ_０、１］および［ｕ_１、υ_１、１］間の角度と定義される。
【００５３】
得られた結果を図５にまとめてある。図５において、本実施形態の結果は、最大フローとして表している。これらデータ集合に対する本実施形態の成績は常に良好である。しかしながら、これらのデータ集合はいずれも、非常に平滑な運動場を特徴とするものであって、この運動場は、運動の不連続部付近のアルゴリズムの挙動を明らかにするものではない。さらに、それにはノイズおよびその他の画素の不一致要素が含まれる。これらは、実画像についてのオプティカルフロー計算の重要な側面であり、本実施形態で特に良好に扱われるものである。
【００５４】
最も驚くべき結果は、本実施形態の方法が他のいずれの方法より数桁も優れた成績を与える「ｓｑｕａｒｅ２」に関するものである。これは、非常に低密度の導関数データが得られている場合であることから、局所的ではなく全体的に平滑化を行うことが有利であることを示すものである。本実施形態が、相関に基づくアルゴリズム（例：Ａｎａｎｄａｎ（前出）；Ｓｉｎｇｈ（前出））より常に良好な成績を与え、他のいかなる方法より大きく劣ることは決してないことが明らかであろう。
２．実画像
実際の条件下での本実施形態の成績を示すため、４種類の実画像についてのフローを調べる（図６（ａ）、図７（ａ）、図８（ａ）および図９（ａ））。これらは、良く知られているルービックキューブ（図６（ａ））、ＮＡＳＡ画像列（図７（ａ））、ハンブルグのタクシー（図８（ａ））およびＳＲＩ樹木（図９（ａ））であり、バロンらの論文（前述）でも検討されている。正しい結果が得られていないため、質的結果のみを示す。
【００５５】
ルービックキューブについての推定フロー場を図６（ｂ）に示してある。このデータ集合は、回転台上で回転するキューブについての特徴を示すものである。フローは、方向および速度のいずれにおいても、回転台およびキューブの運動にそのまま従うことがわかる。フローは、回転台の上面のようなテクスチャのない領域全体で良好に広がっている。さらに、運動の不連続部は良好に保存されている。このフロー場の詳細図が図６（ｃ）にある。図６において、存在する３種類の運動（キューブ、回転台および背景の運動）が正確に再現されている。
【００５６】
ＮＡＳＡ画像列は、カメラズームによって生じる発散フロー場の特徴を示すものである。図７（ａ）に示した画像においてカメラはズームインしている。運動速度は非常に小さく、１画素よりかなり小さいのが普通である。図７（ｂ）に示したように、フローの発散は良好に再現されている。注目すべき点として、炭酸飲料中央部における誤差はほとんどが、反射と少ない運動とが相まって生じたものと考えられる。
【００５７】
ハンブルグのタクシー画像列は、複数の独立した運動の１例である。３台の車が画像列を通じて独立に動いている。得られるフローを図８（ｂ）に示してある。車の動きは良好に再現され、良好に局所化されていることから、運動速度の簡単な閾値処理を行うことで、運動を分割することができる。これは、運動不連続部の正確な再現が必須である場合の例である。
【００５８】
ＳＲＩ樹木の画像列は、水平方向に移動するカメラについての特徴を示すものである。それは、多数の閉塞および低コントラストを特徴とするものである。カメラの動きが普通とは異なることから、運動速度は、場面の奥行きと等価である。従って、図９（ｂ）での結果は、奥行きマップとして示してある。暗い領域は運動が小さいことを示し（大きい奥行き）、明るい領域は運動が大きいことを示している（奥行きが小さい）。結果は、カメラの動きについてのデータを利用し、従って良好な性能を有すると予想される専用の立体アルゴリズムによって得られる結果に非常に近いものである。画像中央にある木の幹に沿って見られるように、奥行きの不連続部は良好に再現されている。他方、注目すべき点として、地表面の平面性が良好に保存されている。それは、高レベルの平滑化を行ないながら、しかもシャープな不連続部を再現することが可能であることを示している。従って、確率的枠組みでオプティカルフローを推定する新規な方法が提供される。簡単なノイズモデルを用いて、画像導関数の固有の不正確さを明瞭に考慮することで、全フローの確率モデルが得られている。フローをそれの角度−速度成分に分離することで、全フローが２段階で計算され、各段階は線形のクリーク電位を用いたＭＲＦのＭＡＰ推定に基づくものである。これらの推定値は最適のものであり、グラフ全体にわたる最大フローの計算によって効果的に得られる。再現されるフロー場は高密度であり、シャープな運動不連続部を保持している。注意深く確率モデルを作成することで、オプティカルフロー推定の問題に固有の大幅な誤差に対して高レベルの堅牢性を得ることができると考えられる。
【００５９】
以上、ある種の利用分野で使用するためのオプティカルフローの推定方法について説明・図示したが、本明細書に添付の請求の範囲のみによって限定される本発明の精神および広義の内容から逸脱しない限りにおいて、変更および修正が可能であることは、当業者には明らかであろう。
【００６０】
【発明の効果】
上記で説明したように本発明によれば、下記のような効果を得ることができる。
（１）複数画像間のオプティカルフローを有効かつ正確に推定することができる。
（２）輝度一定の仮定の制約を行ないながら、画像導関数の測定における誤差を適切にモデル化することができる。
（３）モデル化された環境下で、オプティカルフローに対する画像の全体的な最適解を効率良く得ることができる。
【図面の簡単な説明】
【図１】オプティカルフロー推定に対する輝度制約を示す図であり（図１（ａ））、および、図１（ａ）に図示した輝度制約に相当する従来の条件的分布Ｐ（θ｜Ｉ_ｄ ^０）を示す図（図１（ｂ））である。
【図２】図２は、３つの異なる画像テクスチャを代表する３つの異なる画像導関数を有する３つの異なる画像導関数についての法線フローおよびオプティカルフローの確率分布、すなわち、局所画像変化の程度を描いた図である。
【図３】図３は、図２で用いた３つの異なる画像導関数についてのオプティカルフローの方向および速度の確率分布を描いた図である。
【図４】本発明の方法の全体を示すフローチャートである。
【図５】本発明の方法について結果と比較した、各種合成データ集合についての各種試験アルゴリズムの結果を示す棒グラフである。
【図６】本発明についての試験を行うのに使用される回転台上で回転するキューブの連続画像のうちの１個の画像、（図６（ａ））、および図６（ａ）の画像を含む連続画像を用いて、本発明の方法によって推定されるオプティカルフロー場を示した図（図６（ｂ））、図６（ｂ）に示したオプティカルフロー場の拡大図である（図６（ｃ））。
【図７】本発明についての試験を行うのに使用される、炭酸飲料缶と各種取り合わせた対象物の連続画像中の１画像である（図７（ａ））、および図７（ａ）の画像を含む連続画像を用いて、本発明の方法によって推定したオプティカルフロー場を示した図（図７（ｂ））である。
【図８】本発明についての試験を行うのに使用される、独立に運動する複数の車の連続画像中の１画像を示した図（図８（ａ））、および図８（ａ）の画像を含む連続画像を用いて、本発明の方法によって推定されるオプティカルフロー場を示す図（図８（ｂ））である。
【図９】本発明についての試験を行うのに使用される、カメラが画像を横切って水平方向に移動する、樹木の連続画像中の１画像を示す図（図９（ａ））、および図９（ａ）の画像を含む連続画像を用いて、本発明の方法によって推定されるオプティカルフロー場を示す奥行きマップを示す図（図９（ｂ））である。
【符号の説明】
４００時系列連続画像
４０２、４０４、４０６ステップ
４０８運動方向
４１０、４１２、４１４ステップ
４１６運動速度

Claims

複数の画像間でのオプティカルフローを推定するオプティカルフロー推定方法であって、
（ａ）複数の画像の空間時間導関数を用い、

［式中、Ｓは全画素集合を示し、Ｎ _i は画素ｉに隣接する全画素集合を示し、β ₁ は負ではない任意の平滑化定数を示し、θ _i は画素ｉの方向を示し、Ｉ _d は測定画像導関数を示し、Ｐ（θ｜Ｉ _di ）は画像導関数がＩ _di の場合の方向θの条件的確率を示している。］
で示されるコスト関数からエッジ容量関数を誘導することで第１のフローグラフＧ₁を作成するステップと、
前記第１のフローグラフＧ₁中の第１の最大フローについて解を求めることで、それから第１の最小カットを得るステップと、
前記第１の最小カットから運動方向成分を計算するステップとを有する、
運動方向成分を得るステップと、
（ｂ）前記複数の画像の空間時間導関数および前記運動方向成分を用いて、

［式中、Ｓは全画素集合を示し、Ｎ _i は画素ｉに隣接する全画素集合を示し、β ₂ は負ではない任意の平滑化定数を示し、ｍ _i は画素ｉの速度を示し、Ｉ _d は測定画像導関数を示し、Ｐ（ｍ｜θ _si 、Ｉ _di ）は方向が既知のθ _si であって、画像導関数がＩ _di の場合の速度ｍの条件的確率を示す。］
で示されるコスト関数からエッジ容量関数を誘導することで第２のフローグラフＧ₂を作成するステップと、
前記第２のフローグラフＧ₂中の第２の最大フローについて解を求めることで、それから第２の最小カットを得るステップと、
前記第２の最小カットから運動速度成分を計算するステップとを有する、
運動速度成分を得るステップとを有し、
前記運動方向成分および前記運動速度成分とを組み合わせて、複数画像間のオプティカルフローを推定するオプティカルフロー推定方法。
隣接する画素が４個である、請求項１記載のオプティカルフロー推定方法。
条件的確率分布Ｐ（θ｜Ｉ_di）が［Ｐ（θ｜Ｉ_d ⁰）・Ｐ（Ｉ_d ⁰｜Ｉ_d）］（式中、Ｐ（θ｜Ｉ_d ⁰）は運動方向のモデルを示し、Ｐ（Ｉ_d ⁰｜Ｉ_d）は画像導関数の測定における誤差のモデルを示す）であり、さらに

および

である請求項１記載のオプティカルフロー推定方法。
条件的確率分布Ｐ（θ｜Ｉ_di）と輝度一定の仮定とを組み合わせて、条件的確率分布Ｐ（ｍ｜θ_si、Ｉ_di）を得る、請求項１記載のオプティカルフロー推定方法。