JP5048381B2 - 三次元形状復元装置 - Google Patents

Description

本発明は、対象物体を撮影した画像（二次元形状）から、その対象物体の三次元形状を復元する三次元形状復元装置に関する。

従来より、対象物体を撮影した画像系列を入力し、対象物体について予め設定された複数の特徴点について、その特徴点の位置を画像面上で追跡し、追跡により得られた特徴点位置系列から対象物体の三次元形状を復元する手法の一つとして、因子分解法が知られている。

この因子分解法は、特徴点位置系列から観測行列を生成し、その観測行列を、対象物体の形状（特徴点の三次元位置）を表す形状行列と、対象物体の運動を表す運動行列とに因子分解するものである。

但し、因子分解法は、観測行列を構成するデータが全て揃っていないと、有効な結果が得られないため、特徴点の隠れ，誤追跡，フレームアウト等によって、特徴点位置系列（ひいては観測行列を構成するデータ）の欠損が頻繁に生じる実問題への適用は、非常に困難であるという問題があった。

これに対して、観測行列から欠損のあるデータを除去して、因子分解計算を実行する方法や、観測行列から欠損のあるデータを除いた部分観測行列を生成し、その部分観測行列を因子分解することで得られた形状行列および運動行列から欠損データを推定し、その推定値によって欠損データが置換された観測行列を用いて因子分解する方法（例えば、特許文献１参照）が提案されている。
特許第３７１１２０３号公報

しかし、従来の方法では、観測行列を生成する際に、正常なデータと異常なデータとを識別する必要があるが、その識別が難しいため、自動処理が困難であり、十分な精度が得られないという問題や、識別した異常なデータを除去したり、欠損データの推定を繰り返す手順等が必要となり、必要な演算量が増大するという問題があった。

また、従来の因子分解法では、三次元形状の推定精度を向上させるためには、計算に用いる画像系列を多くする（即ち、観測行列の次元を大きくする）必要があり、それに応じて、計算量（計算時間）が指数関数的に増大する。

そして、計算量が増大すると、三次元形状復元装置を車両に搭載する場合等のように、使用できる演算器の計算能力に制限がある場合には、適用することが困難になるという問題もあった。

本発明は、上記問題点を解決するために、三次元形状復元装置において、観測行列（特徴点位置系列）中に信頼性の低いデータが存在していても因子分解の実行を可能とすること、更には、計算時間を低減することを目的とする。

上記目的を達成するためになされた本発明の三次元形状復元装置では、観測行列生成手段が、対象物体を撮影した画像系列を入力し、対象物体について予め設定された複数の特徴点の位置を画像面上で追跡し、その特徴点の位置を時系列的に並べてなる観測行列を生成し、信頼性評価手段が、特徴点毎に、特徴点の位置の信頼性を表す位置信頼性評価値を生成する。

そして、分解手段が、観測行列生成手段にて生成された観測行列を因子分解して、対象物体の三次元形状を表す形状行列および対象物体の運動を表す運動行列を生成する。このとき、分解手段は、位置信頼性評価値の高い特徴点を重視した因子分解を行う。

つまり、本発明の三次元形状復元装置では、ノイズが重畳する等して信頼性の低い特徴点の除去を行わなくても、信頼性（位置信頼性評価値）の高い特徴点での当てはまりが良くなるように因子分解が行われる。

従って、本発明の三次元形状復元装置によれば、特徴点の隠れ，誤追跡，フレームアウト等に基づく信頼性の低いデータや観測行列中に存在していたとしても、これを除去するといった余計な手順を必要とすることなく、因子分解を確実に実行することができ、しかも信頼性の高いデータほど重視されるため、生成される形状行列や運動行列の精度を向上させることができる。

また、本発明において、分解手段は、観測行列をＷ，運動行列をＭ，形状行列をＳとし、行列Ｘの各要素の値を確率的に示した分布をｑ（Ｘ）で表記するものとして、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、尤度Ｐ（Ｗ｜Ｍ，Ｓ）を最大にする運動行列Ｍの試験分布ｑ（Ｍ）および形状行列Ｓの試験分布ｑ（Ｓ）を、形状行列の試験分布ｑ（Ｓ）の分散に信頼性評価値を反映させた上で求めることにより、確率的に因子分解計算を実行するように構成されている。

このように構成された本発明の三次元形状復元装置によれば、運動行列Ｍや形状行列Ｓがベイズ推定によって確率的に求められることになるため、画像面上で追跡した特徴点の位置に重畳されているノイズの影響を、大幅に軽減することができると共に、形状行列の試験分布ｑ（Ｓ）の分散に信頼性評価値が反映されるため、信頼性評価値の高い特徴点を重視した因子分解を実現することができる。

そして、このようにベイズ推定を用いた因子分解を行う場合、分解手段は、請求項２に記載のように、推定手段が、運動行列の試験分布ｑ（Ｍ）を固定して形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、形状行列の試験分布ｑ（Ｓ）を固定して運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行し、終了手段が、予め設定された終了条件が満たされている場合に、推定手段を終了させ、最終的に得られた形状行列の試験分布ｑ（Ｓ）が形状行列の分布であるものとして、その形状行列の分布から形状行列を得るように構成されていてもよい。

この場合、分解手段は、変分ベイズ推定法を実現するものであり、第１のステップがＥステップ、第２のステップがＭステップに相当する。
但し、変分ベイズ推定法を用いても、計算に用いる画像系列の増加（即ち、観測行列の次元の増大）に応じて、計算量（計算時間）が指数関数的に増大することを避けることはできない。

そこで、分解手段は、請求項３に記載のように、推定手段が、観測行列生成手段にて観測行列が生成される毎に、運動行列の試験分布ｑ（Ｍ）を固定して前記形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、形状行列の試験分布ｑ（Ｓ）を固定して運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて実行し、更新手段が、推定手段にて推定された形状行列の試験分布ｑ（Ｓ）に基づいて、出力とする形状行列の分布を更新し、終了手段が、予め設定された終了条件が満たされている場合に、推定手段及び更新手段を終了させ、更新手段によって最終的に得られた形状行列の分布から形状行列を得るように構成されていることが望ましい。

つまり、推定手段では、計算に用いる画像系列を制限して、第１のステップ（Ｅステップ），第２のステップ（Ｍステップ）を１回ずつ実行し、これを繰り返し実行すると共に、推定手段から計算結果が得られる毎に、その計算結果に基づいて形状行列の分布を更新する。これは、変分ベイズ推定法を逐次計算にて実現（いわゆるオンライン化）するものである。

このように、本発明の三次元形状復元装置によれば、使用する画像系列を複数に分割し、その分割した画像系列毎に因子分解を実行するため、画像系列を増やしても、処理量は推定手段での１回の処理を単位として加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。

また、上述したように、分解手段を、変分ベイズ推定法にて因子分解を実行するように構成した場合、請求項４に記載のように、エネルギー算出手段が、推定手段にて試験分布ｑ（Ｓ），ｑ（Ｍ）の推定結果が得られる毎に、試験事後分布ｑ（Ｓ，Ｍ）＝ｑ（Ｓ）ｑ（Ｍ）に関する自由エネルギーを算出し、終了手段は、エネルギー算出手段での算出結果が、予め設定されたエネルギー閾値より大きいことを終了条件とするように構成されていてもよい。

もともと、変分ベイズ推定法では、自由エネルギーが最大となるように運動行列や形状行列を推定しているため、この自由エネルギーを終了条件として用いることにより、所望の精度を有した形状行列を確実に取得することができる。

また、終了手段は、請求項５に記載のように、推定手段にて形状行列の試験分布ｑ（Ｓ）を算出した回数が、予め設定された回数閾値より大きいことを終了条件とするように構成されていてもよい。

この場合、自由エネルギーを算出する必要がないため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。
ところで、信頼性評価手段は、請求項６に記載のように、位置信頼性評価値として観測行列の分布の分散値を用いるように構成されていてもよい。特に、特徴点の追跡をベイズ推定を用いて行っている場合には、特徴点の位置（観測行列の各要素）の分布の分散値が必ず算出されるため、これを利用することにより処理負荷を軽減することができる。

そして、当該三次元形状復元装置は、請求項７に記載のように、自動車に搭載して使用するように構成されていてもよい。
自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも適用することができる。

以下に本発明の実施形態を図面と共に説明する。
［第１実施形態］
図１は、本発明を適用した三次元形状復元装置としての画像処理装置１の全体構成を示すブロック図である。なお、画像処理装置１は、車両に搭載され、ドライバの頭部を撮影した画像から、ドライバのわき見や居眠り等を検出するための情報である頭部姿勢を時系列的に求める処理を実行する装置である。

＜全体構成＞
画像処理装置１は、ドライバシートの前面上部又は下部（例えば、メータ内）に設置され、図１に示すようにドライバの頭部を含む画像を取得するカメラ２と、カメラ２が取得した画像を取り込む画像取込装置３と、初期化指令が入力された場合に、画像取込装置３を介して取り込んだ入力画像ｚ_t （ｔは時刻を表す）から、ドライバの顔の特徴を表す複数（Ｎ_f個）の特徴点の位置（画像面上での二次元座標：以下、特徴点位置という）ｐ⁽ⁿ⁾（ｎ＝１，２，…Ｎ_f）を抽出する初期特徴点抽出部４と、初期特徴点抽出部４が抽出した特徴点位置ｐ⁽ⁿ⁾ のそれぞれについて、その特徴点位置ｐ⁽ⁿ⁾ を中心とする予め設定された大きさの画像をテンプレートｔｐ⁽ⁿ⁾ として切り出して記憶するテンプレート生成・記憶部５と、特徴点の三次元位置（ひいては、頭部の三次元形状）を規定する顔モデル（形状行列）Ｃ_a を供給する顔モデル供給部６と、入力画像ｚ_t ，テンプレートｔｐ⁽ⁿ⁾ ，顔モデルＣ_a に基づいて、時系列ベイズ推定を行うことにより、頭部姿勢を表すアフィンパラメータＡ_t を推定すると共に、その推定の過程で算出される各特徴点の位置の分布を規定するための特徴点推定分布パラメータ（平均値ｘ_t ，分散Ｖ_t ）を、顔モデル供給部６に供給する推定部８と、推定部８が推定したアフィンパラメータＡ_t に従って、頭部姿勢θｘ，θｙ，θｚを算出する頭部姿勢算出部９から構成される。

このうち、初期特徴点抽出部４では、図２（ａ）に示すように、目尻，目頭，口元、鼻（鼻の穴，鼻の中心，鼻の左右端）等を特徴点として抽出する。但し、本実施形態では、左右目尻、左右目頭，左右口元，鼻の中心の７（＝Ｎ_f）個の点を抽出するものとする。また、初期化指令は、ドライバが正面を向き、全ての特徴点が映っている（特徴点の隠れのない）画像をカメラ２で撮影した時に、ドライバの操作によって入力されるものとする。

テンプレート生成・記憶部５が生成するテンプレートｔｐ⁽ⁿ⁾ は、具体的には、図２（ａ）中の点線白枠にて示すようなものとなる。但し、図では、左口元、鼻の中心、左目尻についてのみ示す。

また、顔モデル供給部６から供給される顔モデルＣ_aは、特徴点（図２（ｂ）に示す線図の各頂点）の位置を三次元的に規定するものであり、図３に示すように、画像面における水平方向の座標軸をｕ軸、垂直方向の座標軸をｖ軸、３次元空間における水平方向の座標軸をＸ軸、垂直方向の座標軸をＹ軸、光軸方向（奥行き方向）の座標軸をＺ軸として（１）（２）式で表される。但し、ｓ⁽ⁿ⁾は、初期特徴点抽出部４で抽出されたｎ番目の特徴点（特徴点位置ｐ⁽ⁿ⁾＝（ｕ⁽ⁿ⁾，ｖ⁽ⁿ⁾））の３次元空間内での座標である。

次に、推定部８で推定されるアフィンパラメータＡ_t は、実世界（三次元のＸＹＺ空間）における座標変換を、射影面である画像面（二次元のｕｖ平面）上での座標変換として表す際に用いる変換行列であり、（３）（４）式に示す関係を有する。

但し、ある特徴点の三次元座標を（Ｘ，Ｙ，Ｚ）^T，座標変換後の三次元座標を（Ｘ’，Ｙ’，Ｚ’）^T、これら三次元座標を画像面上に射影した二次元座標、すなわち画像面上での特徴点の座標を（ｕ，ｖ）^T，（ｕ’，ｖ’）^Tとして、Ｘ＝ｕ，Ｙ＝ｖ（すなわちＸ’＝ｕ’，Ｙ’＝ｖ’）であると仮定する。

そして、頭部姿勢算出部９では、頭部姿勢情報として、頭部の３軸角度θｘ，θｙ，θｚを（５）〜（７）式を用いて算出する。

即ち、三次元座標をＸ軸を中心に角度θｘだけ回転させる変換行列をＲｘ，Ｙ軸を中心に角度θｙだけ回転させる変換行列をＲｙ，Ｚ軸を中心に角度θｚだけ回転させる変換行列をＲｚとし、これら変換行列Ｒｔ（＝ＲｘＲｙＲｚ）を作用させることで三次元座標が（Ｘ，Ｙ，Ｚ）^Tから（Ｘ’，Ｙ’，Ｚ’）^Tに移動するものとすると、これらの関係は、（８）（９）式で表される。

また、上記回転に加えて、三次元座標をＸ軸にｔ₁，Ｙ軸にｔ₂，Ｚ軸にｔ₃だけ平行移動させる変換を考慮した場合、その変換式は（１０）式で表される。

そして、（３）式と（１０）式を比較することにより、アフィンパラメータａ₁〜ａ₈と（１０）式に示す変換行列の各要素ｒ₁〜ｒ₆，ｔ₁，ｔ₂との対応関係が得られ、また、（９）式を整理することで（１１）〜（１８）に示す関係式が得られる。

つまり、これら（１１）〜（１８）式から、（５）〜（７）式が導かれるのである。

そして、時刻ｔにおける頭部姿勢を表すアフィンパラメータＡ_tが得られれば、そのアフィンパラメータＡ_tを顔モデルＣ_aに作用させることで、時刻ｔにおける画像に映っていない隠れた特徴点の位置も含めて、時刻ｔにおける全ての特徴点の位置を予測することができる。

また、このようなアフィンパラメータＡ_tには、カメラパラメータ（カメラの焦点距離、設置場所、向き等）及びドライバの頭部姿勢の変化を全て反映させることができる。
＜推定部の構成＞
ここで、図４は、推定部８が扱う状態空間モデルを表すグラフである。

図４に示すように、推定部８は、時刻ｔにおける顔モデルの姿勢を表すアフィンパラメータａ_t を隠れ状態として求める上位層と、時刻ｔにおける特徴点群の位置ｘ_t＝（ｘ_t ⁽¹⁾，ｘ_t ⁽²⁾， …，ｘ_t ^(Nf)）^Tを隠れ状態として求める特徴点毎に設けられたＮ_f個の下位層（図では一つだけ示す）とを備えており、時刻１からｔまでに入力された入力画像系列ｚ_1:tから、これら隠れ状態Ａ_t，ｘ_tを、時系列ベイズ推定により推定するものである。

なお、時系列ベイズ推定は、図５に示すように、状態変数をすべて確率変数として扱い、隠れ状態に関する時刻ｔ−１における推定結果（推定分布）から時刻ｔにおける予測分布を求めると共に、時刻ｔにおける観測データから、検出すべき隠れ状態であるもっともらしさ（尤度）を求め、これら予測分布と尤度を考慮して時刻ｔにおける推定結果（推定分布）を得るという手順を繰り返すことで、逐次的に隠れ状態を推定するものである。

つまり、推定部８では、入力画像系列（観測データ）ｚ_1:tに基づいて、アフィンパラメータ（頭部姿勢の隠れ状態）Ａ_tを推定した事後確率分布（推定分布）ｐ（Ａ_t｜ｚ_1:t）を求めることになり、これを数式で記述すると（１９）（２０）式で表される。

ここで、ｐ（Ａ_t｜ｚ_1:t-1）はアフィンパラメータＡ_tの事前確率分布（予測分布）、ｐ（ｘ_t｜Ａ_t，ｚ_1:t-1）は特徴点群の位置ｘ_tの事前確率分布（予測分布）、ｐ（x_t｜Ａ_t，ｚ_1:t-1），ｐ（ｚ_t｜ｘ_t）は尤度を表す。

そして、（２０）式の部分を下位層がパーティクルフィルタを用いて推定し、（１９）式の部分を上位層がカルマンフィルタを用いて推定する。
なお、あるパラメータｆの確率分布がガウス分布（正規分布）に従う場合、その確率分布は、平均をμ，分散をΣとすると、（２１）式で表すことができる。つまりこの場合、パラメータの確率分布の計算は、実際には、その平均μと分散Σを求めれば十分である。

次に、図６は推定部８の具体的な構成を表すブロック図である。

図６に示すように、推定部８は、^(Nf)特徴点毎に設けられたＮ_f個の追跡器ＢＫ⁽ⁿ⁾からなり、それぞれが一つの特徴点の位置ｘ⁽ⁿ⁾を追跡して、時刻ｔにおけるその特徴点の位置ｘ_t ⁽ⁿ⁾を推定した確率分布をガウス近似した確率分布を規定する特徴点推定分布パラメータ（平均、標準偏差）を生成する特徴点追跡部１０と、各追跡器ＢＫ⁽ⁿ⁾にて算出された特徴点推定分布パラメータに基づき、カルマンフィルタを用いてアフィンパラメータＡ_t及びアフィンパラメータの予測値Ａ_t+1，その分散Ｖ_t+1（以下、予測値等という）を算出するアフィンパラメータ算出部３０と、アフィンパラメータ算出部３０にて算出されるアフィンパラメータの予測値等を保持する予測値保持部４０と、予測値保持部４０にて保持された時刻ｔ−１において算出されたアフィンパラメータの予測値等に基づいて、時刻ｔの特徴点群の位置ｘ_t（＝（ｘ_t ⁽¹⁾, ｘ_t ⁽²⁾, …，ｘ_t ^(Nf)）^T）を予測した確率分布ｐ（ｘ_t｜Ａ_t）を規定する上位予測分布パラメータ（平均値，分散）を算出して、追跡器ＢＫ⁽ⁿ⁾のそれぞれに供給する予測分布パラメータ算出部４１とを備えている。

つまり、アフィンパラメータ算出部３０が上述の上位層、特徴点追跡部１０を構成する各追跡器ＢＫ⁽ⁿ⁾が上述の下位層に相当する。
なお、このような推定部８を実現する具体的な構成は、例えば、特願２００５−３６８１２４号等に詳述されているため、ここではその詳細についての説明を省略し、以下では、特徴点推定分布パラメータの生成に関係する追跡器ＢＫ⁽ⁿ⁾の概要についてのみ説明する。

＜追跡器の概要＞
まず、特徴点追跡部１０を構成する追跡器ＢＫ⁽ⁿ⁾は、パーティクルフィルタによって一つの特徴点を追跡するものである。ここでは、そのパーティクルフィルタの動作概要を、図７に示す説明図を参照して説明する。

図７に示すように、パーティクルフィルタでは、推定したい対象の状態（特徴点の位置）の実現値（画像面上の座標）をパーティクルで表し、予測，観測，リサンプルという三つの過程を繰り返し実行することで、時系列的にパーティクルの分布を求める。なお、パーティクルフィルタでは、カルマンフィルタとは異なり、ガウス分布に限らず任意の確率分布を対象とすることができる。

まず、予測過程では、推定したい対象の動きを考慮して、状態空間（ここでは画像面上）での各パーティクルの状態（ここでは画像面上での位置）を遷移させ、推定したい対象の動きに加わるノイズを考慮してパーティクルを確率的に散らすことで，推定したい対象が存在しそうな位置にパーティクル群を配置する。これにより、推定したい対象の状態を予測した確率分布が、パーティクル群の位置と数により離散的かつ近似的に表現されることになる。

次に、観測過程では、各パーティクルの状態が推定したい対象の状態であることの尤度（ここでは特徴点を表すテンプレートとの正規化相関値）を、尤度に応じてパーティクルの重みを算出する。これにより、推定したい対象の状態の確率分布が、重みのあるパーティクルによって表現されることになる。

また、リサンプル過程では、重みの小さいパーティクルを消去すると共に重み大きいパーティクルを重みのない複数のパーティクルに増殖させることにより、推定したい対象の状態の確率分布を、重みのあるパーティクルによる表現から、重みのないパーティクルにより表現に確率的に変換する。

そして、このリサンプル（リサンプル過程で生成）されたパーティクルを用いて上述の予測過程を行うことにより、パーティクルで表された推定したい対象の状態の確率分布が時系列的に求められることになる。

ここで、パーティクルの個数をＮ_p、パーティクルの画像面上での座標をｐ_i＝（ｕ_i，ｖ_i）、観測過程にて尤度に応じて算出されるパーティクルの重みをｗ_i（ｉ＝１，２，…，Ｎ_p）とすると、観測によって得られる推定したい対象の状態（即ち、特徴点の位置）の確率分布は、（２２）（２３）式に示す平均値、及び（２４）式（２５）式に示す分散値によって表される。

つまり、時刻ｔにおけるｎ番目の特徴点の位置を表す確率分布の平均値をｘ_t ⁽ⁿ⁾＝（ｘ⁽ⁿ⁾，ｙ⁽ⁿ⁾）、分散をＶ_t ⁽ⁿ⁾＝（Ｖ_x ⁽ⁿ⁾，Ｖ_y ⁽ⁿ⁾）として、ｘ_t＝（ｘ_t ⁽¹⁾，ｘ_t ⁽²⁾，…，ｘ_t ^(Nf)）、Ｖ_t ＝（Ｖ_t ⁽¹⁾，Ｖ_t ⁽²⁾，…，Ｖ_t ^(Nf)）が、特徴点推定分布パラメータとして、各時刻ｔ毎に顔モデル供給部６に供給されることになる。

＜顔モデル供給部＞
次に、本発明の主要部である顔モデル供給部６について詳述する。
顔モデル供給部６は、図１に示すように、頭部を撮影した複数人分の映像から、それぞれ特徴点を抽出し、その抽出した結果に基づいて設定された平均的な特徴点の配置を表す平均顔モデルＳ_Aを記憶する平均顔モデル記憶部６１と、初期化指令が入力されると起動し、推定部８から逐次供給される特徴点推定分布パラメータｘ_t，Ｖ_tに基づいて、カメラ２に被写体（ドライバ）の特徴点の配置を表す個人顔モデルＳ_Pを、学習，生成する個人顔モデル学習部６２と、個人顔モデル学習部６２からの切替信号に従って、平均顔モデル記憶部６１に記憶された平均顔モデルＳ_A、個人顔モデル学習部６２にて生成された個人顔モデルＳ_Pのいずれかを、顔モデルＣ_aとして選択して推定部８に供給する顔モデル切替部６３とを備えている。

＜学習の概要＞
なお、個人顔モデル学習部６２では、以下に示す手法により学習を実行する。
即ち、求めるべき三次元形状（即ち、個人顔モデルＳ_P）を表す形状行列をＳ，特徴点推定分布パラメータｘ_tの系列に基づいて生成される観測行列をＷ，形状行列Ｓの運動を表す運動行列をＭとすると、これらは、Ｗ＝ＭＳの関係を有する。そして、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、尤度Ｐ（Ｗ｜Ｍ，Ｓ）が最大となる運動行列の試験分布ｑ（Ｍ），形状行列の試験分布ｑ（Ｓ）を求めることにより、運動行列Ｍ及び形状行列Ｓを推定する。

これは、（２６）式に定義された自由エネルギーＦ［ｑ（Ｍ，Ｓ）］を最大化するＭ，Ｓを求めることと等価であり、具体的には、（２７）（２８）式を解けばよい。

但し、この計算は簡単には解けないため、（２９）式が成立すること（即ち、運動行列Ｍと形状行列Ｓとに独立な生起確率を仮定すること）、及び、モデル構造（形状行列Ｓ）と運動構造（運動行列Ｍ）の生起確率は、それぞれ正規分布に従うものと仮定し、いわゆる変分ベイズ法を適用して計算を行う。

これにより、求める分布も正規分布に従うため、解くべき問題を、分布の形状を表す特徴点推定分布パラメータｘ_t，Ｖ_tの最適化というように簡略化できる。

そして、変分ベイズ法のＥステップでは、運動行列の分布ｑ（Ｍ）を固定した上で、自由エネルギーを最大化する計算を行えばよく、結局、（２７）（２８）式を変形した（３０）（３１）式からなる方程式を、ラグランジュの未定乗数法を用いて解くことで得られた（３２）式の計算を実行すればよい。但し、式中で記号に付与されている＾（ハット）は、更新値を表すものとする（以下同様）。

同様に、変分ベイズ法のＭステップでは、形状行列の分布ｑ（Ｓ）を固定した上で、自由エネルギーを最大化する計算を行えばよく、結局、（２７）（２８）式を変形した（３３）（３４）式からなる方程式を、ラグランジュの未定乗数法を用いて解くことで得られた（３５）式の計算を実行すればよい。

＜個人顔モデル学習部での処理＞
ここで、個人顔モデル学習部６２が実行する処理を、図８に示すフローチャートに沿って説明する。なお、本処理は、初期化指令の入力により起動される。

但し、ｎ番目の特徴点の三次元座標ｓ⁽ⁿ⁾、求めたい３次元形状を表す形状行列Ｓ、形状行列Ｓに要素がオール１の行ベクトルを加えた同次形状行列〜Ｓを（３６）〜（３８）式で表すものとする。なお、平均顔モデルＳ_A及び個人顔モデルＳ_Pは、同次形状行列 〜Ｓの形式で表現されているものとする。

また、運動のｉ番目の基底ベクトルをｍ_iとして、運動行列Ｍ，同次運動行列〜Ｍ，Θ、ｄ番目の運動行列の行ベクトルθ_dを（３９）〜（４２）式で表すものとする。

本処理が起動すると、まず、Ｓ１１０では、顔モデル切替部６３が推定部８に供給する顔モデルＣ_aとして、平均顔モデルＳ_Aを選択するように切替信号を出力して、Ｓ１２０に進む。

これにより、推定部８では、映像が入力される時刻ｔ毎に、平均顔モデルＳ_Aを用いたアフィンパラメータＡ_t、及び特徴点推定分布パラメータｘ_t，Ｖ_tの算出が行われる。また、頭部姿勢算出部９では、推定部８にて求められたアフィンパラメータＡ_tに基づいて、頭部姿勢θｘ，θｙ，θｚの算出が行われる。

Ｓ１２０では、運動行列の試験分布ｑ（Ｍ）を初期化してＳ１３０に進む。
なお、形状行列の試験分布ｑ（Ｓ）は、形状行列Ｓとその分散を表すΣ_S（以下、形状行列の分布パラメータという）によって規定され、運動行列の試験分布ｑ（Ｍ）は、運動行列Θとその分散を表すＶΘ（以下、運動行列の分布パラメータという）によって規定される。そして、運動行列の分布パラメータΘ，ＶΘは、この初期化によって、全ての要素が、同じ一定値（例えば、０又は１又は０．５）に設定される。

Ｓ１３０では、予め設定された規定個数Ｔの特徴点推定分布パラメータｘ_t，Ｖ_tの取得を完了したか否かを判断し、完了していなければ、完了するまで待機する。
そして規定個数Ｔの特徴点推定分布パラメータｘ_t，Ｖ_tの取得を完了すると、Ｓ１４０に移行し、その取得した特徴点推定分布パラメータｘ_t，Ｖ_tに基づいて、観測行列Ｗおよび信頼度行列Ｑを生成して、Ｓ１５０へ進む。

なお、取得した特徴点推定分布パラメータｘ_tに基づいて作成され、ｎ（ｎ＝１，２，…，Ｎ_f）番目の特徴点の位置（画像面上での２次元座標）を時系列（ｔ＝１〜Ｔ）に並べてなる特徴点ベクトルｘ_1:T ⁽ⁿ⁾を（４３）式で表すものとして、観測行列Ｗは、（４４）式により表される。即ち、観測行列Ｗは、特徴点ベクトルｘ_1:T ⁽ⁿ⁾を列ベクトルして、全ての特徴点について並べることで構成されている。

また、信頼度行列Ｑは、特徴点推定分布パラメータＶ_tに基づいて作成され、ｄ＝２ｔ−１が時刻ｔにおけるｘ成分、ｄ＝２ｔが時刻ｔにおけるｙ成分を表し、ｎ番目の特徴点のｄ番目（ｄ＝１，２，…，２Ｔ）の観測信頼度をσ_x,d ⁽ⁿ⁾で表すものとして、（４５）式により表される。即ち、信頼度行列Ｑは、観測行列Ｗと同じ次数を有する。

Ｓ１５０では、形状行列の試験分布ｑ（Ｓ）を算出する、変分ベイズＥステップとしての処理を実行して、Ｓ１６０に進む。

このＳ１５０では、具体的には、Ｓ１４０にて得られたＷ，Ｑ、及びＳ１２０にて初期設定又は後述するＳ１６０にて更新された最新の運動行列の分布パラメータΘ，ＶΘに基づいて、特徴点（ｎ＝１〜Ｎ_f）毎に、次の（４６）（４７）式を用いて、ｓ⁽ⁿ⁾，Σ_s ⁽ⁿ⁾を求めることにより、形状行列の分布パラメータＳ，Σ_sを更新する。なお、Ｉ_kはｋ次の単位行列である。

但し、（４６）（４７）式は、（３２）式を、特徴点の三次元座標ｓ⁽ⁿ⁾及び分散Σ_s ⁽ⁿ⁾を用いて書き直した（４８）を解くことにより求められたものであり、具体的には、（４９）〜（５６）式を用いて計算を実行する。

Ｓ１６０では、運動行列の試験分布ｑ（Ｍ）を算出する、変分ベイズＭステップとしての処理を実行して、Ｓ１７０に進む。

このＳ１６０では、具体的には、Ｓ１４０にて得られたＷ，Ｑ、及びＳ１５０にて設定された最新の形状行列の分布パラメータＳ，Σ_Sに基づいて、運動行列の行ベクトルθ_d毎に、次の（５７）（５８）式を用いて、θ_d，Ｖ_dを求めることにより、運動行列の分布パラメータΘ，ＶΘを更新する。但し、期待値計算を記述する場合、（５９）式に示す記述ルールを用いるものとする。

但し、（５７）（５８）式は、（３５）式を、特徴点の三次元座標ｓ⁽ⁿ⁾及び分散Σ_s ⁽ⁿ⁾を用いて書き直した（６０）を解くことにより求められたものであり、具体的には、（６１）（６２）式を用いて計算を実行する。

Ｓ１７０では、Ｓ１５０，Ｓ１６０で算出された形状行列の分布パラメータＳ，Σ_s （即ち、試験分布ｑ（Ｓ））及び運動行列の分布パラメータΘ，ＶΘ（即ち、試験分布ｑ（Ｍ））に基づき、上述の（２６）式を用いて自由エネルギーＦ［ｑ（Ｍ，Ｓ）］を算出してＳ１８０に進む。

但し、自由エネルギーは、上述の処理によって得られた運動行列Ｍおよび形状行列Ｓに基づいて推定される位置と、観測により得られた位置（観測行列Ｗ）とが近いほど大きな値となる。

Ｓ１８０では、Ｓ１７０にて算出された自由エネルギーが、予め設定されたエネルギー閾値より大きいことを終了条件として、その終了条件を充足しているか否かを判断し、充足していなければ、Ｓ１５０に戻って、Ｓ１５０〜Ｓ１７０の処理を繰り返し、一方、終了条件を充足していれば、Ｓ１９０に移行する。

ここでは、エネルギー閾値として、例えば、平均顔モデルＳ_Aを表す形状行列と、上述の処理で得られた運動行列Ｍとを用いて算出した自由エネルギーを用いる。但し、エネルギー閾値は、これに限るものではなく、何らかの固定値であってもよい。

Ｓ１９０では、上述の処理により最終的に得られた形状行列Ｓを個人顔モデルＳ_Pとして出力すると共に、顔モデル切替部６３が推定部８に供給する顔モデルＣ_aとして、個人顔モデルＳ_Pを選択するように切替信号を出力して、本処理を終了する。

これにより、以後、推定部８では、映像が入力される時刻ｔ毎に、個人顔モデルＳ_Pを用いたアフィンパラメータＡ_tの算出が行われ、更に、そのアフィンパラメータＡ_tに基づいて、頭部姿勢算出部９では、頭部姿勢θｘ，θｙ，θｚの算出が行われる。

＜効果＞
以上説明したように、画像処理装置１では、頭部姿勢θｘ，θｙ，θｚを推定する処理を開始した直後は、平均顔モデルＳ_Aを用いてアフィンパラメータＡ_tの算出（顔向き推定）を行い、これと並行して、個人顔モデルＳ_Pの学習を実行し、その個人顔モデルＳ_Pの精度が十分に向上した時点で、平均顔モデルＳ_Aから個人顔モデルＳ_Pへの切替を行うようにされている。

従って、画像処理装置１によれば、ドライバが誰であっても、頭部姿勢θｘ，θｙ，θｚの推定を、処理の初期段階から、安定した精度の推定結果を得ることができ、しかも、平均顔モデルＳ_Aから個人顔モデルＳ_Pへの切替後は、高い精度の推定結果を得ることができる。

即ち、個人顔モデルＳ_Pを用いた場合、平均顔モデルＳ_Aを用いる場合より、推定精度が向上する可能性が高いが、学習に時間がかかり、また、計算に失敗すると平均顔モデルＳ_A以上に誤差が大きくなる可能性があるが、画像処理装置１では、このような問題をいずれも解決することができる。

また、画像処理装置１では、平均顔モデルＳ_Aから個人顔モデルＳ_Pへの切替を、自由エネルギーを用いて判断しているため、モデルの切替後は、推定精度を確実に向上させることができる。

即ち、個人顔モデルＳ_Pは、正解データというものが存在しないため、その精度を定量的に評価することが困難である。但し、観測行列Ｗを運動行列Ｍと形状行列Ｓに因子分解する際に、自由エネルギーが最大となるように因子分解を行っている。つまり、これは、自由エネルギーとは、運動行列Ｍおよび形状行列Ｓから算出される予測値が、観測値（観測行列Ｗ）に当てはまっている度合いを定量的に表しているものであると考えることができる。従って、この自由エネルギーを、顔モデルの精度を評価する基準として使用することができるのである。

また、画像処理装置１では、観測行列Ｗを因子分解する際に、観測行列Ｗを構成する各観測値ｘ_tの信頼性（特徴点の位置を表す確率分布（正規分布）の分散）を表す信頼度行列Ｑを用い、これを形状行列の試験分布ｑ（Ｓ）の分散Σ_s に反映させることにより（（４３）〜（５３）式参照）、信頼性の高い特徴点ほど重視した因子分解を行うようにされている。

従って、画像処理装置１によれば、特徴点の隠れ，誤追跡，フレームアウト等に基づく信頼性の低い観測値ｘ_t が観測行列Ｗ中に存在していたとしても、これを除去するといった余計な手順を必要とすることなく、因子分解を確実に実行することができ、しかも信頼性の高いデータほど重視されるため、生成される形状行列や運動行列の精度を向上させることができる。

また、画像処理装置１によれば、観測行列Ｗを因子分解する際に、変分ベイズ法を適用して、運動行列Ｍの試験分布ｑ（Ｍ），形状行列Ｓの試験分布ｑ（Ｓ）を求めることにより、確率的に因子分解計算を実行するようにされているため、観測値ｘ_tに重畳されているノイズの影響を、大幅に軽減することができる。

［第２実施形態］
次に、第２実施形態について説明する。
なお、本実施形態では、個人顔モデル学習部６２での処理内容の一部が、第１実施形態のものとは異なるだけであるため、この相違する部分を中心に説明する。

＜個人顔モデル学習部での処理＞
図８は、個人顔モデル学習部６２が実行する処理の内容を示すフローチャートである。
第１実施形態の場合と同様に、本処理は、初期化指令の入力により起動され、Ｓ２１０〜Ｓ２５０は、第１実施形態のＳ１１０〜Ｓ１５０と同様であるため説明を省略する。

但し、Ｓ２３０にて使用される規定個数は、第１実施形態の場合より少ない数（例えば、数分の一から数十分の一程度）でよい。
そして、Ｓ２６０では、運動行列の試験分布ｑ（Ｍ）を算出する、変分ベイズＭステップとしての処理を実行して、Ｓ２７０に進む。

このＳ２６０では、具体的には、Ｓ２４０にて得られたＷ，Ｑ、Ｓ２５０にて設定された最新の形状行列の分布パラメータＳ，Σ_S、前サイクルで算出された統計量（後述の（６６）（６７）式参照）に基づいて、運動行列の行ベクトルθ_d毎に、次の（６３）〜（６７）式を用いて、θ_d，Ｖ_dを求めることにより、運動行列の分布パラメータΘ，ＶΘを更新する。

Ｓ２７０では、前サイクルで算出された統計量（後述の（７０）式参照）に基づき、（６８）〜（７０）式を用いて、形状行列Ｓを更新する。

以下、Ｓ２８０〜Ｓ３００は、第１実施形態のＳ１７０〜Ｓ１９０と同様であるため説明を省略する。但し、Ｓ２９０では、終了条件を充足していない場合に、Ｓ２３０に戻る。

＜効果＞
以上説明したように、本実施形態では、観測行列Ｗの次数（即ち、規定個数Ｔ）を少ない数に制限して、Ｅステップ，Ｍステップを１回ずつ実行して運動行列の試験分布ｐ（Ｍ）及び形状行列の試験分布ｐ（Ｓ）を求める処理を繰り返し、その結果を過去のデータほど比率が低下するような係数を乗じて足し合わせていく、いわゆるオンライン化した変分ベイズ推定法を用いて、因子分解を行っている。

従って、本実施形態によれば、使用する画像系列の数を増加させても、それによって増加する処理量は、Ｅステップ，Ｍステップを１回ずつ実行する際の処理量を単位として、加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。

その結果、自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも問題なく適用することができる。
［他の実施形態］
以上本発明の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、様々な態様での実施が可能である。

例えば、上記実施形態では、頭部姿勢を推定するためにアフィンパラメータを用いているが、これは頭部姿勢の推定を行うことができれば、例えば、アフィンパラメータを用いずに、拡張カルマンフィルタを利用して直接頭部推定を行う手法等を採用してもよい。

また、因子分解を終了させる際の終了条件として、平均顔モデルＳ_Aから求めた自由エネルギーをエネルギー閾値として、学習中の個人顔モデルＳ_Pから求めた自由エネルギーがそのエネルギー閾値より大きいことを用いているが、エネルギー閾値は予め設定された固定値であってもよい。

更に、自由エネルギーを用いるのではなく、処理の繰り返し回数が、予め設定された回数閾値より大きいことを終了条件としてもよい。この場合、処理を終了させるか否かの判断のために、自由エネルギーを算出する必要がなくなるため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。

上記実施形態では、ベイズ推定法および変分ベイズ推定法により因子分解を行ったが、因子分解の際に、観測行列Ｗの各要素に信頼度を反映させることができれば、どのような方法を用いてもよく、例えば、特異値分解により因子分解を行ってもよい。

発明が適用された画像処理装置の全体構成を示すブロック図。 特徴点，テンプレートについての説明図。 画像処理装置が使用する座標系を示す説明図。 推定部が扱う状態空間モデルを表すグラフ。 時系列ベイズ推定の動作概要を示す説明図。 推定部の構成を示すブロック図。 パーティクルフィルタの動作概要を示す説明図。 第１実施形態における個人顔モデル学習部の処理内容を示すフローチャート。 第２実施形態における個人顔モデル学習部の処理内容を示すフローチャート。

符号の説明

１…画像処理装置 ２…カメラ ３…画像取込装置 ４…初期特徴点抽出部 ５…テンプレート生成・記憶部 ６…顔モデル供給部 ８…推定部 ９…頭部姿勢算出部 １０…特徴点追跡部 ３０…アフィンパラメータ算出部 ４０…予測値保持部 ４１…予測分布パラメータ算出部 ６１…平均顔モデル記憶部 ６２…個人顔モデル学習部 ６３…顔モデル切替部 ＢＫ⁽ⁿ⁾…追跡器

Claims (7)

1. 対象物体を撮影した画像系列を入力し、前記対象物体について予め設定された複数の特徴点の位置を画像面上で追跡し、該特徴点の位置を時系列的に並べてなる観測行列を生成する観測行列生成手段と、
   該観測行列生成手段にて生成された観測行列を因子分解して、前記対象物体の三次元形状を表す形状行列および前記対象物体の運動を表す運動行列を生成する分解手段と、
   を備えた三次元形状復元装置において、
   前記特徴点毎に、前記特徴点の位置の信頼性を表す位置信頼性評価値を生成する信頼性評価手段を備え、
   前記観測行列をＷ，前記運動行列をＭ，前記形状行列をＳとし、行列Ｘの各要素の値を確率的に示した分布をｑ（Ｘ）で表記するものとして、
   前記分解手段は、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、前記尤度Ｐ（Ｗ｜Ｍ，Ｓ）を最大にする前記運動行列の試験分布ｑ（Ｍ）および前記形状行列の試験分布ｑ（Ｓ）を、前記形状行列の試験分布ｑ（Ｓ）の分散に前記信頼性評価値を反映させた上で求めることで、確率的に因子分解計算を実行することにより、前記位置信頼性評価値の高い特徴点を重視した因子分解を実行することを特徴とする三次元形状復元装置。
2. 前記分解手段は、
   前記運動行列の試験分布ｑ（Ｍ）を固定して前記形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、前記形状行列の試験分布ｑ（Ｓ）を固定して前記運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、
   予め設定された終了条件が満たされている場合に、前記推定手段を終了させる終了手段と、
   からなり、最終的に得られた前記形状行列の試験分布ｑ（Ｓ）が前記形状行列の分布であるものとして、該形状行列の分布から形状行列Ｓを得ることを特徴とする請求項１に記載の三次元形状復元装置。
3. 前記分解手段は、
   前記観測行列生成手段にて観測行列が生成される毎に、前記運動行列の試験分布ｑ（Ｍ）を固定して前記形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、前記形状行列の試験分布ｑ（Ｓ）を固定して前記運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、
   前記推定手段にて推定された前記形状行列の試験分布ｑ（Ｓ）の系列に基づいて、出力とする前記形状行列の分布を更新する更新手段と、
   予め設定された終了条件が満たされている場合に、前記推定手段及び更新手段を終了させる終了手段と、
   からなり、前記更新手段によって最終的に得られた前記形状行列の分布から形状行列Ｓを得ることを特徴とする請求項１に記載の三次元形状復元装置。
4. 前記推定手段にて試験分布ｑ（Ｓ），ｑ（Ｍ）の推定結果が得られる毎に、試験事後分布ｑ（Ｓ，Ｍ）＝ｑ（Ｓ）ｑ（Ｍ）に関する自由エネルギーを算出するエネルギー算出手段を備え、
   前記終了手段は、前記エネルギー算出手段での算出結果が、予め設定されたエネルギー閾値より大きいことを前記終了条件とすることを特徴とする請求項２又は請求項３に記載の三次元形状復元装置。
5. 前記終了手段は、前記推定手段にて前記形状行列を算出した回数が、予め設定された回数閾値より大きいことを前記終了条件とすることを特徴とする請求項２又は請求項３に記載の三次元形状復元装置。
6. 前記信頼性評価手段は、前記位置信頼性評価値として前記特徴点の位置を確率的に表す分布の分散値を用いることを特徴とする請求項１乃至請求項５のいずれかに記載の三次元形状復元装置。
7. 自動車に搭載して使用することを特徴とする請求項１乃至請求項６のいずれかに記載の三次元形状復元装置。