(57)【要約】【課題】高速かつ高精度に画素のラベル値を選択する。【解決手段】Ｗ（ｘｐ−ｘｑ＋ｒ）とＶ（ｘｐ）を含むエネルギー関数を評価する各評価関数ｍ（ｔ）ｐ→ｑ（ｘｑ）のｔ＝０の初期値を設定する手段２０２と、ｈ（ｘｐ）を設定する手段２０３と、Ｗが最小値になるｘｑであるＰｅａｋＸｑを算出する手段２０５と、ｘｐごとにｈとＷ（ｘｐ−ＰｅａｋＸｑ＋ｒ）との集約値を算出する手段２０６と、前記集約値と前記ｍ（ｔ）にｘｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’（ｔ）とする手段２０７と、ＷとＷ（ｘｐ−（ｘｑ＋ｋ）＋ｒ）との差分とｍ’（ｔ）ｐ→ｑ（ｘｑ＋ｋ）との集約値と、ｍ’（ｔ）との小さい値の方を最適値としてｍ（ｔ＋１）とする手段２０７と、ｘｑごとにビリーフ値ｂｑ（ｘｑ）を算出しｂｑを最小にするラベル値ｘｑを選択する手段２１３を具備する。【選択図】図２
2008225778.tif 000002
【特許請求の範囲】
【請求項１】
画像の各画素に対応付けられた各ノードｐ、ｑのそれぞれのラベル値ｘｐ、ｘｑを変数とし、ノード対（ｐ、ｑ）およびｐとｑに関連付けられた定数ｒを用いて定義された関数Ｗ（ｘｐ−ｘｑ＋ｒ）と関数Ｖ（ｘｐ）を含む、最小化するｘｐを求めるエネルギー関数を設定する設定手段と、
各ノード対（ｐ、ｑ）に対応した各評価関数ｍ（ｔ）ｐ→ｑ（ｘｑ）を格納する格納手段と、
前記各評価関数のｔ＝０での値である初期値を設定する設定手段と、
【数１】
2008225778.tif 000027
にしたがってｈ（ｘｐ）を設定する設定手段と、
ｘｐごとにＷ（ｘｐ−ｘｑ＋ｒ）が最適値になるｘｑであるＰｅａｋＸｑを算出する算出手段と、
ｘｐごとにｈ（ｘｐ）とＷ（ｘｐ−ＰｅａｋＸｑ＋ｒ）との和または積である集約値を算出する算出手段と、
前記集約値と前記ｍ（ｔ）ｐ→ｑ（ｘｑ）にｘｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’（ｔ）ｐ→ｑ（ｘｑ）とする更新手段と、
ｘｐごとに、Ｗ（ｘｐ−ｘｑ＋ｒ）とＷ（ｘｐ−（ｘｑ＋ｋ）＋ｒ）との差分とｍ’（ｔ）ｐ→ｑ（ｘｑ＋ｋ）との集約値と、ｍ’（ｔ）ｐ→ｑ（ｘｑ）との小さい値の方を最適値としてｍ（ｔ＋１）ｐ→ｑ（ｘｑ）とする更新手段と、
ｔが繰り返し回数Ｔよりも小さい場合には、ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）を計算させ、最終的にｍ（Ｔ）ｐ→ｑ（ｘｑ）を取得する取得手段と、
ｘｑごとに
【数２】
2008225778.tif 000028
を算出する算出手段と、
ｑごとに、ｂｑ（ｘｑ）を最小にするラベル値ｘｑを選択する選択手段と、を具備することを特徴とする画像処理装置。
【請求項２】
前記Ｗ（ｘｐ−ｘｑ＋ｒ）は、ｓをノード対（ｐ、ｑ）に関連付けられた定数として、
Ｗ（ｘｐ−ｘｑ＋ｒ）＝ｓ‖ｘｐ−ｘｑ＋ｒ‖
であることを特徴とする請求項１に記載の画像処理装置。
【請求項３】
精度値δの初期値を設定する設定手段と、
ノードｐごとに精度値δに対応するラベル値の候補を設定する設定手段と、
ノード対（ｐ、ｑ）に対し、前記候補に対応するＷ（ｘｐ−ｘｑ＋ｒ）を設定関数として設定する設定手段と、
前記設定関数を使用して、前記選択手段がラベル値を選択するように制御する制御手段と、
精度が上がるように精度値δをδ’に変更する変更手段と、
精度値δ’がある値になるまで、精度値δ’を新たな初期値としてラベル値を選択するように制御する制御手段と、を具備することを特徴とする請求項１または請求項２に記載の画像処理装置。
【請求項４】
前記ラベル値は、第１の画像に対する第２の画像の対応点の位置を示すことを特徴とする請求項１から請求項３のいずれか１項に記載の画像処理装置。
【請求項５】
前記ノードごとに、関連付けられるラベル値の数が異なることを特徴とする請求項１から請求項４のいずれか１項に記載の画像処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、予め定義されたエネルギー関数の最小化を効率的に計算する技術に関し、例えば、動画像の動きベクトルの高速な算出に利用できる画像処理装置に関する。
【背景技術】
【０００２】
超解像のように高い精度の動きベクトルを必要とするアプリケーションでは、従来のようにローカルの情報だけを用いて動きベクトルを算出するだけでは不十分である。一方、Belief Propagation(ＢＰ)という手法を用いて画面全体の整合性を考えながら動きベクトルを算出すれば、従来よりも信頼性の高い動き推定が可能である。
【０００３】
しかし、ＢＰは、隣接画素対に対してメッセージ値という値を繰り返し計算する最適化手法で、非常に計算量が多いことが知られている。ＢＰによれば、例えば、動き推定では数時間〜数日かかり、消費メモリは８Ｇバイトを超えてしまう。その対策として、全ての画素について動きベクトルの範囲が同じで、評価関数が特定の形の場合のみ適用可能な高速ＢＰが提案された（例えば、非特許文献１参照）。
【非特許文献１】P. F. Felzenszwalb and D. R. Huttenlocher, ``Efficient Belief Propagation for Early Vision,’’ in Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Volume 1, pp.261-268, 2004.
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、依然として計算量や消費メモリが多く、例えば、計算時間は数十分で、消費メモリはＢＰと同じである。
【０００５】
一方、ブロックマッチングにおいては、粗い動きベクトルを求めその周辺のみを探索する階層的探索や、他のフレームの動きベクトルを基準としてその周辺のみを探索する手法などの高速化手法が用いられている。これらの高速化手法と高速ＢＰを組み合わせることができれば、高速で精度の高い動き推定が可能になる。
【０００６】
しかし、階層的探索などの高速化手法を利用すると、画素ごとに動きベクトルの範囲が異なるために高速ＢＰが適用できなくなる。
【０００７】
ところで、高速化ＢＰは、ここで示した動きベクトルの算出の他に、例えば、ステレオ画像の視差の算出、画像の復元に利用される手法として、多次元ベクトル空間上で予め定められた点の集合のうち、予め定義されたエネルギー関数を最小化する点を選択する手法である。これら動きベクトル、視差、奥行き情報は画素のラベル値で決定される。
【０００８】
この発明は、上述した事情を考慮してなされたものであり、高速かつ高精度に画素のラベル値を選択する画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上述の課題を解決するため、本発明の画像処理装置は、画像の各画素に対応付けられた各ノードｐ、ｑのそれぞれのラベル値ｘｐ、ｘｑを変数とし、ノード対（ｐ、ｑ）およびｐとｑに関連付けられた定数ｒを用いて定義された関数Ｗ（ｘｐ−ｘｑ＋ｒ）と関数Ｖ（ｘｐ）を含む、最小化するｘｐを求めるエネルギー関数を設定する設定手段と、各ノード対（ｐ、ｑ）に対応して各評価関数ｍ（ｔ）ｐ→ｑ（ｘｑ）を格納している格納手段と、前記各評価関数のｔ＝０での値である初期値を設定する設定手段と、
【数３】
2008225778.tif 000003
【００１０】
にしたがってｈ（ｘｐ）を設定する設定手段と、ｘｐごとにＷ（ｘｐ−ｘｑ＋ｒ）が最小値になるｘｑであるＰｅａｋＸｑを算出する算出手段と、ｘｐごとにｈ（ｘｐ）とＷ（ｘｐ−ＰｅａｋＸｑ＋ｒ）との和または積である集約値を算出する算出手段と、前記集約値と前記ｍ（ｔ）ｐ→ｑ（ｘｑ）にｘｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’（ｔ）ｐ→ｑ（ｘｑ）とする更新手段と、ｘｐごとに、Ｗ（ｘｐ−ｘｑ＋ｒ）とＷ（ｘｐ−（ｘｑ＋ｋ）＋ｒ）との差分とｍ’（ｔ）ｐ→ｑ（ｘｑ＋ｋ）との集約値と、ｍ’（ｔ）ｐ→ｑ（ｘｑ）との小さい値の方を最適値としてｍ（ｔ＋１）ｐ→ｑ（ｘｑ）とする更新手段と、ｔが繰り返し回数Ｔよりも小さい場合には、ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）を計算させ、最終的にｍ（Ｔ）ｐ→ｑ（ｘｑ）を取得する取得手段と、
ｘｑごとに
【数４】
2008225778.tif 000004
【００１１】
を算出する算出手段と、ｑごとに、ｂｑ（ｘｑ）を最小にするラベル値ｘｑを選択する選択手段と、を具備することを特徴とする。
【発明の効果】
【００１２】
本発明の画像処理装置によれば、高速かつ高精度に画素のラベル値を選択することが可能になる。
【発明を実施するための最良の形態】
【００１３】
以下、図面を参照しながら本発明の実施形態に係る画像処理装置について詳細に説明する。  
まず、ラベル値の例として動きベクトルを例に挙げ、画素ごとに動きベクトルの範囲が異なっている場合でも、高速ＢＰが適用できるようにメッセージ値の計算が修正できる一例をごく簡単に示す。
【００１４】
一例として、第１画素の動きベクトルｐと第２画素の動きベクトルｑの評価関数（後述するメッセージ値に対応）がｆ（ｑ）＋｜ｐ−ｑ｜という形のときに、各ｑに対して評価関数の最小値を求める計算を考える。なお、ここでは、簡単のため動きベクトルがスカラーであるとする。
【００１５】
高速ＢＰでは、各ｑに対して｜ｐ−ｑ｜＝０のときの評価関数（すなわち、ｆ（ｑ））を予め算出し、その後｜ｐ−ｑ｜が非０の場合を差分計算により逐次求めながら評価関数の最小値を求める。
【００１６】
しかし、高速ＢＰでは、ｑの範囲内では｜ｐ−ｑ｜＝０とならない場合が考慮されていない。例えば、ｐ＝−２，−１，０，１，２、ｐ＝−１，０，１，２，３のようにｐとｑの値の範囲が異なる場合、ｐ＝−２，−１，０，１，２、ｐ＝−３，−２，−１，０，１，２，３のように候補数が異なる場合、ｐ＝−２，−１，０，１，２、ｐ＝−２．５，−１．５，−０．５，０．５，１．５のようにｐとｑの差が小数になる場合は考慮されていない。
【００１７】
一方、本実施形態の画像処理装置による手法（拡張高速ＢＰと呼ぶ）では、ｑの範囲内では｜ｐ−ｑ｜＝０とならない場合（ｑの範囲外で最小値を持つ評価関数の場合）については、ｑの範囲境界における最小値を算出してから、高速ＢＰと同様に差分計算による逐次評価を行う。こうすることで階層的探索と拡張高速ＢＰを併用することができ、高速かつ高精度に動きベクトルを推定できるようになる。
【００１８】
（ＢＰ）  
次に、ＢＰについて説明する。  
多次元ベクトル空間上で予め定められた点の集合のうち、予め定義されたエネルギー関数を最小化する点を選択する技術は、多くの応用を持ち、例えば、動画像の動きベクトルの算出、ステレオ画像の視差の算出、画像の復元に利用されている（上記のP. F. Felzenszwalb and D. R. Huttenlocherの文献参照）。一般的なエネルギー関数の厳密な最小化は極めて困難であることが知られている。その近似解の算出手法の１つとして、ビリーフプロパゲーション（ＢＰ：Belief Propagation）と呼ばれる技術が知られている。ＢＰでは、高々２つの変数から成る項の総和として表現される次のエネルギー関数の近似的な最小化を考える（以下、近似的な最小化を単に最小化と呼ぶ）。  
【数５】
2008225778.tif 000005
【００１９】
ラベルの添え字ｉをノード番号、あるいは単にノードと呼ぶ。ｘｉ＝｛０，１，・・・，Ｌｉ−１｝はラベル、Ｎはエネルギー関数が含む各２ノード項を構成する変数対の集合をあらわす。また、各２ノード項を構成する（つまり、Ｎに属する）ノード対をエッジと呼ぶ。ＢＰでは、エネルギー関数を最小化するために次の繰り返し演算を行う。
【００２０】
（ステップ１）ｔ＝０とする（ｔはメッセージ更新回数を表す）。各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ（ｔ）ｐ→ｑ（ｘｑ）の初期値を与える（あるいは全てのｍ（ｔ）ｐ→ｑ（ｘｑ）を０に初期化する）。
【００２１】
（ステップ２）各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）を次のメッセージ更新式により更新する。  
【数６】
2008225778.tif 000006
【００２２】
ただし、Σsは（ｓ，ｐ）∈Ｎかつｓ≠ｑを満たす全てのｓに対する和を表す。なお、ｍ（ｔ）ｐ→ｑ（ｘｑ）とｍ（ｔ）ｑ→ｐ（ｘｑ）は別のメッセージ値であり、区別する必要がある。
【００２３】
（ステップ３）ｔを１増やす。ｔが予め決めておいた繰り返し回数Ｔより小さいなら、ステップ２に戻る。
【００２４】
（ステップ４）各ｑに対して、次のビリーフと呼ばれる値を求める。  
【数７】
2008225778.tif 000007
【００２５】
ただし、Σｐは（ｐ，ｑ）∈Ｎを満たす全てのｐに対する和を表す。
【００２６】
（ステップ５）各ｑに対してｂｑ（ｘｑ）を最小とするラベルｘｑ＝｛０，１，・・・，Ｌｑ−１｝を選ぶ。
【００２７】
ＢＰでは上記のステップ２でメッセージ更新式の繰り返し演算を行うために多くの計算時間が必要となる。
【００２８】
（高速ＢＰ）  
次に、高速ＢＰについて説明する。高速ＢＰは、メッセージ更新式の繰り返し演算を行うための計算時間を少なくする。
【００２９】
高速ＢＰでは、全てのノードｉに対してラベルを同一の範囲０〜（Ｌ−１）とし、かつ、２ノード項Ｗを以下に示すＰｏｔｔｓや線形コストのいずれかの形に制限し、上記のステップ２の繰り返し演算の高速化を実現する。  
Ｐｏｔｔｓ   Ｗ（ｘｐ，ｘｑ）＝０（if ｘｐ＝ｘｑ），ｄ（otherwise）  
線形コスト   Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ‖ｘｐ−ｘｑ‖，ｄ）  
だだし、‖Ａ‖はＡのノルムを示す。ここで、ｐに依存する項をまとめて  
【数８】
2008225778.tif 000008
【００３０】
で表すと、メッセージ更新式を次の式で書き換えられる。  
【数９】
2008225778.tif 000009
【００３１】
高速ＢＰではこの式をもとにメッセージ更新式の算出を高速化する。
【００３２】
＜ＷがＰｏｔｔｓの場合＞    
２ノード項ＷがＰｏｔｔｓの場合、メッセージ更新式は次の式になる。  
【数１０】
2008225778.tif 000010
【００３３】
したがって、はじめに  
【数１１】
2008225778.tif 000011
【００３４】
を計算しておけば、残りの項はｐに依存しないからｍ（ｔ）ｐ→ｑ（ｘｑ）はｑのラベル数のオーダで計算できる。
【００３５】
＜Ｗが線形コストの場合＞    
２ノード項Ｗが線形コストの場合、はじめに２ノード項がＷ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ‖の場合について計算し、得られた結果とｄとのｍｉｎをＰｏｔｔｓと同じ手法で計算すれば良い。２ノード項がＷ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ‖の場合についての計算では、メッセージ更新式は次の通りである。  
【数１２】
2008225778.tif 000012
【００３６】
メッセージ更新式を高速に計算するために、メッセージ更新は次の３つのステップに従って行う。この計算は「in-place」更新により行う、つまりｍ（ｘｑ）の値を各ステップで上書きしていくことで実現する。  
（高速ＢＰ更新ステップ１）ラベルｘｑ＝０〜（Ｌ−１）に対し、ｍ（ｘｑ）←ｈ（ｘｑ）を計算する。  
（高速ＢＰ更新ステップ２）１〜（Ｌ−１）までラベル番号ｘｑを１ずつ増やしながら以下を実行する。  
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ−１）＋ｓ）  
（高速ＢＰ更新ステップ３）（Ｌ−２）〜０までラベル番号ｘｑを１ずつ減らしながら以下を実行する。
【００３７】
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ＋１）＋ｓ）  
高速ＢＰ更新ステップ３が終了した時点で、各ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）の値がｍ（ｘｑ）に得られる。
【００３８】
以上のように、全てのノードに対してラベルを同一の範囲の整数０〜（Ｌ−１）とし、かつ２ノード項ＷがＰｏｔｔｓや線形コストのいずれかであれば、高速化によりＢＰの計算時間を大幅に短縮できる。
【００３９】
（高速ＢＰの問題点）  
高速ＢＰでは、全てのノードに対するラベル数が同一の範囲の整数であることを前提としている。しかし、全てのノードに対するラベル数が同一の範囲の整数でなければならないと仮定してしまうと、例えば、次のような場合には適用できない。
【００４０】
＜ケース１＞ノードｐのラベルがｘｐ＝｛０，１，２，３，４｝、ノードｑのラベルがｘｑ＝｛２，３，４，５，６｝というように、ノードごとにラベルの範囲が異なる場合。  
＜ケース２＞ノードｐのラベルがｘｐ＝｛０，１，２，３，４｝、ノードｑのラベルがｘｑ＝｛０，１，２，３，４，５，６｝というように、ノードごとにラベルの種類数が異なる場合。  
＜ケース３＞ノードｐのラベルに対応する数値が小数を含む場合。例えば２ノード項をＷ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ‖の形で書いたときに、各ラベルに対応する数値がｘｐ＝｛０．３，１．３，２．３，３．３，４．３｝、ｘｑ＝｛０，１，２，３，４，５，６｝となる場合。
【００４１】
以上に示したケースは多くの応用で起こりうる。例えば、動画像の各画素をノード、動きベクトルをラベルとして最適なラベルを選択する問題を考えると、それぞれのケースは下記の状況で起こる。  
＜ケース１＞まず、動きベクトルを低い精度で求め、次に得られた動きベクトルを中心とした小さな範囲で動きベクトルをより高い精度で求める場合を考える。このとき、２回目の動きベクトルの算出では、ノードごとにラベルの範囲が異なる場合を扱う。
【００４２】
＜ケース２＞動画像の動きベクトル推定では、前のフレームからの動きベクトルの変化は比較的小さいと考えられるから、前のフレームで求めた動きベクトルを中心とした小さな範囲のみを次のフレームの動きベクトル候補とすれば高速化できる。このとき、画像中で動きベクトルの変化（変化とは、例えば、過去２フレームの動きベクトルから求めた加速度）が穏やかな領域では、誤推定の可能性は少ないため動きベクトル候補を少なくすれば余計な計算を減らせるが、動きベクトルの変化が激しい領域では、誤推定の可能性を減らすために動きベクトル候補を多くする必要がある。このとき、ノードごとにラベルの種類数が異なる場合を扱う。
【００４３】
＜ケース３＞ケース２と同様に前のフレームで求めた動きベクトルを中心とした小さな範囲を動きベクトル候補とし、ケース１のように動きベクトルの推定を多段構成で行う場合を考える。前のフレームでは動きベクトルを高い精度で求めているから、次のフレームで低い精度で求めるときには、中心とする動きベクトルは小数を含むことになる。このとき、各ノードのラベルとして小数を扱うことになる。
【００４４】
（拡張高速ＢＰ）  
（用語の定義、基本部分）  
本実施形態では、高々２つの変数（ノード）から成る項の総和として表現されるエネルギー関数の近似的な最小化問題を扱う。本実施形態では近似的な最小化を含めて単に最小化と呼んでいる。エネルギー関数の最小化とは、エネルギー関数を最小にするラベルを割り当てることをいう。１変数（１ノード）から成る項をＶ、２ノード（２ノード）から成る項をＷとすると、エネルギー関数は次の式で表せる。  
【数１３】
2008225778.tif 000013
【００４５】
ただし、ラベルの添え字ｉやｊはノードを表し、ｘｉ＝｛０，１，・・・，Ｌｉ−１｝はラベルを表し、整数であることが多い。Ｎはエネルギー関数が含む各２ノード項を構成する変数対の集合をあらわす。なお、この和は積の形で書かれていても良い。関数が積で与えられる場合の解き方は後述する。解くとは、エネルギー関数を最小化するようにラベルを割り当てることをいう。また便宜上、和と積をあわせて集約と呼ぶ。エネルギー関数の最小化においては、上述したように、繰り返し行われるメッセージ値の更新の計算を効率的に行うことが重要である。なお、上記の符号を反転したエネルギー関数  
【数１４】
2008225778.tif 000014
【００４６】
の最大化問題は、符号を変えれば同じ手法で解ける。また、便宜上、最小値を選択する場合の最小値と最大値を選択する場合の最大値をあわせて最適値と呼ぶことにする。
【００４７】
≪ケース１，２：Ｗ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖、Ｌｉが画素ごとに異なる場合≫  
各２ノード項がＷ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖（この一般形をＷ（ｘｐ−ｘｑ＋ｒ）と書くこともある）で与えられ、ラベル数Ｌｉが画素ごとに異なる場合を考える。ｒｐｑは定数で、２ノード項ごとに異なっていても良いものとする。このＷ、Ｌｉに対するメッセージ値は、上記の高速ＢＰの手法では算出できない。
【００４８】
本実施形態によれば、次の手法でメッセージ値を更新するため、２ノード項が前述のＷ（ｘｐ，ｘｑ）＝ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖で与えられる場合でも算出できる。このための手法について図１から図５を参照して説明する。  
本実施形態の画像処理装置の本質的な手法である更新について、はじめにｒｐｑが整数の場合について説明する。ｒｐｑが整数でない場合について、および本実施形態の画像処理装置全体の流れについては後述する。なお、ラベルｘｐ、ｘｑがそれぞれ０〜３、０〜４でｒｐｑ＝−１の場合の例を図３から図５に示す。図３から図５において、縦方向はｍの値を表すものとする。
【００４９】
まず、本実施形態の画像処理装置の更新を行う更新処理装置部分について図２を参照して簡単に説明する。  
更新処理装置部分は、メッセージ値保持部２０１、メッセージ初期値入力部２０２、ｈ（ｘｐ）値算出部２０３、ｘｑ一時メッセージ値保持部２０４、ＰｅａｋＸｑ算出部２０５、集約値算出部２０６、最適値選択部２０７、最適値出力部２０８、ノード対（ｐ，ｑ）選択部２０９、ｘｐ値入力部２１０、ｘｑ値入力部２１１、ｘｑビリーフ値算出部２１２、出力部２１３、ノードｑ選択部２１４を含んでいる。
【００５０】
メッセージ値保持部２０１は、各２ノード対（ｐ、ｑ）に対し、ｐからｑへのメッセージ値ｍｐ→ｑ（ｘｑ）を保持、すなわち格納している。  
メッセージ初期値入力部２０２は、各２ノード対（ｐ、ｑ）に対し、ｐからｑへのメッセージ値ｍｐ→ｑ（ｘｑ）の初期値を設定する。  
ｈ（ｘｐ）値算出部２０３は、２ノード対（ｐ、ｑ）のノードｑに対するラベル値ｘｐを入力とし、更新後のｍｐ→ｑ（ｘｑ）に含まれる項であってノードｐに依存しノードｑに依存しない項（ｈ（ｘｐ））を、ｍｐ→ｑ（ｘｑ）の初期値を用いて算出する。  
ｘｑ一時メッセージ値保持部２０４は、ｘｑ値入力部２１１で取得されたｘｑごとのメッセージ値を一時的に格納する。  
ＰｅａｋＸｑ算出部２０５は、与えられたｘｐに対し関数Ｗ（ｘｐ−ｘｑ＋ｒ）が最適値をとるｘｑであるＰｅａｋＸｑを算出する。  
集約値算出部２０６は、ラベル値ｘｐを入力として、ｈ（ｘｐ）とＷ（ｘｐ−ＰｅａｋＸｑ＋ｒ）を集約した値を算出する。  
最適値選択部２０７は、集約値算出部２０６が算出した集約値と、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が保持する値から最適値を選択し、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が格納している値を更新する。  
最適値出力部２０８は、ラベル値ｘｑおよび差分定数ｋを入力として、関数Ｗ（ｘｐ−ｘｑ＋ｒ）と関数Ｗ（ｘｐ−（ｘｑ＋ｋ）＋ｒ）との差分をｍｐ→ｑ（ｘｑ＋ｋ）に集約した値と、ｍｐ→ｑ（ｘｑ）とを比較し、その最適値でｍｐ→ｑ（ｘｑ＋ｋ）を更新する。  
ノード対（ｐ，ｑ）選択部２０９は、処理すべきノード対（ｐ、ｑ）を選択する。  
ｘｐ値入力部２１０は、処理すべきノード対（ｐ、ｑ）のノードｐに対するラベル値ｘｐを選択する。  
ｘｑ値入力部２１１は、ノード対（ｐ，ｑ）選択部２０９、ノードｑ選択部２１４で選択されたｑに対応するｘｑ値を取得する。  
ｘｑビリーフ値算出部２１２は、ノードｑ選択部２１４からノードｑ、ｘｑ値入力部２１１からラベル値ｘｑ、メッセージ値保持部２０１からメッセージ値を入力し、ラベル値ｘｑに対するｑへのメッセージ値を集約した値であるビリーフ値を算出する。  
出力部２１３は、ノードｑに対し、各ラベル値ｘｑのビリーフ値を算出し、それらのうち最適値を与えるラベル値をノードｑのラベル値として出力する。  
ノードｑ選択部２１４は、処理すべきノードｑを選択する。
【００５１】
（拡張高速ＢＰ更新ステップ１）まず入力としてエネルギー関数が与えられ（Ｓ１０１）、メッセージ初期値入力部２０２がメッセージ値を初期化（ステップＳ１０２）する。次に、ラベルｘｑ＝｛０，１，・・・，Ｌｑ−１｝に対し、ｍ（ｘｑ）←＋∞とするステップを実行する。このステップは後述の方法を使えばなくても良いためフローチャートでは示していないが、実行する場合は例えばステップＳ１０２とステップＳ１０３の間で実行する。ｍ（ｘｑ）はｘｑ一時メッセージ値保持部２０４で保持され、処理ラベルはｘｑ値入力部２１１で制御される。
【００５２】
（拡張高速ＢＰ更新ステップ２）ノード対（ｐ，ｑ）選択部２０９がノード対を選択し、ｘｐ値入力部２１０がｘｐを取得し、ｈ（ｘｐ）値算出部２０３がｈ（ｘｐ）を算出する（ステップＳ１０３）。ＰｅａｋＸｑ算出部２０５が、ｘｐ値入力部２１０で得られるｘｐに対し、関数ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖が最適値をとるｘｑである下記の式で示されるｙｐ（ｘｑ）（ＰｅａｋＸｑと呼ぶ）を算出する（ステップＳ１０４）。集約値算出部２０６がｈ（ｘｐ）＋ｓ‖ｘｐ−ＰｅａｋＸｑ＋ｒｐｑ‖を計算する（ステップＳ１０５）。最適値選択部２０７がｍｉｎ（ｍ（ｙｐ（ｘｐ）），ｈ（ｘｐ）＋ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖）を計算し、最適値選択部２０７が、この値と、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が保持する値とから最適値を選択し、ｘｑ一時メッセージ値保持部２０４に出力し、ｍｐ→ｑ（ｘｑ）を更新する（ステップＳ１０６）。
【００５３】
すなわち、各ラベルｘｐに対して、次の式を実行する。  
ｍ（ｙｐ（ｘｐ））←ｍｉｎ（ｍ（ｙｐ（ｘｐ）），ｈ（ｘｐ）＋ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖）
なお、ｙｐ（ｘｐ）はラベルｘｐごとにラベルｘｑを選択するためのもので、
【数１５】
2008225778.tif 000015
【００５４】
により求める（なお、この式はｘｐではなくｘｑを決める式である）。図１や図２ではｙｐ（ｘｐ）をＰｅａｋＸｑと記載している。
【００５５】
各ｙｐ（ｘｐ）（＝ＰｅａｋＸｑ）についてｍの値を求めると図３のようになる。縦方向はｍの値の大小を表す。拡張高速ＢＰ更新ステップ２ではｍｉｎを算出する操作を行う。この操作は図４に示すように、ラベルｘｑの範囲の制限により両端では２つ以上のｙｐ（ｘｐ）が同じラベル値となることがあるが、このときｍ（ｘｑ）として最小値を選択する操作を表す。なお、この操作によって、ただ１つのｙｐ（ｘｐ）が対応するｘｑについてはその値がｍ（ｘｑ）として用いられ、対応するｙｐ（ｘｐ）が１つも存在しないｘｑについてはｍ（ｘｑ）に＋∞が設定されることになる。
【００５６】
（拡張高速ＢＰ更新ステップ３）ｘｑ＝｛０，１，・・・，Ｌｑ−１｝の範囲で、ラベル番号ｘｑを１から開始して１ずつ増やしながら以下を実行する。  
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ−１）＋ｓ）  
この操作は、図５に示すように、隣接ラベル値でのｍ（ｘｑ）に差分を足した値と注目ラベル値でのｍ（ｘｑ）のうちの小さいほうをｍ（ｘｑ）として選択する操作を表す。
【００５７】
（拡張高速ＢＰ更新ステップ４）ｘｑ＝｛Ｌｑ−２，・・・，０｝の範囲で、ラベル番号ｘｑをＬｑ−２から開始して１ずつ減らしながら以下を実行する。  
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ＋１）＋ｓ）  
なお、拡張高速ＢＰ更新ステップ３と４の順序は逆でもかまわない。
【００５８】
図１では拡張高速ＢＰ更新ステップ３と４をまとめてＳ１０７としている。Ｓ１０７は拡張高速ＢＰ更新ステップ３や４の処理を一般化して書いたもので、最適値出力部２０８で実行される。実行される処理は次の通りである。  
（ステップＳ１０７）ラベル差分をｋとする。ｋが正ならｋからＬｑ−１まで、ｋが負ならＬｑ−１＋ｋから０までの範囲で、１（ｋ＞０）あるいは−１（ｋ＜０）ずつラベル値をずらしながら、以下を実行する。  
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ＋ｋ）＋ΔＷ）  
ここで、ΔＷはラベルをｋずらすことにより生じるＷの差分を表す。拡張高速ＢＰ更新ステップ３はラベル差分ｋ＝１としてＳ１０７を実行する場合、拡張高速ＢＰ更新ステップ４はラベル差分ｋ＝−１としてＳ１０７を実行する場合に相当する。最小値の選択をＡ３０８で行い、選択された最小値はＡ３０５に送られる。
【００５９】
本実施形態の画像処理装置を最も特徴づけるステップは拡張高速ＢＰ更新ステップ２である。  
はじめにｈ（ｘｑ）をラベルごとに求める必要があるが、このときｍ（ｘｑ）に単に求めたｈ（ｘｑ）を代入するのではなく、まず各２ノード項が最小になるラベル値ｙｑ（ｘｑ）を調べ、そのラベル値に対応するｍ（ｙｑ（ｘｑ））を更新する。このとき、複数のｘｑに対してｙｑ（ｘｑ）が同一となる場合にはその最小値を選択しておけば、拡張高速ＢＰ更新ステップ３と４を終えた時点での更新メッセージ値は高速化を施さない場合と完全に一致する。
【００６０】
なお、拡張高速ＢＰ更新ステップ１でメッセージ値を無限大としておくのは、何もしない拡張高速ＢＰ更新ステップ２でｍ（ｘｑ）の値が割り当てられないラベルについて、拡張高速ＢＰ更新ステップ３と４として高速ＢＰと同じ手法で正しいメッセージ値が得られるようにするためである。拡張高速ＢＰ更新ステップ３と４では、Ｌの値がノードに依存して変化する点が高速ＢＰの手法と異なるが、高速ＢＰでの手法と回路やプログラムを共用できる。  
更新ステップ４が終了した時点で、各ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）の値がｍ（ｘｑ）に得られる。ｍ（ｘｑ）を保持するメモリはｍ（ｔ＋１）ｐ→ｑ（ｘｑ）を保持するメモリと同一のものでも良い。エネルギー関数の符号が反転している場合には最大値を選択する。
【００６１】
次に、この拡張高速ＢＰ更新ステップ１〜４を用いた全体の流れについて説明する。  
（最適化ステップ１）入力としてエネルギー関数が与えられる（Ｓ１０１）。ｔ＝０とする（ｔはメッセージ更新回数を表す）。各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ（ｔ）ｐ→ｑ（ｘｑ）の初期値を与えるか、全て０で初期化する（Ｓ１０２）。ブロック図では、メッセージ初期値入力部２０２が初期値を与え、メッセージ値保持部２０１でそれを保持している。
【００６２】
（最適化ステップ２）各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）を拡張高速ＢＰ更新ステップ１〜４により更新する（Ｓ１０８）。ブロック図では、ノード対（ｐ，ｑ）選択部２０９で処理するノード対を選択している。なお、拡張高速ＢＰステップ２の各ノード対の処理は独立しているため並列に実行できる。
【００６３】
（最適化ステップ３）ｔを１増やす。ｔが予め決めておいた繰り返し回数Ｔより小さいなら、拡張高速ＢＰステップ２に戻る。
【００６４】
（最適化ステップ４）各ｑ、各ラベルに対して、次のビリーフと呼ばれる値を求める（Ｓ１０９）。
【数１６】
2008225778.tif 000016
【００６５】
ただし、Σｐは（ｐ，ｑ）∈Ｎを満たす全てのｐに対する和を表す。ブロック図ではノードｑ選択部２１４でノードｑ、ｘｑ値入力部２１１で処理するラベルを選択し、ｘｑビリーフ値算出部２１２で求め出力部２１３に送っている。なお、ｘｑビリーフ値算出部２１２で求めたビリーフ値を出力部２１３に送って評価すれば、全てのビリーフ値を保持する必要はない。
【００６６】
（最適化ステップ５）各ｑに対してｂｑ（ｘｑ）を最小とするラベルｘｑ＝｛０，１，・・・，Ｌｑ−１｝を選ぶ（Ｓ１０９）。ブロック図では出力部２１３がステップ５を行う手段である。
【００６７】
次に、本実施形態の画像処理装置の全体の構成を図６に示す。本実施形態の画像処理装置は、画像入力部６０１、ＭＰＵ（演算ユニット）６０２、データ用メモリ６０３、プログラム用メモリ６０４、ラベル出力部６０５を含んでいる。なお、図２に示した装置部分は、ＭＰＵ６０２、データ用メモリ６０３、プログラム用メモリ６０４に対応する。
【００６８】
本実施形態を実施すべきプログラムはプログラム用メモリ６０４に保存される。プログラムはＲＯＭに格納しても良いし、他の装置（例えばハードディスク）を接続してプログラム実行時にＲＡＭに格納しても良い。プログラムはＭＰＵ６０２により実行される。プログラムは、例えば次の指示を出すように作成する。
【００６９】
（１）画像入力部６０１からの入力した画像をデータ用メモリ６０３に格納する。  
（２）データ用メモリ６０３にメッセージ値や更新中の一時メッセージ値、ビリーフ値を記録するためのバッファを確保する。  
（３）データ用メモリ６０３上のメッセージ値を初期化する。  
（４）画像とメッセージ値を参照しながら、拡張高速ＢＰ更新ステップ１〜４にしたがって一時メッセージ値を算出する。  
（５）求めた一時メッセージ値を利用してメッセージ値を更新する。  
（６）メッセージ値の更新を所定の回数だけ繰り返す。  
（７）メッセージ値を参照しながら、ビリーフ値を算出する。  
（８）最小のビリーフ値を与えるラベル値を算出する。  
（９）ラベル出力部６０５に結果を出力する。
【００７０】
ここで、（１）と（９）以外の処理は全てＭＰＵ６０２、データ用メモリ６０３、プログラム用メモリ６０４のみで行われる。なお、データ用メモリ６０３とプログラム用メモリ６０４は分離されていなくても良い。例えば、パーソナルコンピュータのようにプログラムとデータのいずれにも利用できるメモリを備えていれば、それを必要なメモリサイズで区切り、プログラムとデータの両方に利用できる。
【００７１】
≪ケース３：コスト差分値が小数の場合≫  
ｒｐｑが整数でない場合についてメッセージ値の更新手法は拡張高速ＢＰ更新ステップ２を次のように変更する。  
（拡張高速ＢＰ更新ステップ２（ｒｐｑが整数でない場合））各ラベルｘｑに対して、まず次の集合を求める。
【数１７】
2008225778.tif 000017
【００７２】
ここで、  
【数１８】
2008225778.tif 000018
【００７３】
はそれぞれｘの小数点以下を切り捨てた値、ｘの小数点以下を切り上げた値を表し、コンピュータ言語でfloor(x)、ceil(x)として知られる関数に対応する。求めた｛ｙｐ（ｘｐ）｝内の各要素をｙｐ（ｘｐ）として、以下を実行する。
【００７４】
ｍ（ｙｐ（ｘｐ））←ｍｉｎ（ｍ（ｙｐ（ｘｐ）），ｈ（ｘｐ）＋ｓ‖ｘｐ−ｙｐ（ｘｐ）＋ｒｐｑ‖）
なお、ｒｐｑが整数の場合は、切り捨て値と切り上げ値が一致するために、上述した手法と同等になる。
【００７５】
（拡張高速ＢＰ更新ステップ１で無限大をセットしない手法）  
先の説明では拡張高速ＢＰ更新ステップ１で全てのラベルに対してｍ（ｘｑ）←＋∞と初期化し、拡張高速ＢＰ更新ステップ３と４ではそれら全てを対象とした。しかし、
【数１９】
2008225778.tif 000019
【００７６】
は、ｘｐ＝｛０，１，・・・，Ｌｐ−２｝に対しては明らかにｙｐ（０）とｙｐ（Ｌｐ−１）の間の値をとる。そこで、ｙｐ（０）よりも小さいラベル値およびｙｐ（Ｌｐ−１）よりも大きいラベル値に対するメッセージ値の計算は拡張高速ＢＰ更新ステップ１、３、４で行わなくても、拡張高速ＢＰ更新ステップ４の後に他のメッセージとの差分を求めれば、結果は同じになる。具体的な手順は次の通りである。  
ｙｐ（Ｌｐ−１）よりも大きいラベル値：ｍ（ｙｐ（Ｌｐ−１）＋１）から順に、ラベル値を１ずつ増やしながらｍ（ｘｑ）←ｍ（ｘｑ−１）＋ｓを実行する。  
ｙｐ（０）よりも小さいラベル値：ｍ（ｙｐ（０）−１）から順に、ラベル値を１ずつ減らしながらｍ（ｘｑ）←ｍ（ｘｑ＋１）＋ｓを実行する。  
（ステレオマッチング）  
エネルギー関数最小化のためのラベル選択問題は様々な産業的応用を持つが、ここではその１つであるステレオマッチングを例として説明する。
【００７７】
入力画像として、２つのカメラの視差が画像の水平方向のずれとしてあらわれる画像（左右のカメラから各１枚）が与えられるものとする。このような画像は、例えば、カメラを正確に並べて配置するか、２つのカメラの配置情報を利用し撮影画像に変換を施すことで得られることが知られている。ステレオマッチングでは、左右の画像の各画素における視差をラベルとし、Ｖに各画素を単独で見た場合の各ラベルの信頼度（値が大きくなるほどｘｉが出現しにくいことを表す）を記述し、Ｎを画素とその近傍の画素の組を全て集めた集合とし、Ｗに近傍画素間の関係（値が大きくなるほど、（ｘｉ，ｘｊ）の組が出現しにくいことを表す）を記述したうえで、Ｅを最小化するように各画素のラベルを選択する。得られた視差とカメラの情報があれば、例えば、それらを用いて各画素の奥行き情報を算出できることが知られている。
【００７８】
本実施形態の画像処理装置を利用することにより、より精度の高いステレオマッチング、および柔軟なステレオマッチングが可能になる。なお、視点の異なる画像が３枚以上ある場合には、例えば１枚の画像を基準として、それ以外の各画像に対してステレオマッチングで対応付けを行って基準画像の各画素の奥行き情報を算出し、基準画像の各画素の奥行き情報を各画像で算出した奥行き情報の平均値とすることで高精度な奥行き情報を算出できる。
【００７９】
（階層的探索）  
ステレオマッチングの２ノード項の例として、ここでは次の式を考える。  
Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ‖ｘｐ−ｘｑ‖，ｄｐｑ）  
この式は、線形コストとＰｏｔｔｓの組み合わせによって解ける形をしている。例えば、０〜３０画素の範囲内で１画素単位の視差を求めるのであれば、必要なラベル数は３１として高速ＢＰの手法を適用すれば良いことになる。ところで、奥行き情報を正確に算出するには視差の精度が高いほど良い。しかし、必要なラベル数は視差の精度に比例して増大するため、例えば０．１画素単位の視差を求めるのであればラベル数は３００以上になる。したがって、高い視差の精度を得たい場合、高速ＢＰの手法を用いた場合は多くのメモリが必要となり、また計算時間も増大する。
【００８０】
一方、本実施形態の手法を用いることで、視差の階層的な算出を実現し、必要なメモリ、計算時間を抑制できる。なお、Ｖ（ｘｐ）はｘｐ以外の変数を含まなければどのような形でもかまわないが、例えば次の形が使える。  
Ｖ（ｘｐ）＝ｍｉｎ（‖Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）−Ｉｒｉｇｈｔ（Ｘｐ−ｘｐ，Ｙｐ）‖，ｃｐ）  
Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）、Ｉｒｉｇｈｔ（Ｘｐ，Ｙｐ）はそれぞれ座標（Ｘｐ，Ｙｐ）における左画像・右画像の輝度値を表す。また、ｃｐは定数である。
【００８１】
以下、画素をノードとして、画素をｐ、視差の精度をδ、視差の最小値をγｐ、最大値をγｐ＋（Ｌｐ−１）δで表す。なお、ここで用いた記号γ（ガンマ）はｒ（アール）とは別の記号である。フローチャートを図７に示す。ステレオマッチングは、次の流れで実現できる。
【００８２】
（ステップ１）（初期化）Ｋ＝１とし、各画素ｐに対して、γｐが視差の最小値、δが初期探索での視差の精度、γｐ＋（Ｌｐ−１）δが視差の最大値となるように各パラメータを初期化する（Ｓ７０１）。例えば、視差を０〜３０の範囲、１画素精度とするのであれば、γｐ＝０、δ＝１、Ｌｐ＝３１とする。ノードごとに候補ラベルを決定する（ステップＳ７０２）。
【００８３】
（ステップ２）２ノード項を次の形として求め（ステップＳ７０３）、本実施形態のメッセージ更新を用いたＢＰ手法（ステップＳ１０１〜ステップＳ１０９）で各画素にラベルを割り当てる（Ｓ７０４）。
【００８４】
Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ‖（ｘｐ−ｒｐ）−（ｘｑ−ｒｑ）‖，ｄｐｑ）
（ステップ３）Ｋ←Ｋ＋１とする。Ｋがしきい値に達していれば終了し、達していなければステップ４に進む。
【００８５】
（ステップ４）δを小さくする（Ｓ７０５）。例えば、δに予め定めた１未満の値（例えば０．５）を掛ける。各ｐの視差の候補は｛γｐ，・・・，γｐ＋（Ｌｐ−１）δ｝であるから、これらの候補が予め決めた基準に合致するようにγｐとラベル数Ｌｐを決める。先の例であれば、ステップ２で求めた視差を中心とした−３δ〜＋３δの範囲となるようにγｐとラベル数Ｌｐを決める。これらの値はｐごとに変えても良い。｛δ，γｐ，Ｌｐ｝の更新後、ステップ２に戻り、先の手順を繰り返す（Ｓ７０６）。
【００８６】
この手法を用いると、ラベル数Ｌｐの増加を抑えながら精度の高い（小さなδでの）視差の推定ができる。なお、スケールを表す数値Ｋがしきい値になるまで繰り返す代わりに、精度δが一定値になるまで繰り返しても良い。その場合、Ｋを計算する必要はない。また、ラベル数Ｌｐは全てのｐに対して同じである必要はない。例えば、画素ｐの視差とその４近傍の画素での視差での最大値と最小値を算出し、γｐ＋（Ｌｐ−１）δが最大値以上、γｐが最小値以下となるように個々のγｐとＬｐを決めれば、物体境界のように視差の変動が激しい位置での推定精度を上げることができる。
【００８７】
なお、先のステップ１では０〜３０の範囲で初期化したが、例えば視差の推定をステレオ動画像に対して行う場合には、例えば、直前の時刻における視差を中心とした一定の範囲としても良い。この場合、２ノード項にあらわれるγｐ−ｒｑが小数とならないように各γｐを量子化しても良いし、先に示した小数を扱う手法を用いても良い。
【００８８】
なお、本発明はより複雑なステレオマッチングでも利用できる。例えば文献「J. Sun et al., “Symmetric stereo matching for occlusion handling,” in Proc. IEEE Conference on Computer Vision and Pattern Recognition, vol. 2, pp.399-406, 2005.」では、ステレオの視差（disparity）と画素の隠れ情報（occlusion）の両方を最適化するために、第１の最適化であるステレオの視差の最適化と第２の最適化である画素の隠れ情報の最適化を、交互に繰り返し（例えば予め定めた回数、あるいはラベルの変化量の合計が別途定めたしきい値以下になるまで）行う。このうちステレオの視差の最適化は、視差の精度を上げると計算量やメモリ利用量が大幅に増える問題があるが、本発明を利用すれば、計算量やメモリ利用量をあまり増やすことなく視差の精度を上げられる。なお、例えば文献「J. Sun et al., “Symmetric stereo matching for occlusion handling,” in Proc. IEEE Conference on Computer Vision and Pattern Recognition, vol. 2, pp.399-406, 2005.」にある「セグメンテーションによるエネルギー関数への制約条件（Segmentation as Soft Constraint）」のように、最適化すべきエネルギーのＷ（Ｘｐ，Ｘｑ）項に影響を与えない改良（この改良に限定はされない）は、本発明に何ら変更を加えることなくそのまま利用できる。なお隠れ情報は隠れているかいないかの２値をとるため、従来の手法で最適化すれば良い。
【００８９】
ここまでは、１次元のノードについての説明であるが、以下２次元以上の場合について説明する。  
（オプティカルフロー、多次元動きベクトル）  
画像処理の分野において、２枚(あるいはそれ以上)の画像中の画素を対応付ける問題は、ステレオマッチング以外にもあらわれる。特に、２次元の画像の対応づけ問題は、動画像に対する動きベクトルの推定問題として様々な応用を持つ（例えば画像圧縮、人物や車両の追跡、超解像）。また、例えばＸ線ＣＴやＭＲＩのような医療用画像は３次元画像であることが多いが、これらの対応付け問題は３次元画像中の画素の対応付け問題になる。
【００９０】
以下、次元を表す番号（１）、（２）を記号の右下に付与する。２次元のノードｑに対するラベルをｘｑ＝（ｘ（１）ｑ，ｘ（２）ｑ）で表す。ラベルの範囲をｘ（１）ｑ＝｛０，１，・・・，Ｌ（１）ｑ−１｝、ｘ（２）ｑ＝｛０，１，・・・，Ｌ（２）ｑ−１｝とする。ノードｑの候補となるラベル数はＬｑ＝Ｌ（１）ｑＬ（２）ｑである。２ノード項として、次の形を考える。  
Ｗ（ｘｐ，ｘｑ）＝ｓ（１）‖ｘ（１）ｐ−ｘ（１）ｑ＋ｒ（１）ｐｑ‖＋ｓ（２）‖ｘ（２）ｐ−ｘ（２）ｑ＋ｒ（２）ｐｑ‖
Ｖ（ｘｐ）はｘｐ以外の変数を含まなければどのような形でもかまわないが、例えば次の形が使える。  
Ｖ（ｘｐ）＝ｍｉｎ（‖Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）−Ｉｒｉｇｈｔ（Ｘｐ−ｘ（１）ｐ，Ｙｐ−ｘ（２）ｐ）‖，ｃｐ）
Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）、Ｉｒｉｇｈｔ（Ｘｐ，Ｙｐ）はそれぞれ座標（Ｘｐ，Ｙｐ）における左画像、右画像の輝度値を表す。また、ｃｐは定数である。
【００９１】
この２ノード項に対する拡張高速ＢＰを実現するには、次の手法でメッセージ更新手法を行えばよい。
【００９２】
（拡張高速ＢＰ更新ステップ１）ラベルｘｑ＝｛（０，０），・・・，（Ｌ（１）ｑ−１，Ｌ（２）ｑ−１）｝に対し、ｍ（ｘｑ）←＋∞とする。
【００９３】
（拡張高速ＢＰ更新ステップ２）各ラベルｘｑに対して、次の式を実行する。
【００９４】
ｍ（ｙ（１）ｐ（ｘｐ），ｙ（２）ｐ（ｘｐ））←ｍｉｎ（ｍ（ｙ（１）ｐ（ｘｐ），ｙ（２）ｐ（ｘｐ）），ｈ（ｘｐ）＋ｓ（１）‖ｘ（１）ｐ−ｙ（１）ｐ（ｘｐ）＋ｒ（１）ｐｑ‖＋ｓ（２）‖ｘ（２）ｐ−ｙ（２）ｐ（ｘｐ）＋ｒ（２）ｐｑ‖）
ただし、ｙの値は
【数２０】
2008225778.tif 000020
【００９５】
により算出する。ｙ（２）ｐ（ｘｐ）についても記号の右下に付与されている（１）を（２）に変えただけで全く同じ手法で算出する。なお、コスト差分値が小数の場合も１次元の場合と同じように、切り捨てと切り上げの２つの値それぞれについてｍの更新式を実行すれば良い。
【００９６】
（拡張高速ＢＰ更新ステップ３）次の擬似コードで示した処理を実行する。以下、For A in {…} から Endfor までは{…}の要素をAに順に代入してForとEndforで囲まれた部分を実行することを表す。  
For  ｘ（１）ｑ  in  ｛０，・・・，Ｌ（１）ｑ−１｝  
For  ｘ（２）ｑ  in  ｛０，・・・，Ｌ（２）ｑ−１｝  
ｍ（（ｘ（１）ｑ，ｘ（２）ｑ））←ｍｉｎ（ｍ（（ｘ（１）ｑ，ｘ（２）ｑ）），ｍ（（ｘ（１）ｑ−１，ｘ（２）ｑ））＋ｓ（１））  
ｍ（（ｘ（１）ｑ，ｘ（２）ｑ））←ｍｉｎ（ｍ（（ｘ（１）ｑ，ｘ（２）ｑ）），ｍ（（ｘ（１）ｑ，ｘ（２）ｑ−１））＋ｓ（２））  
ただし、ｍ（（ｘ（１）ｑ，−１）＝＋∞、ｍ（（−１，ｘ（２）ｑ）＝＋∞とする。  
Endfor  
Endfor
（拡張高速ＢＰ更新ステップ４）次の擬似コードで示した処理を実行する。  
For  ｘ（１）ｑ  in  ｛Ｌ（１）ｑ−１，・・・，０｝  
For  ｘ（２）ｑ  in  ｛Ｌ（２）ｑ−１，・・・，０｝  
ｍ（（ｘ（１）ｑ，ｘ（２）ｑ））←ｍｉｎ（ｍ（（ｘ（１）ｑ，ｘ（２）ｑ）），ｍ（（ｘ（１）ｑ＋１，ｘ（２）ｑ））＋ｓ（１））  
ｍ（（ｘ（１）ｑ，ｘ（２）ｑ））←ｍｉｎ（ｍ（（ｘ（１）ｑ，ｘ（２）ｑ）），ｍ（（ｘ（１）ｑ，ｘ（２）ｑ＋１））＋ｓ（２））  
ただし、ｍ（（ｘ（１）ｑ，Ｌ（２）ｑ）＝＋∞、ｍ（（Ｌ（１）ｑ，ｘ（２）ｑ）＝＋∞とする。  
Endfor  
Endfor
もし２ノード項が、
Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ（１）‖ｘ（１）ｐ−ｒ（１）ｑ＋ｒ（１）ｐｑ‖＋ｓ（２）‖ｘ（２）ｐ−ｒ（２）ｑ＋ｒ（２）ｐｑ‖，ｄｐｑ）
のようにPottsと組み合わせた形をしているのであれば、拡張高速ＢＰ更新ステップ４の後で上記＜ＷがＰｏｔｔｓの場合＞で述べたPottsに対するメッセージ更新を行えば良い。
【００９７】
３次元以上の場合であっても、各画素で全ての次元に対する直前あるいは直後の画素との差分を評価すれば、同様にメッセージの更新が可能である。  
（多次元動きベクトルの利用シーン）  
先に１次元の対応付け問題であるステレオマッチングを例として、精度を上げるとラベル数が増加する問題について説明した。２次元の対応付け問題であるオプティカルフロー推定や３次元以上の対応付け問題であっても同じ問題は起こる。しかも、候補ラベル数は次元数に対し指数的に増えるため、２次元以上の対応付け問題では、ステレオマッチングと比べても、精度を上げることによる計算量、メモリ消費量の増大量はさらに大きい。この対策もステレオマッチングと同様で、本実施形態の手法を用いてフローの精度を徐々に上げながらエネルギー最小化を繰り返せば良い。
【００９８】
以下、オプティカルフロー推定問題でフローの精度を上げる手法について説明する。  
以下、画素をノードとして、画素をｐ、フローの精度を（δ（１），δ（２））、２次元の各軸におけるフローの最小値をそれぞれγ（１）ｐ，γ（２）ｐ、フローの最大値をそれぞれγ（１）ｐ＋（Ｌ（１）ｐ−１）δ（１）、γ（２）ｐ＋（Ｌ（２）ｐ−１）δ（２）で表す。オプティカルフロー推定の流れを次に示す。
【００９９】
（ステップ１）（初期化）ｋ＝１とし、各画素ｐに対して、フローの精度（δ（１），δ（２））、各軸におけるフローの最小値γ（１）ｐ，γ（２）ｐおよび最大値γ（１）ｐ＋（Ｌ（１）ｐ−１）δ（１）、γ（２）ｐ＋（Ｌ（２）ｐ−１）δ（２）を初期化する。
【０１００】
（ステップ２）２ノード項を次の形として、本実施形態の拡張高速ＢＰメッセージ更新を用いて各画素にラベルを割り当てる。  
Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ（１）‖ｘ（１）ｐ−ｘ（１）ｑ＋ｒ（１）ｐｑ‖＋ｓ（２）‖ｘ（２）ｐ−ｘ（２）ｑ＋ｒ（２）ｐｑ‖，ｄｐｑ）
（ステップ３）ｋ←ｋ＋１とする。ｋがしきい値に達していれば終了し、達していなければステップ４に進む。
【０１０１】
（ステップ４）フローの精度（δ（１），δ（２））を小さくする。例えば、（δ（１），δ（２））に予め定めた１未満の値（例えば０．５）を掛ける。各ｐにおけるフローの最小値γ（１）ｐ，γ（２）ｐおよび最大値γ（１）ｐ＋（Ｌ（１）ｐ−１）δ（１）、γ（２）ｐ＋（Ｌ（２）ｐ−１）δ（２）が予め決めた基準に合致するように（例えば、ステップ２で求めたフローを中心とし、（δ（１），δ（２））の範囲となるように）γｐと各軸でのラベル数Ｌ（１）ｐ、Ｌ（２）ｐを決める。これらの値はｐごとに変えても良い。パラメータの更新後はステップ２に戻る。
【０１０２】
以上の流れでフロー推定を行えば、計算時間やメモリ消費量の増大を抑えながら高い精度でフローを推定できる。
【０１０３】
なお、ステレオマッチングやフロー推定で示したエネルギー関数の形は一例であり、特にこの形に限定されるものではない。ここでは一例として、エネルギー最小化の前にセグメンテーション（例えば、“D. Comaniciu et al., Mean Shift: A Robust Approach Toward Feature Space Analysis, IEEE Transaction on Pattern Analysis and Machine Intelligence, vol.24, no.5, May 2002”に記載の手法）により画像を位置が近く色の似ている領域に分け、その領域単位で動きを求め、それをエネルギー最小化により修正するフロー推定について述べる（ステレオマッチングでも同様のことが行える）。
【０１０４】
例えば、まず領域の動きをアファイン変換（位置の線形変換と平行移動を組み合わせた６パラメータの変換）で近似して領域単位で動きを推定する（ステレオマッチングの場合は領域単位で視差を推定する）。この推定は、細かい部分の影響を受けにくいため領域全体でみると推定結果は高い信頼性を持つが、個々の画素で見ると必ずしも信頼できないという性質を持つ。そこで、各領域の動きから画素（＝ノード）ｐの動きを求めて得られた動きを（ｘｓｅｇ（１），ｘｓｅｇ（２））として、エネルギー関数のＶを、
Ｖ（ｘｐ）＝ｍｉｎ（‖Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）−Ｉｒｉｇｈｔ（Ｘｐ−ｘ（１）ｐ，Ｙｐ−ｘ（２）ｐ）‖，ｃ１ｐ）＋ｍｉｎ（‖（ｘ（１）ｐ，ｘ（２）ｐ）−（ｘｓｅｇ（１），ｘｓｅｇ（２））‖，ｃ２ｐ）
の形とし（Ｗは前述のものが使える）、エネルギー関数の最小化を行えば良い。なお、他の記号は先の説明と同様で、Ｉｌｅｆｔ（Ｘｐ，Ｙｐ）、Ｉｒｉｇｈｔ（Ｘｐ，Ｙｐ）はそれぞれ座標（Ｘｐ，Ｙｐ）における左画像、右画像の輝度値を表し、ｃ１ｐ、ｃ２ｐは定数である（ステレオマッチングの場合はｘｐ（２）、ｘｓｅｇ（２）に係る項を除去し、１次元になおせば良い）。他にも、例えばＶとして座標（Ｘｐ，Ｙｐ）の輝度誤差を用いる代わりに、座標（Ｘｐ，Ｙｐ）の近傍の輝度誤差の合計を用いることもできる。
【０１０５】
また、文献「J. Sun et al., “Symmetric stereo matching for occlusion handling,” in Proc. IEEE Conference on Computer Vision and Pattern Recognition, vol. 2, pp.399-406, 2005.」でステレオマッチングにおいて視差とオクルージョンを交互に最適化した方法と同じように、フロー推定についてもフローとオクルージョンを交互に最適化することで、オクルージョンを考慮したフロー推定が可能になる。この場合でも、フローの最適化において本発明を利用すればあまりコストを増やさずに高い精度を達成できる。
【０１０６】
Ｖは画素ごとに見た場合の対応位置の選択基準を与えており、これらの例のように別のモデルを与えることで画素ごとの推定精度を上げていくことが可能である。Ｖを改良した場合でも、本実施形態の手法はそのまま適用できる。
【０１０７】
（斥力）  
ここまでの例では２ノード項として、ラベルの差がｒｐｑとなれば最小となる関数を考えてきた。この関数は、２ノード間の関係が「もしノードｐのラベルがｘｐであれば、ノードｑのラベルはｘｐ＋ｒｐｑに近いことが望ましい」ことを定式化したものである。ところで、応用によっては、２ノード間の関係が「もしノードｐのラベルがｘｐであれば、ノードｑのラベルはｘｐ＋ｒｐｑから遠いことが望ましい」という関係の定式化が必要である。その場合にあらわれる２ノード項として、例えば以下の関数がある。ｄｐｑは定数で、ノードｐ、ｑごとに別の値でも良い。
【０１０８】
Ｗ（ｘｐ，ｘｑ）＝ｄｐｑ−ｓ‖ｘｐ−ｘｑ＋ｒｐｑ‖
この２ノード項を用いたエネルギー関数のメッセージ更新は、次の手順で高速に行える。
【０１０９】
（更新ステップ１）両端のメッセージ値をｍ（０）←＋∞、ｍ（Ｌｑ−１）←＋∞で初期化する。
【０１１０】
（更新ステップ２）各ラベルｘｐに対して、まず次の位置を求める。
【数２１】
2008225778.tif 000021
【０１１１】
これを利用して、各ラベルｘｐに対して両端のメッセージ値を次の式で更新する。  
ｍ（０）←ｍｉｎ（ｍ（０），ｈ（ｙｐ（ｘｐ））＋ｄｐｑ−ｓ‖−ｒｐｑ＋ｙｐ（ｘｐ）‖）  
ｍ（Ｌｑ−１）←ｍｉｎ（ｍ（Ｌｑ−１），ｈ（ｙｐ（ｘｐ））＋ｄｐｑ−ｓ‖Ｌｑ−１−ｒｐｑ＋ｙｐ（ｘｐ）‖）
（更新ステップ３）ｘｑ＝｛１，・・・，Ｌｑ−１｝の範囲で、ラベル番号ｘｑを１から開始して１ずつ増やしながら以下を実行する。  
ｍ（ｘｑ）←ｍ（ｘｑ−１）＋ｓ
（更新ステップ４）ｘｑ＝｛Ｌｑ−２，・・・，０｝の範囲で、ラベル番号ｘｑをＬｑ−２から開始して１ずつ減らしながら以下を実行する。
【０１１２】
ｍ（ｘｑ）←ｍｉｎ（ｍ（ｘｑ），ｍ（ｘｑ＋１）＋ｓ）
（イメージ復元）  
本実施形態の手法は、ノイズが加わった画像のノイズを除去するためにも使える。ノードを画素、ラベルを輝度値とすれば、ノイズ除去問題のエネルギー関数の各項は例えば、
Ｖ（ｘｐ）＝‖ｘｐ−Ｉｐ‖、Ｗ（ｘｐ，ｘｑ）＝ｍｉｎ（ｓ‖ｘｐ−ｘｑ‖，ｄｐｑ）
で表せる。ここでＩｐは入力画像の輝度を表し、ラベルｘｐはノイズ除去後の輝度を表す。輝度は８ビットで表現した場合でも２５６種類あり、例えば１２ビットで表現すれば４０９６種類にもなるため、本実施形態の手法を利用した効率的な最小化は効果的である。このエネルギー関数を先に説明したステレオマッチングと同じ手法で最小化でき、選択されたラベルは各画素におけるノイズ除去後の輝度になる。
【０１１３】
（確率積の利用）  
主に画像処理であらわれるマルコフ・ランダム・フィールド（ＭＲＦ）やコンディショナル・ランダム・フィールド（ＣＲＦ）と呼ばれる確率モデルでは、次の関数の最大化と等価な処理を行うことがよくある。  
【数２２】
2008225778.tif 000022
【０１１４】
この関数の最大化問題は、負の対数−ｌｏｇＰ（ｘ１，ｘ２，・・・，ｘｎ）の最小化と等価である。−ｌｏｇＰ（ｘ１，ｘ２，・・・，ｘｎ）はＢＰで扱うエネルギー関数の形であるから、ＭＲＦやＣＲＦをモデルとした確率の最適ラベル選択問題は本実施形態の手法により解ける。
【０１１５】
（ＴＲＷ／Ｗｅｉｇｈｔｅｄ ＢＰの利用）  
ＢＰに似た最適化アルゴリズムとして、２ノード項Ｗ（ｘｐ，ｘｑ）そのものではなく重み付けした値を使ってメッセージ更新を行うアルゴリズムや、２ノード項Ｗ（ｘｐ，ｘｑ）およびメッセージに重み付けを行い、さらに逆方向のメッセージ値も用いてメッセージ更新を行うアルゴリズム（ツリー・リウエイテッド・マックスプロダクト法、以下ＴＲＷ）が知られているが、本実施形態の手法はこれらのアルゴリズムに対してもそのまま適用できる。以下、ＴＲＷを例として本実施形態の手法の適用手法を説明する。
【０１１６】
最小化すべきエネルギー関数として、ＢＰと同じ次の関数を考える。  
【数２３】
2008225778.tif 000023
【０１１７】
ＴＲＷではＢＰとは違い、次の更新式を用いてメッセージ更新を行う。  
【数２４】
2008225778.tif 000024
【０１１８】
ただし、Σsは（ｓ，ｐ）∈Ｎかつｓ≠ｑを満たす全てのｓに対する和を表す。ρｐｑは２ノードｐ、ｑ間の重みをあらわす定数で、もしρｐｑ＝１ならＢＰと等しくなる。この更新式は、  
【数２５】
2008225778.tif 000025
【０１１９】
として、２ノード項としてＷ’（ｘｐ，ｘｑ）＝（１／ρｐｑ）Ｗ（ｘｐ，ｘｑ）を考えれば、ＢＰで用いる更新式と全く同じ形になる。したがって、各２ノード項に対しＷ’（ｘｐ，ｘｑ）＝（１／ρｐｑ）Ｗ（ｘｐ，ｘｑ）を算出するステップを最初に追加すれば、本実施形態の手法をそのまま適用できる。
【０１２０】
（ｓｅｑｕｅｎｔｉａｌ ＢＰ／ＴＲＷの利用）  
メッセージ更新では、ｔ＋１番目のメッセージ値の算出にｔ番目のメッセージ値を用いている。具体的には、  
【数２６】
2008225778.tif 000026
【０１２１】
の算出にｔ番目のメッセージ値を用いている。したがって、ｔ番目、ｔ＋１番目のメッセージ値を保持する必要があり、必要なメモリは全メッセージ値の保持に必要なメモリの２倍になる。更新後のメッセージ値ｍ（ｔ＋１）ｐ→ｑ（ｘｑ）をｍ（ｔ）ｐ→ｑ（ｘｑ）の代わりに用いれば、メッセージ値を上書きするだけでよいため全メッセージ値１つ分のメモリで動作でき、必要なメモリを半分に減らせる。この手法はシーケンシャルＢＰとして知られており、本実施形態の手法でもそのまま利用できる。
【０１２２】
（Ｅｆｆｉｃｉｅｎｔ Ｍｕｌｔｉｓｃａｌｅ ＢＰとの併用）  
高速ＢＰで示されているように、２×２画素を１つのノードとみなしてＢＰを適用し、その結果を各画素のメッセージ値の初期値として使うことで、画像に対するＢＰの繰り返し回数を減らせる。この手法は本実施形態の手法と違いラベル数を減らせないためにメモリ消費量は減らせないが計算時間は減らせる。この手法はマルチスケールＢＰと呼ばれ、本実施形態の手法とも併用できる。本実施形態の手法を用いたラベル数の削減とマルチスケールＢＰを併用すれば、マルチスケールＢＰによる計算時間の削減と本実施形態の手法による計算時間、メモリ消費量の削減の両方の効果が得られる。
【０１２３】
なお、本実施形態で正の無限大を用いた例がいくつかあるが、これらは無限大でない適当な値を使用しても良い（なるべく大きい値であることが望ましい）。
【０１２４】
以上の実施形態の画像処理装置により、階層的探索と拡張高速ＢＰを併用することができるので、高速かつ高精度に画素のラベル値を選択することが可能になる。例えばステレオマッチングやオプティカルフロー推定の階層化のように、従来の高速ＢＰでは扱えなかった２ノード項を持つエネルギー関数に対し、エネルギー関数を高速に最小化できるようになる。典型的な例では、計算時間は数分程度、消費メモリも〜１Ｇ程度に改善される。
【０１２５】
また、本実施形態によれば、２ノード項ＷがＰｏｔｔｓであらわせず、ノードごとにラベルの範囲や種類数が異なる場合や、ラベルが小数を含む場合であってもＢＰを高速に実行できるため、エネルギー最小化を解くことを必要とする場合でも大幅な高速化が実現できる。
【０１２６】
なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。
【図面の簡単な説明】
【０１２７】
【図１】実施形態の画像処理装置の動作の一例を示す図。
【図２】実施形態の画像処理装置の更新処理装置部分のブロック図。
【図３】メッセージ値更新の具体的な一例を示す図。
【図４】メッセージ値更新の具体的な一例での拡張高速ＢＰ更新ステップ２を示す図。
【図５】メッセージ値更新の具体的な一例での拡張高速ＢＰ更新ステップ３を示す図。
【図６】実施形態の画像処理装置のブロック図。
【図７】ステレオマッチングを行う場合の動作の一例を示すフローチャート。
【符号の説明】
【０１２８】
２０１・・・メッセージ値保持部、２０２・・・メッセージ初期値入力部、２０３・・・ｈ（ｘｐ）値算出部、２０４・・・ｘｑ一時メッセージ値保持部、２０５・・・ＰｅａｋＸｑ算出部、２０６・・・集約値算出部、２０７・・・最適値選択部、２０８・・・最適値出力部、２０９・・・ノード対（ｐ，ｑ）選択部、２１０・・・ｘｐ値入力部、２１１・・・ｘｑ値入力部、２１２・・・ｘｑビリーフ値算出部、２１３・・・出力部、２１４・・・ノードｑ選択部、６０１・・・画像入力部、６０３・・・データ用メモリ、６０４・・・プログラム用メモリ、６０５・・・ラベル出力部。
【図１】
2008225778.tif 000029
【図２】
2008225778.tif 000030
【図３】
2008225778.tif 000031
【図４】
2008225778.tif 000032
【図５】
2008225778.tif 000033
【図６】
2008225778.tif 000034
【図７】
2008225778.tif 000035