JP2008225778A

JP2008225778A - 画像処理装置

Info

Publication number: JP2008225778A
Application number: JP2007062360A
Authority: JP
Inventors: Hidenori Takeshima; 秀則竹島; Toshimitsu Kaneko; 敏充金子
Original assignee: Toshiba Corp
Current assignee: Toshiba Corp
Priority date: 2007-03-12
Filing date: 2007-03-12
Publication date: 2008-09-25

Abstract

【課題】高速かつ高精度に画素のラベル値を選択する。
【解決手段】Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）とＶ（ｘ_ｐ）を含むエネルギー関数を評価する各評価関数ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）のｔ＝０の初期値を設定する手段２０２と、ｈ（ｘ_ｐ）を設定する手段２０３と、Ｗが最小値になるｘ_ｑであるＰｅａｋＸｑを算出する手段２０５と、ｘ_ｐごとにｈとＷ（ｘ_ｐ−ＰｅａｋＸｑ＋ｒ）との集約値を算出する手段２０６と、前記集約値と前記ｍ^（ｔ）にｘ_ｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’^（ｔ）とする手段２０７と、ＷとＷ（ｘ_ｐ−（ｘ_ｑ＋ｋ）＋ｒ）との差分とｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ＋ｋ）との集約値と、ｍ’^（ｔ）との小さい値の方を最適値としてｍ^{（ｔ＋１）}とする手段２０７と、ｘ_ｑごとにビリーフ値ｂ_ｑ（ｘ_ｑ）を算出しｂ_ｑを最小にするラベル値ｘ_ｑを選択する手段２１３を具備する。
【選択図】図２

Description

本発明は、予め定義されたエネルギー関数の最小化を効率的に計算する技術に関し、例えば、動画像の動きベクトルの高速な算出に利用できる画像処理装置に関する。

超解像のように高い精度の動きベクトルを必要とするアプリケーションでは、従来のようにローカルの情報だけを用いて動きベクトルを算出するだけでは不十分である。一方、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.

しかし、依然として計算量や消費メモリが多く、例えば、計算時間は数十分で、消費メモリはＢＰと同じである。

一方、ブロックマッチングにおいては、粗い動きベクトルを求めその周辺のみを探索する階層的探索や、他のフレームの動きベクトルを基準としてその周辺のみを探索する手法などの高速化手法が用いられている。これらの高速化手法と高速ＢＰを組み合わせることができれば、高速で精度の高い動き推定が可能になる。

しかし、階層的探索などの高速化手法を利用すると、画素ごとに動きベクトルの範囲が異なるために高速ＢＰが適用できなくなる。

ところで、高速化ＢＰは、ここで示した動きベクトルの算出の他に、例えば、ステレオ画像の視差の算出、画像の復元に利用される手法として、多次元ベクトル空間上で予め定められた点の集合のうち、予め定義されたエネルギー関数を最小化する点を選択する手法である。これら動きベクトル、視差、奥行き情報は画素のラベル値で決定される。

この発明は、上述した事情を考慮してなされたものであり、高速かつ高精度に画素のラベル値を選択する画像処理装置を提供することを目的とする。

上述の課題を解決するため、本発明の画像処理装置は、画像の各画素に対応付けられた各ノードｐ、ｑのそれぞれのラベル値ｘ_ｐ、ｘ_ｑを変数とし、ノード対（ｐ、ｑ）およびｐとｑに関連付けられた定数ｒを用いて定義された関数Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）と関数Ｖ（ｘ_ｐ）を含む、最小化するｘ_ｐを求めるエネルギー関数を設定する設定手段と、各ノード対（ｐ、ｑ）に対応して各評価関数ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）を格納している格納手段と、前記各評価関数のｔ＝０での値である初期値を設定する設定手段と、

にしたがってｈ（ｘ_ｐ）を設定する設定手段と、ｘ_ｐごとにＷ（ｘ_ｐ−ｘ_ｑ＋ｒ）が最小値になるｘ_ｑであるＰｅａｋＸｑを算出する算出手段と、ｘ_ｐごとにｈ（ｘ_ｐ）とＷ（ｘ_ｐ−ＰｅａｋＸｑ＋ｒ）との和または積である集約値を算出する算出手段と、前記集約値と前記ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）にｘ_ｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ）とする更新手段と、ｘ_ｐごとに、Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）とＷ（ｘ_ｐ−（ｘ_ｑ＋ｋ）＋ｒ）との差分とｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ＋ｋ）との集約値と、ｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ）との小さい値の方を最適値としてｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）とする更新手段と、ｔが繰り返し回数Ｔよりも小さい場合には、ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）を計算させ、最終的にｍ^（Ｔ） _ｐ→ｑ（ｘ_ｑ）を取得する取得手段と、
ｘ_ｑごとに

を算出する算出手段と、ｑごとに、ｂ_ｑ（ｘ_ｑ）を最小にするラベル値ｘ_ｑを選択する選択手段と、を具備することを特徴とする。

本発明の画像処理装置によれば、高速かつ高精度に画素のラベル値を選択することが可能になる。

以下、図面を参照しながら本発明の実施形態に係る画像処理装置について詳細に説明する。
まず、ラベル値の例として動きベクトルを例に挙げ、画素ごとに動きベクトルの範囲が異なっている場合でも、高速ＢＰが適用できるようにメッセージ値の計算が修正できる一例をごく簡単に示す。

一例として、第１画素の動きベクトルｐと第２画素の動きベクトルｑの評価関数（後述するメッセージ値に対応）がｆ（ｑ）＋｜ｐ−ｑ｜という形のときに、各ｑに対して評価関数の最小値を求める計算を考える。なお、ここでは、簡単のため動きベクトルがスカラーであるとする。

高速ＢＰでは、各ｑに対して｜ｐ−ｑ｜＝０のときの評価関数（すなわち、ｆ（ｑ））を予め算出し、その後｜ｐ−ｑ｜が非０の場合を差分計算により逐次求めながら評価関数の最小値を求める。

しかし、高速ＢＰでは、ｑの範囲内では｜ｐ−ｑ｜＝０とならない場合が考慮されていない。例えば、ｐ＝−２，−１，０，１，２、ｐ＝−１，０，１，２，３のようにｐとｑの値の範囲が異なる場合、ｐ＝−２，−１，０，１，２、ｐ＝−３，−２，−１，０，１，２，３のように候補数が異なる場合、ｐ＝−２，−１，０，１，２、ｐ＝−２．５，−１．５，−０．５，０．５，１．５のようにｐとｑの差が小数になる場合は考慮されていない。

一方、本実施形態の画像処理装置による手法（拡張高速ＢＰと呼ぶ）では、ｑの範囲内では｜ｐ−ｑ｜＝０とならない場合（ｑの範囲外で最小値を持つ評価関数の場合）については、ｑの範囲境界における最小値を算出してから、高速ＢＰと同様に差分計算による逐次評価を行う。こうすることで階層的探索と拡張高速ＢＰを併用することができ、高速かつ高精度に動きベクトルを推定できるようになる。

（ＢＰ）
次に、ＢＰについて説明する。
多次元ベクトル空間上で予め定められた点の集合のうち、予め定義されたエネルギー関数を最小化する点を選択する技術は、多くの応用を持ち、例えば、動画像の動きベクトルの算出、ステレオ画像の視差の算出、画像の復元に利用されている（上記のP. F. Felzenszwalb and D. R. Huttenlocherの文献参照）。一般的なエネルギー関数の厳密な最小化は極めて困難であることが知られている。その近似解の算出手法の１つとして、ビリーフプロパゲーション（ＢＰ：Belief Propagation）と呼ばれる技術が知られている。ＢＰでは、高々２つの変数から成る項の総和として表現される次のエネルギー関数の近似的な最小化を考える（以下、近似的な最小化を単に最小化と呼ぶ）。

ラベルの添え字ｉをノード番号、あるいは単にノードと呼ぶ。ｘ_ｉ＝｛０，１，・・・，Ｌ_ｉ−１｝はラベル、Ｎはエネルギー関数が含む各２ノード項を構成する変数対の集合をあらわす。また、各２ノード項を構成する（つまり、Ｎに属する）ノード対をエッジと呼ぶ。ＢＰでは、エネルギー関数を最小化するために次の繰り返し演算を行う。

（ステップ１）ｔ＝０とする（ｔはメッセージ更新回数を表す）。各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）の初期値を与える（あるいは全てのｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）を０に初期化する）。

（ステップ２）各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）を次のメッセージ更新式により更新する。

ただし、Σ_sは（ｓ，ｐ）∈Ｎかつｓ≠ｑを満たす全てのｓに対する和を表す。なお、ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）とｍ^（ｔ） _ｑ→ｐ（ｘ_ｑ）は別のメッセージ値であり、区別する必要がある。

（ステップ３）ｔを１増やす。ｔが予め決めておいた繰り返し回数Ｔより小さいなら、ステップ２に戻る。

（ステップ４）各ｑに対して、次のビリーフと呼ばれる値を求める。

ただし、Σ_ｐは（ｐ，ｑ）∈Ｎを満たす全てのｐに対する和を表す。

（ステップ５）各ｑに対してｂ_ｑ（ｘ_ｑ）を最小とするラベルｘ_ｑ＝｛０，１，・・・，Ｌ_ｑ−１｝を選ぶ。

ＢＰでは上記のステップ２でメッセージ更新式の繰り返し演算を行うために多くの計算時間が必要となる。

（高速ＢＰ）
次に、高速ＢＰについて説明する。高速ＢＰは、メッセージ更新式の繰り返し演算を行うための計算時間を少なくする。

高速ＢＰでは、全てのノードｉに対してラベルを同一の範囲０〜（Ｌ−１）とし、かつ、２ノード項Ｗを以下に示すＰｏｔｔｓや線形コストのいずれかの形に制限し、上記のステップ２の繰り返し演算の高速化を実現する。
ＰｏｔｔｓＷ（ｘ_ｐ，ｘ_ｑ）＝０（if ｘ_ｐ＝ｘ_ｑ），ｄ（otherwise）
線形コストＷ（ｘ_ｐ，ｘ_ｑ）＝ｍｉｎ（ｓ‖ｘ_ｐ−ｘ_ｑ‖，ｄ）
だだし、‖Ａ‖はＡのノルムを示す。ここで、ｐに依存する項をまとめて

で表すと、メッセージ更新式を次の式で書き換えられる。

高速ＢＰではこの式をもとにメッセージ更新式の算出を高速化する。

＜ＷがＰｏｔｔｓの場合＞
２ノード項ＷがＰｏｔｔｓの場合、メッセージ更新式は次の式になる。

したがって、はじめに

を計算しておけば、残りの項はｐに依存しないからｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）はｑのラベル数のオーダで計算できる。

＜Ｗが線形コストの場合＞
２ノード項Ｗが線形コストの場合、はじめに２ノード項がＷ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ‖の場合について計算し、得られた結果とｄとのｍｉｎをＰｏｔｔｓと同じ手法で計算すれば良い。２ノード項がＷ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ‖の場合についての計算では、メッセージ更新式は次の通りである。

メッセージ更新式を高速に計算するために、メッセージ更新は次の３つのステップに従って行う。この計算は「in-place」更新により行う、つまりｍ（ｘ_ｑ）の値を各ステップで上書きしていくことで実現する。
（高速ＢＰ更新ステップ１）ラベルｘ_ｑ＝０〜（Ｌ−１）に対し、ｍ（ｘ_ｑ）←ｈ（ｘ_ｑ）を計算する。
（高速ＢＰ更新ステップ２）１〜（Ｌ−１）までラベル番号ｘ_ｑを１ずつ増やしながら以下を実行する。
ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ−１）＋ｓ）
（高速ＢＰ更新ステップ３）（Ｌ−２）〜０までラベル番号ｘ_ｑを１ずつ減らしながら以下を実行する。

ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ＋１）＋ｓ）
高速ＢＰ更新ステップ３が終了した時点で、各ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）の値がｍ（ｘ_ｑ）に得られる。

以上のように、全てのノードに対してラベルを同一の範囲の整数０〜（Ｌ−１）とし、かつ２ノード項ＷがＰｏｔｔｓや線形コストのいずれかであれば、高速化によりＢＰの計算時間を大幅に短縮できる。

（高速ＢＰの問題点）
高速ＢＰでは、全てのノードに対するラベル数が同一の範囲の整数であることを前提としている。しかし、全てのノードに対するラベル数が同一の範囲の整数でなければならないと仮定してしまうと、例えば、次のような場合には適用できない。

＜ケース１＞ノードｐのラベルがｘ_ｐ＝｛０，１，２，３，４｝、ノードｑのラベルがｘ_ｑ＝｛２，３，４，５，６｝というように、ノードごとにラベルの範囲が異なる場合。
＜ケース２＞ノードｐのラベルがｘ_ｐ＝｛０，１，２，３，４｝、ノードｑのラベルがｘ_ｑ＝｛０，１，２，３，４，５，６｝というように、ノードごとにラベルの種類数が異なる場合。
＜ケース３＞ノードｐのラベルに対応する数値が小数を含む場合。例えば２ノード項をＷ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ‖の形で書いたときに、各ラベルに対応する数値がｘ_ｐ＝｛０．３，１．３，２．３，３．３，４．３｝、ｘ_ｑ＝｛０，１，２，３，４，５，６｝となる場合。

以上に示したケースは多くの応用で起こりうる。例えば、動画像の各画素をノード、動きベクトルをラベルとして最適なラベルを選択する問題を考えると、それぞれのケースは下記の状況で起こる。
＜ケース１＞まず、動きベクトルを低い精度で求め、次に得られた動きベクトルを中心とした小さな範囲で動きベクトルをより高い精度で求める場合を考える。このとき、２回目の動きベクトルの算出では、ノードごとにラベルの範囲が異なる場合を扱う。

＜ケース２＞動画像の動きベクトル推定では、前のフレームからの動きベクトルの変化は比較的小さいと考えられるから、前のフレームで求めた動きベクトルを中心とした小さな範囲のみを次のフレームの動きベクトル候補とすれば高速化できる。このとき、画像中で動きベクトルの変化（変化とは、例えば、過去２フレームの動きベクトルから求めた加速度）が穏やかな領域では、誤推定の可能性は少ないため動きベクトル候補を少なくすれば余計な計算を減らせるが、動きベクトルの変化が激しい領域では、誤推定の可能性を減らすために動きベクトル候補を多くする必要がある。このとき、ノードごとにラベルの種類数が異なる場合を扱う。

＜ケース３＞ケース２と同様に前のフレームで求めた動きベクトルを中心とした小さな範囲を動きベクトル候補とし、ケース１のように動きベクトルの推定を多段構成で行う場合を考える。前のフレームでは動きベクトルを高い精度で求めているから、次のフレームで低い精度で求めるときには、中心とする動きベクトルは小数を含むことになる。このとき、各ノードのラベルとして小数を扱うことになる。

（拡張高速ＢＰ）
（用語の定義、基本部分）
本実施形態では、高々２つの変数（ノード）から成る項の総和として表現されるエネルギー関数の近似的な最小化問題を扱う。本実施形態では近似的な最小化を含めて単に最小化と呼んでいる。エネルギー関数の最小化とは、エネルギー関数を最小にするラベルを割り当てることをいう。１変数（１ノード）から成る項をＶ、２ノード（２ノード）から成る項をＷとすると、エネルギー関数は次の式で表せる。

ただし、ラベルの添え字ｉやｊはノードを表し、ｘ_ｉ＝｛０，１，・・・，Ｌ_ｉ−１｝はラベルを表し、整数であることが多い。Ｎはエネルギー関数が含む各２ノード項を構成する変数対の集合をあらわす。なお、この和は積の形で書かれていても良い。関数が積で与えられる場合の解き方は後述する。解くとは、エネルギー関数を最小化するようにラベルを割り当てることをいう。また便宜上、和と積をあわせて集約と呼ぶ。エネルギー関数の最小化においては、上述したように、繰り返し行われるメッセージ値の更新の計算を効率的に行うことが重要である。なお、上記の符号を反転したエネルギー関数

の最大化問題は、符号を変えれば同じ手法で解ける。また、便宜上、最小値を選択する場合の最小値と最大値を選択する場合の最大値をあわせて最適値と呼ぶことにする。

≪ケース１，２：Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖、Ｌ_ｉが画素ごとに異なる場合≫
各２ノード項がＷ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖（この一般形をＷ（ｘ_ｐ−ｘ_ｑ＋ｒ）と書くこともある）で与えられ、ラベル数Ｌ_ｉが画素ごとに異なる場合を考える。ｒ_ｐｑは定数で、２ノード項ごとに異なっていても良いものとする。このＷ、Ｌ_ｉに対するメッセージ値は、上記の高速ＢＰの手法では算出できない。

本実施形態によれば、次の手法でメッセージ値を更新するため、２ノード項が前述のＷ（ｘ_ｐ，ｘ_ｑ）＝ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖で与えられる場合でも算出できる。このための手法について図１から図５を参照して説明する。
本実施形態の画像処理装置の本質的な手法である更新について、はじめにｒ_ｐｑが整数の場合について説明する。ｒ_ｐｑが整数でない場合について、および本実施形態の画像処理装置全体の流れについては後述する。なお、ラベルｘ_ｐ、ｘ_ｑがそれぞれ０〜３、０〜４でｒ_ｐｑ＝−１の場合の例を図３から図５に示す。図３から図５において、縦方向はｍの値を表すものとする。

まず、本実施形態の画像処理装置の更新を行う更新処理装置部分について図２を参照して簡単に説明する。
更新処理装置部分は、メッセージ値保持部２０１、メッセージ初期値入力部２０２、ｈ（ｘｐ）値算出部２０３、ｘｑ一時メッセージ値保持部２０４、ＰｅａｋＸｑ算出部２０５、集約値算出部２０６、最適値選択部２０７、最適値出力部２０８、ノード対（ｐ，ｑ）選択部２０９、ｘｐ値入力部２１０、ｘｑ値入力部２１１、ｘｑビリーフ値算出部２１２、出力部２１３、ノードｑ選択部２１４を含んでいる。

メッセージ値保持部２０１は、各２ノード対（ｐ、ｑ）に対し、ｐからｑへのメッセージ値ｍ_ｐ→ｑ（ｘ_ｑ）を保持、すなわち格納している。
メッセージ初期値入力部２０２は、各２ノード対（ｐ、ｑ）に対し、ｐからｑへのメッセージ値ｍ_ｐ→ｑ（ｘ_ｑ）の初期値を設定する。
ｈ（ｘｐ）値算出部２０３は、２ノード対（ｐ、ｑ）のノードｑに対するラベル値ｘ_ｐを入力とし、更新後のｍ_ｐ→ｑ（ｘ_ｑ）に含まれる項であってノードｐに依存しノードｑに依存しない項（ｈ（ｘ_ｐ））を、ｍ_ｐ→ｑ（ｘ_ｑ）の初期値を用いて算出する。
ｘｑ一時メッセージ値保持部２０４は、ｘｑ値入力部２１１で取得されたｘ_ｑごとのメッセージ値を一時的に格納する。
ＰｅａｋＸｑ算出部２０５は、与えられたｘ_ｐに対し関数Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）が最適値をとるｘ_ｑであるＰｅａｋＸｑを算出する。
集約値算出部２０６は、ラベル値ｘ_ｐを入力として、ｈ（ｘ_ｐ）とＷ（ｘ_ｐ−ＰｅａｋＸｑ＋ｒ）を集約した値を算出する。
最適値選択部２０７は、集約値算出部２０６が算出した集約値と、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が保持する値から最適値を選択し、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が格納している値を更新する。
最適値出力部２０８は、ラベル値ｘ_ｑおよび差分定数ｋを入力として、関数Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）と関数Ｗ（ｘ_ｐ−（ｘ_ｑ＋ｋ）＋ｒ）との差分をｍ_ｐ→ｑ（ｘ_ｑ＋ｋ）に集約した値と、ｍ_ｐ→ｑ（ｘ_ｑ）とを比較し、その最適値でｍ_ｐ→ｑ（ｘ_ｑ＋ｋ）を更新する。
ノード対（ｐ，ｑ）選択部２０９は、処理すべきノード対（ｐ、ｑ）を選択する。
ｘｐ値入力部２１０は、処理すべきノード対（ｐ、ｑ）のノードｐに対するラベル値ｘ_ｐを選択する。
ｘｑ値入力部２１１は、ノード対（ｐ，ｑ）選択部２０９、ノードｑ選択部２１４で選択されたｑに対応するｘ_ｑ値を取得する。
ｘｑビリーフ値算出部２１２は、ノードｑ選択部２１４からノードｑ、ｘｑ値入力部２１１からラベル値ｘ_ｑ、メッセージ値保持部２０１からメッセージ値を入力し、ラベル値ｘ_ｑに対するｑへのメッセージ値を集約した値であるビリーフ値を算出する。
出力部２１３は、ノードｑに対し、各ラベル値ｘ_ｑのビリーフ値を算出し、それらのうち最適値を与えるラベル値をノードｑのラベル値として出力する。
ノードｑ選択部２１４は、処理すべきノードｑを選択する。

（拡張高速ＢＰ更新ステップ１）まず入力としてエネルギー関数が与えられ（Ｓ１０１）、メッセージ初期値入力部２０２がメッセージ値を初期化（ステップＳ１０２）する。次に、ラベルｘ_ｑ＝｛０，１，・・・，Ｌ_ｑ−１｝に対し、ｍ（ｘ_ｑ）←＋∞とするステップを実行する。このステップは後述の方法を使えばなくても良いためフローチャートでは示していないが、実行する場合は例えばステップＳ１０２とステップＳ１０３の間で実行する。ｍ（ｘ_ｑ）はｘｑ一時メッセージ値保持部２０４で保持され、処理ラベルはｘｑ値入力部２１１で制御される。

（拡張高速ＢＰ更新ステップ２）ノード対（ｐ，ｑ）選択部２０９がノード対を選択し、ｘｐ値入力部２１０がｘ_ｐを取得し、ｈ（ｘｐ）値算出部２０３がｈ（ｘ_ｐ）を算出する（ステップＳ１０３）。ＰｅａｋＸｑ算出部２０５が、ｘｐ値入力部２１０で得られるｘ_ｐに対し、関数ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖が最適値をとるｘ_ｑである下記の式で示されるｙ_ｐ（ｘ_ｑ）（ＰｅａｋＸｑと呼ぶ）を算出する（ステップＳ１０４）。集約値算出部２０６がｈ（ｘ_ｐ）＋ｓ‖ｘ_ｐ−ＰｅａｋＸｑ＋ｒ_ｐｑ‖を計算する（ステップＳ１０５）。最適値選択部２０７がｍｉｎ（ｍ（ｙ_ｐ（ｘ_ｐ）），ｈ（ｘ_ｐ）＋ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖）を計算し、最適値選択部２０７が、この値と、ＰｅａｋＸｑに対するｘｑ一時メッセージ値保持部２０４が保持する値とから最適値を選択し、ｘｑ一時メッセージ値保持部２０４に出力し、ｍ_ｐ→ｑ（ｘ_ｑ）を更新する（ステップＳ１０６）。

すなわち、各ラベルｘ_ｐに対して、次の式を実行する。
ｍ（ｙ_ｐ（ｘ_ｐ））←ｍｉｎ（ｍ（ｙ_ｐ（ｘ_ｐ）），ｈ（ｘ_ｐ）＋ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖）
なお、ｙ_ｐ（ｘ_ｐ）はラベルｘ_ｐごとにラベルｘ_ｑを選択するためのもので、

により求める（なお、この式はｘ_ｐではなくｘ_ｑを決める式である）。図１や図２ではｙ_ｐ（ｘ_ｐ）をＰｅａｋＸｑと記載している。

各ｙ_ｐ（ｘ_ｐ）（＝ＰｅａｋＸｑ）についてｍの値を求めると図３のようになる。縦方向はｍの値の大小を表す。拡張高速ＢＰ更新ステップ２ではｍｉｎを算出する操作を行う。この操作は図４に示すように、ラベルｘ_ｑの範囲の制限により両端では２つ以上のｙ_ｐ（ｘ_ｐ）が同じラベル値となることがあるが、このときｍ（ｘ_ｑ）として最小値を選択する操作を表す。なお、この操作によって、ただ１つのｙ_ｐ（ｘ_ｐ）が対応するｘ_ｑについてはその値がｍ（ｘ_ｑ）として用いられ、対応するｙ_ｐ（ｘ_ｐ）が１つも存在しないｘ_ｑについてはｍ（ｘ_ｑ）に＋∞が設定されることになる。

（拡張高速ＢＰ更新ステップ３）ｘ_ｑ＝｛０，１，・・・，Ｌ_ｑ−１｝の範囲で、ラベル番号ｘ_ｑを１から開始して１ずつ増やしながら以下を実行する。
ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ−１）＋ｓ）
この操作は、図５に示すように、隣接ラベル値でのｍ（ｘ_ｑ）に差分を足した値と注目ラベル値でのｍ（ｘ_ｑ）のうちの小さいほうをｍ（ｘ_ｑ）として選択する操作を表す。

（拡張高速ＢＰ更新ステップ４）ｘ_ｑ＝｛Ｌ_ｑ−２，・・・，０｝の範囲で、ラベル番号ｘ_ｑをＬ_ｑ−２から開始して１ずつ減らしながら以下を実行する。
ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ＋１）＋ｓ）
なお、拡張高速ＢＰ更新ステップ３と４の順序は逆でもかまわない。

図１では拡張高速ＢＰ更新ステップ３と４をまとめてＳ１０７としている。Ｓ１０７は拡張高速ＢＰ更新ステップ３や４の処理を一般化して書いたもので、最適値出力部２０８で実行される。実行される処理は次の通りである。
（ステップＳ１０７）ラベル差分をｋとする。ｋが正ならｋからＬ_ｑ−１まで、ｋが負ならＬ_ｑ−１＋ｋから０までの範囲で、１（ｋ＞０）あるいは−１（ｋ＜０）ずつラベル値をずらしながら、以下を実行する。
ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ＋ｋ）＋ΔＷ）
ここで、ΔＷはラベルをｋずらすことにより生じるＷの差分を表す。拡張高速ＢＰ更新ステップ３はラベル差分ｋ＝１としてＳ１０７を実行する場合、拡張高速ＢＰ更新ステップ４はラベル差分ｋ＝−１としてＳ１０７を実行する場合に相当する。最小値の選択をＡ３０８で行い、選択された最小値はＡ３０５に送られる。

本実施形態の画像処理装置を最も特徴づけるステップは拡張高速ＢＰ更新ステップ２である。
はじめにｈ（ｘ_ｑ）をラベルごとに求める必要があるが、このときｍ（ｘ_ｑ）に単に求めたｈ（ｘ_ｑ）を代入するのではなく、まず各２ノード項が最小になるラベル値ｙ_ｑ（ｘ_ｑ）を調べ、そのラベル値に対応するｍ（ｙ_ｑ（ｘ_ｑ））を更新する。このとき、複数のｘ_ｑに対してｙ_ｑ（ｘ_ｑ）が同一となる場合にはその最小値を選択しておけば、拡張高速ＢＰ更新ステップ３と４を終えた時点での更新メッセージ値は高速化を施さない場合と完全に一致する。

なお、拡張高速ＢＰ更新ステップ１でメッセージ値を無限大としておくのは、何もしない拡張高速ＢＰ更新ステップ２でｍ（ｘ_ｑ）の値が割り当てられないラベルについて、拡張高速ＢＰ更新ステップ３と４として高速ＢＰと同じ手法で正しいメッセージ値が得られるようにするためである。拡張高速ＢＰ更新ステップ３と４では、Ｌの値がノードに依存して変化する点が高速ＢＰの手法と異なるが、高速ＢＰでの手法と回路やプログラムを共用できる。
更新ステップ４が終了した時点で、各ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）の値がｍ（ｘ_ｑ）に得られる。ｍ（ｘ_ｑ）を保持するメモリはｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）を保持するメモリと同一のものでも良い。エネルギー関数の符号が反転している場合には最大値を選択する。

次に、この拡張高速ＢＰ更新ステップ１〜４を用いた全体の流れについて説明する。
（最適化ステップ１）入力としてエネルギー関数が与えられる（Ｓ１０１）。ｔ＝０とする（ｔはメッセージ更新回数を表す）。各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）の初期値を与えるか、全て０で初期化する（Ｓ１０２）。ブロック図では、メッセージ初期値入力部２０２が初期値を与え、メッセージ値保持部２０１でそれを保持している。

（最適化ステップ２）各エッジ（ｐ，ｑ）∈Ｎに対し、ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）を拡張高速ＢＰ更新ステップ１〜４により更新する（Ｓ１０８）。ブロック図では、ノード対（ｐ，ｑ）選択部２０９で処理するノード対を選択している。なお、拡張高速ＢＰステップ２の各ノード対の処理は独立しているため並列に実行できる。

（最適化ステップ３）ｔを１増やす。ｔが予め決めておいた繰り返し回数Ｔより小さいなら、拡張高速ＢＰステップ２に戻る。

（最適化ステップ４）各ｑ、各ラベルに対して、次のビリーフと呼ばれる値を求める（Ｓ１０９）。

ただし、Σ_ｐは（ｐ，ｑ）∈Ｎを満たす全てのｐに対する和を表す。ブロック図ではノードｑ選択部２１４でノードｑ、ｘｑ値入力部２１１で処理するラベルを選択し、ｘｑビリーフ値算出部２１２で求め出力部２１３に送っている。なお、ｘｑビリーフ値算出部２１２で求めたビリーフ値を出力部２１３に送って評価すれば、全てのビリーフ値を保持する必要はない。

（最適化ステップ５）各ｑに対してｂ_ｑ（ｘ_ｑ）を最小とするラベルｘ_ｑ＝｛０，１，・・・，Ｌ_ｑ−１｝を選ぶ（Ｓ１０９）。ブロック図では出力部２１３がステップ５を行う手段である。

次に、本実施形態の画像処理装置の全体の構成を図６に示す。本実施形態の画像処理装置は、画像入力部６０１、ＭＰＵ（演算ユニット）６０２、データ用メモリ６０３、プログラム用メモリ６０４、ラベル出力部６０５を含んでいる。なお、図２に示した装置部分は、ＭＰＵ６０２、データ用メモリ６０３、プログラム用メモリ６０４に対応する。

本実施形態を実施すべきプログラムはプログラム用メモリ６０４に保存される。プログラムはＲＯＭに格納しても良いし、他の装置（例えばハードディスク）を接続してプログラム実行時にＲＡＭに格納しても良い。プログラムはＭＰＵ６０２により実行される。プログラムは、例えば次の指示を出すように作成する。

（１）画像入力部６０１からの入力した画像をデータ用メモリ６０３に格納する。
（２）データ用メモリ６０３にメッセージ値や更新中の一時メッセージ値、ビリーフ値を記録するためのバッファを確保する。
（３）データ用メモリ６０３上のメッセージ値を初期化する。
（４）画像とメッセージ値を参照しながら、拡張高速ＢＰ更新ステップ１〜４にしたがって一時メッセージ値を算出する。
（５）求めた一時メッセージ値を利用してメッセージ値を更新する。
（６）メッセージ値の更新を所定の回数だけ繰り返す。
（７）メッセージ値を参照しながら、ビリーフ値を算出する。
（８）最小のビリーフ値を与えるラベル値を算出する。
（９）ラベル出力部６０５に結果を出力する。

ここで、（１）と（９）以外の処理は全てＭＰＵ６０２、データ用メモリ６０３、プログラム用メモリ６０４のみで行われる。なお、データ用メモリ６０３とプログラム用メモリ６０４は分離されていなくても良い。例えば、パーソナルコンピュータのようにプログラムとデータのいずれにも利用できるメモリを備えていれば、それを必要なメモリサイズで区切り、プログラムとデータの両方に利用できる。

≪ケース３：コスト差分値が小数の場合≫
ｒ_ｐｑが整数でない場合についてメッセージ値の更新手法は拡張高速ＢＰ更新ステップ２を次のように変更する。
（拡張高速ＢＰ更新ステップ２（ｒ_ｐｑが整数でない場合））各ラベルｘ_ｑに対して、まず次の集合を求める。

ここで、

はそれぞれｘの小数点以下を切り捨てた値、ｘの小数点以下を切り上げた値を表し、コンピュータ言語でfloor(x)、ceil(x)として知られる関数に対応する。求めた｛ｙ_ｐ（ｘ_ｐ）｝内の各要素をｙ_ｐ（ｘ_ｐ）として、以下を実行する。

ｍ（ｙ_ｐ（ｘ_ｐ））←ｍｉｎ（ｍ（ｙ_ｐ（ｘ_ｐ）），ｈ（ｘ_ｐ）＋ｓ‖ｘ_ｐ−ｙ_ｐ（ｘ_ｐ）＋ｒ_ｐｑ‖）
なお、ｒ_ｐｑが整数の場合は、切り捨て値と切り上げ値が一致するために、上述した手法と同等になる。

（拡張高速ＢＰ更新ステップ１で無限大をセットしない手法）
先の説明では拡張高速ＢＰ更新ステップ１で全てのラベルに対してｍ（ｘ_ｑ）←＋∞と初期化し、拡張高速ＢＰ更新ステップ３と４ではそれら全てを対象とした。しかし、

は、ｘ_ｐ＝｛０，１，・・・，Ｌ_ｐ−２｝に対しては明らかにｙ_ｐ（０）とｙ_ｐ（Ｌ_ｐ−１）の間の値をとる。そこで、ｙ_ｐ（０）よりも小さいラベル値およびｙ_ｐ（Ｌ_ｐ−１）よりも大きいラベル値に対するメッセージ値の計算は拡張高速ＢＰ更新ステップ１、３、４で行わなくても、拡張高速ＢＰ更新ステップ４の後に他のメッセージとの差分を求めれば、結果は同じになる。具体的な手順は次の通りである。
ｙ_ｐ（Ｌ_ｐ−１）よりも大きいラベル値：ｍ（ｙ_ｐ（Ｌ_ｐ−１）＋１）から順に、ラベル値を１ずつ増やしながらｍ（ｘ_ｑ）←ｍ（ｘ_ｑ−１）＋ｓを実行する。
ｙ_ｐ（０）よりも小さいラベル値：ｍ（ｙ_ｐ（０）−１）から順に、ラベル値を１ずつ減らしながらｍ（ｘ_ｑ）←ｍ（ｘ_ｑ＋１）＋ｓを実行する。
（ステレオマッチング）
エネルギー関数最小化のためのラベル選択問題は様々な産業的応用を持つが、ここではその１つであるステレオマッチングを例として説明する。

入力画像として、２つのカメラの視差が画像の水平方向のずれとしてあらわれる画像（左右のカメラから各１枚）が与えられるものとする。このような画像は、例えば、カメラを正確に並べて配置するか、２つのカメラの配置情報を利用し撮影画像に変換を施すことで得られることが知られている。ステレオマッチングでは、左右の画像の各画素における視差をラベルとし、Ｖに各画素を単独で見た場合の各ラベルの信頼度（値が大きくなるほどｘ_ｉが出現しにくいことを表す）を記述し、Ｎを画素とその近傍の画素の組を全て集めた集合とし、Ｗに近傍画素間の関係（値が大きくなるほど、（ｘ_ｉ，ｘ_ｊ）の組が出現しにくいことを表す）を記述したうえで、Ｅを最小化するように各画素のラベルを選択する。得られた視差とカメラの情報があれば、例えば、それらを用いて各画素の奥行き情報を算出できることが知られている。

本実施形態の画像処理装置を利用することにより、より精度の高いステレオマッチング、および柔軟なステレオマッチングが可能になる。なお、視点の異なる画像が３枚以上ある場合には、例えば１枚の画像を基準として、それ以外の各画像に対してステレオマッチングで対応付けを行って基準画像の各画素の奥行き情報を算出し、基準画像の各画素の奥行き情報を各画像で算出した奥行き情報の平均値とすることで高精度な奥行き情報を算出できる。

（階層的探索）
ステレオマッチングの２ノード項の例として、ここでは次の式を考える。
Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｍｉｎ（ｓ‖ｘ_ｐ−ｘ_ｑ‖，ｄ_ｐｑ）
この式は、線形コストとＰｏｔｔｓの組み合わせによって解ける形をしている。例えば、０〜３０画素の範囲内で１画素単位の視差を求めるのであれば、必要なラベル数は３１として高速ＢＰの手法を適用すれば良いことになる。ところで、奥行き情報を正確に算出するには視差の精度が高いほど良い。しかし、必要なラベル数は視差の精度に比例して増大するため、例えば０．１画素単位の視差を求めるのであればラベル数は３００以上になる。したがって、高い視差の精度を得たい場合、高速ＢＰの手法を用いた場合は多くのメモリが必要となり、また計算時間も増大する。

一方、本実施形態の手法を用いることで、視差の階層的な算出を実現し、必要なメモリ、計算時間を抑制できる。なお、Ｖ（ｘ_ｐ）はｘ_ｐ以外の変数を含まなければどのような形でもかまわないが、例えば次の形が使える。
Ｖ（ｘ_ｐ）＝ｍｉｎ（‖Ｉ_ｌｅｆｔ（Ｘ_ｐ，Ｙ_ｐ）−Ｉ_{ｒｉｇｈｔ}（Ｘ_ｐ−ｘ_ｐ，Ｙ_ｐ）‖，ｃ_ｐ）
Ｉ_ｌｅｆｔ（Ｘ_ｐ，Ｙ_ｐ）、Ｉ_{ｒｉｇｈｔ}（Ｘ_ｐ，Ｙ_ｐ）はそれぞれ座標（Ｘ_ｐ，Ｙ_ｐ）における左画像・右画像の輝度値を表す。また、ｃ_ｐは定数である。

以下、画素をノードとして、画素をｐ、視差の精度をδ、視差の最小値をγ_ｐ、最大値をγ_ｐ＋（Ｌ_ｐ−１）δで表す。なお、ここで用いた記号γ（ガンマ）はｒ（アール）とは別の記号である。フローチャートを図７に示す。ステレオマッチングは、次の流れで実現できる。

（ステップ１）（初期化）Ｋ＝１とし、各画素ｐに対して、γ_ｐが視差の最小値、δが初期探索での視差の精度、γ_ｐ＋（Ｌ_ｐ−１）δが視差の最大値となるように各パラメータを初期化する（Ｓ７０１）。例えば、視差を０〜３０の範囲、１画素精度とするのであれば、γ_ｐ＝０、δ＝１、Ｌ_ｐ＝３１とする。ノードごとに候補ラベルを決定する（ステップＳ７０２）。

（ステップ２）２ノード項を次の形として求め（ステップＳ７０３）、本実施形態のメッセージ更新を用いたＢＰ手法（ステップＳ１０１〜ステップＳ１０９）で各画素にラベルを割り当てる（Ｓ７０４）。

Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｍｉｎ（ｓ‖（ｘ_ｐ−ｒ_ｐ）−（ｘ_ｑ−ｒ_ｑ）‖，ｄ_ｐｑ）
（ステップ３）Ｋ←Ｋ＋１とする。Ｋがしきい値に達していれば終了し、達していなければステップ４に進む。

（ステップ４）δを小さくする（Ｓ７０５）。例えば、δに予め定めた１未満の値（例えば０．５）を掛ける。各ｐの視差の候補は｛γ_ｐ，・・・，γ_ｐ＋（Ｌ_ｐ−１）δ｝であるから、これらの候補が予め決めた基準に合致するようにγ_ｐとラベル数Ｌ_ｐを決める。先の例であれば、ステップ２で求めた視差を中心とした−３δ〜＋３δの範囲となるようにγ_ｐとラベル数Ｌ_ｐを決める。これらの値はｐごとに変えても良い。｛δ，γ_ｐ，Ｌ_ｐ｝の更新後、ステップ２に戻り、先の手順を繰り返す（Ｓ７０６）。

この手法を用いると、ラベル数Ｌ_ｐの増加を抑えながら精度の高い（小さなδでの）視差の推定ができる。なお、スケールを表す数値Ｋがしきい値になるまで繰り返す代わりに、精度δが一定値になるまで繰り返しても良い。その場合、Ｋを計算する必要はない。また、ラベル数Ｌ_ｐは全てのｐに対して同じである必要はない。例えば、画素ｐの視差とその４近傍の画素での視差での最大値と最小値を算出し、γ_ｐ＋（Ｌ_ｐ−１）δが最大値以上、γ_ｐが最小値以下となるように個々のγ_ｐとＬ_ｐを決めれば、物体境界のように視差の変動が激しい位置での推定精度を上げることができる。

なお、先のステップ１では０〜３０の範囲で初期化したが、例えば視差の推定をステレオ動画像に対して行う場合には、例えば、直前の時刻における視差を中心とした一定の範囲としても良い。この場合、２ノード項にあらわれるγ_ｐ−ｒ_ｑが小数とならないように各γ_ｐを量子化しても良いし、先に示した小数を扱う手法を用いても良い。

なお、本発明はより複雑なステレオマッチングでも利用できる。例えば文献「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）」のように、最適化すべきエネルギーのＷ（Ｘｐ，Ｘｑ）項に影響を与えない改良（この改良に限定はされない）は、本発明に何ら変更を加えることなくそのまま利用できる。なお隠れ情報は隠れているかいないかの２値をとるため、従来の手法で最適化すれば良い。

ここまでは、１次元のノードについての説明であるが、以下２次元以上の場合について説明する。
（オプティカルフロー、多次元動きベクトル）
画像処理の分野において、２枚(あるいはそれ以上)の画像中の画素を対応付ける問題は、ステレオマッチング以外にもあらわれる。特に、２次元の画像の対応づけ問題は、動画像に対する動きベクトルの推定問題として様々な応用を持つ（例えば画像圧縮、人物や車両の追跡、超解像）。また、例えばＸ線ＣＴやＭＲＩのような医療用画像は３次元画像であることが多いが、これらの対応付け問題は３次元画像中の画素の対応付け問題になる。

以下、次元を表す番号（１）、（２）を記号の右下に付与する。２次元のノードｑに対するラベルをｘ_ｑ＝（ｘ_（１）ｑ，ｘ_（２）ｑ）で表す。ラベルの範囲をｘ_（１）ｑ＝｛０，１，・・・，Ｌ_（１）ｑ−１｝、ｘ_（２）ｑ＝｛０，１，・・・，Ｌ_（２）ｑ−１｝とする。ノードｑの候補となるラベル数はＬ_ｑ＝Ｌ_（１）ｑＬ_（２）ｑである。２ノード項として、次の形を考える。
Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｓ_（１）‖ｘ_（１）ｐ−ｘ_（１）ｑ＋ｒ_{（１）ｐｑ}‖＋ｓ_（２）‖ｘ_（２）ｐ−ｘ_（２）ｑ＋ｒ_{（２）ｐｑ}‖
Ｖ（ｘ_ｐ）はｘ_ｐ以外の変数を含まなければどのような形でもかまわないが、例えば次の形が使える。
Ｖ（ｘ_ｐ）＝ｍｉｎ（‖Ｉ_ｌｅｆｔ（Ｘ_ｐ，Ｙ_ｐ）−Ｉ_{ｒｉｇｈｔ}（Ｘ_ｐ−ｘ_（１）ｐ，Ｙ_ｐ−ｘ_（２）ｐ）‖，ｃ_ｐ）
Ｉ_ｌｅｆｔ（Ｘ_ｐ，Ｙ_ｐ）、Ｉ_{ｒｉｇｈｔ}（Ｘ_ｐ，Ｙ_ｐ）はそれぞれ座標（Ｘ_ｐ，Ｙ_ｐ）における左画像、右画像の輝度値を表す。また、ｃ_ｐは定数である。

この２ノード項に対する拡張高速ＢＰを実現するには、次の手法でメッセージ更新手法を行えばよい。

（拡張高速ＢＰ更新ステップ１）ラベルｘ_ｑ＝｛（０，０），・・・，（Ｌ_（１）ｑ−１，Ｌ_（２）ｑ−１）｝に対し、ｍ（ｘ_ｑ）←＋∞とする。

（拡張高速ＢＰ更新ステップ２）各ラベルｘ_ｑに対して、次の式を実行する。

ｍ（ｙ_（１）ｐ（ｘ_ｐ），ｙ_（２）ｐ（ｘ_ｐ））←ｍｉｎ（ｍ（ｙ_（１）ｐ（ｘ_ｐ），ｙ_（２）ｐ（ｘ_ｐ）），ｈ（ｘ_ｐ）＋ｓ_（１）‖ｘ_（１）ｐ−ｙ_（１）ｐ（ｘ_ｐ）＋ｒ_{（１）ｐｑ}‖＋ｓ_（２）‖ｘ_（２）ｐ−ｙ_（２）ｐ（ｘ_ｐ）＋ｒ_{（２）ｐｑ}‖）
ただし、ｙの値は

により算出する。ｙ_（２）ｐ（ｘ_ｐ）についても記号の右下に付与されている（１）を（２）に変えただけで全く同じ手法で算出する。なお、コスト差分値が小数の場合も１次元の場合と同じように、切り捨てと切り上げの２つの値それぞれについてｍの更新式を実行すれば良い。

（拡張高速ＢＰ更新ステップ３）次の擬似コードで示した処理を実行する。以下、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.」でステレオマッチングにおいて視差とオクルージョンを交互に最適化した方法と同じように、フロー推定についてもフローとオクルージョンを交互に最適化することで、オクルージョンを考慮したフロー推定が可能になる。この場合でも、フローの最適化において本発明を利用すればあまりコストを増やさずに高い精度を達成できる。

Ｖは画素ごとに見た場合の対応位置の選択基準を与えており、これらの例のように別のモデルを与えることで画素ごとの推定精度を上げていくことが可能である。Ｖを改良した場合でも、本実施形態の手法はそのまま適用できる。

（斥力）
ここまでの例では２ノード項として、ラベルの差がｒ_ｐｑとなれば最小となる関数を考えてきた。この関数は、２ノード間の関係が「もしノードｐのラベルがｘ_ｐであれば、ノードｑのラベルはｘ_ｐ＋ｒ_ｐｑに近いことが望ましい」ことを定式化したものである。ところで、応用によっては、２ノード間の関係が「もしノードｐのラベルがｘ_ｐであれば、ノードｑのラベルはｘ_ｐ＋ｒ_ｐｑから遠いことが望ましい」という関係の定式化が必要である。その場合にあらわれる２ノード項として、例えば以下の関数がある。ｄ_ｐｑは定数で、ノードｐ、ｑごとに別の値でも良い。

Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｄ_ｐｑ−ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ_ｐｑ‖
この２ノード項を用いたエネルギー関数のメッセージ更新は、次の手順で高速に行える。

（更新ステップ１）両端のメッセージ値をｍ（０）←＋∞、ｍ（Ｌ_ｑ−１）←＋∞で初期化する。

（更新ステップ２）各ラベルｘ_ｐに対して、まず次の位置を求める。

これを利用して、各ラベルｘ_ｐに対して両端のメッセージ値を次の式で更新する。
ｍ（０）←ｍｉｎ（ｍ（０），ｈ（ｙ_ｐ（ｘ_ｐ））＋ｄ_ｐｑ−ｓ‖−ｒ_ｐｑ＋ｙ_ｐ（ｘ_ｐ）‖）
ｍ（Ｌ_ｑ−１）←ｍｉｎ（ｍ（Ｌ_ｑ−１），ｈ（ｙ_ｐ（ｘ_ｐ））＋ｄ_ｐｑ−ｓ‖Ｌ_ｑ−１−ｒ_ｐｑ＋ｙ_ｐ（ｘ_ｐ）‖）
（更新ステップ３）ｘ_ｑ＝｛１，・・・，Ｌ_ｑ−１｝の範囲で、ラベル番号ｘ_ｑを１から開始して１ずつ増やしながら以下を実行する。
ｍ（ｘ_ｑ）←ｍ（ｘ_ｑ−１）＋ｓ
（更新ステップ４）ｘ_ｑ＝｛Ｌ_ｑ−２，・・・，０｝の範囲で、ラベル番号ｘ_ｑをＬ_ｑ−２から開始して１ずつ減らしながら以下を実行する。

ｍ（ｘ_ｑ）←ｍｉｎ（ｍ（ｘ_ｑ），ｍ（ｘ_ｑ＋１）＋ｓ）
（イメージ復元）
本実施形態の手法は、ノイズが加わった画像のノイズを除去するためにも使える。ノードを画素、ラベルを輝度値とすれば、ノイズ除去問題のエネルギー関数の各項は例えば、
Ｖ（ｘ_ｐ）＝‖ｘ_ｐ−Ｉ_ｐ‖、Ｗ（ｘ_ｐ，ｘ_ｑ）＝ｍｉｎ（ｓ‖ｘ_ｐ−ｘ_ｑ‖，ｄ_ｐｑ）
で表せる。ここでＩ_ｐは入力画像の輝度を表し、ラベルｘ_ｐはノイズ除去後の輝度を表す。輝度は８ビットで表現した場合でも２５６種類あり、例えば１２ビットで表現すれば４０９６種類にもなるため、本実施形態の手法を利用した効率的な最小化は効果的である。このエネルギー関数を先に説明したステレオマッチングと同じ手法で最小化でき、選択されたラベルは各画素におけるノイズ除去後の輝度になる。

（確率積の利用）
主に画像処理であらわれるマルコフ・ランダム・フィールド（ＭＲＦ）やコンディショナル・ランダム・フィールド（ＣＲＦ）と呼ばれる確率モデルでは、次の関数の最大化と等価な処理を行うことがよくある。

この関数の最大化問題は、負の対数−ｌｏｇＰ（ｘ_１，ｘ_２，・・・，ｘ_ｎ）の最小化と等価である。−ｌｏｇＰ（ｘ_１，ｘ_２，・・・，ｘ_ｎ）はＢＰで扱うエネルギー関数の形であるから、ＭＲＦやＣＲＦをモデルとした確率の最適ラベル選択問題は本実施形態の手法により解ける。

（ＴＲＷ／ＷｅｉｇｈｔｅｄＢＰの利用）
ＢＰに似た最適化アルゴリズムとして、２ノード項Ｗ（ｘ_ｐ，ｘ_ｑ）そのものではなく重み付けした値を使ってメッセージ更新を行うアルゴリズムや、２ノード項Ｗ（ｘ_ｐ，ｘ_ｑ）およびメッセージに重み付けを行い、さらに逆方向のメッセージ値も用いてメッセージ更新を行うアルゴリズム（ツリー・リウエイテッド・マックスプロダクト法、以下ＴＲＷ）が知られているが、本実施形態の手法はこれらのアルゴリズムに対してもそのまま適用できる。以下、ＴＲＷを例として本実施形態の手法の適用手法を説明する。

最小化すべきエネルギー関数として、ＢＰと同じ次の関数を考える。

ＴＲＷではＢＰとは違い、次の更新式を用いてメッセージ更新を行う。

ただし、Σ_sは（ｓ，ｐ）∈Ｎかつｓ≠ｑを満たす全てのｓに対する和を表す。ρ_ｐｑは２ノードｐ、ｑ間の重みをあらわす定数で、もしρ_ｐｑ＝１ならＢＰと等しくなる。この更新式は、

として、２ノード項としてＷ’（ｘ_ｐ，ｘ_ｑ）＝（１／ρ_ｐｑ）Ｗ（ｘ_ｐ，ｘ_ｑ）を考えれば、ＢＰで用いる更新式と全く同じ形になる。したがって、各２ノード項に対しＷ’（ｘ_ｐ，ｘ_ｑ）＝（１／ρ_ｐｑ）Ｗ（ｘ_ｐ，ｘ_ｑ）を算出するステップを最初に追加すれば、本実施形態の手法をそのまま適用できる。

（ｓｅｑｕｅｎｔｉａｌＢＰ／ＴＲＷの利用）
メッセージ更新では、ｔ＋１番目のメッセージ値の算出にｔ番目のメッセージ値を用いている。具体的には、

の算出にｔ番目のメッセージ値を用いている。したがって、ｔ番目、ｔ＋１番目のメッセージ値を保持する必要があり、必要なメモリは全メッセージ値の保持に必要なメモリの２倍になる。更新後のメッセージ値ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）をｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）の代わりに用いれば、メッセージ値を上書きするだけでよいため全メッセージ値１つ分のメモリで動作でき、必要なメモリを半分に減らせる。この手法はシーケンシャルＢＰとして知られており、本実施形態の手法でもそのまま利用できる。

（ＥｆｆｉｃｉｅｎｔＭｕｌｔｉｓｃａｌｅＢＰとの併用）
高速ＢＰで示されているように、２×２画素を１つのノードとみなしてＢＰを適用し、その結果を各画素のメッセージ値の初期値として使うことで、画像に対するＢＰの繰り返し回数を減らせる。この手法は本実施形態の手法と違いラベル数を減らせないためにメモリ消費量は減らせないが計算時間は減らせる。この手法はマルチスケールＢＰと呼ばれ、本実施形態の手法とも併用できる。本実施形態の手法を用いたラベル数の削減とマルチスケールＢＰを併用すれば、マルチスケールＢＰによる計算時間の削減と本実施形態の手法による計算時間、メモリ消費量の削減の両方の効果が得られる。

なお、本実施形態で正の無限大を用いた例がいくつかあるが、これらは無限大でない適当な値を使用しても良い（なるべく大きい値であることが望ましい）。

以上の実施形態の画像処理装置により、階層的探索と拡張高速ＢＰを併用することができるので、高速かつ高精度に画素のラベル値を選択することが可能になる。例えばステレオマッチングやオプティカルフロー推定の階層化のように、従来の高速ＢＰでは扱えなかった２ノード項を持つエネルギー関数に対し、エネルギー関数を高速に最小化できるようになる。典型的な例では、計算時間は数分程度、消費メモリも〜１Ｇ程度に改善される。

また、本実施形態によれば、２ノード項ＷがＰｏｔｔｓであらわせず、ノードごとにラベルの範囲や種類数が異なる場合や、ラベルが小数を含む場合であってもＢＰを高速に実行できるため、エネルギー最小化を解くことを必要とする場合でも大幅な高速化が実現できる。

なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。

実施形態の画像処理装置の動作の一例を示す図。実施形態の画像処理装置の更新処理装置部分のブロック図。メッセージ値更新の具体的な一例を示す図。メッセージ値更新の具体的な一例での拡張高速ＢＰ更新ステップ２を示す図。メッセージ値更新の具体的な一例での拡張高速ＢＰ更新ステップ３を示す図。実施形態の画像処理装置のブロック図。ステレオマッチングを行う場合の動作の一例を示すフローチャート。

符号の説明

２０１・・・メッセージ値保持部、２０２・・・メッセージ初期値入力部、２０３・・・ｈ（ｘｐ）値算出部、２０４・・・ｘｑ一時メッセージ値保持部、２０５・・・ＰｅａｋＸｑ算出部、２０６・・・集約値算出部、２０７・・・最適値選択部、２０８・・・最適値出力部、２０９・・・ノード対（ｐ，ｑ）選択部、２１０・・・ｘｐ値入力部、２１１・・・ｘｑ値入力部、２１２・・・ｘｑビリーフ値算出部、２１３・・・出力部、２１４・・・ノードｑ選択部、６０１・・・画像入力部、６０３・・・データ用メモリ、６０４・・・プログラム用メモリ、６０５・・・ラベル出力部。

Claims

画像の各画素に対応付けられた各ノードｐ、ｑのそれぞれのラベル値ｘ_ｐ、ｘ_ｑを変数とし、ノード対（ｐ、ｑ）およびｐとｑに関連付けられた定数ｒを用いて定義された関数Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）と関数Ｖ（ｘ_ｐ）を含む、最小化するｘ_ｐを求めるエネルギー関数を設定する設定手段と、
各ノード対（ｐ、ｑ）に対応した各評価関数ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）を格納する格納手段と、
前記各評価関数のｔ＝０での値である初期値を設定する設定手段と、

にしたがってｈ（ｘ_ｐ）を設定する設定手段と、
ｘ_ｐごとにＷ（ｘ_ｐ−ｘ_ｑ＋ｒ）が最適値になるｘ_ｑであるＰｅａｋＸｑを算出する算出手段と、
ｘ_ｐごとにｈ（ｘ_ｐ）とＷ（ｘ_ｐ−ＰｅａｋＸｑ＋ｒ）との和または積である集約値を算出する算出手段と、
前記集約値と前記ｍ^（ｔ） _ｐ→ｑ（ｘ_ｑ）にｘ_ｑ＝ＰｅａｋＸｑを代入した値とを比較し、小さい値の方を最適値として選択してｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ）とする更新手段と、
ｘ_ｐごとに、Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）とＷ（ｘ_ｐ−（ｘ_ｑ＋ｋ）＋ｒ）との差分とｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ＋ｋ）との集約値と、ｍ’^（ｔ） _ｐ→ｑ（ｘ_ｑ）との小さい値の方を最適値としてｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）とする更新手段と、
ｔが繰り返し回数Ｔよりも小さい場合には、ｍ^{（ｔ＋１）} _ｐ→ｑ（ｘ_ｑ）を計算させ、最終的にｍ^（Ｔ） _ｐ→ｑ（ｘ_ｑ）を取得する取得手段と、
ｘ_ｑごとに

を算出する算出手段と、
ｑごとに、ｂ_ｑ（ｘ_ｑ）を最小にするラベル値ｘ_ｑを選択する選択手段と、を具備することを特徴とする画像処理装置。
前記Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）は、ｓをノード対（ｐ、ｑ）に関連付けられた定数として、
Ｗ（ｘ_ｐ−ｘ_ｑ＋ｒ）＝ｓ‖ｘ_ｐ−ｘ_ｑ＋ｒ‖
であることを特徴とする請求項１に記載の画像処理装置。
精度値δの初期値を設定する設定手段と、
ノードｐごとに精度値δに対応するラベル値の候補を設定する設定手段と、
ノード対（ｐ、ｑ）に対し、前記候補に対応するＷ（ｘ_ｐ−ｘ_ｑ＋ｒ）を設定関数として設定する設定手段と、
前記設定関数を使用して、前記選択手段がラベル値を選択するように制御する制御手段と、
精度が上がるように精度値δをδ’に変更する変更手段と、
精度値δ’がある値になるまで、精度値δ’を新たな初期値としてラベル値を選択するように制御する制御手段と、を具備することを特徴とする請求項１または請求項２に記載の画像処理装置。
前記ラベル値は、第１の画像に対する第２の画像の対応点の位置を示すことを特徴とする請求項１から請求項３のいずれか１項に記載の画像処理装置。
前記ノードごとに、関連付けられるラベル値の数が異なることを特徴とする請求項１から請求項４のいずれか１項に記載の画像処理装置。