JP3349632B2 - ３次元計測方法及び装置 - Google Patents

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、３次元空間を移動
する物体の動画像からその物体の３次元情報を求めるた
めの３次元計測方法及び装置に関する。

【０００２】３次元空間内を動く物体をカメラで撮影
し、撮影された結果から物体の３次元情報（運動と形
状）を求めることは、ＦＡ（ファクトリー・オートメー
ション）のロボットの目などにとって重要な課題であ
る。また、最近のマルチメディア技術の発達にともなっ
て、人間の動きに適切に反応することで人間と機械との
コミュニケーションを行うシステムが提案されている
が、そのようなシステムを現実のものとするためにも３
次元情報を求めることは重要な課題である。

【０００３】３次元情報を得るための一般的な３次元計
測装置として、レンジセンサが用いられている。しか
し、レンジセンサのような特別な機材は高価であるの
で、廉価なシステムを構築するためのにはＴＶカメラな
どのように通常に用いられているカメラを利用すること
が好ましい。

【０００４】カメラを複数台用いることによって、それ
らの間の視差情報に基づいて３次元情報を算出すること
ができるが、その場合にはカメラを複数台配置しそれら
の位置関係を固定しておかなければならず、操作が煩雑
である。したがって、通常のカメラを１台だけ用いて３
次元情報を算出することが最も望まれるところである。

【０００５】

【従来の技術】従来より、１台のＴＶカメラにより撮影
された動画像に基づいて３次元情報を得る方法が提案さ
れている。１台のＴＶカメラにより撮影される動画像
は、映画のコマのような単位時間当たり一定の枚数の静
止画像からなる。

【０００６】図１０は動画像から３次元情報を算出する
流れを示す図、図１１は特徴点の対応付けに基づく３次
元情報算出の流れを示す図、図１２はオプティカルフロ
ーに基づく３次元情報算出の流れを示す図である。

【０００７】動画像から３次元情報を算出するには、大
きく分類して、図１１に示すような隣接する画像間の特
徴点の対応に基づく方式と、図１２に示すようなオプテ
ィカルフローに基づく方式の２つの方式がある。

【０００８】図１１に示す方式において、特徴点とは、
物体の一定箇所を画像上に投影した点のことである。こ
の方式では、隣接する画像間で特徴点の対応付けを行う
のであるが、その操作又は処理は容易ではない。

【０００９】図１２に示す方式において、オプティカル
フローとは、投影された物体像の２次元的な動きのこと
である。この方式では、動画像からオプティカルフロー
を算出し、算出されたオプティカルフローから３次元情
報を算出するという、２段階の処理に分解して３次元情
報の算出が行われる。

【００１０】従来において、動画像からオプティカルフ
ローを算出する種々の方法が提案されている。例えば、
「時空間に関する滑らかさ制約を用いたオプティカルフ
ローの算出」（遠藤，鳥生，吉田、電子情報通信学会Ｄ
−II，vol.J74-D-II，no.12,pp.1678-1685 ，1991．）
には、勾配法にフローは時間的にも空間的にも滑らかで
あるという制約条件を加えることでオプティカルフロー
を高精度に算出することの可能な方法が提案されてい
る。また、特開平４−１１５３７５号公報には、オプテ
ィカルフローを算出し易い物体境界付近と算出し難い物
体内部とに対して互いに別処理を行うことでオプティカ
ルフロー算出の精度を高める方法が提案されている。さ
らに、特開平６−１７６１４８号公報には、前の時刻に
おけるフロー算出結果を予測値として用いることでフロ
ー算出の精度を高める方法が提案されている。また、オ
プティカルフローを算出する専用のハードウェアも作成
されている（「動画像処理プロセッサＩＳＨＴＡＲによ
るリアルタイム・オプティカル・フロー抽出」佐々木，
塩原、日本ロボット学会誌，vol.13，no.3，pp.331-33
4，1995．）。

【００１１】他方、オプティカルフローから３次元情報
を求める方法も、これまでにいくつか提案されている。
オプティカルフローは、３次元空間内を運動する物体が
投影されて作られる２次元的な動きである。３次元空間
内を運動する物体をカメラの投影面に投影する際に、元
の奥行き情報が欠落するので、与えられたオプティカル
フローから元の３次元情報を算出するためには何らかの
仮定や条件が必要となる。ここでは、対象となる物体を
１つの剛体とし、３次元空間内を運動する１つの剛体を
カメラで撮影する場合について考える。剛体とは３次元
的な形が変化しないものをいう。複数の剛体が存在する
場合には、前処理で剛体を分離することにより、ここで
考える場合に帰着させることができるので、この仮定は
実用上十分である。複数の剛体を分離する方法として、
例えば特開平４−１４８２８３号公報に記載の方法があ
る。

【００１２】この仮定の下では、求めるべき３次元情報
は３次元空間（Ｘ，Ｙ，Ｚ）内における剛体の運動と形
状である。剛体の運動は、剛体の並進運動（ｕ_X 、
ｕ_Y 、ｕ _Z ）と回転運動（ω_X 、ω_Y 、ω_Z ）とを表す
６つのパラメータによって表される。剛体の形状は、剛
体の各観測点までの距離、例えば各観測点の座標で表さ
れる。剛体の各観測点までの距離は、観測点の個数に等
しい個数の未知数で表される。しかし、オプティカルフ
ローから３次元情報を算出する過程において、遠くにあ
って速く動く剛体と近くにあってゆっくり動く剛体とを
区別することができないため、並進速度の絶対的な大き
さを求めることができない。

【００１３】そのため、例えば並進運動が単位速度で行
われていると仮定することによって、つまり例えば、 ｕ_X ² ＋ｕ_Y ² ＋ｕ_Z ² ＝１ と仮定することによって、剛体の２つのパラメータによ
って残りの１つのパラメータが導かれるので、運動を表
すパラメータは５つに減少する。

【００１４】初期の頃に提案された３次元情報の算出方
法では、正確なオプティカルフローが与えられた場合に
は正しい３次元情報を求めることができるが、オプティ
カルフローに雑音が少し乗っただけで３次元情報が求め
られなくなる欠点があった。そのため、現在では雑音に
対する頑強性を増す研究が行われている。その際に、与
えられるオプティカルフローは、実際の３次元物体が投
影されて作られる真のオプティカルフローに白色雑音が
加わったものであるとする統計的なモデルが用いられて
きた（図１３を参照）。

【００１５】なお、雑音の原因としては、カメラのレン
ズの歪み、ビデオ信号をデジタル化する際の量子化誤
差、物体の移動による輝度（濃度）の変化などである。
図１３に示すモデルは、動画像からオプティカルフロー
を算出する従来の方法で求められたオプティカルフロー
には必ずしも厳密な意味では当てはまらないが、近似的
には当てはまると考えられている。したがって、図１３
に示すモデルを仮定して、オプティカルフローから３次
元情報を算出する改良された方法を提案することは、技
術的にも実際的にも意義のあることである。

【００１６】図１３に示すモデルの下では、算出される
３次元情報（以下「推定量」という）の偏り（真の値か
らの系統的なずれ）や分散（算出される値のばらつき）
を計算することができるので、偏りがなくかつ分散が小
さい推定量を与える３次元情報の算出方法ほど精度の良
い算出方法であるという基準が得られる。ただし、オプ
ティカルフローからの３次元情報算出においては、各観
測点における剛体までの距離を精度良く算出することは
原理的に不可能であるので、分散を小さくできるのは剛
体の運動を表す５つのパラメータに関してだけである。
剛体までの距離の算出精度を改善するためには、他の方
法を併用する必要がある。本明細書においてはこれには
触れないが、例えば、「フローの時系列と運動方程式を
用いた３次元復元手法」（為清 彰，鳥生 隆，遠藤利
生，人工知能学会全国大会，pp.601-604，1992．）及び
特開平５−３４６３０８号公報には、時系列で得られた
情報を統合する方法が提案されている。また、特開平５
−１８７８４４号公報には、色情報等を併合することで
精度を高める方法が提案されている。

【００１７】偏りがなく分散が小さい推定量を与える推
定方法ほど良い推定方法であるという基準に基づいて、
オプティカルフローからの３次元情報の算出方法の改良
が行われてきた。現在では、あまり分散は小さくならな
いが少ない演算量（行列の固有値計算）で３次元情報を
求めることができる線形法と、分散が小さくなる最尤推
定法（オプティカルフローの誤差に関する最小２乗法）
が良く用いられている。

【００１８】

【発明が解決しようとする課題】多くの統計的な問題で
は、最尤推定法は、最適な推定量つまり偏りがなく分散
が最小の推定量を与えるので、オプティカルフローから
の３次元情報の算出問題でも最適になるのではないかと
予想されていたが、本発明者らによって偏りがなくかつ
より分散が小さくなる別の推定方法が提案された（「 A
superior estimator to the maximum likelihood est
imator on 3-D motion estimation from noisy optical
flow 」，T.Endoh, T.Toriu and N.Tagawa, IEICE Tra
ns. Info.& Syst., vol.E77-D, no.11, 1994. 以下
「提案文献」という）。

【００１９】しかし、この提案文献による推定方法の優
位性は、オプティカルフローに加わる雑音の分散があま
り小さくない場合にしか成り立たず、オプティカルフロ
ーに加わる雑音の大きさによらず最尤推定法よりも優れ
ている推定方法が存在するかどうかは未解決の課題であ
った。

【００２０】なお、従来からベイズ推定法という特定の
事前知識を前提とする推定方法が用いられてきたが、こ
れはその事前知識（例えば物体の形状が特定の形である
など）が正しい場合には最尤推定法よりも優れているも
のの、事前知識が間違っている場合には一般には推定量
が偏ってしまい、良い推定方法とはならないという問題
があった。

【００２１】本発明は、上述の問題に鑑みてなされたも
ので、オプティカルフローに加わる雑音の大きさにかか
わらず最尤推定法よりも精度良く安定に３次元情報を算
出することを可能とする３次元計測方法及び装置を提供
することを目的とする。

【００２２】

【課題を解決するための手段】請求項１の発明に係る方
法は、３次元空間を移動する物体の動画像からオプティ
カルフローを求め、求めたオプティカルフローから前記
物体の３次元情報を求めるための３次元計測方法であっ
て、１つのオプティカルフローに基づいて初期推定を行
い、３次元情報の初期推定量を求める第１のステップ
と、前記第１のステップで求められた初期推定量に基づ
いて推定関数の重みを定める第２のステップと、前記１
つのオプティカルフローに基づいて、前記第２のステッ
プで定められた前記重みを用いて前記推定関数が小さく
なるように推定を行い、３次元情報の推定量を求める第
３のステップと、前記第３のステップで求められた推定
量を３次元情報として出力する第４のステップと、を有
する。

【００２３】請求項２の発明に係る方法は、前記第１の
ステップにおいて、最尤推定法を用いてオプティカルフ
ローの２乗誤差を最小にするように前記初期推定を行
う。請求項３の発明に係る方法は、前記第３のステップ
において、仮の３次元情報に基づいて仮のオプティカル
フローを求め、求められた仮のオプティカルフローと前
記１つのオプティカルフローとについての誤差及び前記
重みを含む推定関数の値を求め、前記推定関数の値が小
さくなるように前記仮の３次元情報を修正し、これを再
帰的に繰り返すことによって前記推定を行う。

【００２４】請求項４の発明に係る方法では、前記推定
関数は、前記仮のオプティカルフローと前記１つのオプ
ティカルフローとの２乗誤差をパラメータで微分して得
られる項に、前記仮のオプティカルフローと前記１つの
オプティカルフローとの誤差を前記重みで平均して得ら
れる修正項を付け加えたものである。

【００２５】請求項５の発明に係る方法は、前記第２の
ステップにおいて、前記第１のステップで求められた初
期推定量における物体までの距離の平均と分散とを用い
て前記重みを定める方法である。

【００２６】請求項６の発明に係る装置は、１つのオプ
ティカルフローに基づいて初期推定を行い、３次元情報
の初期推定量を求める初期推定部と、前記初期推定部で
求められた初期推定量に基づいて重みを算出する重み算
出部と、前記１つのオプティカルフローに基づいて且つ
前記重み算出部で算出された前記重みを用いて推定を行
い、３次元情報の推定量を求める重み付き推定部と、前
記重み付き推定部で求められた推定量を３次元情報とし
て出力する出力部と、を有して構成される。

【００２７】請求項７の発明に係る装置は、入力される
オプティカルフローを格納する格納部と、格納部から読
み出された１つのオプティカルフローに基づいて初期推
定を行い、３次元情報の初期推定量を求める初期推定部
と、前記初期推定部で求められた初期推定量に基づいて
重みを算出する重み算出部と、前記１つのオプティカル
フローに基づいて且つ前記重み算出部で算出された前記
重みを用いて推定を行い、３次元情報の推定量を求める
重み付き推定部と、前記重み付き推定部で求められた推
定量を３次元情報として出力する出力部と、を有して構
成される。

【００２８】第１のステップは、例えば線形法又は最尤
推定法などによって行われる。推定関数の重みＣを定め
る第２のステップは、例えば後述の（２９）式における
特定の重みＣ^* を定めるステップである。３次元情報の
推定量を求める第３のステップは、Ｃ^* を用いて後述の
（１８）式で示される推定関数Ｋ（θ）＝０を解いて並
進運動と回転運動とを表すパラメータθを求めるステッ
プである。第４のステップは求めたθを３次元情報とし
て出力するステップである。

【００２９】重みＣは、例えば物体の各観測点までの距
離ｒの平均μとその分散Σに依存して定められる。分散
Σは物体の幅（奥行き長さ）の２乗平均に対応する。し
たがって、重みＣは物体の幅に応じて定められるともい
える。

【００３０】第３のステップを実行する重み付き推定部
は、従来の最尤推定法において２乗誤差を最小にすると
ころを、推定関数Ｋ（θ）の値を０とするように置き換
えたものといえる。推定関数Ｋ（θ）は、例えば（１
８）式で表されるように、オプティカルフローの２乗誤
差を微分したものに、オプティカルフローの２乗誤差を
重みＣで平均した修正項を付け加えたものである。重み
Ｃを０とすると従来の最尤推定法と同じになる。重み付
き推定部は、２乗誤差を最小として得られる推定量を、
物体の幅に応じて分散が最小になるように重みＣを定め
て算出し直すものであるといえる。

【００３１】

【発明の実施の形態】

〔３次元計測装置の構成〕図１は本発明に係る３次元計
測装置１のブロック図、図２は初期推定部１４における
処理手順を示す図、図３は重み付き推定部１６における
処理手順を示す図、図４は３次元空間を移動する物体の
投影モデルを示す図、図５は３次元計測装置１の処理の
概要を示すフローチャート、図６はオプティカルフロー
ＰＦの例を示す図、図７は２乗誤差ＥＲを算出する過程
を示す図、図８は３次元情報ＤＴの例を示す図である。
本実施形態において、対象とする物体は３次元的な形が
変化しない剛体であるとする。

【００３２】３次元計測装置１は、画像格納部１１、オ
プティカルフロー算出部１２、オプティカルフロー格納
部１３、初期推定部１４、重み算出部１５、重み付き推
定部１６、及び出力部１７などからなる。

【００３３】画像格納部１１は、ビデオカメラＶＣによ
って撮影された動画像ＦＭを入力して格納する。動画像
ＦＭは、例えば１秒当たり３０フレームの静止画像の集
合である。しかし、動画像ＦＭはこれに限られず、例え
ば２フレームのみの静止画像であってもよい。

【００３４】オプティカルフロー算出部１２は、入力さ
れた動画像ＦＭに基づいてオプティカルフローＰＦを算
出する。オプティカルフロー算出部１２の構成は従来よ
り公知のものが用いられる。つまり、従来の技術の項で
述べた公知文献に記載された方法などによって、動画像
ＦＭからオプティカルフローＰＦが算出される。算出さ
れたオプティカルフローＰＦは、オプティカルフロー格
納部１３に格納される。なお、他の処理装置によって作
成されたオプティカルフローＰＦを入力し、それをオプ
ティカルフロー格納部１３に格納してもよい。

【００３５】オプティカルフローＰＦは、時間的に隣接
するフレーム（画像）間において、物体の各観測点に対
応する点のＸ方向及びＹ方向の変化量（又は移動量）に
よって表される。図６の例では、各点のＸ方向の変化量
はＸ成分ＰＦｘによって、Ｙ方向の変化量はＹ成分ＰＦ
ｙによって表される。Ｘ成分ＰＦｘ及びＹ成分ＰＦｙに
は、３行３列の各点の変化量が示されている。ここで
は、これらの各点は、画像の各画素に一致するものとす
る。

【００３６】図６において、例えば、第１行第１列の画
素は、Ｘ方向に「３」、Ｙ方向に「０」の変化量を持
つ。第１行第２列の画素は、Ｘ方向に「２」、Ｙ方向に
「−１」の変化量を持つ。

【００３７】初期推定部１４は、オプティカルフローＰ
Ｆに基づいて初期推定を行い、３次元情報の初期推定量
ＤＴｉを算出する。初期推定部１４における初期推定の
処理は、例えば従来より公知の方法である最尤推定法に
よって行われる。最尤推定法については後述する。

【００３８】重み算出部１５は、初期推定部１４で算出
された初期推定量ＤＴｉに基づいて、重みＣを算出す
る。重みＣは、例えば物体の各観測点までの距離ｒの平
均μとその分散Σに依存して定められる。各観測点まで
の距離をｒ_1,ｒ_2,…ｒ_n とし、観測点数（画素数）をｎ
とすると、その平均μは、（ｒ₁ ＋ｒ₂ ＋…＋ｒ_n ）／
ｎであり、分散Σは、〔（ｒ₁ −μ）² ＋（ｒ₂ −μ）
² ＋…＋（ｒ_n −μ）²〕／ｎである。分散Σは物体の
幅（奥行き長さ）の２乗平均に対応する。したがって、
重みＣは物体の幅に応じて定められるともいえる。

【００３９】重み付き推定部１６は、オプティカルフロ
ーＰＦに基づいて、且つ重み算出部１５で算出された重
みＣを用いて推定を行い、３次元情報の推定量ＤＴｐを
算出する。出力部１７は、重み付き推定部１６で算出さ
れた推定量ＤＴｐを３次元情報ＤＴとして出力する。

【００４０】３次元計測装置１におけるオプティカルフ
ローＰＦから３次元情報ＤＴを出力する処理の流れの概
略を纏めると、次のようになる。すなわち、図５に示す
ように、初期推定部１４がオプティカルフロー格納部１
３に格納されたオプティカルフローＰＦを読み込み（＃
１）、３次元情報の初期推定量ＤＴｉを算出する（＃
２）。算出された初期推定量ＤＴｉに基づいて、重み算
出部１５が重みＣを算出する（＃３）。算出された重み
Ｃを用い、重み付き推定部１６が３次元情報の推定量Ｄ
Ｔｐを算出し（＃４）、出力部１７がその推定量ＤＴｐ
を３次元情報ＤＴとして出力する（＃５）。 〔３次元情報の算出〕次に、オプティカルフローＰＦに
基づいて３次元情報の推定量ＤＴｐを算出する方法につ
いて詳しく説明する。

【００４１】オプティカルフローＰＦが算出された後に
おいて、そのオプティカルフローＰＦに対して、まず従
来の方法による初期推定によって物体までの距離の平均
（距離の逆数の重み付き平均）と分散とを求め、それら
から定まる重みＣを用いて、オプティカルフローＰＦか
ら３次元情報ＤＴを算出し直す。重みＣを適切に設定す
ることにより、最尤推定法よりも算出精度が良くなる。
ベイズ推定法と異なり、初期推定で得られた推定結果が
仮に間違っていたとしても、偏りが生じないため、安定
して３次元情報ＤＴを算出することができる。 〔最尤推定法〕まず、初期推定部１４での処理、つまり
最尤推定法に基づく最尤推定量の計算手順を図２を参照
して説明する。

【００４２】最尤推定量とは、与えられたオプティカル
フローＰＦと、３次元空間内を運動する物体の投影像と
の２乗誤差が最小になるような物体の３次元情報のこと
である。

【００４３】図２において、適当な仮の３次元情報ＤＴ
ｔが初期値として与えられる（＃２１）。仮の３次元情
報ＤＴｔとしては、例えば、直前のフレームについて算
出された３次元情報、又は総てのパラメータが０である
３次元情報などが用いられる。

【００４４】与えられた３次元情報ＤＴｔについて、そ
のような３次元情報ＤＴｔを持つ物体の投影結果である
仮のオプティカルフローＰＦｔを計算する（＃２２）。
仮のオプティカルフローＰＦｔとオプティカルフロー格
納部１３に格納されているオプティカルフローＰＦとの
２乗誤差を計算する（＃２３）。具体的には、例えば図
７に示すように、仮のオプティカルフローＰＦｔとオプ
ティカルフローＰＦとについて、Ｘ成分ＰＦｘ及びＹ成
分ＰＦｙのそれぞれについて各画素毎の変化量の差を求
め、それぞれの差を２乗したものを合計し、その結果を
２乗誤差ＥＲとする。

【００４５】２乗誤差ＥＲを評価し、２乗誤差ＥＲが十
分小さければ、つまり例えば所定のしきい値よりも小さ
ければ（＃２４でイエス）、そのときに３次元情報ＤＴ
ｔとして格納されている情報を最尤推定量として出力す
る（＃２６）。

【００４６】２乗誤差ＥＲが小さくなければ（＃２４で
ノー）、２乗誤差ＥＲが小さくなるように、そのときに
格納されている３次元情報ＤＴｔを微小変化させ、３次
元情報ＤＴｔを更新する（＃２５）。そのための手法と
して、例えば３次元情報ＤＴｔをランダムに複数回にわ
たり微小変化させ、その中で２乗誤差ＥＲが最も小さく
なるものを選択する。更新された３次元情報ＤＴｔにつ
いて、２乗誤差ＥＲが十分に小さくなるまで、さらに更
新を繰り返す。

【００４７】統計学の定理によると、観測点が十分に多
く且つ推定すべき母数（未知数）が観測点に比べて十分
に少なければ、一般に、最尤推定量は最適になる、つま
り最尤推定量は偏りがなく分散が最小となる。しかし、
オプティカルフローＰＦから３次元情報ＤＴを算出する
場合に、物体の運動は並進運動と回転運動とを表す５つ
のパラメータで表現できるものの、物体の形状は観測点
毎の物体までの距離で表されるため、未知数の個数は
「観測点数＋５」となり、観測点数に比べて十分少ない
とは言えず、したがってこの性質は成り立たない。

【００４８】しかし、仮に物体の形状が事前に分かって
いれば、推定すべき母数は物体の運動を表す５つのパラ
メータだけになるため、最尤推定量が最適になる。形状
がいくつかのパラメータで表されること、例えば、半径
が未知の球であるとか、２次式の曲面であるとかを知っ
ている場合も同様である。さらに、物体の形状が少数の
パラメータで表現される特定の確率分布にしたがって変
化すると知っている場合も同様である。物体形状が特定
の確率分布にしたがって変化すると知っている場合の最
尤推定量をベイズ推定量と呼称する。ベイズ推定量は、
物体の形状に関する事前知識が正しい限り、常に、何の
事前知識も持たない最尤推定量よりも小さい分散を持つ
が、事前知識が間違っている場合には偏りが大きくなる
欠点がある。 〔重み付き推定〕最尤推定量よりも優れた推定量を算出
するために、本実施形態における重み付き推定部１６に
おいては、最尤推定量が２乗誤差ＥＲを最小にするとこ
ろを、以下に述べる推定関数Ｋ（θ）に置き換える。つ
まり、推定関数Ｋ（θ）の値が「０」となるような推定
量（３次元情報）ＤＴｐを出力することで、新しい推定
量を定義する。

【００４９】ここに用いる推定関数Ｋ（θ）は、具体的
にはオプティカルフローＰＦの２乗誤差ＥＲをパラメー
タで微分した式に、オプティカルフローＰＦの２乗誤差
ＥＲを可変な重みＣで平均した修正項を付け加えたもの
である。この場合に、２乗誤差ＥＲの微分は式の上で予
め行っておくことができるため、推定関数Ｋ（θ）の計
算は２乗誤差ＥＲの計算と同程度の演算量で行うことが
できる。重み付き推定部１６における推定量ＤＴｐの計
算手順は図３に示されている。

【００５０】図３において、適当な仮の３次元情報ＤＴ
ｔが初期値として与えられる（＃３１）。仮の３次元情
報ＤＴｔとしては、例えば、初期推定部１４で算出され
た初期推定量ＤＴｉ、直前のフレームについて算出され
た３次元情報、又は総てのパラメータが「０」である３
次元情報などが用いられる。

【００５１】与えられた３次元情報ＤＴｔについて、そ
のような３次元情報ＤＴｔを持つ物体の投影結果である
仮のオプティカルフローＰＦｔを計算する（＃３２）。
仮のオプティカルフローＰＦｔとオプティカルフロー格
納部１３に格納されているオプティカルフローＰＦとに
よって定まる推定関数Ｋ（θ）の値を計算する（＃３
３）。

【００５２】推定関数Ｋ（θ）の値が十分小さければ、
つまり例えば所定のしきい値よりも小さければ（＃３４
でイエス）、そのときに３次元情報ＤＴｔとして格納さ
れている情報を推定量ＤＴｐとして出力する（＃３
６）。

【００５３】推定関数Ｋ（θ）の値が小さくなければ
（＃３４でノー）、推定関数Ｋ（θ）の値が小さくなる
ように、そのときに格納されている３次元情報ＤＴｔを
微小変化させ、３次元情報ＤＴｔを更新する（＃３
５）。そのための手法として、例えば３次元情報ＤＴｔ
をランダムに複数回にわたり微小変化させ、その中で推
定関数Ｋ（θ）の値が最も小さくなるものを選択する。
更新された３次元情報ＤＴｔについて、推定関数Ｋ
（θ）の値が十分に小さくなるまでさらに更新を繰り返
す。

【００５４】さて、２乗誤差ＥＲをパラメータで微分し
た式を０とするような３次元情報と、２乗誤差ＥＲを最
小にする３次元情報とは一致するので、重みＣが０の場
合には、＃３６の推定量ＤＴｐは＃２６の最尤推定量と
一致する。したがって、初期推定において３次元情報を
推定し、それに応じて重みＣを適切に定めれば、得られ
る推定量ＤＴｐが最尤推定量よりも小さい分散を持つこ
とが期待される。実際、ある重みＣに対してそれから得
られる推定量の分散が、他の重みＣに対する推定量の中
では最小になることを示すことができる。この最適な重
みＣは、物体までの距離の逆数の重み付き平均と分散と
によって表される。オプティカルフローＰＦからの３次
元情報算出においては、物体までの距離自体は精度良く
算出できないが、その平均及び分散は精度良く算出する
ことができるので、この最適な重みＣは初期推定から決
定可能である。したがって、推定された最適な重みＣを
持つ推定関数Ｋ（θ）が「０」となるような３次元情報
ＤＴを計算することによって、最尤推定量よりも分散が
小さくなる推定量ＤＴｐが算出できる。

【００５５】さらに、この推定関数Ｋ（θ）に関して
は、どんな重みＣに対してもそれから定まる推定量ＤＴ
ｐが不偏になることを示せるので、初期推定で多少の誤
差が入っても安定して３次元情報ＤＴを算出することが
可能となる。このような推定関数Ｋ（θ）と最適な重み
Ｃの詳細は後で述べる。 〔推定関数〕次に、最尤推定量よりも分散が小さくなる
推定量ＤＴｐを導くための推定関数Ｋ（θ）の構成と、
その最適な重みＣの計算方法について説明する。説明
は、投影モデルと記号の定義、雑音の統計的モデルと最
尤推定量、推定関数Ｋ（θ）の定義と最適な重みＣの導
出について、この順で行う。 〔投影モデル〕３次元空間において１つの物体が運動し
ており、それがある投影面に中心投影されているものと
する。この投影像から元の物体の運動を推定する問題を
考える。物体上の１つの観測点の３次元位置ベクトルを
Ｘとする。物体の並進速度ベクトルをｕ、回転速度ベク
トルをωとすると、物体の運動方程式は次の（１）式で
示される。

【００５６】Ｘ’＝ｕ＋ω×Ｘ ……（１） 但し、Ｘ’は時間微分、×は外積を表す 観測者が原点にいる座標系を考え、物体をある平面に中
心投影して観測するものとする。投影面上の観測点の３
次元位置ベクトル〔投影面がＺ＝１の平面の場合なら
ば、（ｘ，ｙ，１）^t 〕をη、その点における投影面の
外向き単位法線ベクトル〔投影面がＺ＝１の平面の場合
ならば（０，０，１）^t 〕をζとする（｜ζ｜＝１）。
上添字ｔは転置行列を表す。但し、投影面の原点からの
距離（焦点距離）は１に正規化されていると仮定する。
この投影関係は次の（２）式で表される。

【００５７】 η＝νＸ ν＝１／（ζ・Ｘ） ……（２） 但し、・は内積を表す νは法線方向に測った観測点Ｘまでの距離の逆数であ
る。投影面がＺ＝１の平面の場合の投影モデルが図４に
示されている。

【００５８】ηは観測点が移動するにつれて時間ととも
に変化するが、η・ζ＝１が常に成り立つことに注意す
る。したがって、ベクトル三重積の公式である次の
（３）式、 （ｘ×ｙ）×ｚ＝（ｘ・ｚ）ｙ−（ｙ・ｚ）ｘ ……（３） より、ｐを任意のベクトルとして、次の（４）式が成り
立つ。

【００５９】 ｐ−（ｐ・ζ）η＝ζ×（ｐ×η） ……（４） （２）式を（１）式に代入し、ｐ＝νｕ＋ω×ηに対し
て（４）式を適用することにより、次の（５）式で表さ
れる復元方程式を得る。

【００６０】 ξ＝η’＝ζ×［ν（ｕ×η）＋（ω×η）×η］ ……（５） 上の（５）式の左辺ξは、投影面上の点の時間微分η’
であり、観測されるオプティカルフローを表す。右辺
は、物体の運動パラメータｕ，ωと距離の逆数νの関数
になる。また、投影面上の観測点は複数（ｎ個）存在す
るので、それらをｋの添字（ｋ＝１〜ｎ）で区別し、ｋ
番目の観測点の３次元位置ベクトルをη_k、法線ベクト
ルをζ_k 、その点での物体までの距離の逆数をν_k 、そ
の点でのオプティカルフローをξ_k と表す。したがっ
て、ｋ番目の点における復元方程式は、次の（６）式で
表される。

【００６１】 ξ_k ＝ζ_k ×［ν_k （ｕ×η_k ）＋（ω×η_k ）×η_k ］ ……（６） 並進速度ｕと物体までの距離の逆数ν_k との間には定数
倍の不定性があるので、ｕ≠０の場合には｜ｕ｜＝１と
置いてこの不定性を除く。ｕ＝０の場合は別途これを判
定するものとし、以下｜ｕ｜＝１とする。この仮定によ
り、ｕの各成分は独立ではなくなるので、単位球面｜ｕ
｜＝１を表す２次元ベクトルをφとし、φを推定すべき
未知数と見なす。φとωを合わせてθと書く。θ＝（φ
^t ，ω^t）^t である。したがって、オプティカルフロー
からの３次元情報算出問題における未知数の総数は、物
体の並進速度の角度φを指定するために２個、回転速度
ωに３個、各観測点における物体までの距離の逆数ν_k
がｎ個あり、合計「ｎ＋５」個である。観測量はオプテ
ィカルフローξ_k であり、各観測点毎に投影面の法線ベ
クトルζ_k に直交する２成分が与えられるので、合計
「２ｎ」である。したがって、正確なオプティカルフロ
ーＰＦが与えられた場合には、一般には観測点が５点以
上あれば解は定まる。 〔雑音の統計的モデルと最尤推定量〕推定量の偏り及び
分散について議論するためには、オプティカルフローＰ
Ｆに加わる雑音のモデルを定める必要がある。多くの研
究では、観測点毎、成分毎に独立に平均０の正規分布に
したがう雑音がオプティカルフローＰＦに対して加法的
に加わるモデルを採用している。ここでもこのモデルを
採用する。すなわち、ｋ番目の観測点におけるオプティ
カルフローξ_k には、以下の仮定にしたがう雑音が加わ
るとする。 観測点毎の独立性 各観測点毎に独立である。 成分毎の独立性 投影面の法線ベクトルζ_k に直交する２成分には、独立
に平均０分散σ² の雑音が加わるとし、ζ_k に平行な成
分には雑音が加わらないとする（または分散０の雑音が
加わるとする）。この仮定は、雑音によりオプティカル
フローが投影面から飛び出さないことを意味する。

【００６２】ここで、Ｅ［ξ_k ］＝ξ _k ⁽⁰⁾ とする。上
の仮定より、次の（７）式で示す等式が成り立つ。ただ
しＥ［Ｘ］は確率変数Ｘの期待値を表す。 ξ_k ・ζ_k ＝０ ……（７） ξ _k ⁽⁰⁾ ＝ζ_k ×［ν_k （ｕ×η_k ） ＋（ω×η_k ）×η_k ］ ……（８） Ｖ［ξ_k ］＝σ² Ｎ_k 但し、Ｎ_k ＝Ｉ₃ −ζ_k ζ _k ^t ……（９） 但し、Ｖ［ξ_k ］はξ_k の共分散行列を表す。特に投影
面がＺ＝１の平面の場合には、Ｖ［ξ_k ］＝σ² diag
（１，１，０）になる。

【００６３】雑音が含まれているオプティカルフローか
ら３次元情報を算出する場合に、オプティカルフローの
平均からのへだたりを各観測点毎に２乗して加えたもの
をＪ ^(MLE) とし、これを最小にするようなｕ，ω，ν_K
を推定値とする推定方法が考えられる。上で述べた雑音
のモデルの下では、Ｊ^(MLE) は確率密度関数ｆ_k と次の
（１０）式の関係がある。

【００６４】 Π_k=1,n ｆ_k ( ξ_k ） ＝ const・ exp［−Ｊ^(MLE) ( ν_1,ν_2,…，ν_n ，ｕ，ω） ／2 σ² ］……（１０） したがって、Ｊ^(MLE) を最小にする３次元情報は、確率
密度関数の値を最大にする、すなわち最尤推定量とな
る。ここで、Ｊ^(MLE) は次の（１１）式で表される。

【００６５】 Ｊ^(MLE) ( ν_1,ν_2,…，ν_n ，ｕ，ω） ＝Σ_k=1,n ｜ζ_k ×｛（ν_k ｕ＋ω×η_k ）×η_k ｝−ξ_k ｜² ……（１１） 以下においては、Ｊ^(MLE) をν_k に対して最小化するこ
とで運動パラメータのみの関数に変形することを試み
る。

【００６６】Ｊ^(MLE) をν_k に対して偏微分して０と置
くと、次の（１２）式が得られる。 ν_k ＝［ξ_k ・ｕ_k −［ζ_k ×｛（ω×η_k ）×η_k ｝］・ｕ_k ］ ／｜ζ_k ×（ｕ×η_k ）｜² ……（１２） となる。ただし、ｕ_k = ζ_k ×（ｕ×η_k ）である。つ
まり、物体までの距離の逆数ν_k の推定値はその点での
オプティカルフローξ_k の観測値と運動パラメータｕ，
ωの推定値だけで定まる。したがって、観測点数ｎが無
限大になり運動パラメータの推定値が真の値に収束する
ような場合でも、各ν_k の推定値にはσ程度の誤差が含
まれることが分かる。このように、オプティカルフロー
から任意形状の物体の３次元情報を算出する問題では、
物体形状の算出精度は原理的に高めることができない。
そこで、以下では観測点数の増加によって算出精度を上
げることが可能な運動パラメータｕ，ωの推定方法の検
討を行うこととし、物体形状を表すν_k に関しては例え
ば時系列データの統合などの別の手段でその算出精度を
上げるという方針を選択する。ここでは後者の問題は扱
わないが、時系列データとして与えられたオプティカル
フローを用いて精度良く物体形状を推定する方法とし
て、例えば運動方程式を利用するものがある（前掲「フ
ローの時系列と運動方程式を用いた３次元復元手
法」）。

【００６７】（１２）式を（１１）式に代入すると、距
離の逆数ν_k を消去した評価関数が次の（１３）式のよ
うに導かれる。 Ｊ^(MLE) ( ｕ，ω） ＝Σ_k=1,n ｜（ｕ×η_k ）・｛（ξ_k −（ω×η_k ）｝｜² ／｜ζ_k ×（ｕ×η_k ）｜² ……（１３） この（１３）式を、真のフローξ _k ⁽⁰⁾ に加わる雑音Ｘ
_k を使って書き直すと次の（１４）式の形に表すことが
できる。

【００６８】 Ｊ^(MLE) ( θ) ＝Σ_k=1,n ｜ａ^(k)t( θ) Ｘ_k ＋ｂ^(k) （θ）｜² ……（１４） 但し、ａ^(k) は単位ベクトル、ｂ^(k) はスカラであり、
いずれも運動パラメータθの関数である。雑音がないと
きの真の運動パラメータは誤差０になるので、ｂ
^(k) （θ⁽⁰⁾ ）＝０が成り立つ。但し、（θ⁽⁰⁾ ）は真
の値である。Ｘ_k は平均０、共分散行列σ² Ｎ_k の正規
分布に独立にしたがう。また、ａ^(k)t( θ) Ｎ _k ａ^(k)
(θ) ＝１も成り立つ。 〔推定関数の定義と最適な重みの導出〕一般に、確率変
数Ｘ_k の観測値からｍ次元の母数ベクトルθを推定する
問題を考える。オプティカルフローからの３次元情報算
出問題の場合には、ｍ＝５である。その最尤推定量θ
^(MLE) が最小にする評価関数Ｊ^(MLE) が次式で表される
とする。

【００６９】 Ｊ^(MLE) （θ）＝Σ_k=1,n ｜ａ^(k)t( θ) Ｘ_k ＋ｂ^(k) （θ）｜² ……（１５） 但し、Ｘ_k は平均０、共分散行列σ² Ｎ_k の正規分布に
独立にしたがうとする。また、θの真の値をθ⁽⁰⁾ と
し、ｂ^(k) （θ⁽⁰⁾ ）＝０とする。ここで、 Ｊ_k （θ）＝ａ^(k)t( θ) Ｘ_k ＋ｂ^(k) （θ） ……（１６） と置き、Ｊ（θ）＝（Ｊ₁ （θ）_, Ｊ₂ （θ）_, …Ｊ_n
（θ））^t とする。すると、 Ｊ^(MLE) （θ) ＝｜Ｊ（θ）｜² ……（１７） と表すことができる。

【００７０】さて、ここでｍ次元の推定関数Ｋ（θ）を
次の（１８）式で定義する。

【００７１】

【数１】

【００７２】但し、下付きのθはθによる微分を表す。
Ｃはｍ×ｌの行列、Ｓはｎ×ｌの行列とする。これらは
後から定める行列である。Ｓの列ベクトルは線形独立と
する。この推定関数Ｋ（θ）を０に等しいと置いて得ら
れるθに関する方程式の解を推定量θ^(C) とする。これ
は、確率変数Ｘ_k の関数であるから、やはり確率変数と
なる。Ｅ［Ｋ（θ⁽⁰⁾ ）］＝０であるから、θ^(C) は漸
近的に不偏となる（前掲の提案文献を参照）。もし、Ｃ
＝０であるならば、

【００７３】

【数２】

【００７４】であるから、θ^(C) は最尤推定量Ｊ^(MLE)
に一致する。次に、推定量θ^(C) の漸近分散を計算す
る。漸近分散とは観測点の個数ｎが大きいときに良く一
致する分散の近似であり、計算が簡単であるという利点
を持つ。Ｖ［θ^(C) ］について、次の（１９）式が成り
立つ（前掲の提案文献を参照）。

【００７５】 Ｖ［θ^(C) ］＝σ² （Ｂ^t Ｂ＋ＣＳ^t Ｂ）^-1 ［σ² Σ_k=1,n Ａ_k ＋（Ｂ^t ＋ＣＳ^t ）（Ｂ＋ＳＣ^t ）］ （Ｂ^t Ｂ＋Ｂ^t ＳＣ^t ）^-1 ＝σ² （Ｂ^t Ｂ）^-1＋σ² （Ｂ^t Ｂ＋ＣＳ^t Ｂ）^-1 （σ² Σ_k=1,n Ａ_k ＋ＣＳ^t Ｑ_B ＳＣ^t ） （Ｂ^t Ｂ＋Ｂ^t ＳＣ^t ）^-1 ……（１９） 但し、

【００７６】

【数３】

【００７７】また、Ｑ_B ＝Ｉ_n −Ｂ（Ｂ^t Ｂ）^-1Ｂ^t で
ある。右辺第２項は半正定値になるので、この漸近共分
散行列は任意の不偏推定量の下界を表すクラメル・ラオ
の下界σ²(Ｂ^t Ｂ）^-1よりも小さくはならないことが分
かる。ただし、２つの行列の大小関係は差が半正定値行
列になるかどうかで定めるものとする。

【００７８】ここで、従来方式を用いた初期推定で、Ｂ
^t Ｂ、Ｂ^t Ｓ、Σ_k Ａ_k 、σ² の値が精度良く推定でき
たと仮定する。これらはｎによらない個数の要素からな
るので、観測点数ｎに対してＯ_p （１／ｎ^1/2 ）の精度
で推定できることが期待される。ただし、Ｏ_p はｎが大
きいときの確率変数の大きさの程度を表す記号で、Ｙ＝
Ｏ_p （１／ｎ^1/2 ）とはｎ^1/2 ｜Ｙ｜がＮより大きくな
る確率はＮとｎが大きくなるにつれて０に近づくことを
意味する。

【００７９】すると、以下で示すようにＶ［θ^(C) ］を
最も小さくするような重みＣを求めることが可能にな
る。まず適当に記号を置き換えて、 Ｖ［θ^(C) ］＝σ² （Ｂ^t Ｂ）^-1 ＋σ² （Ｂ^t Ｂ）^-1( Ｉ_m ＋ＣＨ^t ）^-1（Ａ＋ＣＬＣ^t ） （Ｉ_m ＋ＨＣ^t ）^-1（Ｂ^t Ｂ）^-1 ……（２２） と表す。ただし、Ａ＝σ² Σ_k Ａ_k 、Ｈ＝（Ｂ^t Ｂ）^-1
Ｂ^t Ｓ、Ｌ＝Ｓ^t Ｑ_B Ｓである。

【００８０】これより、Ｖ［θ^(C) ］のＣに関する最小
化問題は、一般にｍ×ｍの半正定値対称行列Ａ、ｌ×ｌ
の半正定値対称行列Ｌとｍ×ｌの行列Ｈが与えられたと
き、ｍ×ｌの行列Ｘに対して、 Ｆ（Ｘ）＝（Ｉ_m ＋ＸＨ_t ）^-1（Ａ＋ＸＬＸ^t ）（Ｉ_m ＋ＨＸ^t ）^-1 ……（２３） を最小化する問題に帰着される。この問題を解くため
に、まずＦ（Ｘ）をＸの（ｉ，ｊ）成分で微分したＦ
_i,j （Ｘ）を０と置いて、最小値の候補を求めることを
考える。微分の公式（Ｙ^-1）' ＝−Ｙ^-1Ｙ' Ｙ^-1より、 Ｆ_i,j （Ｘ）＝（Ｉ_m ＋ＸＨ^t ）^-1（Ｅ_i,j Ｙ^t ＋ＹＥ _i,j ^t ） （Ｉ_m ＋ＨＸ^t ）^-1 ……（２４） と表される。ただし、Ｙ＝ＸＬ−（Ａ＋ＸＬＸ^t ）（Ｉ
_m ＋ＨＸ^t ）^-1Ｈであり、Ｅ_i,j は（ｉ，ｊ）成分のみ
１の行列である。すべての（ｉ，ｊ）の組に対してＥ
_i,j Ｙ^t ＋ＹＥ _i,j ^t ＝０を満たすＹは０に限られるの
で、微分が０となるＸはＸＬ−（Ａ＋ＸＬＸ^t ）（Ｉ_m
＋ＨＸ^t ）^-1Ｈ＝０を満たす。ここで、Ｚ＝（Ｉ_m ＋Ｘ
Ｈ^t ）^-1（Ａ＋ＸＬＸ^t ）と変数変換すればＺ＝Ａが解
になる。したがって、Ｌが正則ならば Ｘ^* ＝ＡＨＬ^-1 ……（２５） が解である。ここではＬは正則と仮定しておく。Ｌ＝Ｓ
^t Ｑ_B Ｓが正則という条件は、元のＶ［θ^(C) ］のＣに
関する最小化問題においてはＳの列ベクトルの線形結合
がＢの列ベクトルの線形結合では表せないという条件と
等価になる。このときのＦの値は Ｆ（Ｘ^* ）＝（Ｉ_m ＋ＡＨＬ^-1Ｈ^t ）^-1（Ａ＋ＡＨＬ^-1Ｈ^t Ａ） （Ｉ_m ＋ＨＬ^-1Ｈ^t Ａ）^-1 ＝Ａ−ＡＨ（Ｌ＋Ｈ^t ＡＨ）^-1Ｈ^t Ａ ……（２６） と表される。

【００８１】次に、このＦ（Ｘ^* ）がＦ（Ｘ）の最小値
を与えることを証明する。Ｌ’＝ＨＬ^-1Ｈ^t と置くと、 Ｆ（Ｘ^* ＋Ｘ）−Ｆ（Ｘ^* ） ＝（Ｉ_m ＋ＡＬ’＋ＸＨ^t ）^-1 （Ａ＋ＸＬＸ^t ＋ＡＬ’Ａ＋ＡＨＸ^t ＋ＸＨ^t Ａ） （Ｉ_m ＋Ｌ’A ＋ＨＸ^t ）^-1 −（Ｉ_m ＋ＡＬ’）^-1（Ａ＋ＡＬ’Ａ）（Ｉ_m ＋Ｌ’Ａ）^-1 ＝（Ｉ_m ＋ＡＬ’+ ＸＨ^t ）^-1Ｘ ［Ｌ−Ｈ^t （Ｉ_m ＋ＡＬ’）^-1（Ａ＋ＡＬ’Ａ） （Ｉ_m ＋Ｌ’Ａ）^-1Ｈ］Ｘ^t （Ｉ_m ＋Ｌ’Ａ＋ＨＸ^t ）^-1 ……（２７） と計算できる。ここで、 Ｌ−Ｈ^t （Ｉ_m ＋ＡＬ’）^-1（Ａ＋ＡＬ’Ａ）（Ｉ_m ＋Ｌ’Ａ）^-1Ｈ ＝Ｌ−Ｈ^t Ｆ（Ｘ^* ）Ｈ ＝Ｌ−Ｈ^t ［Ａ−ＡＨ（Ｌ＋Ｈ^t ＡＨ）^-1Ｈ^t Ａ］Ｈ ＝Ｌ（Ｌ＋Ｈ^t ＡＨ）^-1Ｌ≧０ ……（２８） を用いると、Ｆ（Ｘ^* ＋Ｘ）−Ｆ（Ｘ^* ）≧０がしたが
う。つまり左辺は半正定値行列になる。

【００８２】元のＶ［θ^(C) ］のＣに関する最小化問題
において、Ｖ［θ^(C) ］の最小値を与えるＣ＝Ｃ^* とそ
の最小値を表すと、 Ｃ^* ＝σ² （Σ_k Ａ_k ）（Ｂ^t Ｂ）^-1Ｂ^t Ｓ（Ｓ^t Ｑ_B Ｓ）^-1 ……（２９） Ｖ［θ^(C) ］ ＝σ² （Ｂ^t Ｂ）^-1＋σ⁴ （Ｂ^t Ｂ）^-1Σ_k Ａ_k （Ｂ^t Ｂ）^-1 −σ⁶ （Ｂ^t Ｂ）^-1（Σ_k Ａ_k ）（Ｂ^t Ｂ）^-1Ｂ^t Ｓ ［Ｓ^t Ｑ_B Ｓ＋σ² Ｓ^t Ｂ（Ｂ^t Ｂ）^-1（Σ_k Ａ_k ） （Ｂ^t Ｂ）^-1Ｂ^t Ｓ］^-1Ｓ^t Ｂ（Ｂ^t Ｂ）^-1 （Σ_k Ａ_k ）（Ｂ^t Ｂ）^-1 ……（３０） となる。最尤推定量の漸近共分散行列は上式の第１項と
第２項の和に等しいので、Ｃ＝Ｃ^* のときの推定量θ
^(C) は最尤推定量よりも常に小さい漸近共分散行列を持
つことが分かる。

【００８３】また、ｌ≧ｍで且つＬ’＝ＨＬ^-1Ｈ^t の逆
行列が存在する場合には、 Ｆ（Ｘ^* ）＝Ｌ’^-1−Ｌ’^-1（Ａ＋Ｌ’^-1）^-1Ｌ’^-1 ……（３１） が成り立つ。Ｐ_s ＝Ｓ（Ｓ^t Ｓ）^-1Ｓ^t に対して次の関
係式、 （Ｂ^t Ｂ）^-1＋｛Ｂ^t Ｓ（Ｓ^t Ｑ_B Ｓ) ^-1Ｓ^t Ｂ｝^-1 ＝（Ｂ^t Ｐ_s Ｂ）^-1 ……（３２） が成り立つことを用いると、ｌ≧ｍにおいて、 Ｖ［θ^(C) ］ ＝σ² （Ｂ^t Ｐ_s Ｂ）^-1−σ² ｛（Ｂ^t Ｐ_s Ｂ）^-1−（Ｂ^t Ｂ）^-1｝ ｛（Ｂ^t Ｐ_s Ｂ）^-1−（Ｂ^t Ｂ）^-1＋σ² （Ｂ^t Ｂ）^-1 Σ_k Ａ_k （Ｂ^t Ｂ）^-1｝^-1 ｛（Ｂ^t Ｐ_s Ｂ）^-1−（Ｂ^t Ｂ）^-1｝ ……（３３） が成り立つ。オプティカルフローに加わる雑音が余り小
さくないときに、最尤推定量よりも分散が小さくなる一
般化された重みＣ＝Ｐ_F を持つ評価関数Ｊ^(C) を最小に
する推定量θ^(C) の漸近共分散行列は、σ² （Ｂ^t Ｐ_F
Ｂ）^-1で近似できるので、ここで定式化した推定量でＳ
＝Ｆとした推定量θ^(C) はそれと同じか又はより小さい
漸近共分散行列を持つことが分かる。ただし、Ｆは座標
η_k の２次式からなるｎ×６（ｌ＝６）の行列である。

【００８４】σ、Ｂ^t Ｂ、Ｂ^t Ｓ、Σ_k Ａ_k の初期推定
の方法に関してまとめておく。最尤推定法等で運動パラ
メータを推定すると相対誤差Ｏ_p （１／ｎ^1/2 ）の精度
で算出できると考えられる。これから距離の逆数ν_k を
最尤推定法で推定する。すなわち、ｕ＾などを運動パラ
メータの推定値として、 ν＾_k ＝（ξ_k −ζ_k ×ω＾_k ）・ｕ＾_k ／｜ｕ＾_k ｜² ……（３４） とする。但し、ω_k ＝（ω×η_k ）×η_k であり、ω＾
_k はその推定値である。このとき、 ν＾_k ＝ν_k ⁽⁰⁾ ＋（ｕ_k ⁽⁰⁾ ・δξ_k ／｜ｕ＾_k ｜² ） ＋Ｏ_p （１／ｎ^1/2 ） ……（３５） が成り立つ。ここで、δξ_k ＝ξ_k −ξ_k ⁽⁰⁾ とした。
Ｂ^t Ｓはδξ_k の１次式になるのでＯ_p （１／ｎ^1/2 ）
の相対誤差で推定できる。Σ_k Ａ_k は運動パラメータの
みの式でありν_k ⁽⁰⁾ を含まないので、Ｏ_p （１／ｎ
^1/2 ）の相対誤差で計算できる。σも同程度の精度で推
定できる。Ｂ^t Ｂはδξ_k の２次式になるため偏るが、
その偏りは推定できるのでやはり同程度の精度で推定可
能である。したがって、重みＣも同程度の精度で計算可
能である。故に漸近共分散行列の主要項（１／ｎの係
数）は真のＢ^t Ｂなどを知っている場合と同じ値になる
ので、初期推定で得られた値を基に計算された推定量も
常に最尤推定量より小さい共分散行列を持つことが分か
る。

【００８５】上述の説明において、（１８）式の右辺の
第１項は、（１７）式の右辺を微分したものであり、第
２項は、Ｊ（θ）に対して重みＣを係数とした修正項で
ある。つまり、（１８）式で示される推定関数Ｋ（θ）
は、（１７）式の右辺で示される２乗誤差ＥＲの微分
と、オプティカルフローの各画素の値を重みＣに応じて
混ぜ合わせたものとからなる。画素の値の混ぜ合わせ方
は、ＣＳ^t による。Ｃはｍ×ｌの行列であり、Ｓはｎ×
ｌの行列である。ｍは未知数の個数、例えば「５」であ
る。ｌは例えば５〜１０の値である。ｎは画素数であ
り、例えば６４０×４８０である。Ｃ及びＳの各行列に
おいて、ｌ個の重み列は互いに異なる数列である。この
推定関数Ｋ（θ）の値を０又は最小とするようなθつま
りｕとωが求められる。そのために、重みＣとして、
（２９）式に示される特別な重みであるＣ^* が用いられ
る。（２９）式において、Ｂ^t Ｂは距離の分散に関連
し、Ｂ^t Ｓは距離の平均に関連する。ここでＢは、（２
０）式に示されているように距離を並べたものであると
いえる。求められたθすなわち推定量ＤＴｐは、図８に
その例を示すように５つのパラメータ値であり、形状情
報ＤＴｆとともに３次元情報ＤＴとして出力される。

【００８６】したがって、本発明における重みを定める
第２のステップは、（２９）式におけるＣ^* を定めるス
テップに対応し、３次元情報の推定量を求める第３のス
テップは、Ｃ^* を用いて（１８）式で示される推定関数
Ｋ（θ）＝０を解いてθを求めるステップに対応し、第
４のステップは求めたθを出力するステップに対応す
る。なお、本発明の第１のステップは、例えば線形法又
は最尤推定法などによって行われる。

【００８７】なお、上述した３次元計測装置１は、適当
なプログラムをロードした制御用のコンピュータ、パー
ソナルコンピュータ、又はワークステーションなどを用
いて実現することができる。 〔実験〕上述した３次元計測装置１による推定方法の有
用性を検証するため、簡単な計算機実験を行った。

【００８８】図９は、オプティカルフローに加わる雑音
の標準偏差（分散の平方根）を横軸とし、算出される並
進速度の標準偏差を縦軸としてこれらの関係を示したも
のである。図９において、umleは最尤推定量、upolw は
本発明者が提案文献において先に提案した推定量、umod
は本実施形態の推定関数Ｋ（θ）に基づく推定量を表
す。オプティカルフローに加わる雑音の大きさによら
ず、本実施形態によるumodは、umle及びupolw よりも常
に小さい分散を持つことが分かる。

【００８９】

【発明の効果】請求項１乃至請求項７の発明によると、
オプティカルフローに加わる雑音の大きさにかかわらず
最尤推定法よりも分散が小さくなり、精度の良い３次元
情報を安定に算出することができる。

【００９０】また、ベイズ推定法のように特定の事前知
識を前提とするものではないので、推定量が偏ってしま
うということがなく、精度のよい３次元情報を安定して
算出できる。

【００９１】このように本発明によると、物体の運動を
精確に計測することができるので、例えば本発明をロボ
ットの目などに適用した場合にロボットが物体を掴み損
ねる可能性が小さくなるなど、種々の制御を正確に行う
ことができる。

【図面の簡単な説明】

【図１】本発明に係る３次元計測装置のブロック図であ
る。

【図２】初期推定部における処理手順を示す図である。

【図３】重み付き推定部における処理手順を示す図であ
る。

【図４】３次元空間を移動する物体の投影モデルを示す
図である。

【図５】３次元計測装置の処理の概要を示すフローチャ
ートである。

【図６】オプティカルフローの例を示す図である。

【図７】２乗誤差を算出する過程を示す図である。

【図８】３次元情報の例を示す図である。

【図９】オプティカルフローに加わる雑音と並進速度の
標準偏差との関係を示すグラフである。

【図１０】動画像から３次元情報を算出する流れを示す
図である。

【図１１】特徴点の対応付けに基づく３次元情報算出の
流れを示す図である。

【図１２】オプティカルフローに基づく３次元情報算出
の流れを示す図である。

【図１３】オプティカルフローに加わる雑音の統計的な
モデルを示す図である。

【符号の説明】

１ ３次元計測装置 １３ オプティカルフロー格納部（格納部） １４ 初期推定部 １５ 重み算出部 １６ 重み付き推定部 １７ 出力部 Ｋ（θ） 推定関数 Ｃ 重み ＤＴｉ 初期推定量 ＤＴｐ 推定量 ＰＦ オプティカルフロー ＰＦｔ 仮のオプティカルフロー ＤＴ ３次元情報 ＤＴｔ 仮の３次元情報 ＥＲ ２乗誤差

Claims (7)

(57)【特許請求の範囲】
1. 【請求項１】３次元空間を移動する物体の動画像からオ
   プティカルフローを求め、求めたオプティカルフローか
   ら前記物体の３次元情報を求めるための３次元計測方法
   であって、１つの オプティカルフローに基づいて初期推定を行い、
   ３次元情報の初期推定量を求める第１のステップと、 前記第１のステップで求められた初期推定量に基づいて
   推定関数の重みを定める第２のステップと、 前記１つのオプティカルフローに基づいて、前記第２の
   ステップで定められた前記重みを用いて前記推定関数が
   小さくなるように推定を行い、３次元情報の推定量を求
   める第３のステップと、 前記第３のステップで求められた推定量を３次元情報と
   して出力する第４のステップと、 を有してなることを特徴とする３次元計測方法。
2. 【請求項２】前記第１のステップにおいて、 最尤推定法を用いてオプティカルフローの２乗誤差を最
   小にするように前記初期推定を行う、 請求項１記載の３次元計測方法。
3. 【請求項３】前記第３のステップにおいて、 仮の３次元情報に基づいて仮のオプティカルフローを求
   め、求められた仮のオプティカルフローと前記１つのオ
   プティカルフローとについての誤差及び前記重みを含む
   推定関数の値を求め、前記推定関数の値が小さくなるよ
   うに前記仮の３次元情報を修正し、これを再帰的に繰り
   返すことによって前記推定を行う、 請求項１又は請求項２記載の３次元計測方法。
4. 【請求項４】前記推定関数は、 前記仮のオプティカルフローと前記１つのオプティカル
   フローとの２乗誤差をパラメータで微分して得られる項
   に、前記仮のオプティカルフローと前記１つのオプティ
   カルフローとの誤差を前記重みで平均して得られる修正
   項を付け加えたものである、 請求項３記載の３次元計測方法。
5. 【請求項５】前記第２のステップにおいて、 前記第１のステップで求められた初期推定量における物
   体までの距離の平均と分散とを用いて前記重みを定め
   る、 請求項１乃至請求項４のいずれかに記載の３次元計測方
   法。
6. 【請求項６】３次元空間を移動する物体の動画像からオ
   プティカルフローを求め、求めたオプティカルフローか
   ら前記物体の３次元情報を求めるための３次元計測装置
   であって、１つ のオプティカルフローに基づいて初期推定を行い、
   ３次元情報の初期推定量を求める初期推定部と、 前記初期推定部で求められた初期推定量に基づいて重み
   を算出する重み算出部と、 前記１つのオプティカルフローに基づいて且つ前記重み
   算出部で算出された前記重みを用いて推定を行い、３次
   元情報の推定量を求める重み付き推定部と、 前記重み付き推定部で求められた推定量を３次元情報と
   して出力する出力部と、 を有してなることを特徴とする３次元計測装置。
7. 【請求項７】オプティカルフローから物体の３次元情報
   を求めるための３次元計測装置であって、 入力されるオプティカルフローを格納する格納部と、 格納部から読み出された１つのオプティカルフローに基
   づいて初期推定を行い、３次元情報の初期推定量を求め
   る初期推定部と、 前記初期推定部で求められた初期推定量に基づいて重み
   を算出する重み算出部と、 前記１つのオプティカルフローに基づいて且つ前記重み
   算出部で算出された前記重みを用いて推定を行い、３次
   元情報の推定量を求める重み付き推定部と、 前記重み付き推定部で求められた推定量を３次元情報と
   して出力する出力部と、 を有してなることを特徴とする３次元計測装置。