三次元形状復元装置
【課題】三次元形状復元装置において、観測データ行列（特徴点位置系列）中に複数の物体に由来するデータが混在していても因子分解の実行を可能とする。
【解決手段】複数の物体に由来する特徴点の位置を表す観測データ行列ｘを因子分解するにあたり、物体の形状を表す形状行列をｚ、特徴点のそれぞれが属する物体を指定する指定行列をｈとし、また、物体ｕの運動行列Ｗe|u ，物体ｕに由来する特徴点の測定の正確さτ|u，特徴点がそれぞれの物体ｕ（ｕ＝１，２，…，Ｕ）に属している確率（混合比）ｇをまとめたものをモデルパラメータをθとして、変分ベイズ推定法を適用し、Ｅステップでは、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を求め、Ｍステップでは、第２試験分布ｑ（ｚ，ｈ）を固定して第１試験分布ｑ（θ）を求めることで、確率的に因子分解計算を実行する。
【解決手段】複数の物体に由来する特徴点の位置を表す観測データ行列ｘを因子分解するにあたり、物体の形状を表す形状行列をｚ、特徴点のそれぞれが属する物体を指定する指定行列をｈとし、また、物体ｕの運動行列Ｗe|u ，物体ｕに由来する特徴点の測定の正確さτ|u，特徴点がそれぞれの物体ｕ（ｕ＝１，２，…，Ｕ）に属している確率（混合比）ｇをまとめたものをモデルパラメータをθとして、変分ベイズ推定法を適用し、Ｅステップでは、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を求め、Ｍステップでは、第２試験分布ｑ（ｚ，ｈ）を固定して第１試験分布ｑ（θ）を求めることで、確率的に因子分解計算を実行する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物体を撮影した画像（二次元形状）系列から、その物体の三次元形状を復元する三次元形状復元装置に関する。
【背景技術】
【０００２】
従来より、物体を撮影した画像系列を入力し、その物体について予め設定された複数の特徴点について、その特徴点の位置を画像面上で追跡し、追跡により得られた特徴点位置系列から物体の三次元形状を復元する装置が知られている。
【０００３】
この主の三次元復元装置の一つとして、特徴点位置系列から観測データ行列を生成し、その観測データ行列を、物体の形状（特徴点の三次元位置）を表す形状行列と、物体の運動を表す運動行列とに因子分解する手法を用いたものがある（例えば、特許文献１〜３参照）。
【０００４】
但し、因子分解法は、観測データ行列を構成するデータが全て揃っていないと、有効な結果が得られないため、特徴点の隠れ，誤追跡，フレームアウト等によって、特徴点位置系列（ひいては観測データ行列を構成するデータ）の欠損が頻繁に生じる実問題への適用は、非常に困難であるという問題があった。
【０００５】
これに対して出願人は、データを統計的に扱い、データの確率分布の分散をデータの信頼度として考慮した計算を行うことによって、データの欠損に対処する手法を提案している（特許文献４参照）。
【０００６】
しかし、この手法でも、形状を推定すべき物体が人物の顔である場合、サングラスやスカーフ等の装飾品によって特徴点が隠れた状態が継続すると、検出精度が劣化してしまう。
【０００７】
このため、予め固定的に設定された特徴点（以下「固定特徴点」という）だけでなく、その都度検出される予め規定された特徴（例えば、撮像された物体のエッジやコーナー等）を有する地点（以下「任意特徴点」という）を検出し、この任意特徴点も利用することが検討されている。
【特許文献１】特開２０００−１９４８６３号公報
【特許文献２】特開２００３−１４１５５２号公報
【特許文献３】特開２００４−３４７３３２号公報
【特許文献４】特開２００８−２６１６６２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかし、任意特徴点は、物体上に限らず、物体の背景や物体以外の物体上であっても、その特徴を有していれば抽出されてしまう。
【０００９】
しかし、特許文献１〜４に記載の従来装置では、いずれも、観測データ行列を構成する特徴点が、単一の対象物に由来していることを前提とするものであるため、特徴点が複数の対象物に由来している場合には、形状行列（ひいては物体の形状）を正しく求めることができないという問題があった。
【００１０】
なお、各特徴点が、どの物体に由来するかを識別し、その識別された特徴点を用いて物体毎に作成した観測データ行列を、それぞれ因子分解する手法が考えられる。しかし、特徴点の識別を、別の画像処理等によって行う必要があり、処理量が増大するだけでなく、その識別が困難であり、十分な精度が得られないという問題があった。
【００１１】
また、観測データ行列を因子分解して形状行列を得る手法において、三次元形状の推定精度を向上させるには、計算に用いる画像系列を増やす（即ち、観測データ行列の次元を大きくする）必要がある。
【００１２】
しかし、因子分解の際に必要となる計算量は、通常、観測データ行列の行数×列数の３乗に比例することが知られており、観測データ行列の次元を大きくすることによって爆発的に増大する。
【００１３】
そして、計算量が増大すると、三次元形状復元装置を車両に搭載する場合等のように、使用できる演算器の計算能力に制限がある場合には、適用することが困難になるという問題もあった。
【００１４】
本発明は、上記問題点を解決するために、三次元形状復元装置において、観測データ行列（特徴点位置系列）中に複数の物体に由来するデータが混在していても因子分解の実行を可能とすること、更には、計算時間を低減することを目的とする。
【課題を解決するための手段】
【００１５】
上記目的を達成するためになされた本発明の三次元形状復元装置では、特徴点追跡手段が、画像系列を入力し、入力された各画像から予め設定された特徴を有する複数の特徴点を抽出すると共に、特徴点のそれぞれについて該特徴点の位置を画像面上で追跡する。
【００１６】
また、観測データ行列生成手段は、特徴点抽出追跡手段での抽出，追跡結果に従って、特徴点の位置を表す観測データ行列を生成する。そして、分解手段が、観測データ行列生成手段にて生成された観測データ行列を因子分解して、画像に写された物体の三次元形状を表す形状行列および前記物体の運動を表す運動行列を生成する。
【００１７】
このとき、分解手段は、観測データ行列をｘ、形状行列をｚ、特徴点のそれぞれが属する物体を指定する指定行列をｈ、運動行列，特徴点の測定の正確さ，特徴点がそれぞれの物体に属している確率をまとめたモデルパラメータをθとして、観測データ行列ｘから事後確率Ｐ（θ，ｚ，ｈ｜ｘ）を計算することで、θ，ｚ，ｈを求める。
【００１８】
このように構成された本発明の三次元形状復元装置によれば、観測データ行列ｘが複数の物体に由来する特徴点で構成された混合データ行列であっても、因子分解を実行することができ、指定行列ｈによって分類される特徴点毎、即ち、特徴点が由来する物体毎に、形状行列や運動行列を求めることができる。
【００１９】
なお、本発明において、分解手段は、請求項２に記載のように、事後確率Ｐ（θ，ｚ，ｈ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）が、第１試験分布ｑ（θ）と第２試験分布ｑ（ｈ，ｚ）とに因子化されるものとして、事後確率Ｐ（θ，ｚ，ｈ｜ｘ）と試験分布ｑ（θ，ｈ，ｚ）とのＫＬ距離を最小にする第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）を求めることで、確率的に因子分解計算を実行するように構成されていてもよい。
【００２０】
このように構成された本発明の三次元形状復元装置によれば、形状行列ｚを求めるのと同時に、指定行列ｈもベイズ推定によって確率的に求められるため、因子分解以外の特別な処理を行うことなく、簡単に特徴点を物体毎に識別することができる。
【００２１】
そして、このようにベイズ推定を用いた因子分解を行う場合、請求項３に記載のように、観測データ行列生成手段は、予め設定されたフレーム数分の画像に基づいて観測データ行列ｘを生成し、分解手段は、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、予め設定された終了条件が満たされている場合に、前記推定手段を終了させる終了手段とからなり、最終的に得られた前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得るように構成されていてもよい。
【００２２】
この場合、分解手段は、変分ベイズ推定法を実現するものであり、第１のステップがＥステップ、第２のステップがＭステップに相当する。
【００２３】
但し、変分ベイズ推定法を用いても、計算に用いる画像系列の増加（即ち、観測データ行列の次元の増大）に応じて、計算量（計算時間）が指数関数的に増大することを避けることはできない。
【００２４】
そこで、請求項４に記載のように、観測データ行列生成手段は、１フレーム分の画像毎に前記観測データ行列ｘを生成し、分解手段は、観測データ行列生成手段にて観測データ行列ｘが生成される毎に、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、第２試験分布ｑ（ｚ，ｈ）を固定して第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、予め設定された終了条件が満たされている場合に、推定手段及び更新手段を終了させる終了手段とからなり、推定手段によって逐次的に得られる第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得るように構成されていていることが望ましい。
【００２５】
つまり、推定手段では、計算に用いる画像系列を制限して、第１のステップ（Ｅステップ），第２のステップ（Ｍステップ）を１回ずつ実行し、これを繰り返し実行すると共に、推定手段から計算結果が得られる毎に、その計算結果に基づいて形状行列の分布を更新する。これは、変分ベイズ推定法を逐次計算にて実現（いわゆるオンライン化）するものである。
【００２６】
このように、本発明の三次元形状復元装置によれば、使用する画像系列を複数に分割し、その分割した画像系列毎に因子分解を実行するため、画像系列を増やしても、処理量は推定手段での１回の処理を単位として加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【００２７】
そして、当該三次元形状復元装置は、請求項５に記載のように、自動車に搭載して使用するように構成されていてもよい。
【００２８】
自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも適用することができる。
【発明を実施するための最良の形態】
【００２９】
以下に本発明の実施形態を図面と共に説明する。
【００３０】
［第１実施形態］
図１は、本発明を適用した三次元形状復元装置としての画像処理装置１の全体構成を示すブロック図である。なお、画像処理装置１は、車両に搭載され、ドライバの頭部を撮影した画像から、ドライバのわき見や居眠り等を検出するための情報である頭部姿勢を時系列的に求める処理を実行する装置である。
【００３１】
＜全体構成＞
図１に示すように、画像処理装置１は、ドライバシートの前面上部又は下部（例えば、メータ内）に設置され、ドライバの頭部を含む画像を取得するカメラ２と、カメラ２が取得した画像を取り込む画像取込装置３と、画像取込装置３を介して取り込んだ入力画像ｓf （ｆはフレームを識別する識別子）から、予め設定された特徴を有する複数の特徴点の位置（画像面上での二次元座標）ｘ(p) （ｐ＝１，２，…Ｐ）を抽出すると共に、その抽出された特徴点を追跡して、その特徴点の位置の分布を規定するための特徴点推定分布パラメータ（平均値ｘf ，分散Ｖf ）を求める特徴点追跡部４とを備えている。
【００３２】
また、画像処理装置１は、特徴点の三次元位置（ひいては、頭部の三次元形状）を規定する顔モデル（形状行列）Ｃa を供給する顔モデル供給部５と、特徴点追跡部４から供給される特徴点推定分布パラメータ、顔モデル供給部５から供給される顔モデルＣa に基づき、カルマンフィルタを用いて、ドライバの頭部姿勢を表すアフィンパラメータＡf を推定する推定部６と、推定部６が推定したアフィンパラメータＡf に従って、頭部姿勢θｘ，θｙ，θｚを算出する頭部姿勢算出部７から構成される。
【００３３】
＜特徴点追跡部＞
特徴点追跡部４は、目尻，目頭，口元、鼻（鼻の穴，鼻の中心，鼻の左右端）等（図２（ａ）参照）、ドライバの顔の特徴を表す複数の特徴点（以下「固定特徴点」ともいう）の位置を抽出，追跡するする固定特徴点追跡部４１と、画像中の追跡し易い特徴（例えば、撮像された物体のエッジやコーナ等）を有する部位（以下「任意特徴点」ともいう）の位置を抽出，追跡する任意特徴点追跡部４３とからなる。
【００３４】
そして、固定特徴点追跡部４１は、予め用意されたテンプレートを用いて入力画像ｓf とのテンプレートマッチングによって、固定特徴点の抽出，追跡を行う。なお、テンプレートは、具体的には、図２（ａ）中の点線白枠にて示すようなものとなる。但し、図では、左口元、鼻の中心、左目尻についてのみ示す。
【００３５】
一方、任意特徴点追跡部４３は、ＳＩＦＴ（Scale Invariant Feature Transform） やＫＬＴ(Kanade-Lucas-Tomasi) といった周知の特徴点検出器を用いて任意特徴点の抽出，追跡を行う。
【００３６】
ここで、図３は、任意特徴点追跡部４３により抽出される任意特徴点の例を示す説明図である。但し、（ａ）が入力画像、（ｂ）が固定特徴点追跡部４１によって抽出される固定特徴点、（ｃ）が任意特徴点追跡部４３によって抽出される任意特徴点を示す。
【００３７】
図３に示すように、任意特徴点は、ドライバ頭部からだけでなく、背景（車室）からも抽出されることになる。但し、ドライバがサングラス等を装着して固定特徴点（目尻，目頭等）が隠れたとしても、サングラスの輪郭等から多数の任意特徴点が抽出されることになる。
【００３８】
つまり、特徴点追跡部４では、常時追跡対象となる固定特徴点以外に、任意特徴点が適宜、追跡対象として追加されることになる。
【００３９】
但し、任意特徴点追跡部４は、任意特徴点の中で固定特徴点と一致するものを追跡対象から除外すると共に、後述する顔モデル供給部５からドライバ頭部に位置する特徴点（以下「頭部特徴点」という）の情報が得られた後は、その情報に基づいて、頭部特徴点以外の任意特徴点を追跡対象から除外するように構成されている。
【００４０】
また、固定特徴点及び任意特徴点の追跡には、いずれも、パーティクルフィルタが用いられており、その処理に必要な事後確率ｐ（ｘf｜Ａf）は、推定部６から供給を受けるように構成されている。
【００４１】
＜顔モデル＞
顔モデル供給部５から供給される顔モデルＣa は、特徴点追跡部４にて追跡される特徴点のうち、ドライバの顔上に位置するものの位置を三次元的に規定するものであり、図４に示すように、画像面における水平方向の座標軸をｕ軸、垂直方向の座標軸をｖ軸、三次元空間における水平方向の座標軸をＸ軸、垂直方向の座標軸をＹ軸、光軸方向（奥行き方向）の座標軸をＺ軸として（１）（２）式で表される。但し、ｚ(p) は、特徴点追跡部４で抽出されたｐ番目の特徴点（特徴点位置ｘ(p) ＝（ｕ(p)，ｖ(p)）の３次元空間内での座標である。
【００４２】
【数１】
＜推定部＞
推定部６で推定されるアフィンパラメータＡf は、実世界（三次元のＸＹＺ空間）における座標変換を、射影面である画像面（二次元のｕｖ平面）上での座標変換として表す際に用いる変換行列であり、（３）（４）式に示す関係を有する。
【００４３】
但し、ある特徴点の三次元座標を（Ｘ，Ｙ，Ｚ）T ，座標変換後の三次元座標を（Ｘ’，Ｙ’，Ｚ’）T 、これら三次元座標を画像面上に射影した二次元座標、すなわち画像面上での特徴点の座標を（ｕ，ｖ）T，（ｕ’，ｖ’）Tとして、Ｘ＝ｕ，Ｙ＝ｖ（すなわちＸ’＝ｕ’，Ｙ’＝ｖ’）であると仮定する。
【００４４】
【数２】
＜頭部姿勢算出部＞
頭部姿勢算出部７では、頭部姿勢情報として、頭部の３軸角度θｘ，θｙ，θｚを（５）〜（７）式を用いて算出する。
【００４５】
【数３】
そして、フレームｆにおける頭部姿勢を表すアフィンパラメータＡf が得られれば、そのアフィンパラメータＡf を顔モデルＣa に作用させることで、フレームｆにおける画像に映っていない隠れた特徴点の位置も含めて、フレームｆにおける全ての特徴点の位置を予測することができる。
【００４６】
また、このようなアフィンパラメータＡf には、カメラパラメータ（カメラの焦点距離、設置場所、向き等）及びドライバの頭部姿勢の変化を全て反映させることができる。
【００４７】
＜特徴点追跡部及び推定部での処理の概要＞
ここで、図５は、特徴点追跡部４及び推定部６が扱う状態空間モデルを表すグラフである。
【００４８】
図５に示すように、フレームｆにおける顔モデルの姿勢を表すアフィンパラメータＡf を隠れ状態として求める上位層（推定部６）と、フレームｆにおける特徴点群の位置ｘf ＝（ｘf(1)，ｘf(2)，…，ｘf(P)）T を隠れ状態として求めるＰ個の追跡器からなる下位層（特徴点追跡部４）を備えており、フレーム１からｆまでに入力された入力画像系列ｓ1:f から、これら隠れ状態Ａf ，ｘf を、時系列ベイズ推定により推定するものである。
【００４９】
なお、時系列ベイズ推定は、図６に示すように、状態変数をすべて確率変数として扱い、隠れ状態に関するフレームｆ−１における推定結果（推定分布）からフレームｆにおける予測分布を求めると共に、フレームｆにおける観測データから、検出すべき隠れ状態であるもっともらしさ（尤度）を求め、これら予測分布と尤度を考慮してフレームｆにおける推定結果（推定分布）を得るという手順を繰り返すことで、逐次的に隠れ状態を推定するものである。
【００５０】
つまり、入力画像系列（観測データ）ｓ1:f に基づいて、アフィンパラメータ（頭部姿勢の隠れ状態）Ａf を推定した事後確率（推定分布）ｐ（Ａf｜ｓ1:f）を求めることになり、これを数式で記述すると（８）（９）式で表される。
【００５１】
【数４】
ここで、ｐ（Ａf｜ｓ1:f-1）はアフィンパラメータＡf の事前確率分布（予測分布）、ｐ（ｘf｜Ａf，ｓ1:f-1 ）は特徴点群の位置ｘf の事前確率分布（予測分布）、ｐ（ｘf ｜Ａf ，ｓ1:f-1 ），ｐ（ｓf ｜ｘf ）は尤度を表す。
【００５２】
そして、（８）式の部分を下位層（特徴点追跡部４）がパーティクルフィルタを用いて推定し、（９）式の部分を上位層（推定部６）がカルマンフィルタを用いて推定する。
【００５３】
なお、あるパラメータの確率分布がガウス分布（正規分布）に従う場合、そのパラメータの確率分布の計算は、実際には、その平均ｘf と分散Ｖf を求めれば十分である。
【００５４】
なお、このような推定部６の具体的な構成は、例えば、特願２００５−３６８１２４号等に、また、特徴点追跡部４の具体的な構成は、例えば、特願２００７−１０２９１４号等に詳述されているため、ここではその詳細についての説明を省略する。
【００５５】
＜顔モデル供給部＞
次に、本発明の主要部である顔モデル供給部５について詳述する。
【００５６】
顔モデル供給部５は、図１に示すように、頭部を撮影した複数人分の映像から、それぞれ固定特徴点を抽出し、その抽出した結果に基づいて設定された平均的な固定特徴点の配置を表す平均顔モデルを記憶する平均顔モデル記憶部５１と、外部からの起動指令によって起動し、特徴点追跡部４から逐次供給される特徴点推定分布パラメータｘf に基づいて、カメラ２の被写体であるドライバ頭部に由来する特徴点（固定特徴点及び任意特徴点）の配置を表す個人顔モデルを、学習，生成する個人顔モデル学習部５２と、個人顔モデル学習部５２からの切替信号に従って、平均顔モデル記憶部５１に記憶された平均顔モデル、個人顔モデル学習部５２にて生成された個人顔モデルのいずれかを、顔モデルＣa として選択して推定部６に供給する顔モデル切替部５３とを備えている。
【００５７】
なお、個人顔モデル学習部５２が本発明の三次元形状復元装置に相当し、複数の物体（例えば、ドライバ，車両，車両の窓に映る車外の風景など）に由来する特徴点ｘf が、それぞれどの物体に由来するものであるかを示す情報と、特徴点が由来する各物体（特にドライバの頭部）の形状とを同時に学習するものである。
【００５８】
＜学習の概要＞
個人顔モデル学習部５２では、以下に示す手法により学習を行う。
【００５９】
＜モデル＞
まず、学習に適用するモデルについて説明する。
【００６０】
ここでは、物体が複数存在する場合の観測データ（特徴点群）の生成過程を考えると、生成過程が複数あって、各過程が確率的に選択されるMixtures of PPCA（Probabilistic Principal Component Analysis;確率的主成分分析）によってモデル化するものとする。
【００６１】
まず、カメラによる画像の撮影を、三次元空間から二次元平面への線形写像（アフィンカメラ）と仮定する。撮像される物体の形状は時不変の剛体だとして、Ｐ個の特徴点の三次元位置を表す行列ｚを（１０）式、フレームｆの画像上での各特徴点の位置を表す行列ｘf を（１１）式で表すものとする。但し、ｐ番目の特徴点の三次元位置をｚ(p) 、フレームｆにおけるｐ番目の特徴点の二次元位置をｘf(p)とする。
【００６２】
【数５】
アフィンカメラの下で特徴点の二次元位置ｘf(p)の生成過程は、観測ノイズεが加わって（１２）式で表すことができる。
【００６３】
【数６】
但し、Ｗf は射影行列、μf は並進成分を表す。
【００６４】
そして、Ｐ個の特徴点をフレーム数Ｆに渡って観測（特徴点追跡部４から取得）したとすると、（１３）式に示す観測データ行列ｘが得られる。
【００６５】
【数７】
但し、Ｗe は（１４）式で定義され、以下では運動行列と称する。また、ｚe は（１５）式で定義され、以下では形状行列と称する。
【００６６】
【数８】
画像に複数の物体が撮影されている場合は、それぞれの物体に対して運動行列Ｗe が存在するため、物体ｕ（ｕ＝１，２，…，Ｕ）に対応するものをＷe/u で表すものとする。但し、画像上では実際には各特徴点がどの物体に由来するのかわからないため、各物体の特徴点が混在する観測データ行列ｘが測定されることになる。
【００６７】
また、物体ｕに関して測定された特徴点の数をＰu で表すものとすると、特徴点の総数Ｐは、（１６）式で表される。
【００６８】
【数９】
ここで、特徴点ｐがどの物体に由来するかを指定する指示変数ｈ(p) を（１７）式で表すものとする。この指示変数は、仮定する物体の総数と同じＵ個ｈu(p)の要素を持つ。
【００６９】
【数１０】
そして、特徴点ｐについての指示変数ｈ(p) は、特徴点ｐが物体ｕに由来するものであるということを表現する場合、ｕ番目の要素のみが１（ｈu(p)＝１）となり、他の要素は全て０となる。また、指示変数ｈ(p) は、特徴点の数Ｐだけ存在するため、これをまとめたものを指定行列ｈと呼び、（１８）式で表すものとする。
【００７０】
【数１１】
また、全特徴点の内からランダムに１点を選択した時に、その選択した特徴点が物体ｕに属している確率をｇu として、これを全ての物体についてまとめたものを混合比ｇと呼び、（１９）式で表すものとする。
【００７１】
【数１２】
なお、各物体に関して測定されている特徴点数（Ｐ1，…，ＰU）が既知である場合、ｇu はＰ個の特徴点の中からランダムに選択して特定のＰu 個の特徴点のうちのどれかを選ぶ確率なのでｇu＝Ｐu／Ｐとなる。但し、Ｐu は未知なので、ｇも推定する必要がある。
【００７２】
また、物体ｕに関する測定に混入するノイズεは期待値０、分散τ|u-1の正規分布に従うものとする。
【００７３】
特徴点ｐが物体ｕに属すると仮定すると、特徴点ｐのノイズの分布ｐ（εp ）は、（２０）式で表されることになる。但し、ｄ次元ベクトルとして表された確率変数ｘが期待値ｍ、共分散Σの正規分布に従うことをＮ（ｘ｜ｍ，Σ）で示し、その内容は（２１）式で定義されるものとする。
【００７４】
【数１３】
以上の設定から、形状行列ｚ、指示行列ｈ、混合比ｇ、物体毎のアフィンパラメータＷe|u （ｕ＝１，２，…，Ｕ），各物体の特徴点位置の測定の正確さ（ノイズの分散の逆数）τ|uが推定（学習）の対象となる。
【００７５】
これらの関係を図示したグラフィカルモデルを図７に示す。
【００７６】
なお、以下では、式の記述を簡略化するために、（２２）（２３）式に示すモデルパラメータθ，Ｊを用いるものとする。
【００７７】
【数１４】
図７に示すグラフィカルモデルの下では、特徴点（観測データ）ｘが得られれば、形状行列ｚとモデルパラメータθは事後確率ｐ（θ，ｚ，ｈ｜ｘ）を計算することで推定することができる。
【００７８】
また、事後分布（事後確率）ｐ（θ，ｚ，ｈ｜ｘ）はベイズ推定により、（２４）式のように記述される。
【００７９】
【数１５】
この式を計算するために、形状行列ｚ，モデルパラメータθ（即ち、Ｗe|u ，τ|u，ｇ）の事前分布ｐ0 を（２５）〜（２８）式によって設定する。但し、αe|u は、４次元ベクトルであり、ここではαe|u ＝［１，１，１，１］Tが用いられるが、何等かの分布を設けてもよい。
【００８０】
【数１６】
但し、ｗe,i|u はＷe|u のｉ行目に対応した列ベクトルである。また、ガンマ分布Ｇａｍ（ｘ｜ａ，ｂ）は（２９）式、ディリクレ分布Ｄｉｒ（ｘ｜γx ）は（３０）式で定義される。
【００８１】
【数１７】
なお、（２４）式を解析的に計算することは不可能なため、ここでは変分推定法を利用して事後分布を近似的に求める。
【００８２】
＜変分推定法＞
変分推定法では、事後分布ｐ（θ，ｈ，ｚ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）を計算する。但し、ｑ（θ，ｈ，ｚ）＝ｑ（θ）ｑ（ｈ，ｚ）と因子化されると仮定して、試験分布ｑ（θ，ｈ，ｚ）と真の事後分布ｐ（θ，ｈ，ｚ｜ｘ）との間のＫＬ（Kullback-Leibler）疑似距離を最小化するように、試験分布ｑ（θ）ｑ（ｈ，ｚ）を調節する。
【００８３】
ところで、ＫＬ疑似距離を最小にする試験分布を求めるということは、自由エネルギーを最大にする試験分布を求めることに等しいことが知られている。つまり、変分推定法では、この自由エネルギーを最大にする試験分布が最も妥当な推定解であるということになる。
【００８４】
本件の場合、自由エネルギーＦE は（３１）式で表される。また、この場合、自由エネルギーＦE を最大化する試験分布ｑ（θ，ｈ，ｚ）は、事前の同時分布ｐ0（Ｊ） に基づいて、（３２）式及び（３３）式に示す条件を満たすことが知られている。
【００８５】
【数１８】
但し、<>ｑは、分布ｑに基づく期待値を求める操作を意味する。
【００８６】
この条件式（３２）（３３）を展開して得られる式は、推定対象であるθ，ｈ，ｚが互いに依存した形式となっているため、明示的な解を得ることができない。しかし、条件式（３２）（３３）を繰り返し適用することは、自由エネルギーＦE を対象とした勾配法を行うことに等しいため、近似的な解を得ることができる。
【００８７】
具体的には、変分ベイズ法のＥステップでは、試験分布ｑ（θ）を固定した上で、試験分布ｑ（ｈ，ｚ）を求め、変分ベイズ法のＭステップでは、試験分布ｑ（ｈ，ｚ）を固定した上で、試験分布ｑ（θ）を求めるものとする。
【００８８】
ここで、Ｅステップで求める各試験分布を（３４）〜（３６）式のようにパラメータ付けする。
【００８９】
【数１９】
また、Ｍステップで求める各試験分布を（３７）〜（４０）式のようにパラメータ付けする。
【００９０】
【数２０】
（２５）〜（２８）式と、以下に示す（４１）式に基づいて、（３２）（３３）式を展開する。
【００９１】
【数２１】
すると、Ｅステップで求める各試験分布のパラメータには、次の（４２）〜（４６）式に示す条件が課されることになる。但し、Σz|u は物体ｕの推定構造の共分散行列（信頼度）、ｚ|u(p) は特徴点ｐが物体ｕ由来であると仮定したときの推定座標、ｃu(p)は特徴点ｐが物体ｕ由来であるという仮定の信頼度、ｃu は試験分布ｑ（ｈ(p) ）の正規化項を表す。
【００９２】
【数２２】
同様に、Ｍステップで求める各試験分布のパラメータには、次の（４７）〜（５２）式に示す条件が課されることになる。但し、γguの比は混合比の期待値（推定値）、ＭWe|uはＷe|u の期待値（推定値）、Σｗe|u はＷe|u の推定の各行の共分散行列（信頼度）、ｍτ|u はτ|uの期待値（推定値）、γτ|u はτ|uの推定の分散（信頼度）ｍτ|u／γτ|uの構成要素を表す。
【００９３】
【数２３】
なお、Ｅステップで計算される全てのパラメータをまとめた分布パラメータΘS は（５３）式で表され、Ｍステップで計算される全てのパラメータをまとめた分布パラメータΘM は（５４）式で表される。
【００９４】
【数２４】
また、ｈ，ｚに依存する十分統計量ｒを（５５）式で表すものとする。具体的には、（４７）（４８）（４９）（５１）（５２）式のそれぞれにおいて、<<>>fcで囲われた項が十分統計量となる。
【００９５】
【数２５】
＜個人顔モデル学習部での処理＞
ここで、個人顔モデル学習部５２が実行する処理を、図８に示すフローチャートに沿って説明する。なお、本処理は、外部からの起動指令の入力により起動される。
【００９６】
本処理が起動すると、まず、Ｓ１１０では、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、平均顔モデル記憶部５１に記憶された平均顔モデルを選択するように切替信号を出力して、Ｓ１２０に進む。
【００９７】
これにより、推定部６では、画像（フレーム）が入力される毎に、平均顔モデルを用いたアフィンパラメータＡf 、及び特徴点推定分布パラメータｘf，Ｖfの算出が行われる。また、頭部姿勢算出部７では、推定部６にて求められたアフィンパラメータＡf に基づいて、頭部姿勢θｘ，θｙ，θｚの算出が行われる。
【００９８】
Ｓ１２０では、モデルパラメータθの試験分布である第１試験分布ｑ（θ）を初期化してＳ１３０に進む。
【００９９】
具体的には、（５４）式に示した分布パラメータΘM 及び、（５５）式に示した十分統計量ｒ初期化され、例えば、全ての要素が、同じ一定値（例えば、０又は１又は０．５）に設定される。
【０１００】
Ｓ１３０では、特徴点追跡部４から予め設定された規定フレーム数Ｆ個分の特徴点パラメータｘf （（１１）式参照）の取得を完了したか否かを判断し、完了していなければ、完了するまで待機する。
【０１０１】
そして規定フレーム数Ｆ分の特徴点パラメータｘf の取得を完了すると、これらをまとめた観測データ行列ｘ（（１３）式参照）を生成してＳ１４０に移行し、設定されている分布パラメータΘM に基づいて、第２試験分布ｑ（ｚ，ｈ）の算出、即ち、変分ベイズのＥステップとしての処理を実行して、Ｓ１５０へ進む。
【０１０２】
具体的には、（５３）式に示す分布パラメータΘS を、（４２）〜（４６）式に従って算出する。
【０１０３】
Ｓ１５０では、観測データ行列ｘ及びＳ１４０での第２試験分布ｑ（ｚ，ｈ）の算出結果に従って、十分統計量ｒを更新し、続くＳ１６０では、第１試験分布ｑ（θ）の算出、即ち、変分ベイズのＭステップとしての処理を実行して、Ｓ１７０に進む。
【０１０４】
具体的には、（５４）式に示す分布パラメータΘM を、Ｓ１５０にて算出された十分統計量ｒを用い、（４７）〜（５２）式に従って算出する。
【０１０５】
Ｓ１７０では、Ｓ１４０，Ｓ１６０での算出結果に基づき、上述の（３１）式を用いて自由エネルギーＦE を算出してＳ１９０に進む。
【０１０６】
但し、自由エネルギーＦE は、上述の処理によって得られた第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）から求められる運動行列Ｗおよび形状行列ｚに基づいて推定される位置と、観測により得られた位置（観測データ行列ｘ）とが近いほど大きな値となる。
【０１０７】
Ｓ１８０では、Ｓ１７０にて算出された自由エネルギーＦE が、予め設定されたエネルギー閾値より大きいことを終了条件として、その終了条件を充足しているか否かを判断し、充足していなければ、Ｓ１４０に戻って、Ｓ１４０〜Ｓ１７０の処理を繰り返し、一方、終了条件を充足していれば、Ｓ１９０に移行する。
【０１０８】
ここでは、エネルギー閾値として、例えば、平均顔モデルを表す形状行列と、上述の処理で得られた運動行列Ｗとを用いて算出した自由エネルギー、またはその自由エネルギーより大きな値に設定された固定値を用いる。
【０１０９】
Ｓ１９０では、各特徴点ｐ（ｐ＝１，２，…，Ｐ）について、その特徴点が由来する物体ｕを、（５６）式に示す判別式を用いて判別し、その判別結果からドライバの頭部に由来する特徴点（以下「頭部特徴点」という）のみを抽出して特徴点追跡部４に通知する。
【０１１０】
【数２６】
これにより、特徴点追跡部４では、以後、通知された頭部特徴点に基づいて、任意特徴点の中からドライバの頭部に由来するもののみを選択して追跡することが可能となる。
【０１１１】
Ｓ２００では、第２試験分布ｑ（ｈ，ｚ）から求められた形状行列のうち、Ｓ１９０にて得られた頭部特徴点に基づく部分のみを抽出した部分形状行列を生成し、その部分形状行列を個人顔モデルとして出力すると共に、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、個人顔モデルを選択するように切替信号を出力して、本処理を終了する。
【０１１２】
これにより、以後、推定部６では、画像（フレーム）が入力される毎に、個人顔モデルを用いたアフィンパラメータＡf の算出が行われ、更に、そのアフィンパラメータＡf に基づいて、頭部姿勢算出部７では、頭部姿勢θｘ，θｙ，θｚの算出が行われる。
【０１１３】
＜効果＞
以上説明したように、画像処理装置１では、頭部姿勢θｘ，θｙ，θｚを推定する処理を開始した直後は、固定特徴点に基づく平均顔モデルを用いてアフィンパラメータＡf の算出（顔向き推定）を行い、これと並行して、任意特徴点も利用した個人顔モデルの学習を実行し、その個人顔モデルの精度が十分に向上した時点で、平均顔モデルから個人顔モデルへの切替を行うようにされている。
【０１１４】
従って、画像処理装置１によれば、ドライバが誰であっても、また、ドライバが身につける物品（サングラス，マスク，帽子，スカーフ等）によって、固定特徴点を抽出することができない場合であっても、頭部姿勢θｘ，θｙ，θｚの推定を、処理の初期段階から、安定した精度の推定結果を得ることができ、しかも、平均顔モデルから個人顔モデルへの切替後は、高い精度の推定結果を得ることができる。
【０１１５】
また、画像処理装置１では、観測データ行列ｘを因子分解する際に、指定行列ｈを導入して、各特徴点がどの物体に由来するかを同時に推定するようにされている。
【０１１６】
従って、画像処理装置１によれば、観測データ行列ｘに、ドライバ頭部に由来する特徴点（頭部特徴点）だけでなく、ドライバの頭部以外の物体（例えば、車両や窓に映る車外の景色等）に由来する特徴点（非頭部特徴点）が含まれていたとしても、即ち、観測データ行列ｘが混合データ行列であったとしても、その観測データ行列ｘを因子分解した結果から、物体の形状を個別に推定することができ、ひいてはドライバ頭部の形状を精度よく推定することができる。
【０１１７】
また、画像処理装置１では、個人顔モデル学習部５２にて得られた頭部特徴点の情報を、特徴点追跡部４にて供給するようにされているため、特徴点追跡部４では、頭部特徴点の情報に基づいて任意特徴点を選別することができ、その選別結果に基づき、例えば、非頭部特徴点の追跡を中止することによって、不要な追跡処理が実行されることを防止することができる。
【０１１８】
［第２実施形態］
次に、第２実施形態について説明する。
【０１１９】
なお、本実施形態では、個人顔モデル学習部５２での処理内容の一部が、第１実施形態のものとは異なるだけであるため、この相違する部分を中心に説明する。
【０１２０】
詳しくは、第１実施形態では、規定フレーム数Ｆ分の特徴点データが得られた後に、因子分解を一括して行う、いわゆるバッチ手法を用いているが、本実施形態では、１フレーム分の特徴点データを取得する毎に、逐次的にθ，ｈ，ｚを推定する、所謂、オンライン推定法を実行するものである。
【０１２１】
＜オンライン推定法の概要＞
オンライン推定法では、（３１）式に示す自由エネルギーＦE を、（５７）式に示すように、オンライン推定に適した形式のオンライン自由エネルギーＦEλに置き換えて使用する。
【０１２２】
但し、忘却係数をλ（但し、０＜λ＜１）、現フレームをｆc 、フレームｆまでの観測に基づいた試験分布を（５８）式、十分統計量を（５９）式、忘却係数の影響を正規化するための係数を（６０）式で表すものとする。
【０１２３】
【数２７】
（３１）式と（５７）式を比較すると、オンライン自由エネルギーＦEλの十分統計量ｒを、（６１）式に示す十分統計量ｒfのｆ＝１，２，…，ｆｃを通した減衰付き和で置き換えたものとなっている。
【０１２４】
そして、十分統計量ｒfの減衰係数付き和は（６２）式に示すように逐次的に計算することができ、また、正規化項も（６３）式に示すように逐次的に計算することができる。
【０１２５】
【数２８】
＜個人顔モデル学習部での処理＞
図９は、個人顔モデル学習部５２が実行する処理の内容を示すフローチャートである。
【０１２６】
第１実施形態の場合と同様に、本処理は、初期化指令の入力により起動され、Ｓ３１０〜Ｓ３４０は、第１実施形態のＳ１１０〜Ｓ１４０と同様であるため説明を省略する。
【０１２７】
但し、Ｓ１３０にて使用される規定個数は、第１実施形態の場合より少ない数（例えば、数分の一から数十分の一程度）でよく、ここでは１とされている。
【０１２８】
そして、Ｓ３５０では、観測データ行列ｘf 及びＳ３４０での第２試験分布ｑ（ｚ，ｈ）の算出結果に基づき、（６２）（６３）式を用いて、十分統計量<<ｒ>>fcを更新し、続くＳ３６０では、第１試験分布ｑ（θ）の算出、即ち、変分ベイズのＭステップとしての処理を実行して、Ｓ１７０に進む。
【０１２９】
具体的には、（５４）式に示す分布パラメータΘM を、Ｓ３５０にて算出された十分統計量<<ｒ>>fcを、（４７）〜（５２）式の十分統計量の項（<<>>ｆｃで囲われた項）に代入して算出する。
【０１３０】
Ｓ３７０では、Ｓ１９０と同様に頭部特徴点を抽出して特徴点追跡部４に通知し、続くＳ３８０では、（５７）式に従ってオンライン自由エネルギーＦEλを算出し、Ｓ３９０では、Ｓ１８０と同様に、その算出したオンライン自由エネルギーＦEλが、予め設定されたエネルギー閾値より大きいことを終了条件として、その終了条件を充足しているか否かを判断する。
【０１３１】
そして、終了条件を充足していなければ、Ｓ３３０に戻って、Ｓ３３０〜Ｓ３８０の処理を繰り返し、一方、終了条件を充足していれば、Ｓ４００に移行する。
【０１３２】
Ｓ４００では、Ｓ２００と同様に、頭部特徴点に基づく部分形状行列を個人顔モデルとして出力すると共に、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、個人顔モデルを選択するように切替信号を出力して、本処理を終了する。
【０１３３】
＜効果＞
以上説明したように、本実施形態では、１フレーム毎に、Ｅステップ，Ｍステップを１回ずつ実行して第１試験分布ｐ（θ）及び第２試験分布ｐ（ｈ，ｚ）を求める処理を繰り返し、その結果を過去のデータほど比率が低下するように減衰係数ηfcを乗じて足し合わせていく、いわゆるオンライン化した変分ベイズ推定法を用いて、因子分解を行っている。
【０１３４】
従って、本実施形態によれば、使用する画像系列の数（フレーム数）を増加させても、それによって増加する処理量は、Ｅステップ，Ｍステップを１回ずつ実行する際の処理量を単位として、加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【０１３５】
その結果、自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも問題なく適用することができる。
【０１３６】
更に本実施形態では、Ｍステップで求める第１試験分布を規定するパラメータの算出式中に現れる十分統計量<<ｒ>>fc 及び減衰係数ηfcを、前フレームでの算出結果を利用した簡易な式（（６２）（６３）式参照）によって逐次的求めるようにされているため、因子分解に要する処理負荷を、更に削減することができる。
【０１３７】
また、本実施形態では、十分統計量<<ｒ>>fcが求められれば観測データ行列ｘfcは不要であり、特徴点の増減に対応するために、過去の観測データ行列を記憶しておく必要がないため、装置の記憶容量も削減することができる。
【０１３８】
［実験］
図１０に示すような二つの物体（四角柱，三角柱）を適当に動かして画像に投影したものを入力画像系列として、第２実施形態で説明した処理を実行する画像処理装置１に処理を行わせて形状行列ｚを求めた。なお、特徴点としては、各物体の頂点を用いた。
【０１３９】
図１１は、処理の結果得られた形状行列ｚから推定し二つの物体の形状の推定誤差（ＲＭＳ誤差）の推移を観測フレーム数に対して示したものであり、（ａ）が四角柱についての結果、（ｂ）が三角柱についての結果である。
【０１４０】
図１１からは、フレーム数が少ないと推定精度が低いが、フレーム数の増加に共に推定精度が向上していること、及び、二つの物体を正しく分離できていることがわかる。
【０１４１】
［他の実施形態］
以上本発明の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、様々な態様での実施が可能である。
【０１４２】
例えば、上記実施形態では、頭部姿勢を推定するためにアフィンパラメータを用いているが、これは頭部姿勢の推定を行うことができれば、例えば、アフィンパラメータを用いずに、拡張カルマンフィルタを利用して直接頭部推定を行う手法等を採用してもよい。
【０１４３】
また、因子分解を終了させる際の終了条件として、平均顔モデルから求めた自由エネルギーに基づいて設定されたエネルギー閾値を用いているが、例えば、処理の繰り返し回数が、予め設定された回数閾値より大きいことを終了条件としてもよい。この場合、処理を終了させるか否かの判断のために、自由エネルギーを算出する必要がなくなるため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。
【図面の簡単な説明】
【０１４４】
【図１】発明が適用された画像処理装置の全体構成を示すブロック図。
【図２】固定特徴点の説明図。
【図３】任意特徴点の説明図。
【図４】画像処理装置が使用する座標系を示す説明図。
【図５】特徴点追跡部及び推定部が扱う状態空間モデルを表すグラフ。
【図６】時系列ベイズ推定の動作概要を示す説明図。
【図７】個人顔モデル学習部が扱うグラフィックモデルを示す説明図。
【図８】第１実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【図９】第２実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【図１０】実験に用いた物体の形状を示す説明図。
【図１１】物体の形状と実験により得られた形状との誤差を示すグラフ。
【符号の説明】
【０１４５】
１…画像処理装置 ２…カメラ ３…画像取込装置 ４…特徴点追跡部 ５…顔モデル供給部 ６…推定部 ７…頭部姿勢算出部 ４１…固定特徴点追跡部 ４３…任意特徴点追跡部 ５１…平均顔モデル記憶部 ５２…個人顔モデル学習部 ５３…顔モデル切替部
【技術分野】
【０００１】
本発明は、物体を撮影した画像（二次元形状）系列から、その物体の三次元形状を復元する三次元形状復元装置に関する。
【背景技術】
【０００２】
従来より、物体を撮影した画像系列を入力し、その物体について予め設定された複数の特徴点について、その特徴点の位置を画像面上で追跡し、追跡により得られた特徴点位置系列から物体の三次元形状を復元する装置が知られている。
【０００３】
この主の三次元復元装置の一つとして、特徴点位置系列から観測データ行列を生成し、その観測データ行列を、物体の形状（特徴点の三次元位置）を表す形状行列と、物体の運動を表す運動行列とに因子分解する手法を用いたものがある（例えば、特許文献１〜３参照）。
【０００４】
但し、因子分解法は、観測データ行列を構成するデータが全て揃っていないと、有効な結果が得られないため、特徴点の隠れ，誤追跡，フレームアウト等によって、特徴点位置系列（ひいては観測データ行列を構成するデータ）の欠損が頻繁に生じる実問題への適用は、非常に困難であるという問題があった。
【０００５】
これに対して出願人は、データを統計的に扱い、データの確率分布の分散をデータの信頼度として考慮した計算を行うことによって、データの欠損に対処する手法を提案している（特許文献４参照）。
【０００６】
しかし、この手法でも、形状を推定すべき物体が人物の顔である場合、サングラスやスカーフ等の装飾品によって特徴点が隠れた状態が継続すると、検出精度が劣化してしまう。
【０００７】
このため、予め固定的に設定された特徴点（以下「固定特徴点」という）だけでなく、その都度検出される予め規定された特徴（例えば、撮像された物体のエッジやコーナー等）を有する地点（以下「任意特徴点」という）を検出し、この任意特徴点も利用することが検討されている。
【特許文献１】特開２０００−１９４８６３号公報
【特許文献２】特開２００３−１４１５５２号公報
【特許文献３】特開２００４−３４７３３２号公報
【特許文献４】特開２００８−２６１６６２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかし、任意特徴点は、物体上に限らず、物体の背景や物体以外の物体上であっても、その特徴を有していれば抽出されてしまう。
【０００９】
しかし、特許文献１〜４に記載の従来装置では、いずれも、観測データ行列を構成する特徴点が、単一の対象物に由来していることを前提とするものであるため、特徴点が複数の対象物に由来している場合には、形状行列（ひいては物体の形状）を正しく求めることができないという問題があった。
【００１０】
なお、各特徴点が、どの物体に由来するかを識別し、その識別された特徴点を用いて物体毎に作成した観測データ行列を、それぞれ因子分解する手法が考えられる。しかし、特徴点の識別を、別の画像処理等によって行う必要があり、処理量が増大するだけでなく、その識別が困難であり、十分な精度が得られないという問題があった。
【００１１】
また、観測データ行列を因子分解して形状行列を得る手法において、三次元形状の推定精度を向上させるには、計算に用いる画像系列を増やす（即ち、観測データ行列の次元を大きくする）必要がある。
【００１２】
しかし、因子分解の際に必要となる計算量は、通常、観測データ行列の行数×列数の３乗に比例することが知られており、観測データ行列の次元を大きくすることによって爆発的に増大する。
【００１３】
そして、計算量が増大すると、三次元形状復元装置を車両に搭載する場合等のように、使用できる演算器の計算能力に制限がある場合には、適用することが困難になるという問題もあった。
【００１４】
本発明は、上記問題点を解決するために、三次元形状復元装置において、観測データ行列（特徴点位置系列）中に複数の物体に由来するデータが混在していても因子分解の実行を可能とすること、更には、計算時間を低減することを目的とする。
【課題を解決するための手段】
【００１５】
上記目的を達成するためになされた本発明の三次元形状復元装置では、特徴点追跡手段が、画像系列を入力し、入力された各画像から予め設定された特徴を有する複数の特徴点を抽出すると共に、特徴点のそれぞれについて該特徴点の位置を画像面上で追跡する。
【００１６】
また、観測データ行列生成手段は、特徴点抽出追跡手段での抽出，追跡結果に従って、特徴点の位置を表す観測データ行列を生成する。そして、分解手段が、観測データ行列生成手段にて生成された観測データ行列を因子分解して、画像に写された物体の三次元形状を表す形状行列および前記物体の運動を表す運動行列を生成する。
【００１７】
このとき、分解手段は、観測データ行列をｘ、形状行列をｚ、特徴点のそれぞれが属する物体を指定する指定行列をｈ、運動行列，特徴点の測定の正確さ，特徴点がそれぞれの物体に属している確率をまとめたモデルパラメータをθとして、観測データ行列ｘから事後確率Ｐ（θ，ｚ，ｈ｜ｘ）を計算することで、θ，ｚ，ｈを求める。
【００１８】
このように構成された本発明の三次元形状復元装置によれば、観測データ行列ｘが複数の物体に由来する特徴点で構成された混合データ行列であっても、因子分解を実行することができ、指定行列ｈによって分類される特徴点毎、即ち、特徴点が由来する物体毎に、形状行列や運動行列を求めることができる。
【００１９】
なお、本発明において、分解手段は、請求項２に記載のように、事後確率Ｐ（θ，ｚ，ｈ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）が、第１試験分布ｑ（θ）と第２試験分布ｑ（ｈ，ｚ）とに因子化されるものとして、事後確率Ｐ（θ，ｚ，ｈ｜ｘ）と試験分布ｑ（θ，ｈ，ｚ）とのＫＬ距離を最小にする第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）を求めることで、確率的に因子分解計算を実行するように構成されていてもよい。
【００２０】
このように構成された本発明の三次元形状復元装置によれば、形状行列ｚを求めるのと同時に、指定行列ｈもベイズ推定によって確率的に求められるため、因子分解以外の特別な処理を行うことなく、簡単に特徴点を物体毎に識別することができる。
【００２１】
そして、このようにベイズ推定を用いた因子分解を行う場合、請求項３に記載のように、観測データ行列生成手段は、予め設定されたフレーム数分の画像に基づいて観測データ行列ｘを生成し、分解手段は、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、予め設定された終了条件が満たされている場合に、前記推定手段を終了させる終了手段とからなり、最終的に得られた前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得るように構成されていてもよい。
【００２２】
この場合、分解手段は、変分ベイズ推定法を実現するものであり、第１のステップがＥステップ、第２のステップがＭステップに相当する。
【００２３】
但し、変分ベイズ推定法を用いても、計算に用いる画像系列の増加（即ち、観測データ行列の次元の増大）に応じて、計算量（計算時間）が指数関数的に増大することを避けることはできない。
【００２４】
そこで、請求項４に記載のように、観測データ行列生成手段は、１フレーム分の画像毎に前記観測データ行列ｘを生成し、分解手段は、観測データ行列生成手段にて観測データ行列ｘが生成される毎に、第１試験分布ｑ（θ）を固定して第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、第２試験分布ｑ（ｚ，ｈ）を固定して第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、予め設定された終了条件が満たされている場合に、推定手段及び更新手段を終了させる終了手段とからなり、推定手段によって逐次的に得られる第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得るように構成されていていることが望ましい。
【００２５】
つまり、推定手段では、計算に用いる画像系列を制限して、第１のステップ（Ｅステップ），第２のステップ（Ｍステップ）を１回ずつ実行し、これを繰り返し実行すると共に、推定手段から計算結果が得られる毎に、その計算結果に基づいて形状行列の分布を更新する。これは、変分ベイズ推定法を逐次計算にて実現（いわゆるオンライン化）するものである。
【００２６】
このように、本発明の三次元形状復元装置によれば、使用する画像系列を複数に分割し、その分割した画像系列毎に因子分解を実行するため、画像系列を増やしても、処理量は推定手段での１回の処理を単位として加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【００２７】
そして、当該三次元形状復元装置は、請求項５に記載のように、自動車に搭載して使用するように構成されていてもよい。
【００２８】
自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも適用することができる。
【発明を実施するための最良の形態】
【００２９】
以下に本発明の実施形態を図面と共に説明する。
【００３０】
［第１実施形態］
図１は、本発明を適用した三次元形状復元装置としての画像処理装置１の全体構成を示すブロック図である。なお、画像処理装置１は、車両に搭載され、ドライバの頭部を撮影した画像から、ドライバのわき見や居眠り等を検出するための情報である頭部姿勢を時系列的に求める処理を実行する装置である。
【００３１】
＜全体構成＞
図１に示すように、画像処理装置１は、ドライバシートの前面上部又は下部（例えば、メータ内）に設置され、ドライバの頭部を含む画像を取得するカメラ２と、カメラ２が取得した画像を取り込む画像取込装置３と、画像取込装置３を介して取り込んだ入力画像ｓf （ｆはフレームを識別する識別子）から、予め設定された特徴を有する複数の特徴点の位置（画像面上での二次元座標）ｘ(p) （ｐ＝１，２，…Ｐ）を抽出すると共に、その抽出された特徴点を追跡して、その特徴点の位置の分布を規定するための特徴点推定分布パラメータ（平均値ｘf ，分散Ｖf ）を求める特徴点追跡部４とを備えている。
【００３２】
また、画像処理装置１は、特徴点の三次元位置（ひいては、頭部の三次元形状）を規定する顔モデル（形状行列）Ｃa を供給する顔モデル供給部５と、特徴点追跡部４から供給される特徴点推定分布パラメータ、顔モデル供給部５から供給される顔モデルＣa に基づき、カルマンフィルタを用いて、ドライバの頭部姿勢を表すアフィンパラメータＡf を推定する推定部６と、推定部６が推定したアフィンパラメータＡf に従って、頭部姿勢θｘ，θｙ，θｚを算出する頭部姿勢算出部７から構成される。
【００３３】
＜特徴点追跡部＞
特徴点追跡部４は、目尻，目頭，口元、鼻（鼻の穴，鼻の中心，鼻の左右端）等（図２（ａ）参照）、ドライバの顔の特徴を表す複数の特徴点（以下「固定特徴点」ともいう）の位置を抽出，追跡するする固定特徴点追跡部４１と、画像中の追跡し易い特徴（例えば、撮像された物体のエッジやコーナ等）を有する部位（以下「任意特徴点」ともいう）の位置を抽出，追跡する任意特徴点追跡部４３とからなる。
【００３４】
そして、固定特徴点追跡部４１は、予め用意されたテンプレートを用いて入力画像ｓf とのテンプレートマッチングによって、固定特徴点の抽出，追跡を行う。なお、テンプレートは、具体的には、図２（ａ）中の点線白枠にて示すようなものとなる。但し、図では、左口元、鼻の中心、左目尻についてのみ示す。
【００３５】
一方、任意特徴点追跡部４３は、ＳＩＦＴ（Scale Invariant Feature Transform） やＫＬＴ(Kanade-Lucas-Tomasi) といった周知の特徴点検出器を用いて任意特徴点の抽出，追跡を行う。
【００３６】
ここで、図３は、任意特徴点追跡部４３により抽出される任意特徴点の例を示す説明図である。但し、（ａ）が入力画像、（ｂ）が固定特徴点追跡部４１によって抽出される固定特徴点、（ｃ）が任意特徴点追跡部４３によって抽出される任意特徴点を示す。
【００３７】
図３に示すように、任意特徴点は、ドライバ頭部からだけでなく、背景（車室）からも抽出されることになる。但し、ドライバがサングラス等を装着して固定特徴点（目尻，目頭等）が隠れたとしても、サングラスの輪郭等から多数の任意特徴点が抽出されることになる。
【００３８】
つまり、特徴点追跡部４では、常時追跡対象となる固定特徴点以外に、任意特徴点が適宜、追跡対象として追加されることになる。
【００３９】
但し、任意特徴点追跡部４は、任意特徴点の中で固定特徴点と一致するものを追跡対象から除外すると共に、後述する顔モデル供給部５からドライバ頭部に位置する特徴点（以下「頭部特徴点」という）の情報が得られた後は、その情報に基づいて、頭部特徴点以外の任意特徴点を追跡対象から除外するように構成されている。
【００４０】
また、固定特徴点及び任意特徴点の追跡には、いずれも、パーティクルフィルタが用いられており、その処理に必要な事後確率ｐ（ｘf｜Ａf）は、推定部６から供給を受けるように構成されている。
【００４１】
＜顔モデル＞
顔モデル供給部５から供給される顔モデルＣa は、特徴点追跡部４にて追跡される特徴点のうち、ドライバの顔上に位置するものの位置を三次元的に規定するものであり、図４に示すように、画像面における水平方向の座標軸をｕ軸、垂直方向の座標軸をｖ軸、三次元空間における水平方向の座標軸をＸ軸、垂直方向の座標軸をＹ軸、光軸方向（奥行き方向）の座標軸をＺ軸として（１）（２）式で表される。但し、ｚ(p) は、特徴点追跡部４で抽出されたｐ番目の特徴点（特徴点位置ｘ(p) ＝（ｕ(p)，ｖ(p)）の３次元空間内での座標である。
【００４２】
【数１】
＜推定部＞
推定部６で推定されるアフィンパラメータＡf は、実世界（三次元のＸＹＺ空間）における座標変換を、射影面である画像面（二次元のｕｖ平面）上での座標変換として表す際に用いる変換行列であり、（３）（４）式に示す関係を有する。
【００４３】
但し、ある特徴点の三次元座標を（Ｘ，Ｙ，Ｚ）T ，座標変換後の三次元座標を（Ｘ’，Ｙ’，Ｚ’）T 、これら三次元座標を画像面上に射影した二次元座標、すなわち画像面上での特徴点の座標を（ｕ，ｖ）T，（ｕ’，ｖ’）Tとして、Ｘ＝ｕ，Ｙ＝ｖ（すなわちＸ’＝ｕ’，Ｙ’＝ｖ’）であると仮定する。
【００４４】
【数２】
＜頭部姿勢算出部＞
頭部姿勢算出部７では、頭部姿勢情報として、頭部の３軸角度θｘ，θｙ，θｚを（５）〜（７）式を用いて算出する。
【００４５】
【数３】
そして、フレームｆにおける頭部姿勢を表すアフィンパラメータＡf が得られれば、そのアフィンパラメータＡf を顔モデルＣa に作用させることで、フレームｆにおける画像に映っていない隠れた特徴点の位置も含めて、フレームｆにおける全ての特徴点の位置を予測することができる。
【００４６】
また、このようなアフィンパラメータＡf には、カメラパラメータ（カメラの焦点距離、設置場所、向き等）及びドライバの頭部姿勢の変化を全て反映させることができる。
【００４７】
＜特徴点追跡部及び推定部での処理の概要＞
ここで、図５は、特徴点追跡部４及び推定部６が扱う状態空間モデルを表すグラフである。
【００４８】
図５に示すように、フレームｆにおける顔モデルの姿勢を表すアフィンパラメータＡf を隠れ状態として求める上位層（推定部６）と、フレームｆにおける特徴点群の位置ｘf ＝（ｘf(1)，ｘf(2)，…，ｘf(P)）T を隠れ状態として求めるＰ個の追跡器からなる下位層（特徴点追跡部４）を備えており、フレーム１からｆまでに入力された入力画像系列ｓ1:f から、これら隠れ状態Ａf ，ｘf を、時系列ベイズ推定により推定するものである。
【００４９】
なお、時系列ベイズ推定は、図６に示すように、状態変数をすべて確率変数として扱い、隠れ状態に関するフレームｆ−１における推定結果（推定分布）からフレームｆにおける予測分布を求めると共に、フレームｆにおける観測データから、検出すべき隠れ状態であるもっともらしさ（尤度）を求め、これら予測分布と尤度を考慮してフレームｆにおける推定結果（推定分布）を得るという手順を繰り返すことで、逐次的に隠れ状態を推定するものである。
【００５０】
つまり、入力画像系列（観測データ）ｓ1:f に基づいて、アフィンパラメータ（頭部姿勢の隠れ状態）Ａf を推定した事後確率（推定分布）ｐ（Ａf｜ｓ1:f）を求めることになり、これを数式で記述すると（８）（９）式で表される。
【００５１】
【数４】
ここで、ｐ（Ａf｜ｓ1:f-1）はアフィンパラメータＡf の事前確率分布（予測分布）、ｐ（ｘf｜Ａf，ｓ1:f-1 ）は特徴点群の位置ｘf の事前確率分布（予測分布）、ｐ（ｘf ｜Ａf ，ｓ1:f-1 ），ｐ（ｓf ｜ｘf ）は尤度を表す。
【００５２】
そして、（８）式の部分を下位層（特徴点追跡部４）がパーティクルフィルタを用いて推定し、（９）式の部分を上位層（推定部６）がカルマンフィルタを用いて推定する。
【００５３】
なお、あるパラメータの確率分布がガウス分布（正規分布）に従う場合、そのパラメータの確率分布の計算は、実際には、その平均ｘf と分散Ｖf を求めれば十分である。
【００５４】
なお、このような推定部６の具体的な構成は、例えば、特願２００５−３６８１２４号等に、また、特徴点追跡部４の具体的な構成は、例えば、特願２００７−１０２９１４号等に詳述されているため、ここではその詳細についての説明を省略する。
【００５５】
＜顔モデル供給部＞
次に、本発明の主要部である顔モデル供給部５について詳述する。
【００５６】
顔モデル供給部５は、図１に示すように、頭部を撮影した複数人分の映像から、それぞれ固定特徴点を抽出し、その抽出した結果に基づいて設定された平均的な固定特徴点の配置を表す平均顔モデルを記憶する平均顔モデル記憶部５１と、外部からの起動指令によって起動し、特徴点追跡部４から逐次供給される特徴点推定分布パラメータｘf に基づいて、カメラ２の被写体であるドライバ頭部に由来する特徴点（固定特徴点及び任意特徴点）の配置を表す個人顔モデルを、学習，生成する個人顔モデル学習部５２と、個人顔モデル学習部５２からの切替信号に従って、平均顔モデル記憶部５１に記憶された平均顔モデル、個人顔モデル学習部５２にて生成された個人顔モデルのいずれかを、顔モデルＣa として選択して推定部６に供給する顔モデル切替部５３とを備えている。
【００５７】
なお、個人顔モデル学習部５２が本発明の三次元形状復元装置に相当し、複数の物体（例えば、ドライバ，車両，車両の窓に映る車外の風景など）に由来する特徴点ｘf が、それぞれどの物体に由来するものであるかを示す情報と、特徴点が由来する各物体（特にドライバの頭部）の形状とを同時に学習するものである。
【００５８】
＜学習の概要＞
個人顔モデル学習部５２では、以下に示す手法により学習を行う。
【００５９】
＜モデル＞
まず、学習に適用するモデルについて説明する。
【００６０】
ここでは、物体が複数存在する場合の観測データ（特徴点群）の生成過程を考えると、生成過程が複数あって、各過程が確率的に選択されるMixtures of PPCA（Probabilistic Principal Component Analysis;確率的主成分分析）によってモデル化するものとする。
【００６１】
まず、カメラによる画像の撮影を、三次元空間から二次元平面への線形写像（アフィンカメラ）と仮定する。撮像される物体の形状は時不変の剛体だとして、Ｐ個の特徴点の三次元位置を表す行列ｚを（１０）式、フレームｆの画像上での各特徴点の位置を表す行列ｘf を（１１）式で表すものとする。但し、ｐ番目の特徴点の三次元位置をｚ(p) 、フレームｆにおけるｐ番目の特徴点の二次元位置をｘf(p)とする。
【００６２】
【数５】
アフィンカメラの下で特徴点の二次元位置ｘf(p)の生成過程は、観測ノイズεが加わって（１２）式で表すことができる。
【００６３】
【数６】
但し、Ｗf は射影行列、μf は並進成分を表す。
【００６４】
そして、Ｐ個の特徴点をフレーム数Ｆに渡って観測（特徴点追跡部４から取得）したとすると、（１３）式に示す観測データ行列ｘが得られる。
【００６５】
【数７】
但し、Ｗe は（１４）式で定義され、以下では運動行列と称する。また、ｚe は（１５）式で定義され、以下では形状行列と称する。
【００６６】
【数８】
画像に複数の物体が撮影されている場合は、それぞれの物体に対して運動行列Ｗe が存在するため、物体ｕ（ｕ＝１，２，…，Ｕ）に対応するものをＷe/u で表すものとする。但し、画像上では実際には各特徴点がどの物体に由来するのかわからないため、各物体の特徴点が混在する観測データ行列ｘが測定されることになる。
【００６７】
また、物体ｕに関して測定された特徴点の数をＰu で表すものとすると、特徴点の総数Ｐは、（１６）式で表される。
【００６８】
【数９】
ここで、特徴点ｐがどの物体に由来するかを指定する指示変数ｈ(p) を（１７）式で表すものとする。この指示変数は、仮定する物体の総数と同じＵ個ｈu(p)の要素を持つ。
【００６９】
【数１０】
そして、特徴点ｐについての指示変数ｈ(p) は、特徴点ｐが物体ｕに由来するものであるということを表現する場合、ｕ番目の要素のみが１（ｈu(p)＝１）となり、他の要素は全て０となる。また、指示変数ｈ(p) は、特徴点の数Ｐだけ存在するため、これをまとめたものを指定行列ｈと呼び、（１８）式で表すものとする。
【００７０】
【数１１】
また、全特徴点の内からランダムに１点を選択した時に、その選択した特徴点が物体ｕに属している確率をｇu として、これを全ての物体についてまとめたものを混合比ｇと呼び、（１９）式で表すものとする。
【００７１】
【数１２】
なお、各物体に関して測定されている特徴点数（Ｐ1，…，ＰU）が既知である場合、ｇu はＰ個の特徴点の中からランダムに選択して特定のＰu 個の特徴点のうちのどれかを選ぶ確率なのでｇu＝Ｐu／Ｐとなる。但し、Ｐu は未知なので、ｇも推定する必要がある。
【００７２】
また、物体ｕに関する測定に混入するノイズεは期待値０、分散τ|u-1の正規分布に従うものとする。
【００７３】
特徴点ｐが物体ｕに属すると仮定すると、特徴点ｐのノイズの分布ｐ（εp ）は、（２０）式で表されることになる。但し、ｄ次元ベクトルとして表された確率変数ｘが期待値ｍ、共分散Σの正規分布に従うことをＮ（ｘ｜ｍ，Σ）で示し、その内容は（２１）式で定義されるものとする。
【００７４】
【数１３】
以上の設定から、形状行列ｚ、指示行列ｈ、混合比ｇ、物体毎のアフィンパラメータＷe|u （ｕ＝１，２，…，Ｕ），各物体の特徴点位置の測定の正確さ（ノイズの分散の逆数）τ|uが推定（学習）の対象となる。
【００７５】
これらの関係を図示したグラフィカルモデルを図７に示す。
【００７６】
なお、以下では、式の記述を簡略化するために、（２２）（２３）式に示すモデルパラメータθ，Ｊを用いるものとする。
【００７７】
【数１４】
図７に示すグラフィカルモデルの下では、特徴点（観測データ）ｘが得られれば、形状行列ｚとモデルパラメータθは事後確率ｐ（θ，ｚ，ｈ｜ｘ）を計算することで推定することができる。
【００７８】
また、事後分布（事後確率）ｐ（θ，ｚ，ｈ｜ｘ）はベイズ推定により、（２４）式のように記述される。
【００７９】
【数１５】
この式を計算するために、形状行列ｚ，モデルパラメータθ（即ち、Ｗe|u ，τ|u，ｇ）の事前分布ｐ0 を（２５）〜（２８）式によって設定する。但し、αe|u は、４次元ベクトルであり、ここではαe|u ＝［１，１，１，１］Tが用いられるが、何等かの分布を設けてもよい。
【００８０】
【数１６】
但し、ｗe,i|u はＷe|u のｉ行目に対応した列ベクトルである。また、ガンマ分布Ｇａｍ（ｘ｜ａ，ｂ）は（２９）式、ディリクレ分布Ｄｉｒ（ｘ｜γx ）は（３０）式で定義される。
【００８１】
【数１７】
なお、（２４）式を解析的に計算することは不可能なため、ここでは変分推定法を利用して事後分布を近似的に求める。
【００８２】
＜変分推定法＞
変分推定法では、事後分布ｐ（θ，ｈ，ｚ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）を計算する。但し、ｑ（θ，ｈ，ｚ）＝ｑ（θ）ｑ（ｈ，ｚ）と因子化されると仮定して、試験分布ｑ（θ，ｈ，ｚ）と真の事後分布ｐ（θ，ｈ，ｚ｜ｘ）との間のＫＬ（Kullback-Leibler）疑似距離を最小化するように、試験分布ｑ（θ）ｑ（ｈ，ｚ）を調節する。
【００８３】
ところで、ＫＬ疑似距離を最小にする試験分布を求めるということは、自由エネルギーを最大にする試験分布を求めることに等しいことが知られている。つまり、変分推定法では、この自由エネルギーを最大にする試験分布が最も妥当な推定解であるということになる。
【００８４】
本件の場合、自由エネルギーＦE は（３１）式で表される。また、この場合、自由エネルギーＦE を最大化する試験分布ｑ（θ，ｈ，ｚ）は、事前の同時分布ｐ0（Ｊ） に基づいて、（３２）式及び（３３）式に示す条件を満たすことが知られている。
【００８５】
【数１８】
但し、<>ｑは、分布ｑに基づく期待値を求める操作を意味する。
【００８６】
この条件式（３２）（３３）を展開して得られる式は、推定対象であるθ，ｈ，ｚが互いに依存した形式となっているため、明示的な解を得ることができない。しかし、条件式（３２）（３３）を繰り返し適用することは、自由エネルギーＦE を対象とした勾配法を行うことに等しいため、近似的な解を得ることができる。
【００８７】
具体的には、変分ベイズ法のＥステップでは、試験分布ｑ（θ）を固定した上で、試験分布ｑ（ｈ，ｚ）を求め、変分ベイズ法のＭステップでは、試験分布ｑ（ｈ，ｚ）を固定した上で、試験分布ｑ（θ）を求めるものとする。
【００８８】
ここで、Ｅステップで求める各試験分布を（３４）〜（３６）式のようにパラメータ付けする。
【００８９】
【数１９】
また、Ｍステップで求める各試験分布を（３７）〜（４０）式のようにパラメータ付けする。
【００９０】
【数２０】
（２５）〜（２８）式と、以下に示す（４１）式に基づいて、（３２）（３３）式を展開する。
【００９１】
【数２１】
すると、Ｅステップで求める各試験分布のパラメータには、次の（４２）〜（４６）式に示す条件が課されることになる。但し、Σz|u は物体ｕの推定構造の共分散行列（信頼度）、ｚ|u(p) は特徴点ｐが物体ｕ由来であると仮定したときの推定座標、ｃu(p)は特徴点ｐが物体ｕ由来であるという仮定の信頼度、ｃu は試験分布ｑ（ｈ(p) ）の正規化項を表す。
【００９２】
【数２２】
同様に、Ｍステップで求める各試験分布のパラメータには、次の（４７）〜（５２）式に示す条件が課されることになる。但し、γguの比は混合比の期待値（推定値）、ＭWe|uはＷe|u の期待値（推定値）、Σｗe|u はＷe|u の推定の各行の共分散行列（信頼度）、ｍτ|u はτ|uの期待値（推定値）、γτ|u はτ|uの推定の分散（信頼度）ｍτ|u／γτ|uの構成要素を表す。
【００９３】
【数２３】
なお、Ｅステップで計算される全てのパラメータをまとめた分布パラメータΘS は（５３）式で表され、Ｍステップで計算される全てのパラメータをまとめた分布パラメータΘM は（５４）式で表される。
【００９４】
【数２４】
また、ｈ，ｚに依存する十分統計量ｒを（５５）式で表すものとする。具体的には、（４７）（４８）（４９）（５１）（５２）式のそれぞれにおいて、<<>>fcで囲われた項が十分統計量となる。
【００９５】
【数２５】
＜個人顔モデル学習部での処理＞
ここで、個人顔モデル学習部５２が実行する処理を、図８に示すフローチャートに沿って説明する。なお、本処理は、外部からの起動指令の入力により起動される。
【００９６】
本処理が起動すると、まず、Ｓ１１０では、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、平均顔モデル記憶部５１に記憶された平均顔モデルを選択するように切替信号を出力して、Ｓ１２０に進む。
【００９７】
これにより、推定部６では、画像（フレーム）が入力される毎に、平均顔モデルを用いたアフィンパラメータＡf 、及び特徴点推定分布パラメータｘf，Ｖfの算出が行われる。また、頭部姿勢算出部７では、推定部６にて求められたアフィンパラメータＡf に基づいて、頭部姿勢θｘ，θｙ，θｚの算出が行われる。
【００９８】
Ｓ１２０では、モデルパラメータθの試験分布である第１試験分布ｑ（θ）を初期化してＳ１３０に進む。
【００９９】
具体的には、（５４）式に示した分布パラメータΘM 及び、（５５）式に示した十分統計量ｒ初期化され、例えば、全ての要素が、同じ一定値（例えば、０又は１又は０．５）に設定される。
【０１００】
Ｓ１３０では、特徴点追跡部４から予め設定された規定フレーム数Ｆ個分の特徴点パラメータｘf （（１１）式参照）の取得を完了したか否かを判断し、完了していなければ、完了するまで待機する。
【０１０１】
そして規定フレーム数Ｆ分の特徴点パラメータｘf の取得を完了すると、これらをまとめた観測データ行列ｘ（（１３）式参照）を生成してＳ１４０に移行し、設定されている分布パラメータΘM に基づいて、第２試験分布ｑ（ｚ，ｈ）の算出、即ち、変分ベイズのＥステップとしての処理を実行して、Ｓ１５０へ進む。
【０１０２】
具体的には、（５３）式に示す分布パラメータΘS を、（４２）〜（４６）式に従って算出する。
【０１０３】
Ｓ１５０では、観測データ行列ｘ及びＳ１４０での第２試験分布ｑ（ｚ，ｈ）の算出結果に従って、十分統計量ｒを更新し、続くＳ１６０では、第１試験分布ｑ（θ）の算出、即ち、変分ベイズのＭステップとしての処理を実行して、Ｓ１７０に進む。
【０１０４】
具体的には、（５４）式に示す分布パラメータΘM を、Ｓ１５０にて算出された十分統計量ｒを用い、（４７）〜（５２）式に従って算出する。
【０１０５】
Ｓ１７０では、Ｓ１４０，Ｓ１６０での算出結果に基づき、上述の（３１）式を用いて自由エネルギーＦE を算出してＳ１９０に進む。
【０１０６】
但し、自由エネルギーＦE は、上述の処理によって得られた第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）から求められる運動行列Ｗおよび形状行列ｚに基づいて推定される位置と、観測により得られた位置（観測データ行列ｘ）とが近いほど大きな値となる。
【０１０７】
Ｓ１８０では、Ｓ１７０にて算出された自由エネルギーＦE が、予め設定されたエネルギー閾値より大きいことを終了条件として、その終了条件を充足しているか否かを判断し、充足していなければ、Ｓ１４０に戻って、Ｓ１４０〜Ｓ１７０の処理を繰り返し、一方、終了条件を充足していれば、Ｓ１９０に移行する。
【０１０８】
ここでは、エネルギー閾値として、例えば、平均顔モデルを表す形状行列と、上述の処理で得られた運動行列Ｗとを用いて算出した自由エネルギー、またはその自由エネルギーより大きな値に設定された固定値を用いる。
【０１０９】
Ｓ１９０では、各特徴点ｐ（ｐ＝１，２，…，Ｐ）について、その特徴点が由来する物体ｕを、（５６）式に示す判別式を用いて判別し、その判別結果からドライバの頭部に由来する特徴点（以下「頭部特徴点」という）のみを抽出して特徴点追跡部４に通知する。
【０１１０】
【数２６】
これにより、特徴点追跡部４では、以後、通知された頭部特徴点に基づいて、任意特徴点の中からドライバの頭部に由来するもののみを選択して追跡することが可能となる。
【０１１１】
Ｓ２００では、第２試験分布ｑ（ｈ，ｚ）から求められた形状行列のうち、Ｓ１９０にて得られた頭部特徴点に基づく部分のみを抽出した部分形状行列を生成し、その部分形状行列を個人顔モデルとして出力すると共に、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、個人顔モデルを選択するように切替信号を出力して、本処理を終了する。
【０１１２】
これにより、以後、推定部６では、画像（フレーム）が入力される毎に、個人顔モデルを用いたアフィンパラメータＡf の算出が行われ、更に、そのアフィンパラメータＡf に基づいて、頭部姿勢算出部７では、頭部姿勢θｘ，θｙ，θｚの算出が行われる。
【０１１３】
＜効果＞
以上説明したように、画像処理装置１では、頭部姿勢θｘ，θｙ，θｚを推定する処理を開始した直後は、固定特徴点に基づく平均顔モデルを用いてアフィンパラメータＡf の算出（顔向き推定）を行い、これと並行して、任意特徴点も利用した個人顔モデルの学習を実行し、その個人顔モデルの精度が十分に向上した時点で、平均顔モデルから個人顔モデルへの切替を行うようにされている。
【０１１４】
従って、画像処理装置１によれば、ドライバが誰であっても、また、ドライバが身につける物品（サングラス，マスク，帽子，スカーフ等）によって、固定特徴点を抽出することができない場合であっても、頭部姿勢θｘ，θｙ，θｚの推定を、処理の初期段階から、安定した精度の推定結果を得ることができ、しかも、平均顔モデルから個人顔モデルへの切替後は、高い精度の推定結果を得ることができる。
【０１１５】
また、画像処理装置１では、観測データ行列ｘを因子分解する際に、指定行列ｈを導入して、各特徴点がどの物体に由来するかを同時に推定するようにされている。
【０１１６】
従って、画像処理装置１によれば、観測データ行列ｘに、ドライバ頭部に由来する特徴点（頭部特徴点）だけでなく、ドライバの頭部以外の物体（例えば、車両や窓に映る車外の景色等）に由来する特徴点（非頭部特徴点）が含まれていたとしても、即ち、観測データ行列ｘが混合データ行列であったとしても、その観測データ行列ｘを因子分解した結果から、物体の形状を個別に推定することができ、ひいてはドライバ頭部の形状を精度よく推定することができる。
【０１１７】
また、画像処理装置１では、個人顔モデル学習部５２にて得られた頭部特徴点の情報を、特徴点追跡部４にて供給するようにされているため、特徴点追跡部４では、頭部特徴点の情報に基づいて任意特徴点を選別することができ、その選別結果に基づき、例えば、非頭部特徴点の追跡を中止することによって、不要な追跡処理が実行されることを防止することができる。
【０１１８】
［第２実施形態］
次に、第２実施形態について説明する。
【０１１９】
なお、本実施形態では、個人顔モデル学習部５２での処理内容の一部が、第１実施形態のものとは異なるだけであるため、この相違する部分を中心に説明する。
【０１２０】
詳しくは、第１実施形態では、規定フレーム数Ｆ分の特徴点データが得られた後に、因子分解を一括して行う、いわゆるバッチ手法を用いているが、本実施形態では、１フレーム分の特徴点データを取得する毎に、逐次的にθ，ｈ，ｚを推定する、所謂、オンライン推定法を実行するものである。
【０１２１】
＜オンライン推定法の概要＞
オンライン推定法では、（３１）式に示す自由エネルギーＦE を、（５７）式に示すように、オンライン推定に適した形式のオンライン自由エネルギーＦEλに置き換えて使用する。
【０１２２】
但し、忘却係数をλ（但し、０＜λ＜１）、現フレームをｆc 、フレームｆまでの観測に基づいた試験分布を（５８）式、十分統計量を（５９）式、忘却係数の影響を正規化するための係数を（６０）式で表すものとする。
【０１２３】
【数２７】
（３１）式と（５７）式を比較すると、オンライン自由エネルギーＦEλの十分統計量ｒを、（６１）式に示す十分統計量ｒfのｆ＝１，２，…，ｆｃを通した減衰付き和で置き換えたものとなっている。
【０１２４】
そして、十分統計量ｒfの減衰係数付き和は（６２）式に示すように逐次的に計算することができ、また、正規化項も（６３）式に示すように逐次的に計算することができる。
【０１２５】
【数２８】
＜個人顔モデル学習部での処理＞
図９は、個人顔モデル学習部５２が実行する処理の内容を示すフローチャートである。
【０１２６】
第１実施形態の場合と同様に、本処理は、初期化指令の入力により起動され、Ｓ３１０〜Ｓ３４０は、第１実施形態のＳ１１０〜Ｓ１４０と同様であるため説明を省略する。
【０１２７】
但し、Ｓ１３０にて使用される規定個数は、第１実施形態の場合より少ない数（例えば、数分の一から数十分の一程度）でよく、ここでは１とされている。
【０１２８】
そして、Ｓ３５０では、観測データ行列ｘf 及びＳ３４０での第２試験分布ｑ（ｚ，ｈ）の算出結果に基づき、（６２）（６３）式を用いて、十分統計量<<ｒ>>fcを更新し、続くＳ３６０では、第１試験分布ｑ（θ）の算出、即ち、変分ベイズのＭステップとしての処理を実行して、Ｓ１７０に進む。
【０１２９】
具体的には、（５４）式に示す分布パラメータΘM を、Ｓ３５０にて算出された十分統計量<<ｒ>>fcを、（４７）〜（５２）式の十分統計量の項（<<>>ｆｃで囲われた項）に代入して算出する。
【０１３０】
Ｓ３７０では、Ｓ１９０と同様に頭部特徴点を抽出して特徴点追跡部４に通知し、続くＳ３８０では、（５７）式に従ってオンライン自由エネルギーＦEλを算出し、Ｓ３９０では、Ｓ１８０と同様に、その算出したオンライン自由エネルギーＦEλが、予め設定されたエネルギー閾値より大きいことを終了条件として、その終了条件を充足しているか否かを判断する。
【０１３１】
そして、終了条件を充足していなければ、Ｓ３３０に戻って、Ｓ３３０〜Ｓ３８０の処理を繰り返し、一方、終了条件を充足していれば、Ｓ４００に移行する。
【０１３２】
Ｓ４００では、Ｓ２００と同様に、頭部特徴点に基づく部分形状行列を個人顔モデルとして出力すると共に、顔モデル切替部５３が推定部６に供給する顔モデルＣa として、個人顔モデルを選択するように切替信号を出力して、本処理を終了する。
【０１３３】
＜効果＞
以上説明したように、本実施形態では、１フレーム毎に、Ｅステップ，Ｍステップを１回ずつ実行して第１試験分布ｐ（θ）及び第２試験分布ｐ（ｈ，ｚ）を求める処理を繰り返し、その結果を過去のデータほど比率が低下するように減衰係数ηfcを乗じて足し合わせていく、いわゆるオンライン化した変分ベイズ推定法を用いて、因子分解を行っている。
【０１３４】
従って、本実施形態によれば、使用する画像系列の数（フレーム数）を増加させても、それによって増加する処理量は、Ｅステップ，Ｍステップを１回ずつ実行する際の処理量を単位として、加算的に増大するだけで、指数関数的に増大することがないため、因子分解に要する処理負荷（ひいては処理に要する時間）を大幅に削減することができる。
【０１３５】
その結果、自動車に搭載される演算器は一般的に計算能力が低いが、そのような場合でも問題なく適用することができる。
【０１３６】
更に本実施形態では、Ｍステップで求める第１試験分布を規定するパラメータの算出式中に現れる十分統計量<<ｒ>>fc 及び減衰係数ηfcを、前フレームでの算出結果を利用した簡易な式（（６２）（６３）式参照）によって逐次的求めるようにされているため、因子分解に要する処理負荷を、更に削減することができる。
【０１３７】
また、本実施形態では、十分統計量<<ｒ>>fcが求められれば観測データ行列ｘfcは不要であり、特徴点の増減に対応するために、過去の観測データ行列を記憶しておく必要がないため、装置の記憶容量も削減することができる。
【０１３８】
［実験］
図１０に示すような二つの物体（四角柱，三角柱）を適当に動かして画像に投影したものを入力画像系列として、第２実施形態で説明した処理を実行する画像処理装置１に処理を行わせて形状行列ｚを求めた。なお、特徴点としては、各物体の頂点を用いた。
【０１３９】
図１１は、処理の結果得られた形状行列ｚから推定し二つの物体の形状の推定誤差（ＲＭＳ誤差）の推移を観測フレーム数に対して示したものであり、（ａ）が四角柱についての結果、（ｂ）が三角柱についての結果である。
【０１４０】
図１１からは、フレーム数が少ないと推定精度が低いが、フレーム数の増加に共に推定精度が向上していること、及び、二つの物体を正しく分離できていることがわかる。
【０１４１】
［他の実施形態］
以上本発明の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、様々な態様での実施が可能である。
【０１４２】
例えば、上記実施形態では、頭部姿勢を推定するためにアフィンパラメータを用いているが、これは頭部姿勢の推定を行うことができれば、例えば、アフィンパラメータを用いずに、拡張カルマンフィルタを利用して直接頭部推定を行う手法等を採用してもよい。
【０１４３】
また、因子分解を終了させる際の終了条件として、平均顔モデルから求めた自由エネルギーに基づいて設定されたエネルギー閾値を用いているが、例えば、処理の繰り返し回数が、予め設定された回数閾値より大きいことを終了条件としてもよい。この場合、処理を終了させるか否かの判断のために、自由エネルギーを算出する必要がなくなるため、当該装置での処理負荷を、より一層軽減することができ、また、常に一定時間内で処理を終了させることができる。
【図面の簡単な説明】
【０１４４】
【図１】発明が適用された画像処理装置の全体構成を示すブロック図。
【図２】固定特徴点の説明図。
【図３】任意特徴点の説明図。
【図４】画像処理装置が使用する座標系を示す説明図。
【図５】特徴点追跡部及び推定部が扱う状態空間モデルを表すグラフ。
【図６】時系列ベイズ推定の動作概要を示す説明図。
【図７】個人顔モデル学習部が扱うグラフィックモデルを示す説明図。
【図８】第１実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【図９】第２実施形態における個人顔モデル学習部の処理内容を示すフローチャート。
【図１０】実験に用いた物体の形状を示す説明図。
【図１１】物体の形状と実験により得られた形状との誤差を示すグラフ。
【符号の説明】
【０１４５】
１…画像処理装置 ２…カメラ ３…画像取込装置 ４…特徴点追跡部 ５…顔モデル供給部 ６…推定部 ７…頭部姿勢算出部 ４１…固定特徴点追跡部 ４３…任意特徴点追跡部 ５１…平均顔モデル記憶部 ５２…個人顔モデル学習部 ５３…顔モデル切替部
【特許請求の範囲】
【請求項１】
画像系列を入力し、入力された各画像から予め設定された特徴を有する複数の特徴点を抽出すると共に、前記特徴点のそれぞれについて該特徴点の位置を画像面上で追跡する特徴点抽出追跡手段と、
該特徴点抽出追跡手段での抽出，追跡結果に従って、前記特徴点の位置を表す観測データ行列を生成する観測データ行列生成手段と、
該観測データ行列生成手段にて生成された観測データ行列を因子分解して、前記画像に写された物体の三次元形状を表す形状行列および該物体の運動を表す運動行列を生成する分解手段と、
を備えた三次元形状復元装置において、
前記分解手段は、
前記観測データ行列をｘ、前記形状行列をｚ、前記特徴点のそれぞれが属する物体を指定する指定行列をｈ、前記運動行列，前記特徴点の測定の正確さ，前記特徴点がそれぞれの物体に属している確率をまとめたモデルパラメータをθとして、前記観測データ行列ｘから事後確率Ｐ（θ，ｚ，ｈ｜ｘ）を計算することで、θ，ｚ，ｈを求めることを特徴とする三次元形状復元装置。
【請求項２】
前記分解手段は、
前記事後確率Ｐ（θ，ｚ，ｈ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）が、第１試験分布ｑ（θ）と第２試験分布ｑ（ｈ，ｚ）とに因子化されるものとして、前記事後確率Ｐ（θ，ｚ，ｈ｜ｘ）と前記試験分布ｑ（θ，ｈ，ｚ）とのＫＬ距離を最小にする前記第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）を求めることで、確率的に因子分解計算を実行することを特徴とする請求項１に記載の三次元形状復元装置。
【請求項３】
前記観測データ行列生成手段は、予め設定されたフレーム数分の画像に基づいて前記観測データ行列ｘを生成し、
前記分解手段は、
前記第１試験分布ｑ（θ）を固定して前記第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、
予め設定された終了条件が満たされている場合に、前記推定手段を終了させる終了手段と、
からなり、最終的に得られた前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得ることを特徴とする請求項２に記載の三次元形状復元装置。
【請求項４】
前記観測データ行列生成手段は、１フレーム分の画像毎に前記観測データ行列ｘを生成し、
前記分解手段は、
前記観測データ行列生成手段にて前記観測データ行列ｘが生成される毎に、前記第１試験分布ｑ（θ）を固定して前記第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、
予め設定された終了条件が満たされている場合に、前記推定手段及び更新手段を終了させる終了手段と、
からなり、前記推定手段によって逐次的に得られる前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得ることを特徴とする請求項２に記載の三次元形状復元装置。
【請求項５】
自動車に搭載して使用することを特徴とする請求項１乃至請求項４のいずれかに記載の三次元形状復元装置。
【請求項１】
画像系列を入力し、入力された各画像から予め設定された特徴を有する複数の特徴点を抽出すると共に、前記特徴点のそれぞれについて該特徴点の位置を画像面上で追跡する特徴点抽出追跡手段と、
該特徴点抽出追跡手段での抽出，追跡結果に従って、前記特徴点の位置を表す観測データ行列を生成する観測データ行列生成手段と、
該観測データ行列生成手段にて生成された観測データ行列を因子分解して、前記画像に写された物体の三次元形状を表す形状行列および該物体の運動を表す運動行列を生成する分解手段と、
を備えた三次元形状復元装置において、
前記分解手段は、
前記観測データ行列をｘ、前記形状行列をｚ、前記特徴点のそれぞれが属する物体を指定する指定行列をｈ、前記運動行列，前記特徴点の測定の正確さ，前記特徴点がそれぞれの物体に属している確率をまとめたモデルパラメータをθとして、前記観測データ行列ｘから事後確率Ｐ（θ，ｚ，ｈ｜ｘ）を計算することで、θ，ｚ，ｈを求めることを特徴とする三次元形状復元装置。
【請求項２】
前記分解手段は、
前記事後確率Ｐ（θ，ｚ，ｈ｜ｘ）の近似関数である試験分布ｑ（θ，ｈ，ｚ）が、第１試験分布ｑ（θ）と第２試験分布ｑ（ｈ，ｚ）とに因子化されるものとして、前記事後確率Ｐ（θ，ｚ，ｈ｜ｘ）と前記試験分布ｑ（θ，ｈ，ｚ）とのＫＬ距離を最小にする前記第１試験分布ｑ（θ）及び第２試験分布ｑ（ｈ，ｚ）を求めることで、確率的に因子分解計算を実行することを特徴とする請求項１に記載の三次元形状復元装置。
【請求項３】
前記観測データ行列生成手段は、予め設定されたフレーム数分の画像に基づいて前記観測データ行列ｘを生成し、
前記分解手段は、
前記第１試験分布ｑ（θ）を固定して前記第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて交互に繰り返し実行する推定手段と、
予め設定された終了条件が満たされている場合に、前記推定手段を終了させる終了手段と、
からなり、最終的に得られた前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得ることを特徴とする請求項２に記載の三次元形状復元装置。
【請求項４】
前記観測データ行列生成手段は、１フレーム分の画像毎に前記観測データ行列ｘを生成し、
前記分解手段は、
前記観測データ行列生成手段にて前記観測データ行列ｘが生成される毎に、前記第１試験分布ｑ（θ）を固定して前記第２試験分布ｑ（ｚ，ｈ）を推定する第１のステップ、及び、前記第２試験分布ｑ（ｚ，ｈ）を固定して前記第１試験分布ｑ（θ）を推定する第２のステップを、互いの推定結果を用いて実行する推定手段と、
予め設定された終了条件が満たされている場合に、前記推定手段及び更新手段を終了させる終了手段と、
からなり、前記推定手段によって逐次的に得られる前記第２試験分布ｑ（ｚ，ｈ）から形状行列ｚ及び指定行列ｈを得ることを特徴とする請求項２に記載の三次元形状復元装置。
【請求項５】
自動車に搭載して使用することを特徴とする請求項１乃至請求項４のいずれかに記載の三次元形状復元装置。
【図１】
【図４】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図５】
【図６】
【図４】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図５】
【図６】
【公開番号】特開２０１０−１２３０７０（Ｐ２０１０−１２３０７０Ａ）
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願番号】特願２００８−２９８４７８（Ｐ２００８−２９８４７８）
【出願日】平成２０年１１月２１日（２００８．１１．２１）
【出願人】（５０４１４３４４１）国立大学法人 奈良先端科学技術大学院大学 (226)
【出願人】（０００００４２６０）株式会社デンソー (27,639)
【Ｆターム（参考）】
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願日】平成２０年１１月２１日（２００８．１１．２１）
【出願人】（５０４１４３４４１）国立大学法人 奈良先端科学技術大学院大学 (226)
【出願人】（０００００４２６０）株式会社デンソー (27,639)
【Ｆターム（参考）】
[ Back to top ]