(57)【要約】【課題】処理開始直後から安定した推定結果が得られ、且つ、推定精度を向上させること可能な推定装置を提供する。【解決手段】顔モデル供給部６は、初期化指令が入力され個人顔モデル学習部６２での処理が起動すると、平均顔モデル記憶部６１に記憶された平均顔モデルＳAが顔モデル切替部６３を介して頭部姿勢を求めるためのアフィンパラメータＡtを推定する推定部８に供給する。また、個人顔モデル学習部６２では、推定部８にて画像上の特徴点を追跡した結果（観測値）を取得して、個人顔モデルＳPの学習を実行し、その個人顔モデルＳPの自由エネルギーが、平均顔モデルＳAの自由エネルギーより大きくなった時点で、個人顔モデルの学習を終了し、推定部８に供給する顔モデルＣaを、平均顔モデルＳAから個人顔モデルＳPへ切り替える。【選択図】図１
2008262288.tif 000002
【特許請求の範囲】
【請求項１】
形状に個体差のある対象物体について、該対象物体の三次元形状を表す形状モデルを設定する形状モデル設定手段と、
前記対象物体の一つを特定対象物体として、該特定対象物体を撮影した画像系列を入力し、前記対象物体上に存在する複数の特徴点を抽出する特徴点抽出手段と、
前記形状モデル設定手段にて設定された形状モデル、及び前記特徴点抽出手段にて抽出された特徴点に基づいて、前記特定対象物体の向きを推定する推定手段と、
を備え、
前記形状モデル生成手段は、
前記対象物体の平均的な形状を表す平均モデルを記憶する記憶手段と、
前記特徴点抽出手段での抽出結果に基づいて、前記特定対象物体の形状を表す個別モデルを学習する学習手段と、
当該装置の起動時には前記記憶手段に記憶された平均モデルを前記形状モデルとして設定し、前記学習手段にて学習された個別モデルが前記平均モデルより信頼性が高くなったことを表す切替条件が成立すると、前記平均モデルの代わりに前記個別モデルを前記形状モデルとして設定するモデル切替手段と、
からなることを特徴とする推定装置。
【請求項２】
前記学習手段は、
前記特徴点抽出手段にて抽出された特徴点の位置を時系列的に並べてなる観測行列を生成する観測行列生成手段と、
前記観測行列生成手段にて生成された観測行列を因子分解して、前記特定対象物体の三次元形状を表す形状行列および前記対象物体の運動を表す運動行列を生成する分解手段と、
からなり、前記形状行列を前記個別モデルとすることを特徴とする請求項１に記載の推定装置。
【請求項３】
前記観測行列をＷ，前記運動行列をＭ，前記形状行列をＳとし、行列Ｘの各要素の値を確率的に示した分布をｑ（Ｘ）で表記するものとして、
前記分解手段は、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、前記尤度Ｐ（Ｗ｜Ｍ，Ｓ）を最大にする前記運動行列の試験分布ｑ（Ｍ）および前記形状行列の試験分布ｑ（Ｓ）を求めることにより、確率的に因子分解計算を実行することを特徴とする請求項２に記載の推定装置。
【請求項４】
前記分解手段は、
前記運動行列の試験分布ｑ（Ｍ）を固定して前記形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、前記形状行列の試験分布ｑ（Ｓ）を固定して前記運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、
予め設定された終了条件が成立している場合に、前記推定手段を終了させる終了手段と、
からなり、最終的に得られた前記形状行列の試験分布ｑ（Ｓ）が前記形状行列の分布であるものとして、該形状行列の分布から形状行列Ｓを取得し、
前記モデル切替手段は、前記終了条件を前記切替条件として使用することを特徴とする請求項３に記載の推定装置。
【請求項５】
前記分解手段は、
前記観測行列生成手段にて観測行列が生成される毎に、前記運動行列の試験分布ｑ（Ｍ）を固定して前記形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、前記形状行列の試験分布ｑ（Ｓ）を固定して前記運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、
前記推定手段にて推定された前記形状行列の試験分布ｑ（Ｓ）の系列に基づいて、出力とする前記形状行列の分布を更新する更新手段と、
予め設定された終了条件が満たされている場合に、前記推定手段及び更新手段を終了させる終了手段と、
からなり、前記更新手段によって最終的に得られた前記形状行列の分布から形状行列Ｓを取得し、
前記モデル切替手段は、前記終了条件を前記切替条件として使用することを特徴とする請求項３に記載の推定装置。
【請求項６】
前記推定手段にて試験分布ｑ（Ｓ），ｑ（Ｍ）の推定結果が得られる毎に、試験事後分布ｑ（Ｓ，Ｍ）＝ｑ（Ｓ）ｑ（Ｍ）に関する自由エネルギーを算出するエネルギー算出手段を備え、
前記終了手段は、前記エネルギー算出手段での算出結果が、予め設定されたエネルギー閾値より大きいことを前記終了条件とすることを特徴とする請求項４又は請求項５に記載の推定装置。
【請求項７】
前記終了手段は、前記推定手段にて前記形状行列の試験分布ｑ（Ｓ）を算出した回数が、予め設定された回数閾値より大きいことを前記終了条件とすることを特徴とする請求項４又は請求項５に記載の推定装置。
【請求項８】
前記対象物体とは人間の頭部であることを特徴とする請求項１乃至請求項７のいずれかに記載の推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象物体を撮影した画像（二次元形状）から、その対象物体の三次元空間内での向きを推定する推定装置に関する。
【背景技術】
【０００２】
従来より、予め記憶されている人間の頭部（以下、対象物体ともいう）の三次元形状（以下、形状モデルともいう）と、対象物体を撮影した画像から抽出された複数の特徴点系列とに基づいて、撮影した対象物体の三次元空間内での向きを推定する装置が知られている。
【０００３】
ところで、このような装置では、形状モデルを生成する際にモデルとなった対象物体と、撮影された対象物体とが同一である場合には、その向きを精度良く推定することが可能であるが、その両者が異なる場合には、推定精度が低下してしまうという問題があった。
【０００４】
これに対して、向きの推定を開始する前に、撮影対象となる対象物体の形状モデル（個別モデル）を個別に学習し、その学習した個別モデルを用いて対象物体の向きを推定する方法（例えば、特許文献１参照）や、複数の対象物体についての平均的な形状を表す形状モデル（平均モデル）を予め作成しておき、その平均モデルを用いて対象物体の向きを推定する方法等が知られている。
【特許文献１】特開２００３−１４１５５２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、個別モデルを用いた手法では、どの対象物体に対しても推定精度が向上する可能性が高いが、学習に時間を要するため、対象物体の向きの推定を直ちに開始することができないという問題や、学習に失敗すると、却って推定精度が大幅に低下してしまうという問題があった。
【０００６】
また、平均モデルを用いた手法では、モデルの学習を必要としないため、対象物体の向きの推定を直ちに開始することができるが、平均モデルが示す三次元形状と対象物体の三次元形状とが一致することはないため、推定結果にはある程度の誤差が常に含まれることになり、推定精度の向上を望むことができないという問題があった。
【０００７】
本発明は、上記問題点を解決するために、処理開始直後から安定した推定結果が得られ、且つ、推定精度を向上させること可能な推定装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するためになされた本発明の推定装置では、形状モデル設定手段が、形状に個体差のある対象物体について、その対象物体の三次元形状を表す形状モデルを設定し、特徴点抽出手段が、対象物体の一つを特定対象物体として、その特定対象物体を撮影した画像系列を入力し、特定対象物体上の複数の特徴点を抽出する。
【０００９】
即ち、対象物体の形状に個体差があるということは、対象物体毎に、複数の特徴点間の位置関係に個体差があることになる。
そして、推定手段が、形状モデル設定手段にて設定された形状モデル、及び特徴点抽出手段にて抽出された特徴点に基づいて、特定対象物体の向きを推定する。
【００１０】
但し、形状モデル生成手段では、記憶手段に、対象物体の平均的な形状を表す平均モデルが記憶されており、学習手段が、特徴点抽出手段での抽出結果に基づいて、特定対象物体の形状を表す個別モデルを学習し、モデル切替手段が、当該装置の起動時には記憶手段に記憶された平均モデルを前記形状モデルとして設定し、学習手段にて学習された個別モデルが平均モデルより信頼性が高くなったことを表す切替条件が成立すると、平均モデルの代わりに個別モデルを形状モデルとして設定する。
【００１１】
つまり、本発明の推定装置では、当該装置の起動時には、平均モデルが形状モデルに設定され、対象物体の向きの推定を、起動直後から安定した精度にて実行する。これと並行して、学習手段が個別モデルの学習を実行し、個別モデルが平均モデルより信頼性が高くなった（切替条件が成立した）後は、個別モデルが形状モデルに設定され、平均モデルを用いる場合よりも高い推定精度にて、対象物体の向きを推定する。
【００１２】
従って、本発明の推定装置によれば、処理開始直後には平均モデルを用いることで安定した推定結果を得ることができ、且つ、個別モデルに切り替わった後は、高い推定結果を得ることができる。
【００１３】
なお、学習手段は、例えば、請求項２に記載のように、観測行列生成手段が、特徴点抽出手段にて抽出された特徴点の位置を時系列的に並べてなる観測行列を生成し、分解手段が、観測行列生成手段にて生成された観測行列を因子分解して、特定対象物体の三次元形状を表す形状行列および特定対象物体の運動を表す運動行列を生成し、その生成された形状行列を個別モデルとするように構成されていてもよい。
【００１４】
この場合、請求項３に記載のように、観測行列をＷ，運動行列をＭ，形状行列をＳとし、行列Ｘの各要素の値を確率的に示した分布をｑ（Ｘ）で表記するものとして、分解手段は、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、尤度Ｐ（Ｗ｜Ｍ，Ｓ）を最大にする運動行列の試験分布ｑ（Ｍ）および前記形状行列の試験分布ｑ（Ｓ）を求めることにより、確率的に因子分解計算を実行することが望ましい。
【００１５】
このように構成された本発明の推定装置によれば、運動行列Ｍや形状行列Ｓがベイズ推定によって確率的に求められることになるため、画像面上で追跡した特徴点の位置に重畳されているノイズの影響を、大幅に軽減することができる。
【００１６】
そして、このようにベイズ推定を用いた因子分解を行う場合、分解手段は、請求項４に記載のように、推定手段が、運動行列の試験分布ｑ（Ｍ）を固定して形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、形状行列の試験分布ｑ（Ｓ）を固定して運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行し、終了手段が、予め設定された終了条件が満たされている場合に、推定手段を終了させ、最終的に得られた形状行列の試験分布ｑ（Ｓ）が形状行列の分布であるものとして、その形状行列の分布から形状行列を取得するように構成されていてもよい。そして、この場合、モデル切替手段は、終了手段での終了条件を切替条件として使用すればよい。
【００１７】
なお、分解手段は、変分ベイズ推定法を実現するものであり、第１のステップがＥステップ、第２のステップがＭステップに相当する。
また、ベイズ推定を用いた因子分解を行う場合、分解手段は、請求項５に記載のように、推定手段が、観測行列生成手段にて観測行列が生成される毎に、運動行列の試験分布ｑ（Ｍ）を固定して形状行列の試験分布ｑ（Ｓ）を推定する第１のステップ、及び、形状行列の試験分布ｑ（Ｓ）を固定して運動行列の試験分布ｑ（Ｍ）を推定する第２のステップを、互いの推定結果を用いて実行し、更新手段が、推定手段にて推定された形状行列の試験分布ｑ（Ｓ）に基づいて、出力とする形状行列の分布を更新し、終了手段が、予め設定された終了条件が満たされている場合に、推定手段及び更新手段を終了させ、更新手段によって最終的に得られた形状行列の分布から形状行列を得るように構成されていてもよい。そして、この場合、モデル切替手段は、終了手段での終了条件を切替条件として使用すればよい。
【００１８】
つまり、推定手段では、計算に用いる画像系列を制限して、第１のステップ（Ｅステップ），第２のステップ（Ｍステップ）を１回ずつ実行し、これを繰り返し実行すると共に、推定手段から計算結果が得られる毎に、その計算結果に基づいて形状行列の分布を更新する。これは、変分ベイズ推定法を逐次計算にて実現（いわゆるオンライン化）するものである。
【００１９】
このように構成された本発明の推定装置によれば、形状行列の学習に使用する画像系列を複数に分割し、その分割した画像系列毎に因子分解を実行するため、画像系列を増やしても、処理量は推定手段での１回の処理を単位として加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【００２０】
また、上述したように、分解手段を、変分ベイズ推定法にて因子分解を実行するように構成した場合、請求項６に記載のように、エネルギー算出手段が、推定手段にて試験分布ｑ（Ｓ），ｑ（Ｍ）の推定結果が得られる毎に、試験事後分布ｑ（Ｓ，Ｍ）＝ｑ（Ｓ）ｑ（Ｍ）に関する自由エネルギーを算出し、終了手段は、エネルギー算出手段での算出結果が、予め設定されたエネルギー閾値より大きいことを終了条件とするように構成されていてもよい。
【００２１】
もともと、変分ベイズ推定法では、自由エネルギーが最大となるように運動行列や形状行列を推定しているため、この自由エネルギーを終了条件として用いることにより、所望の精度を有した形状行列を確実に取得することができる。
【００２２】
また、終了手段は、請求項７に記載のように、推定手段にて形状行列の試験分布ｑ（Ｓ）を算出した回数が、予め設定された回数閾値より大きいことを終了条件とするように構成されていてもよい。
【００２３】
この場合、自由エネルギーを算出する必要がないため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。
ところで、対象物体としては、例えば、請求項８に記載のように人間の頭部が考えられる。特に、その人間が車両のドライバーである場合、当該推定装置での推定結果を、車両の運転をアシストする制御や、車両の走行安全性を向上させる制御等に好適に用いることができる。
【発明を実施するための最良の形態】
【００２４】
以下に本発明の実施形態を図面と共に説明する。
［第１実施形態］
図１は、本発明を適用した推定装置としての画像処理装置１の全体構成を示すブロック図である。なお、画像処理装置１は、車両に搭載され、ドライバの頭部を撮影した画像から、ドライバのわき見や居眠り等を検出するための情報である頭部姿勢を時系列的に求める処理を実行する装置である。
【００２５】
＜全体構成＞
画像処理装置１は、ドライバシートの前面上部又は下部（例えば、メータ内）に設置され、図１に示すように、ドライバの頭部を含む画像を取得するカメラ２と、カメラ２が取得した画像を取り込む画像取込装置３と、初期化指令が入力された場合に、画像取込装置３を介して取り込んだ入力画像ｚt （ｔは時刻を表す）から、ドライバの顔の特徴を表す複数（Ｎf個）の特徴点の位置（画像面上での二次元座標：以下、特徴点位置という）ｐ(n) （ｎ＝１，２，…Ｎf）を抽出する初期特徴点抽出部４と、初期特徴点抽出部４が抽出した特徴点位置ｐ(n) のそれぞれについて、その特徴点位置ｐ(n) を中心とする予め設定された大きさの画像をテンプレートｔｐ(n) として切り出して記憶するテンプレート生成・記憶部５と、特徴点の三次元位置（ひいては、頭部の三次元形状）を規定する顔モデル（形状行列）Ｃa を供給する顔モデル供給部６と、入力画像ｚt ，テンプレートｔｐ(n) ，顔モデルＣa に基づいて、時系列ベイズ推定を行うことにより、頭部姿勢を表すアフィンパラメータＡt を推定すると共に、その推定の過程で算出される各特徴点の位置の分布を規定するための特徴点推定分布パラメータ（平均値ｘt ，分散Ｖt ）を、顔モデル供給部６に供給する推定部８と、推定部８が推定したアフィンパラメータＡt に従って、頭部姿勢θｘ，θｙ，θｚを算出する頭部姿勢算出部９から構成される。
【００２６】
このうち、初期特徴点抽出部４では、図２（ａ）に示すように、目尻，目頭，口元、鼻（鼻の穴，鼻の中心，鼻の左右端）等を特徴点として抽出する。但し、本実施形態では、左右目尻、左右目頭，左右口元，鼻の中心の７（＝Ｎf）個の点を抽出するものとする。また、初期化指令は、ドライバが正面を向き、全ての特徴点が映っている（特徴点の隠れのない）画像をカメラ２で撮影した時に、ドライバの操作によって入力されるものとする。
【００２７】
テンプレート生成・記憶部５が生成するテンプレートｔｐ(n) は、具体的には、図２（ａ）中の点線白枠にて示すようなものとなる。但し、図では、左口元、鼻の中心、左目尻についてのみ示す。
【００２８】
また、顔モデル供給部６から供給される顔モデルＣaは、特徴点（図２（ｂ）に示す線図の各頂点）の位置を三次元的に規定するものであり、図３に示すように、画像面における水平方向の座標軸をｕ軸、垂直方向の座標軸をｖ軸、３次元空間における水平方向の座標軸をＸ軸、垂直方向の座標軸をＹ軸、光軸方向（奥行き方向）の座標軸をＺ軸として（１）（２）式で表される。但し、ｓ(n) は、初期特徴点抽出部４で抽出されたｎ番目の特徴点（特徴点位置ｐ(n) ＝（ｕ(n) ，ｖ(n) ））の３次元空間内での座標である。
【００２９】
【数１】
2008262288.tif 000003
次に、推定部８で推定されるアフィンパラメータＡt は、実世界（三次元のＸＹＺ空間）における座標変換を、射影面である画像面（二次元のｕｖ平面）上での座標変換として表す際に用いる変換行列であり、（３）（４）式に示す関係を有する。
【００３０】
但し、ある特徴点の三次元座標を（Ｘ，Ｙ，Ｚ）T ，座標変換後の三次元座標を（Ｘ’，Ｙ’，Ｚ’）T 、これら三次元座標を画像面上に射影した二次元座標、すなわち画像面上での特徴点の座標を（ｕ，ｖ）T ，（ｕ’，ｖ’）T として、Ｘ＝ｕ，Ｙ＝ｖ（すなわちＸ’＝ｕ’，Ｙ’＝ｖ’）であると仮定する。
【００３１】
【数２】
2008262288.tif 000004
そして、頭部姿勢算出部９では、頭部姿勢情報として、頭部の３軸角度θｘ，θｙ，θｚを（５）〜（７）式を用いて算出する。
【００３２】
【数３】
2008262288.tif 000005
即ち、三次元座標をＸ軸を中心に角度θｘだけ回転させる変換行列をＲｘ，Ｙ軸を中心に角度θｙだけ回転させる変換行列をＲｙ，Ｚ軸を中心に角度θｚだけ回転させる変換行列をＲｚとし、これら変換行列Ｒｔ（＝ＲｘＲｙＲｚ）を作用させることで三次元座標が（Ｘ，Ｙ，Ｚ）T から（Ｘ’，Ｙ’，Ｚ’）T に移動するものとすると、これらの関係は、（８）（９）式で表される。
【００３３】
【数４】
2008262288.tif 000006
また、上記回転に加えて、三次元座標をＸ軸にｔ1 ，Ｙ軸にｔ2 ，Ｚ軸にｔ3 だけ平行移動させる変換を考慮した場合、その変換式は（１０）式で表される。
【００３４】
【数５】
2008262288.tif 000007
そして、（３）式と（１０）式を比較することにより、アフィンパラメータａ1 〜ａ8 と（１０）式に示す変換行列の各要素ｒ1 〜ｒ6 ，ｔ1 ，ｔ2 との対応関係が得られ、また、（９）式を整理することで（１１）〜（１８）に示す関係式が得られる。
【００３５】
【数６】
2008262288.tif 000008
つまり、これら（１１）〜（１８）式から、（５）〜（７）式が導かれるのである。
【００３６】
そして、時刻ｔにおける頭部姿勢を表すアフィンパラメータＡt が得られれば、そのアフィンパラメータＡt を顔モデルＣa に作用させることで、時刻ｔにおける画像に映っていない隠れた特徴点の位置も含めて、時刻ｔにおける全ての特徴点の位置を予測することができる。
【００３７】
また、このようなアフィンパラメータＡt には、カメラパラメータ（カメラの焦点距離、設置場所、向き等）及びドライバの頭部姿勢の変化を全て反映させることができる。
＜推定部の構成＞
ここで、図４は、推定部８が扱う状態空間モデルを表すグラフである。
【００３８】
図４に示すように、推定部８は、時刻ｔにおける顔モデルの姿勢を表すアフィンパラメータａt を隠れ状態として求める上位層と、時刻ｔにおける特徴点群の位置ｘt＝（ｘt(1)，ｘt(2)， …，ｘt(Nf)）T を隠れ状態として求める特徴点毎に設けられたＮf個の下位層（図では一つだけ示す）とを備えており、時刻ｔまでに入力された入力画像系列ｚ1:t から、これら隠れ状態Ａt ，ｘtを、時系列ベイズ推定により推定するものである。
【００３９】
なお、時系列ベイズ推定は、図５に示すように、状態変数をすべて確率分布として扱い、隠れ状態に関する時刻ｔ−１における推定結果（推定分布）から時刻ｔにおける予測分布を求めると共に、時刻ｔにおける観測データから、検出すべき隠れ状態であるもっともらしさ（尤度）を求め、これら予測分布と尤度を考慮して時刻ｔにおける推定結果（推定分布）を得るという手順を繰り返すことで、逐次的に隠れ状態を推定するものである。
【００４０】
つまり、推定部８では、入力画像系列（観測データ）ｚ1:t に基づいて、アフィンパラメータ（頭部姿勢の隠れ状態）Ａt を推定した事後確率分布（推定分布）ｐ（Ａt ｜ｚ1:t ）を求めることになり、これを数式で記述すると（１９）（２０）式で表される。
【００４１】
【数７】
2008262288.tif 000009
ここで、ｐ（Ａt ｜ｚ1:t-1 ）はアフィンパラメータＡt の事前確率分布（予測分布）、ｐ（ｘt｜Ａt ，ｚ1:t-1 ）は特徴点群の位置ｘtの事前確率分布（予測分布）、ｐ（ｘt ｜Ａt ，ｚ1:t-1 ），ｐ（ｚt ｜ｘt）は尤度を表す。
【００４２】
そして、（２０）式の部分を下位層がパーティクルフィルタを用いて推定し、（１９）式の部分を上位層がカルマンフィルタを用いて推定する。
なお、あるパラメータｆの確率分布がガウス分布（正規分布）に従う場合、その確率分布は、平均をμ，分散をΣとすると、（２１）式で表すことができる。つまりこの場合、パラメータの確率分布の計算は、実際には、その平均μと分散Σを求めれば十分である。
【００４３】
【数８】
2008262288.tif 000010
次に、図６は推定部８の具体的な構成を表すブロック図である。
【００４４】
図６に示すように、推定部８は、特徴点毎に設けられたＮf個の追跡器ＢＫ(n) からなり、それぞれが一つの特徴点の位置ｘ(n) を追跡して、時刻ｔにおけるその特徴点の位置ｘt(n) を推定した確率分布をガウス近似した確率分布を規定する特徴点推定分布パラメータ（平均、標準偏差）を生成する特徴点追跡部１０と、各追跡器ＢＫ(n) にて算出された特徴点推定分布パラメータに基づき、カルマンフィルタを用いてアフィンパラメータＡt 及びアフィンパラメータの予測値Ａt+1  ，その分散Ｖt+1 （以下、予測値等という）を算出するアフィンパラメータ算出部３０と、アフィンパラメータ算出部３０にて算出されるアフィンパラメータの予測値等を保持する予測値保持部４０と、予測値保持部４０にて保持された時刻ｔ−１において算出されたアフィンパラメータの予測値等に基づいて、時刻ｔの特徴点群の位置ｘt（＝（ｘt(1), ｘt(2), …，ｘt(Nf)）T ）を予測した確率分布ｐ（ｘt｜Ａt ）を規定する上位予測分布パラメータ（平均値，分散）を算出して、追跡器ＢＫ(n) のそれぞれに供給する予測分布パラメータ算出部４１とを備えている。
【００４５】
つまり、アフィンパラメータ算出部３０が上述の上位層、特徴点追跡部１０を構成する各追跡器ＢＫ(n) が上述の下位層に相当する。
なお、このような推定部８を実現する具体的な構成は、例えば、特願２００５−３６８１２４号等に詳述されているため、ここではその詳細についての説明を省略し、以下では、特徴点推定分布パラメータの生成に関係する追跡器ＢＫ(n) の概要についてのみ説明する。
【００４６】
＜追跡器の概要＞
まず、特徴点追跡部１０を構成する追跡器ＢＫ(n) は、パーティクルフィルタによって一つの特徴点を追跡するものである。ここでは、そのパーティクルフィルタの動作概要を、図７に示す説明図を参照して説明する。
【００４７】
図７に示すように、パーティクルフィルタでは、推定したい対象の状態（特徴点の位置）の実現値（画像面上の座標）をパーティクルで表し、予測，観測，リサンプルという三つの過程を繰り返し実行することで、時系列的にパーティクルの分布を求める。なお、パーティクルフィルタでは、カルマンフィルタとは異なり、ガウス分布に限らず任意の確率分布を対象とすることができる。
【００４８】
まず、予測過程では、推定したい対象の動きを考慮して、状態空間（ここでは画像面上）での各パーティクルの状態（ここでは画像面上での位置）を遷移させ、推定したい対象の動きに加わるノイズを考慮してパーティクルを確率的に散らすことで，推定したい対象が存在しそうな位置にパーティクル群を配置する。これにより、推定したい対象の状態を予測した確率分布が、パーティクル群の位置と数により離散的かつ近似的に表現されることになる。
【００４９】
次に、観測過程では、各パーティクルの状態が推定したい対象の状態であることの尤度（ここでは特徴点を表すテンプレートとの正規化相関値）を、尤度に応じてパーティクルの重みを算出する。これにより、推定したい対象の状態の確率分布が、重みのあるパーティクルによって表現されることになる。
【００５０】
また、リサンプル過程では、重みの小さいパーティクルを消去すると共に重み大きいパーティクルを重みのない複数のパーティクルに増殖させることにより、推定したい対象の状態の確率分布を、重みのあるパーティクルによる表現から、重みのないパーティクルにより表現に確率的に変換する。
【００５１】
そして、このリサンプル（リサンプル過程で生成）されたパーティクルを用いて上述の予測過程を行うことにより、パーティクルで表された推定したい対象の状態の確率分布が時系列的に求められることになる。
【００５２】
ここで、パーティクルの個数をＮp 、パーティクルの画像面上での座標をｐi ＝（ｕi ，ｖi ）、観測過程にて尤度に応じて算出されるパーティクルの重みをｗi （ｉ＝１，２，…，Ｎp ）とすると、観測によって得られる推定したい対象の状態（即ち、特徴点の位置）の確率分布は、（２２）（２３）式に示す平均値、及び（２４）式（２５）式に示す分散値によって表される。
【００５３】
【数９】
2008262288.tif 000011
つまり、時刻ｔにおけるｎ番目の特徴点の位置を表す確率分布の平均値をｘt(n)＝（ｘ(n) ，ｙ(n) ）、分散をＶt(n)＝（Ｖx(n)，Ｖy(n)）として、ｘt ＝（ｘt(1)，ｘt(2)，…，ｘt(Nf) ）、Ｖt ＝（Ｖt(1)，Ｖt(2)，…，Ｖt(Nf) ）が、特徴点推定分布パラメータとして、各時刻ｔ毎に顔モデル供給部６に供給されることになる。
【００５４】
＜顔モデル供給部＞
次に、本発明の主要部である顔モデル供給部６について詳述する。
顔モデル供給部６は、図１に示すように、頭部を撮影した複数人分の映像から、それぞれ特徴点を抽出し、その抽出した結果に基づいて設定された平均的な特徴点の配置を表す平均顔モデルＳAを記憶する平均顔モデル記憶部６１と、初期化指令が入力されると起動し、推定部８から逐次供給される特徴点推定分布パラメータｘt ，Ｖt に基づいて、カメラ２に被写体（ドライバ）の特徴点の配置を表す個人顔モデルＳPを、学習，生成する個人顔モデル学習部６２と、個人顔モデル学習部６２からの切替信号に従って、平均顔モデル記憶部６１に記憶された平均顔モデルＳA、個人顔モデル学習部６２にて生成された個人顔モデルＳPのいずれかを、顔モデルＣaとして選択して推定部８に供給する顔モデル切替部６３とを備えている。
【００５５】
＜学習の概要＞
なお、個人顔モデル学習部６２では、以下に示す手法により学習を実行する。
即ち、求めるべき三次元形状（即ち、個人顔モデルＳP）を表す形状行列をＳ，特徴点推定分布パラメータｘt の系列に基づいて生成される観測行列をＷ，形状行列Ｓの運動を表す運動行列をＭとすると、これらは、Ｗ＝ＭＳの関係を有する。そして、事後確率Ｐ（Ｍ，Ｓ｜Ｗ）と尤度Ｐ（Ｗ｜Ｍ，Ｓ）との関係を規定するベイズ推定の式に基づき、尤度Ｐ（Ｗ｜Ｍ，Ｓ）が最大となる運動行列の試験分布ｑ（Ｍ），形状行列の試験分布ｑ（Ｓ）を求めることにより、運動行列Ｍ及び形状行列Ｓを推定する。
【００５６】
これは、（２６）式に定義された自由エネルギーＦ［ｑ（Ｍ，Ｓ）］を最大化するＭ，Ｓを求めることと等価であり、具体的には、（２７）（２８）式を解けばよい。
【００５７】
【数１０】
2008262288.tif 000012
但し、この計算は簡単には解けないため、（２９）式が成立すること（即ち、運動行列Ｍと形状行列Ｓに独立な生起確率を仮定すること）、及び、モデル構造（形状行列Ｓ）と運動構造（運動行列Ｍ）の生起確率は、それぞれ正規分布に従うものと仮定し、いわゆる変分ベイズ法を適用して計算を行う。
【００５８】
【数１１】
2008262288.tif 000013
これにより、求める分布も正規分布に従うため、解くべき問題を、分布の形状を表す特徴点推定分布パラメータｘt ，Ｖt の最適化というように簡略化できる。
【００５９】
そして、変分ベイズ法のＥステップでは、運動行列の分布ｑ（Ｍ）を固定した上で、自由エネルギーを最大化する計算を行えばよく、結局、（２７）（２８）式を変形した（３０）（３１）式からなる方程式を、ラグランジュの未定乗数法を用いて解くことで得られた（３２）式の計算を実行すればよい。但し、式中で記号に付与されている＾（ハット）は、更新値を表すものとする（以下同様）。
【００６０】
【数１２】
2008262288.tif 000014
同様に、変分ベイズ法のＭステップでは、形状行列の分布ｑ（Ｓ）を固定した上で、自由エネルギーを最大化する計算を行えばよく、結局、（２７）（２８）式を変形した（３３）（３４）式からなる方程式を、ラグランジュの未定乗数法を用いて解くことで得られた（３５）式の計算を実行すればよい。
【００６１】
【数１３】
2008262288.tif 000015
＜個人顔モデル学習部での処理＞
ここで、個人顔モデル学習部６２が実行する処理を、図８に示すフローチャートに沿って説明する。なお、本処理は、初期化指令の入力により起動される。
【００６２】
但し、ｎ番目の特徴点の三次元座標ｓ(n) 、求めたい３次元形状を表す形状行列Ｓ、形状行列Ｓに要素がオール１の行ベクトルを加えた同次形状行列 〜Ｓを（３６）〜（３８）式で表すものとする。なお、平均顔モデルＳA及び個人顔モデルＳPは、同次形状行列 〜Ｓの形式で表現されているものとする。
【００６３】
【数１４】
2008262288.tif 000016
また、運動のｉ番目の基底ベクトルをｍi として、運動行列Ｍ，同次運動行列 〜Ｍ，Θ、ｄ番目の運動行列の行ベクトルθd を（３９）〜（４２）式で表すものとする。
【００６４】
【数１５】
2008262288.tif 000017
本処理が起動すると、まず、Ｓ１１０では、顔モデル切替部６３が推定部８に供給する顔モデルＣa として、平均顔モデルＳAを選択するように切替信号を出力して、Ｓ１２０に進む。
【００６５】
これにより、推定部８では、映像が入力される時刻ｔ毎に、平均顔モデルＳAを用いたアフィンパラメータＡt 、及び特徴点推定分布パラメータｘt ，Ｖt の算出が行われる。また、頭部姿勢算出部９では、推定部８にて求められたアフィンパラメータＡt に基づいて、頭部姿勢θｘ，θｙ，θｚの算出が行われる。
【００６６】
Ｓ１２０では、運動行列の試験分布ｑ（Ｍ）を初期化してＳ１３０に進む。
なお、形状行列の試験分布ｑ（Ｓ）は、形状行列Ｓとその分散を表すΣS （以下、形状行列の分布パラメータという）によって規定され、運動行列の試験分布ｑ（Ｍ）は、運動行列Θとその分散を表すＶΘ （以下、運動行列の分布パラメータという）によって規定される。そして、運動行列の分布パラメータΘ，ＶΘ は、この初期化によって、全ての要素が、同じ一定値（例えば、０又は１又は０．５）に設定される。
【００６７】
Ｓ１３０では、予め設定された規定個数Ｔの特徴点推定分布パラメータｘt ，Ｖt の取得を完了したか否かを判断し、完了していなければ、完了するまで待機する。
そして規定個数Ｔの特徴点推定分布パラメータｘt ，Ｖt の取得を完了すると、Ｓ１４０に移行し、その取得した特徴点推定分布パラメータｘt ，Ｖt に基づいて、観測行列Ｗおよび信頼度行列Ｑを生成して、Ｓ１５０へ進む。
【００６８】
なお、取得した特徴点推定分布パラメータｘt に基づいて作成され、ｎ（ｎ＝１，２，…，Ｎf）番目の特徴点の位置（画像面上での２次元座標）を時系列（ｔ＝１〜Ｔ）に並べてなる特徴点ベクトルｘ1:T(n)を（４３）式で表すものとして、観測行列Ｗは、（４４）式により表される。即ち、観測行列Ｗは、特徴点ベクトルｘ1:T(n)を列ベクトルして、全ての特徴点について並べることで構成されている。
【００６９】
【数１６】
2008262288.tif 000018
また、信頼度行列Ｑは、特徴点推定分布パラメータＶt に基づいて作成され、ｄ＝２ｔ−１が時刻ｔにおけるｘ成分、ｄ＝２ｔが時刻ｔにおけるｙ成分を表し、ｎ番目の特徴点のｄ番目（ｄ＝１，２，…，２Ｔ）の観測信頼度をσx,d(n)で表すものとして、（４５）式により表される。即ち、信頼度行列Ｑは、観測行列Ｗと同じ次数を有する。
【００７０】
【数１７】
2008262288.tif 000019
Ｓ１５０では、形状行列の試験分布ｑ（Ｓ）を算出する、変分ベイズＥステップとしての処理を実行して、Ｓ１６０に進む。
【００７１】
このＳ１５０では、具体的には、Ｓ１４０にて得られたＷ，Ｑ、及びＳ１２０にて初期設定又は後述するＳ１６０にて更新された最新の運動行列の分布パラメータΘ，ＶΘ に基づいて、特徴点（ｎ＝１〜Ｎf）毎に、次の（４６）（４７）式を用いて、ｓ(n) ，Σs(n)を求めることにより、形状行列の分布パラメータＳ，Σs を更新する。なお、Ｉk はｋ次の単位行列である。
【００７２】
【数１８】
2008262288.tif 000020
但し、（４６）（４７）式は、（３２）式を、特徴点の三次元座標ｓ(n) 及び分散Σs(n)を用いて書き直した（４８）を解くことにより求められたものであり、具体的には、（４９）〜（５６）式を用いて計算を実行する。
【００７３】
【数１９】
2008262288.tif 000021
Ｓ１６０では、運動行列の試験分布ｑ（Ｍ）を算出する、変分ベイズＭステップとしての処理を実行して、Ｓ１７０に進む。
【００７４】
このＳ１６０では、具体的には、Ｓ１４０にて得られたＷ，Ｑ、及びＳ１５０にて設定された最新の形状行列の分布パラメータＳ，ΣS に基づいて、運動行列の行ベクトルθd 毎に、次の（５７）（５８）式を用いて、θd ，Ｖd を求めることにより、運動行列の分布パラメータΘ，ＶΘ を更新する。但し、期待値計算を記述する場合、（５９）式に示す記述ルールを用いるものとする。
【００７５】
【数２０】
2008262288.tif 000022
但し、（５７）（５８）式は、（３５）式を、特徴点の三次元座標ｓ(n) 及び分散Σs(n)を用いて書き直した（６０）を解くことにより求められたものであり、具体的には、（６１）（６２）式を用いて計算を実行する。
【００７６】
【数２１】
2008262288.tif 000023
Ｓ１７０では、Ｓ１５０，Ｓ１６０で算出された形状行列の分布パラメータＳ，Σ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 を求めることにより、運動行列の分布パラメータΘ，ＶΘ を更新する。
【００９４】
【数２２】
2008262288.tif 000024
Ｓ２７０では、前サイクルで算出された統計量（後述の（７０）式参照）に基づき、（６８）〜（７０）式を用いて、形状行列Ｓを更新する。
【００９５】
【数２３】
2008262288.tif 000025
以下、Ｓ２８０〜Ｓ３００は、第１実施形態のＳ１７０〜Ｓ１９０と同様であるため説明を省略する。但し、Ｓ２９０では、終了条件を充足していない場合に、Ｓ２３０に戻る。
【００９６】
＜効果＞
以上説明したように、本実施形態では、観測行列Ｗの次数（即ち、規定個数Ｔ）を少ない数に制限して、Ｅステップ，Ｍステップを１回ずつ実行して運動行列の試験分布ｐ（Ｍ）及び形状行列の試験分布ｐ（Ｓ）を求める処理を繰り返し、その結果を過去のデータほど比率が低下するような係数を乗じて足し合わせていく、いわゆるオンライン化した変分ベイズ推定法を用いて、因子分解を行っている。
【００９７】
従って、本実施形態によれば、使用する画像系列の数を増加させても、それによって増加する処理量は、Ｅステップ，Ｍステップを１回ずつ実行する際の処理量を単位として、加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【００９８】
その結果、自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも問題なく適用することができる。
［他の実施形態］
以上本発明の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、様々な態様での実施が可能である。
【００９９】
例えば、上記実施形態では、頭部姿勢を推定するためにアフィンパラメータを用いているが、これは頭部姿勢の推定を行うことができれば、例えば、アフィンパラメータを用いずに、拡張カルマンフィルタを利用して直接頭部推定を行う手法等を採用してもよい。
また、因子分解を終了させる際の終了条件として、平均顔モデルＳAから求めた自由エネルギーをエネルギー閾値として、学習中の個人顔モデルＳPから求めた自由エネルギーがそのエネルギー閾値より大きいことを用いているが、エネルギー閾値は予め設定された固定値であってもよい。
【０１００】
更に、自由エネルギーを用いるのではなく、処理の繰り返し回数が、予め設定された回数閾値より大きいことを終了条件としてもよい。この場合、処理を終了させるか否かの判断のために、自由エネルギーを算出する必要がなくなるため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。
【０１０１】
上記実施形態では、ベイズ推定法および変分ベイズ推定法により因子分解を行ったが、因子分解の際に、観測行列Ｗの各要素に信頼度を反映させることができれば、どのような方法を用いてもよく、例えば、特異値分解により因子分解を行ってもよい。
【図面の簡単な説明】
【０１０２】
【図１】発明が適用された画像処理装置の全体構成を示すブロック図。
【図２】特徴点，テンプレートについての説明図。
【図３】画像処理装置が使用する座標系を示す説明図。
【図４】推定部が扱う状態空間モデルを表すグラフ。
【図５】時系列ベイズ推定の動作概要を示す説明図。
【図６】推定部の構成を示すブロック図。
【図７】パーティクルフィルタの動作概要を示す説明図。
【図８】第１実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【図９】第２実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【符号の説明】
【０１０３】
１…画像処理装置  ２…カメラ  ３…画像取込装置  ４…初期特徴点抽出部  ５…テンプレート生成・記憶部  ６…顔モデル供給部  ８…推定部  ９…頭部姿勢算出部  １０…特徴点追跡部  ３０…アフィンパラメータ算出部  ４０…予測値保持部  ４１…予測分布パラメータ算出部  ６１…平均顔モデル記憶部  ６２…個人顔モデル学習部  ６３…顔モデル切替部  ＢＫ(n)…追跡器
【図３】
2008262288.tif 000026
【図６】
2008262288.tif 000027
【図８】
2008262288.tif 000028
【図９】
2008262288.tif 000029
【図１】
2008262288.tif 000030
【図２】
2008262288.tif 000031
【図４】
2008262288.tif 000032
【図５】
2008262288.tif 000033
【図７】
2008262288.tif 000034