(57)【特許請求の範囲】
【請求項１】
複数画像の画像列から得られた画像データから、３次元シーン構造及びカメラの動きを復元する方法であって、
前記画像列の基準画像が基準透視点にあるカメラにより撮像され、前記画像列の１以上の相次ぐ画像が前記カメラを平行移動及び／または回転することによる１以上の相次ぐ異なる透視点において撮像され、
前記方法は、
（ａ）前記各相次ぐ画像に対して、前記基準画像に対する、前記基準透視点から前記相次ぐ異なる透視点までの前記カメラの平行移動及び／または回転から導出される画像データ・シフトを決定する段階と、
（ｂ）各々の画像に対する前記画像データ・シフトを含むシフト・データ行列を構成する段階と、
（ｃ）特異値分解（ＳＶＤ）を使用して、前記シフト・データ行列から、一方のランク（階数）３の因子行列が３次元構造に対応する３次元構造行列であり他方のランク３の因子行列が前記カメラの動きに対応する動き行列である、２つのランク３の因子行列を計算する段階と、
（ｄ）１次方程式を解くことにより、前記３次元構造行列から前記３次元構造を復元する段階と、
（ｅ）前記復元された３次元構造を使用して、前記動き行列から前記カメラの動きを復元する段階と、
を有し、
前記シフト・データ行列を構成する段階は、
Ｈ及び
【外１】
0003651590.tif 000075
を計算することを含み、ここでＨは（Ｎtot−３）×Ｎtot行列であり、
【数１】
0003651590.tif 000076
であり、Ｎtot≡２Ｎp＋２ＮL＋ＮXであって、Ｎp，ＮL及びＮXは点、線及び強度の数にそれぞれ等しく、ＣはＣii'≡δii'＋１である（ＮI−１）×（ＮI−１）行列であり、ここでＣii'はすべての画像に対する定数表の値であり、δii'はｉ＝ｉ'ならば１に等しくｉ≠ｉ'ならば０に等しい演算子であって、
【数２】
0003651590.tif 000077
であり、ＨTはＨ行列の転置であって、ここでＨは、ＨTが恒等行列でありかつ
【外２】
0003651590.tif 000078
を零化する(annihilate)ように定義された行列であって、
【数３】
0003651590.tif 000079
であり、ｙ及びｚ成分に対しても同様であり、ωI及びωLは利用者が設定する定数の重みであり、ここで
【数４】
0003651590.tif 000080
であり、ここでＡlは基準画像において前記線と前記カメラの中心を含む平面への単位法線であり、ここで
【外３】
0003651590.tif 000081
はｘ，ｙ及びｚ方向内の単位ベクトルであり、ここで
【数５】
0003651590.tif 000082
であって、ｑ＝(ｘ,ｙ)は前記基準画像内の追跡点の前記位置であり、
３つの点回転フローｒ(1)(ｘ,ｙ)，ｒ(2)(ｘ,ｙ)，ｒ(3)(ｘ,ｙ)は、
【数６】
0003651590.tif 000083
により定義され、ここでΨIx≡−[∇Ｉ・ｒ(1)(ｐ)]，ΨIy≡−[∇Ｉ・ｒ(2)(ｐ)]，ΨIz≡−[∇Ｉ・ｒ(3)(ｐ)]であり、またここでｐはピクセル位置の画像座標を与え、
ここでΛは、各行が異なる画像ｉに対応しかつ［[ＰU・δＡi]T[ＰL・δＡi]T］に等しいＮI×２ＮL行列であり、ＰU及びＰLはそれぞれ上部及び下部方向と呼ぶ２つの方向
【外４】
0003651590.tif 000084
に射影する単位３−ベクトルであり、δＡiは線フローδＡil≡Ａil−Ａ0lであり、
【外５】
0003651590.tif 000085
は前記ｚ方向内の単位ベクトルであり、ここでＳは、各行が異なる画像ｉに対応しかつ［[ｓix]T[ｓiy]T］に等しいＮi×２Ｎp行列であり、ｓim≡ｑim−ｑ0mであってｍ番目の追跡点に対する画像移動を表し、ΔはＮI×ＮX行列であって各行は画像ｉに対応し[ΔＩi]Tに等しく、ΔＩは前記基準画像に対する画像強度の変化であり、Ｉiはｉ番目の強度画像を表し、Ｉin＝Ｉi(ｐn)はＩi内のｎ番目のピクセル位置における画像強度を表し、表記｛Ｖ｝はＶaにより与えられる成分を有するベクトルを表すために使用され、
前記２つのランク３の因子行列を計算する段階は、
ＳＶＤを使用して
【数７】
0003651590.tif 000086
の最良のランク３の因数分解を計算することを含み、Ｍ(3)，Ｓ(3)は、それぞれ動き及び構造に対応するランク３行列であり、
前記３次元構造を復元する段階は、
１次方程式
【数８】
0003651590.tif 000087
を使用して、Ｕ及びΩに関する３Ｎtot線形制約条件を得るために、
【外６】
0003651590.tif 000088
から構造未知数Ｑz，Ｂz及びＺ-1を除外することを含み、ここでＵ及びΩは未知の３×３行列であり、
【外７】
0003651590.tif 000089
はそれぞれ総合平行移動フロー・ベクトル及び総合回転フロー・ベクトルを表し、ＳＶＤを使用してΟ(Ｎtot)計算量でこれらの制約条件を解き、ここで前記総合平行移動フロー・ベクトルは
【数９】
0003651590.tif 000090
により定義され、ｙ及びｚ成分に対しても同様に定義され、
【数１０】
0003651590.tif 000091
であり、ここでＡは、射影中心及び基準画像内の描かれた線を通る第１の平面の法線であり、第２の平面の平面の法線であるＢは、３次元線Ｌ上の任意の点Ｑに対して、Ｂ・Ａ＝０及びＢ・Ｑ＝−１であることを条件とする(require)ことにより定義され、ＢU≡Ｂ・ＰU及びＢLＬ≡Ｂ・ＰLはＢの上部及び下部成分であり、ここで
【数１１】
0003651590.tif 000092
であり、ここでΦIx≡−[Ｚ-1Ｉx]，ΦIy≡−[Ｚ-1Ｉy]，ΦIz≡[Ｚ-1(∇Ｉ・ｐ)]であり、
ここで、Ｑは基準画像内の画像ピクセルに対応する３次元追跡点に対する３次元座標であり、Ｚnはカメラからカメラの光軸に沿ったｎ番目のピクセルに描かれた３次元の点までの深度であり、
Ｕ及びΩを仮定して、
【数１２】
0003651590.tif 000093
から、前記構造未知数Ｑz，Ｂz及びＺ-1を復元することを含み、
前記カメラの動きを復元する段階は、
前記平行移動を復元するために、
【数１３】
0003651590.tif 000094
を使用することと、
【数１４】
0003651590.tif 000095
から、回転ωix，ωiy，ωizを復元することを含み、
ここでＣはＣii'≡δii'＋１である定数（ＮI−１）×（ＮI−１）行列であり、Ｔは前記平行移動を表す、３次元シーン構造及びカメラの動きを復元する方法。
【請求項２】
前記画像データは、点、線及び強度からなる群から選択される１以上のものである請求項１に記載の方法。
【請求項３】
前記画像データ・シフトを決定する段階は、最初にカメラの回転を復元し補償することを含む、請求項１に記載の方法。
【請求項４】
複数画像の画像列から得られた画像データから、３次元シーン構造及びカメラの動きを復元する方法であって、
前記画像列の基準画像が基準透視点にあるカメラにより撮像され、前記画像列の１以上の相次ぐ画像が前記カメラを平行移動及び／または回転することによる１以上の相次ぐ異なる透視点において撮像され、前記方法は、
（ａ）あらかじめ推定された平行移動変位により基準画像を歪ませた(warp)後に、最初にカメラの回転を復元し補償する段階と、
（ｂ）前記相次ぐ各画像に対して、前記基準画像に対する、前記基準透視点から前記相次ぐ異なる透視点までの前記カメラの平行移動及び／または回転から導出される画像データ・シフトを決定する段階と、
（ｃ）各々の画像に対する前記画像データ・シフトを含むシフト・データ行列を構成する段階と、
（ｄ）前記画像シフトにδＡil，ｓim，ΔＩinを乗算し、点に対して分母因子(denominator factor)１−Ｑ-1zＴzを、強度に対して分母因子１−Ｚ-1Ｔzを、線に対して分母因子１−Ｂ・Ｔを含むようにこれらの量を再定義することにより、前記画像シフトを修正する段階であって、ここでｓim≡ｑim−ｑ0mであり、ｓimはｍ番目の追跡点に対する画像変位を表し、ΔＩinは画像強度における変化であり、Ｉiはｉ番目の強度画像を表し、Ｉin＝Ｉi(Ｐn)はＩi内のｎ番目のピクセル位置における画像強度を表し、δＡilは線フローを表す、前記画像シフトを修正する段階と、
（ｅ）ＳＶＤ（特異値分解）を使用して、前記修正されたシフト・データ行列から、一方のランク（階数）３の因子行列が３次元構造に対応する３次元構造行列であり他方のランク３の因子行列が前記カメラの動きに対応する動き行列である、２つのランク３の因子行列を計算する段階と、
（ｆ）１次方程式を解くことにより、前記３次元構造行列から前記３次元構造を復元する段階と、
（ｇ）前記復元された３次元構造を使用して、前記カメラの動き行列から前記カメラの動きを復元する段階と、
を有し、
前記シフト・データ行列を構成する段階は、
Ｈ及び
【外８】
0003651590.tif 000096
を計算することを含み、ここでＨは（Ｎtot−３）×Ｎtot行列であり、
【数１５】
0003651590.tif 000097
であり、Ｎtot≡２Ｎp＋２ＮL＋ＮXであって、Ｎp，ＮL及びＮXは点、線及び強度の数にそれぞれ等しく、ＣはＣii'≡δii'＋１である（ＮI−１）×（ＮI−１）行列であり、ここでＣii'はすべての画像に対する定数表の値であり、δii'はｉ＝ｉ'ならば１に等しくｉ≠ｉ'ならば０に等しい演算子であって、
【数１６】
0003651590.tif 000098
であり、ＨTはＨ行列の転置であって、ここでＨは、ＨTが恒等行列でありかつ
【外９】
0003651590.tif 000099
を零化(annihilate)するように定義された行列であって、
【数１７】
0003651590.tif 000100
であり、ｙ及びｚ成分に対しても同様であり、ωI及びωLは利用者が設定する定数の重みであり、ここで
【数１８】
0003651590.tif 000101
であり、ここでＡlは基準画像において前記線と前記カメラの中心を含む平面への単位法線であり、ここで
【外１０】
0003651590.tif 000102
はｘ，ｙ及びｚ方向内の単位ベクトルであり、ここで
【数１９】
0003651590.tif 000103
であって、ｑ＝(ｘ,ｙ)は前記基準画像内の追跡点の画像位置であり、
３つの点回転フローｒ(1)(ｘ,ｙ)，ｒ(2)(ｘ,ｙ)，ｒ(3)(ｘ,ｙ)は、
【数２０】
0003651590.tif 000104
により定義され、ここでΨIx≡−[∇Ｉ・ｒ(1)(ｐ)]，ΨIy≡−[∇Ｉ・ｒ(2)(ｐ)]，ΨIz≡−[∇Ｉ・ｒ(3)(ｐ)]であり、またここでｐはピクセル位置の画像座標を与え、
ここでΛは、各行が異なる画像ｉに対応しかつ［[ＰU・δＡi]T[ＰL・δＡi]T］に等しいＮI×２ＮL行列であり、ＰU及びＰLはそれぞれ上部及び下部方向と呼ぶ２つの方向
【外１１】
0003651590.tif 000105
に射影する単位３−ベクトルであり、δＡiは線フローδＡil≡Ａil−Ａ0lであり、
【外１２】
0003651590.tif 000106
は前記ｚ方向内の単位ベクトルであり、ここでＳは、各行が異なる画像ｉに対応しかつ［[ｓix]T[ｓiy]T］に等しいＮi×２Ｎp行列であり、ｓim≡ｑim−ｑ0mであってｍ番目の追跡点に対する画像移動を表し、ΔはＮI×ＮX行列であって各行は画像ｉに対応し[ΔＩi]Tに等しく、ΔＩは前記基準画像に対する画像強度の変化であり、Ｉiはｉ番目の強度画像を表し、Ｉin＝Ｉi(ｐn)はＩi内のｎ番目のピクセル位置における画像強度を表し、表記｛Ｖ｝はＶaにより与えられる成分を有するベクトルを表すために使用され、
前記２つのランク３の因子行列を計算する段階は、
ＳＶＤを使用して
【数２１】
0003651590.tif 000107
の最良のランク３の因数分解を計算することを含み、Ｍ(3)，Ｓ(3)は、それぞれ動き及び構造に対応するランク３行列であり、
前記３次元構造を復元する段階は、
１次方程式
【数２２】
0003651590.tif 000108
を使用して、Ｕ及びΩに関する３Ｎtot線形制約条件を得るために、
【外１３】
0003651590.tif 000109
から構造未知数Ｑz，Ｂz及びＺ-1を除外することを含み、ここでＵ及びΩは未知の３×３行列であり、
【外１４】
0003651590.tif 000110
はそれぞれ総合平行移動フロー・ベクトル及び総合回転フロー・ベクトルを表し、ＳＶＤを使用してΟ(Ｎtot)計算量でこれらの制約条件を解き、ここで前記総合平行移動フロー・ベクトルは
【数２３】
0003651590.tif 000111
により定義され、ｙ及びｚ成分に対しても同様に定義され、
【数２４】
0003651590.tif 000112
であり、ここでＡは、射影中心及び基準画像内の描かれた線を通る第１の平面の法線であり、第２の平面の平面の法線であるＢは、３次元線Ｌ上の任意の点Ｑに対して、Ｂ・Ａ＝０及びＢ・Ｑ＝−１であることを条件とする(require)ことにより定義され、ＢU≡Ｂ・ＰU及びＢLＬ≡Ｂ・ＰLはＢの上部及び下部成分であり、ここで
【数２５】
0003651590.tif 000113
であり、ここでΦIx≡−[Ｚ-1Ｉx]，ΦIy≡−[Ｚ-1Ｉy]，ΦIz≡[Ｚ-1(∇Ｉ・ｐ)]であり、
ここで、Ｑは基準画像内の追跡された３次元点に対する３次元座標であり、Ｚnはカメラからカメラの光軸に沿ったｎ番目のピクセルに描かれた３次元の点までの深度であり、
Ｕ及びΩを仮定して、
【数２６】
0003651590.tif 000114
から、前記構造未知数Ｑz，Ｂz及びＺ-1を復元することを含み、
前記カメラの動きを復元する段階は、
前記平行移動を復元するために、
【数２７】
0003651590.tif 000115
を使用することと、
【数２８】
0003651590.tif 000116
から、回転ωix，ωiy，ωizを復元することを含み、
ここでＣはＣii'≡δii'＋１である定数（ＮI−１）×（ＮI−１）行列であり、Ｔは前記平行移動を表す、３次元シーン構造及びカメラの動きを復元する方法。
【請求項５】
前記画像データは、点、線及び強度からなる群から選択される１以上のものである請求項４に記載の方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、３次元構造及びカメラの動き(motion)を復元する方法に関し、特に、点、線から、及び／または画像強度から直接に、構造及び動きデータを復元する線形アルゴリズムに関する。
【０００２】
【従来の技術】
２次元（２Ｄ）画像から導出した情報から３次元（３Ｄ）モデルを描画する技術は、コンピュータ・グラフィクス以前からあり、その根源は、写真測量及びコンピュータ・ビジョンの分野にある。
【０００３】
写真測量は、写真が撮影されるとき、３次元世界は平らな２次元画像の平面に透視法（遠近法）により射影されるという、基本的な考えに基づいている。結果として、個々の点で見られる２次元画像上の特徴物(feature)は、実際には、カメラに始まり無限遠点に延びるある光線(ray)に沿って位置している。２枚の異なる写真で同じ特徴物を見た場合に、特徴物が２つの光線の交点に位置すると制約することにより、実際の位置を決定することができる。この処理は三角測量として知られている。この処理を使用して、少なくとも２つの画像に見られる任意の点を３次元中に配置することができる。同様に、十分な数の点を使用すれば、未知のカメラ位置を決定することも可能である。写真測量及び三角測量の技術は、航空画像から地勢図を作成するような応用に使用されていた。しかしながら、写真測量処理は時間を要し非効率的である。
【０００４】
コンピュータ・ビジョン技術は、ステレオ（立体）画像から３次元シーン構造(3D scene structure)を復元することを含み、２つの画像の間の対応は、反復アルゴリズムによって、その２つの画像から自動的に定められる。３次元シーンを再構成するために、反復アルゴリズムは、２つの画像間での点の一致を検索する。ステレオ画像から得られた３次元シーン構造を使用して、カメラの位置と動きを決定することも可能である。
現在のコンピュータ技術は、運動(motion)ベースの再構成(reconstruction)に関心を集中しており、２次元画像から３次元構造（ジオメトリ）を推測する問題に対するコンピュータテクノロジーの自然な応用である。これは運動からの構造（ＳＦＭ；Structure-from-Motion）として知られている。運動からの構造、すなわち複数の２次元画像から未知の３次元シーンを再構成する問題は、コンピュータ・ビジョンにおいて最も研究された問題の１つである。
【０００５】
現在知られている大多数のＳＦＭアルゴリズムは、事前に計算された特徴対応(feature correspondence)、通常は追跡点(tracked point)からシーンを再構成する。他のアルゴリズムは、対応関係を別に計算する過程を用いずに、画像から強度を再構成する直接法である。従来の直接法は、画像の数が少数に限定されており、シーンに対して通常は平坦であるという強い仮定を必要とし、あるいは反復最適法を使用し、初期推定値を必要とした。
【０００６】
これらの方法は、相補的な長所と短所を有する。通常、画像データのある部分は、対応を判定するために使用できないような低い品質のものである。直接法が、データ全体の冗長性を利用することによってデータの一部の低い品質を補償しようと試みるのに対して、特徴ベースの方法は、追跡が比較的容易な少数の区別しやすい点あるいは線の特徴物を予め選択することによって、この問題に対処している。特徴ベースの方法は、入力データが比較的信頼できるという利点を有するが、利用できる画像情報の大部分を無視しており、３次元シーンのまばらな再構成しかできない。直接法は、入力データの冗長性によって、緻密で正確な３次元再構成ができる可能性を有するが、データのわずかな部分にある大きな誤差によって過度に影響を受けるおそれがある。
【０００７】
追跡線(tracked line)にもとづく方法は、 M. Spetsakis,“A linear Algorithm for Point and Line based Structure From Motion（点及び線にもとづく運動からの構造用の線形アルゴリズム）”, CVGIP, 56:2, 230-241,1992に説明されており、ここでは、３つの画像中の１３本の線に対するオリジナルの線形アルゴリズムが提案されている。最適化による方法は、C. J. Taylor, D. Kriegmann,“Structure and Motion from Line Segment in Multiple Images（複数の画像中の線分からの構造と動き）”, PAMI, 17:11, 1021-1032, 1995に開示されている。さらに、K. Morris and I. Kanade,“A unified factorization algorithm for points, line segments and planes with uncertainty Models（不確定性モデルを有する点、線分及び平面に対する統一因数分解アルゴリズム）”, ICCV, 696-702, 1998は、アフィン・フレームワークにおける線に関する研究について説明している。線及び点に対する射影法は、B. Triggs,“Factorization methods for projective structure and motion（射影構造及び動きに対する因数分解法）”, CVPR, 845-851, 1996に説明されており、少数のフレームから射影深度(projective depth)を計算することを含んでいる。Hartley,“In Defense of the Eight-Point Algorithm（８点アルゴリズムの擁護）”, PAMI, 19, 580-593, 1995において、Hartleyは、３つ以上の画像を追跡した点及び線から再構成する全透視法(full perspective approach)を提案した。
【０００８】
M. Irani,“Multi-Frame Optical Flow Estimation using Subspace Constraints（部分空間制約条件を使用するマルチフレーム・オプティカル・フロー推定）”, ICCV, 626-633, 1999に説明されている方法は、画像強度から直接再構成を行なっている。対応を復元するためのIraniにおける本質的な過程は、B. Lucas and T. Kanade,“An Iterative Image Registration Technique with an Application to Stereo Vision（反復画像位置決め技法と立体視への応用）”, IJCAI, 674-679, 1981に説明されているオプティカル・フロー法のマルチフレーム一般化であり、剛性制約条件(rigidity constraint)ではなく、平滑度制約条件(smoothness constraint)によっている。初期状態で計算ができなかったＤのエントリを書き込むために、Iraniの方法では、Ｄの因数分解を単純に使用している。
【０００９】
輝度不変性方程式(brightness constancy equation)（後述の式(7)）を、Iraniは行列形式ではΔ＝−ＤＩと書いており、ここでＤはシフトｄiを表にしたものであり、Ｉは強度勾配∇Ｉ(ｐn)を含んでいる。Iraniは、Ｄが、（既知の較正を有するカメラに対して）ランク（階数）６を有することを指摘しており、それはΔがランク６を持たねばならないことを意味している。雑音の影響を減少するために、Iraniは観測されたΔをランク６の１つに射影している。次に、行列方程式ＤＩ2＝−Δ2をもたらすオプティカル・フローを復元するために、Iraniは、Lucas-Kanade法のマルチ画像形式を適用しており、ここで、Ｉ2のエントリは「スムージング」ウィンドウにわたって合計された２乗強度勾配(squared intensity gradient)Ｉa Ｉbであり、Δ2のエントリはＩaΔＩの形式を有する。追加されたLucas-Kanade平滑化制約条件によって、画像強度が少なくとも１つの方向で一定であるスムージング・ウィンドウの中を除いて、シフトＤあるいはｄinは、Ｄ＝−Δ2[Ｉ2]+が一般逆行列を表すとして計算することができる。Ｄに対してランク制約条件を使用して、Iraniは、強度が１方向に一定であるウィンドウに対するＤの付加的なエントリを決定している。
【００１０】
【発明が解決しようとする課題】
本発明の目的は、点、線から、及び／または画像強度から直接に、３次元構造やカメラの動きデータを復元する線形アルゴリズムを提供することにある。
【００１１】
【課題を解決するための手段】
本発明は、画像列の基準画像が基準透視点にあるカメラにより撮像され、画像列の１以上の連続した画像がカメラを平行移動及び／または回転することによる１以上の連続した異なる透視点において撮像される、複数画像列から得られた画像データから３次元シーン構造及びカメラの動きを復元する方法を指向している。本発明の方法は、
（ａ）基準画像に対する連続した画像の各々に対して、基準透視点から連続した異なる透視点までのカメラの平行移動及び／または回転から導出される画像データ・シフトを決定する段階と、
（ｂ）各々の画像に対する画像データ・シフトを含むシフト・データ行列を構成する段階と、
（ｃ）ＳＶＤ（特異値分解）を使用して、シフト・データ行列から、一方のランク（階数）３の因子行列が３次元構造に対応する３次元構造行列であり他方のランク３の因子行列がカメラの動きに対応する動き行列である、２つのランク３の因子行列を計算する段階と、
（ｄ）１次方程式を解くことにより、ＳＶＤを用いて３次元構造行列から３次元構造を復元する段階と、
（ｅ）復元された３次元構造を使用して、カメラの動き行列からカメラの動きを復元する段階と、を有する。
【００１２】
本発明による方法は、一般的な運動アルゴリズムであり、連続した射影の各々に対するカメラの位置は、単一の平面の上には位置していない。この方法は、点、線あるいは強度から画像を再構成することができる。
【００１３】
本発明は、特徴にもとづく再構成と直接法とを組み合わせた本質的に線形アルゴリズムである。本発明は、特徴対応が利用できれば、特徴対応を及び／または画像強度を直接に使用することができる。B. K. P. Horn and B. G. Schunck,“Determining Optical Flow（オプティカル・フローの決定）”, AI, 17, 185-203, 1981に説明されているようなオプティカル・フロー方法とは異なり、このアルゴリズムは剛性制約条件を利用し、動きが小さくかつ輝度不変性方程式が有効であると仮定して、原則として平滑化制約条件を必要としない。しかし実際には、時にはこのアルゴリズムにおいても平滑度制約条件を課すことが必要である。
【００１４】
本発明による方法は、点及び線の特徴物の両方から再構成をすることができる。特に人工の環境では直線が豊富に存在するので、これは本発明の重要な側面である。多くの場合、両方の特徴形式を追跡することが可能である。一本の線上に多くのサンプル点を見出せる可能性があるために、線は非常に正確に局所化できる。線は、疑似Ｔ接続点(spurious T-junction)及びアパーチャー効果のような、点を追跡する場合の問題のいくつかの影響を受けない。ここで説明する方法は、入力データの３つの形式のすべてから任意の組み合わせで、あるいは、全体透視図の下の任意の多くの画像列にわたって別々に、再構成することができる最初のものである。
【００１５】
【発明の実施の形態】
本発明による方法のこれらの及び他の特徴、態様及び利点は、下記の説明、添付した特許請求の範囲及び添付図面を参照して良く理解されるであろう。
【００１６】
［定義］
本明細書において、本発明は、ここで定義する次の用語にしたがって説明される。本発明による方法は、較正は既知であり、焦点距離は１(wlog)であると、仮定する。未較正の画像列に対する方法への修正は、後に説明される。
【００１７】
Ｎi個の一連の画像を想定する。０番目の画像を基準画像(reference image)として選択する。Ｒi，Ｔiは、基準画像と画像ｉの間の回転と平行移動を表すものとする。小さい回転を、回転速度ω≡(ωix,ωiy,ωiz)Tによりパラメータ化する。[Ｔi]2≡(ＴixＴiy)Tとする。
【００１８】
３次元の点Ｐを、Ｐ'＝Ｒ（Ｐ−Ｔ）として変換する。Ｎp個の点とＮL個の線が、一連の画像上で追跡されると仮定する。理解しやすいように、追跡点に対しては、強度画像(intensity image)内のピクセル位置に対するものとは異なる表記法を使用する。ｑim≡(ｑx,ｑy)imはｍ番目の追跡点を表し、Ａil≡(α,β,γ)Tilはｉ番目の画像中のｉ番目の線を表すものとする。
【００１９】
【数２９】
0003651590.tif 000002
【００２０】
とする。Ｒ＊ｑは、回転：
【００２１】
【数３０】
0003651590.tif 000003
【００２２】
後のｑから得られた画像点を表すものとする。３つの点回転フロー(point rotational flow) ｒ(1)(ｘ,ｙ)，ｒ(2)(ｘ,ｙ)，ｒ(3)(ｘ,ｙ)を、
【００２３】
【数３１】
0003651590.tif 000004
【００２４】
により定義する。
【００２５】
【外１５】
0003651590.tif 000005
【００２６】
は、３次元単位ベクトル
【００２７】
【数３２】
0003651590.tif 000006
【００２８】
を表す。同様に、
【００２９】
【数３３】
0003651590.tif 000007
【００３０】
である。Ｓim≡ｑim−ｑ0mはｍ番目の追跡点に対する画像移動を表すものとし、Ｑm＝(Ｑx,Ｑy,Ｑz)Tmは追跡点ｑmに対応する３次元シーン点を表すものとする。次に３次元線Ｌを、その線が位置する２つの場所によりパラメータで表す。Ａにより表される第１の平面は、射影の中心及び基準画像内の描かれた線(imaged line)を通過している。第２の平面Ｂに対する法線は、Ｌ上の任意の点Ｑに対してＢ・Ａ＝０及びＢ・Ｑ＝−１を求めることにより定義される。
【００３１】
ｐn≡(ｘn,ｙn)Tを、ｎ番目のピクセル位置の画像座標であるとする。Ｎxをピクセル総数として、１からＮxまでのピクセルに順序をつける。Ｉiはｉ番目の強度画像を表し、Ｉin＝Ｉi(ｐn)はＩi内のｎ番目のピクセル位置における画像強度を表すものとする。Ｐnは、基準画像内のｐnにおいて描かれた３次元の点を表し、Ｉ0の座標系においてＰn≡(Ｘn,Ｘn,Ｘｎ)Tである。ｄinは、３次元の点ＰnのＩ0からＩiまでの画像位置のシフトを表すものとする。
【００３２】
Ｖaは、整数ａにより索引付けされた１組の量であるとする。Ｖaにより与えられる成分を有するベクトルを表すために、表記｛Ｖ｝を使用する。
【００３３】
［アルゴリズムの説明−予備処理］
本発明による方法は、平行移動運動が大き過ぎないことを必要とする（たとえば、｜Ｔ／Ｚmin｜≦１／３であり、カメラの位置は平面内にないこと）。追跡点及び線を仮定すると、式(1)を最小化することにより、回転Ｒiに対して、基準画像と後続の画像のそれぞれの間の回転をほぼ復元できる。
【００３４】
【数３４】
0003651590.tif 000008
【００３５】
ここで、点の測定と線の測定との相対的な確度にしたがって定数を調整するべきである。式(1)を最小化することは、平行移動があまり大きくない限り、正確な回転を与える。回転を復元した後に、回転に対して補償が行われ、以後回転は小さいものとして処理することが可能である。
【００３６】
［点］
点の移動は、参考文献として本明細書に包含される、J. Oliensis, “A Multi-Frame Structure from Motion Algorithm under Perspective Projection（透視投影下の運動アルゴリズムからのマルチフレーム構造）”, IJCV, 34:2/3, 163-192, 1999によって計算される。点の移動は、
【００３７】
【数３５】
0003651590.tif 000009
【００３８】
であり、ここで回転フローｆ(Ｒi,ｑim)≡ｑim−(Ｒi)-1＊ｑimである。この仮定の下で、
【００３９】
【数３６】
0003651590.tif 000010
【００４０】
と近似することができ、ここで最後の３つの項は１次の回転フローを与える。それに対応して、３つの長さ２Ｎpの平行移動フロー・ベクトル、
【００４１】
【数３７】
0003651590.tif 000011
【００４２】
及び長さ２Ｎpの回転フロー
【００４３】
【数３８】
0003651590.tif 000012
【００４４】
を定義する。
【００４５】
次に、シフトｓim ｓimをＮi×２Ｎp行列Ｓの中に収集する、この行列で、各行は、異なる画像ｉに対応し、［[Ｓix]T [Ｓiy]T］に等しい。したがって、
【００４６】
【数３９】
0003651590.tif 000013
【００４７】
となる。
【００４８】
［線］
線測定モデルは次のように説明される。測定された線の雑音は、線の測定された位置と真の位置の間の統合化した画像距離(integrated image distance)に比例していると、合理的に想定することができる。ΘFOVは、視野（ＦＯＶ；Field of View）をラジアンで表したものとする。通常、ΘFOV＜１である。測定された線が真の線とはｚ回転だけ異なるならば、これはΟ((ΘFOV／２)2)の雑音を与える。ｘ−ｙ平面内の軸の周りの回転は、典型的にはｚ回転のものより大きいΟ(ΘFOV)の雑音を与える。これらの推定値は、真の線と測定された線により仕切られた変位領域(displacement region)が、第１の場合には２つの三角形のように見え、第２の場合には１つの四辺形のように見える事実を反映している。本明細書において説明する表記に対して、これは線の測定は、ＡzをＡx，Ａyよりおよそ１／ΘFOV倍だけ正確に決定することを意味する。
【００４９】
線に対する画像フロー(image flow)は次のように説明できる。Ａ≡(Ａ;０)及びＢ≡(Ｂ;−１)を、線Ｌに対する平面法線(plane normal)Ａ，Ｂに対応する長さ４の同次ベクトルであるとする。Ｑ≡(Ｑ;１)をＬ上の点に対応する同次ベクトルであるとする。すると、Ａ・Ｂ＝Ａ・Ｑ＝Ｂ・Ｑ＝０である。回転Ｒ及び平行移動Ｔの後に、必要条件Ａ'・Ｑ'＝Ｂ'・Ｑ'＝０から変換されたＡ'Ｂ'を決定する。
【００５０】
【数４０】
0003651590.tif 000014
【００５１】
を使用して、この必要条件を満足させることができる。しかし、線の新しい画像に対するＡ'はＡ'4＝０を有しなければならないが、Ａ*は必ずしもＡ'4＝０を有さなくてもよい。そこで、Ａ'≡Ａ*−Ａ*4Ｂ*／Ｂ*4と設定する、これは新しい画像線(image line)が
【００５２】
【数４１】
0003651590.tif 000015
【００５３】
により与えられることを意味する。
【００５４】
ここで、線フロー(line flow) δＡil＝Ａil−Ａ0lを考える。δＡR＝Ａi−Ｒ-1Ａi及び
【００５５】
【数４２】
0003651590.tif 000016
【００５６】
と定義すると、δＡ＝δＡT＋δＡRＲである。｜Ｔ・Ｂ｜≪１である、小さい回転及び平行移動に対しては、
【００５７】
【数４３】
0003651590.tif 000017
【００５８】
である。
【００５９】
Ｂ・Ａ＝０は、
【００６０】
【数４４】
0003651590.tif 000018
【００６１】
であることを意味する。各Ａを定義する際の尺度あいまい性(scale ambiguity)を排除するために、基準画像において｜Ａ01｜＝１とし、直線フローが０＝δＡil・Ａol≡（Ａil−Ａ0l）・Ａolを正確に満足することを必要とする。測定誤差がすべての線に対して同様であるべきであるという事実を反映するために、基準画像内ですべてのＡ0を同じ大きさに正規化する。０＝δＡil・Ａolという必要条件は、ここでの線の表現及び小さい運動の仮定と矛盾しないやり方で、「ゲージを決定」するものである。
【００６２】
「ゲージの決定」の後、各δＡilは自由度を２つ含んでいる。δＡilを、それぞれ上部方向及び下部方向と呼ぶ２つの方向、
【００６３】
【外１６】
0003651590.tif 000019
【００６４】
に沿ったその射影により表す。単位３−ベクトルＰlU，ＰlLがこれらの２つの方向に射影するものとする。典型的なΘFOV＜１に対して、
【００６５】
【数４５】
0003651590.tif 000020
【００６６】
であり、かつＡlの上部成分はＡl,zにほぼ等しい。したがって、画像測定は、上部成分を下部成分よりもおおよそ１／ΘFOV倍正確に決定する。点の定義と同様に、線の平行移動フローを
【００６７】
【数４６】
0003651590.tif 000021
【００６８】
と定義する。ここで、これらは長さ２ＮLベクトルである。ＢU≡Ｂ・ＰU及びＢL≡Ｂ・ＰLは、Ｂの上部及び下部成分である。同様に、回転フローを定義する。
【００６９】
【数４７】
0003651590.tif 000022
【００７０】
∧を、各行が異なる画像ｉに対応するＮl×２ＮL行列であり、［[ＰU・δＡi]T[ＰU・δＡi]T］に等しいとすれば、
【００７１】
【数４８】
0003651590.tif 000023
【００７２】
となる。
【００７３】
［強度］
追跡特徴点(tracked feature point)に対するｓimと同様に、同じ議論をｄinすなわちピクセル位置ｐnにおいて描かれた３次元の点に対応する画像シフト(image shift)に適用することができる。したがって、
【００７４】
【数４９】
0003651590.tif 000024
【００７５】
となる。∇Ｉn≡∇Ｉ(Ｐn)≡(Ｉx,Ｉy)Tnとして、ピクセルの離散格子に対してある適切な定義を有する、Ｉ0に対する画像強度の勾配を表すものとする。ΔＩinは、基準画像に対する強度の変化を定義するものとする。その最も単純な定義は、ΔＩin＝Ｉin−Ｉ0nである。輝度不変性方程式は、
ΔＩin＋∇Ｉn・ｄin＝０          (7)
である。
【００７６】
この式及び前述の式は、次式を意味する。
【００７７】
【数５０】
0003651590.tif 000025
【００７８】
長さがＮxである３つの平行移動フロー・ベクトルΦIx≡−[Ｚ-1Ｉx]，ΦIy≡−[Ｚ-1Ｉy]，ΦIz≡−[Ｚ-1(∇Ｉ・ｐ)]、及び、長さがＮxである回転フローΨIx≡−[∇Ｉ・ｒ(1)(ｐ)]，ΨIy≡−[∇Ｉ・ｒ(2)(ｐ)]，ΨIz≡−[∇Ｉ・ｒ(3)(ｐ)]を定義する。
【００７９】
Δは、各行が画像ｉに対応し[ΔＩi]Tに等しいＮI×Ｎx行列であるとすれば、
【００８０】
【数５１】
0003651590.tif 000026
【００８１】
である。
【００８２】
［再構成］
３次元シーンの再構成を次に説明する。点、線及び画像強度に対する全体の回転フロー・ベクトルを
【００８３】
【数５２】
0003651590.tif 000027
【００８４】
により定義し、ｙ及びｚ成分に対しても同様に定義する。ここでωL及びωIは、点、線及び強度データの相対的な雑音度(noisiness)に応じて利用者が設定すべき一定の重みである。
【００８５】
【外１７】
0003651590.tif 000028
【００８６】
は、２Ｎp＋２ＮL＋ＮX≡Ｎtotの長さを有する。この定義から、
【００８７】
【外１８】
0003651590.tif 000029
【００８８】
は測定された量から計算可能であって、既知の値であるかのように扱うことができることを確かめることができる。ハウスホルダー(Householder)行列を使用して、（Ｎtot−３）×Ｎtot行列Ｈを計算し、３つの
【００８９】
【外１９】
0003651590.tif 000030
【００９０】
を零化することができる。Ｈによる計算と乗算は、Ο(Ｎtot)程度の計算量を必要とする。ＮI×Ｎtot行列
【００９１】
【数５３】
0003651590.tif 000031
【００９２】
を定義する。ＣをＣii'≡δii'＋１である定数の（ＮI−１）×（ＮI−１）行列であるとする。特別な処理に対する基準画像を選び出すことによる偏りを無効にするために、Ｃは考慮される。
【００９３】
【数５４】
0003651590.tif 000032
【００９４】
を定義する。
【００９５】
【数５５】
0003651590.tif 000033
【００９６】
により、全体の平行移動フロー・ベクトルを定義する。ｙ及びｚに対しても同様に定義する。以上から、式(2),(6),(8)から、式(9)が得られる。
【００９７】
【数５６】
0003651590.tif 000034
【００９８】
【外２０】
0003651590.tif 000035
【００９９】
は、およそランク３である。
【０１００】
本発明による方法の基本的アルゴリズムは、以下の１）〜５）の過程からなる。
【０１０１】
過程１）  Ｈを定義し、
【０１０２】
【外２１】
0003651590.tif 000036
【０１０３】
を計算する。特異値分割（ＳＶＤ；singular value decomposition）を使用して、
【０１０４】
【数５７】
0003651590.tif 000037
【０１０５】
の最良のランク３因数分解を計算する、ここでＭ(3)，Ｓ(3)はそれぞれ動きと構造に対応するランク３行列である。
【０１０６】
過程２）  Ｓ(3)は、式(10)を満足する。
【０１０７】
【数５８】
0003651590.tif 000038
【０１０８】
ここでＵ及びΩは、未知の３×３行列である。この連立方程式において、Ｕ及びΩにおける１８個の未知数に対する３Ｎtot個の線形制約条件(linear constraint)を得るために、
【０１０９】
【外２２】
0003651590.tif 000039
【０１１０】
から未知数Ｑz，Ｂz及びＺ-1を除外する。ＳＶＤを使用して、Ο(Ｎtot)程度の計算量でこれらの制約条件は解かれる。
【０１１１】
過程３）  Ｕ及びΩを仮定して、
【０１１２】
【数５９】
0003651590.tif 000040
【０１１３】
から構造未知数(structure unknown)Ｑz，Ｂz及びＺ-1を復元する、ここで、Ｑは基準画像内の画像ピクセルに対する３次元座標であり、Ｂはカメラの中心から３次元の線までの最短のベクトルであり、Ｚはカメラの光軸に沿ったカメラから３次元シーンまでの深度である。
【０１１４】
過程４）  Ｕを仮定して、平行移動を復元するために、
【０１１５】
【数６０】
0003651590.tif 000041
【０１１６】
を使用する。
【０１１７】
過程５）  回転ωix，ωiy，ωizを、
【０１１８】
【数６１】
0003651590.tif 000042
【０１１９】
から復元する。
【０１２０】
以上の説明は、偏りの訂正のために重要な一部の段階を除いている。すなわち、１）
【０１２１】
【外２３】
0003651590.tif 000043
【０１２２】
内の上部線成分は、これらの成分の測定におけるより高い精度がもたらされるように、１／ΘFOVに比例した重み付けがなされている。２）ＦＯＶは通常は小さいという事実による偏りを修正するために、上記の過程２）〜４）において、ｘ及びｙ成分はそのままとする一方で、ＴizをＦＯＶに比例して重み付けし、
【０１２３】
【外２４】
0003651590.tif 000044
【０１２４】
をこの逆だけ重み付けしている。３）上記のアルゴリズムの各過程は、より良い結果を与え、小さい運動であるという仮定に対して修正するために、反復することができる。これは、元の特徴点シフト(original feature point shift)に１−Ｚ-1Ｔzを乗算し、線シフト(line shift)に１−Ｂ・Ｔを乗算することを含む。運動及び雑音が小さく、かつカメラの位置が１つの平面上に位置しなければ、正しい再構成にこのアルゴリズムが収束するように、保証されている。
【０１２５】
当然ながら、このアルゴリズムのうち強度に基づく部分に対する結果は、導関数を計算するために使用される技術に決定的に依存する。これをさらに頑強（ロバスト）にするために、粗密モード(coarse-to-fine mode)で再構成をねじり(warp)、再計算することにより、移動ｄinの大きさを反復して小さくすべきである。１つの具体化例は、画像の動きが小さいとした仮定が現在の尺度に対して有効であるか否かを判定するために、異なる画像に対して計算された空間導関数(spatial derivative)の間の一致性を使用して、単一の尺度で画像導関数(image derivative)を単純に計算する。仮定が有効であるピクセル・サイトのみが、動きの計算のために使用される。
【０１２６】
時には、アルゴリズムを実行する前に、画像強度に前処理を行なうことが有用であることが判明している。まずラプラシアン画像(Laplacian image)を計算し、次に、輪郭領域を強調しテクスチャのない画像を抑制するために、スティグマ(stigma)関数により強度を変換する。これは、未処理画像と追跡特徴物(tracked feature)の選択された集合との間の中間の形式を強度画像に与える。
【０１２７】
［線に対する浅浮き彫り(bas-relief)］
浅浮き彫りのあいまい性によって、点特徴に対するベクトル[Ｚ-1]の定数成分を復元することは難しいが、通常すべての他の成分は正確に復元することができる。同様の結果を線に対しても得ることができる。線に対する回転フローはω×Ａである。典型的なΘFOV≪１に対して、｜Ａz｜≪｜Ａx,y｜である。したがって、
【０１２８】
【外２５】
0003651590.tif 000045
【０１２９】
回転は、
【０１３０】
【外２６】
0003651590.tif 000046
【０１３１】
にほぼ比例するフローを与える。前述の式(4)から、ｚ−平行移動Ｔzの効果は、｜Ａz｜及び、回転の効果からの[Ｂz]定数成分(constant component)の一定成分による平行移動フローにより抑制され、この成分を正確に復元することを難しくする。
【０１３２】
次のように、この効果のさらに定量的な分析を得ることが可能である。入力データが線のみからなると仮定すると、
【０１３３】
【数６２】
0003651590.tif 000047
【０１３４】
は、
【０１３５】
【数６３】
0003651590.tif 000048
【０１３６】
を意味し、ここで
【０１３７】
【外２７】
0003651590.tif 000049
【０１３８】
はＳ(3)を零化する射影行列である。小さな重複部分
【０１３９】
【外２８】
0003651590.tif 000050
【０１４０】
を作るそれら[Ｂ]成分は、小さな重複部分
【０１４１】
【外２９】
0003651590.tif 000051
【０１４２】
も作る。標準摂動理論から、雑音がこれらの成分をＢに対する真の解に混合するであろうという結果になり、汚染の量は固有値の大きさに逆比例する。Ｂ≡［[Ｂx]，[Ｂy]，[Ｂz]］と定義する。本発明者らは、いくつかの系列に対してＨBの固有値を計算し、
【０１４３】
【数６４】
0003651590.tif 000052
【０１４４】
であるように定義し、そして、予想したように、その固有ベクトルが[Ｂz]の定数成分にほぼ等しい１つの小さな固有値が存在することを発見した。したがって、この成分を復元する場合の不正確さを推定するために、入力データを使用することができる。
【０１４５】
［強度に対する一般的な動きのアルゴリズム］
単に強度に対する一般的な動きのアルゴリズムは、次の通りである。
【０１４６】
過程０：  最初に、平行移動が０であると仮定し、回転を復元し、すべての画像Ｉ1, ..., ＩNI-1を基準画像Ｉ0に向けてねじる(warp)。この場合、画像移動ｄinを回転されていない画像と呼ぶ。
【０１４７】
過程１：  次にＨ及びΔCH≡Ｃ-1/2ΔＨTを計算し、特異値分割（ＳＶＤ）を使用して最良のランク３因数分解
【０１４８】
【外３０】
0003651590.tif 000053
【０１４９】
を計算する。ここでＭ(3)及びＳ(3)はランク３であり、それぞれ動きと構造に対応する。
【０１５０】
過程２ａ：  さらに、Ｓ(3)はΦ＝ＨTＳ(3)Ｕ＋ΨΩを満たし、ここでＵとΩは未知の３×３行列である。この連立方程式において、Ｕ及びΩにおける１８個の未知数に対する３Ｎp線形制約条件を得るために、
【０１５１】
【外３１】
0003651590.tif 000054
【０１５２】
から未知数Ｚ-1を除く。ＳＶＤを使用することにより、Ο(Ｎp)回の計算でこれらの制約条件は解かれる。次にＵ及びΩを仮定して、構造未知数Ｚ-1を式(5)から復元する。
【０１５３】
過程２ｂ：  次にＵを仮定して、平行移動を復元するために
【０１５４】
【数６５】
0003651590.tif 000055
【０１５５】
及び式(4)を使用する。
【０１５６】
過程３：  次に、回転Ｗは
【０１５７】
【数６６】
0003651590.tif 000056
【０１５８】
から復元することが可能である。
【０１５９】
［アルゴリズムの解析］
次に、上述した強度に対する一般的な動きのアルゴリズムの各過程を詳しく説明する。
【０１６０】
過程０：
J. Oliensis,“Rigorous Bounds for Two-Frame Structure from Motion（運動からの２フレームの構造に対する厳密な限界）”, IUW, 1225-1231, 1994に説明されているように、上述の基本的な一般的な動きのアルゴリズムにおける平行移動の無視は、平行移動が穏やかなものであって｜Ｔitrue｜／Ｚmin≦１／４である場合に、回転を復元する際には小さい誤差をしか生まない。回転を復元するための２つのマルチフレーム技術を説明する。式(4)において平行移動を無視することは、
【０１６１】
【数６７】
0003651590.tif 000057
【０１６２】
を与え、これはＷに対して解くことができる。この技術は、平坦なシーンを再構成することに関し、ここでは基準画像内のノイズを過度に強調することによる偏りを減少させるＣ-1/2を使用していることを除いて、L. Zelnik-Manor and M. Irani,“Multi-Frame Alignment of Planes（平面のマルチフレーム・アライメント）”, CVPR, I:151-156, 1999で説明されている射影アルゴリズムに類似している。例によって、Ｗに対して反復してねじりを加えて解を再度求めることにより、本技術をより大きな回転に拡張することができる。
【０１６３】
１次の回転変位を用いているので、上述の技術は、小さな回転に対して最も適したものである。他の方法は、P. Sturm and B. Triggs,“A factorization based algorithm for multi-image projective structure and motion（マルチイメージ射影構造と動きに対する因数分解にもとづくアルゴリズム）”, ECCV, 709-720, 1996で述べられているマルチフレームSturm- Triggs技術についての論文であるJ. Oliensis and M. Werman,“Structure from Motion from Points, Lines and Intensities（点、線及び強度からの運動からの構造）”, CVPR, 2000に説明されているように、平坦なシーンへの拡張から始まる。このマルチフレーム方法は、基準画像と後続の画像の間の任意の大きさの相同(homography)を復元する。復元された相同に最も近い、正の行列式の直交行列を選択することにより、回転を復元することができる。
【０１６４】
段階２ａ：
式(4)から、Ｓ(3)の列は、Φx，Φy，Φzにより生成された部分空間とほぼ同じ部分空間を生成する。すなわち、
【０１６５】
【数６８】
0003651590.tif 000058
【０１６６】
である。ここでＵは３×３行列である。上述の式(7)は、Ｕ及びＺ-1nに対して直接解くことが可能な、過大に条件付きの(overconstrained)連立１次方程式を与える。しかし、これはΟ(Ｎ3p)程度の計算量である。大きい画像に対しては、代わりに次のΟ(Ｎp)算法を使用する。
【０１６７】
まず、射影行列
【０１６８】
【数６９】
0003651590.tif 000059
【０１６９】
を定義する。式(7)にＨTを乗算すると、
【０１７０】
【数７０】
0003651590.tif 000060
【０１７１】
を与え、これは
【０１７２】
【数７１】
0003651590.tif 000061
【０１７３】
を意味する。
【０１７４】
上述の式(8)からＺ-1nを除き、Ｕ及びΩに対して直接解く。Ｕの列をＵ≡［Ｕ1Ｕ2Ｕ3］により表し、同様にΩの列も表す。Ｕ'3≡ｓ-1Ｕ3及びΩ'3≡ｓ-1Ω3とする。ここで尺度ｓは、画像の中心からの画像点の平均距離に等しい。アルゴリズムの偏りを減少させるために、ｓを導入する。Φx，Φy，Φzの定義から、上述の式(8)は、
【０１７５】
【数７２】
0003651590.tif 000062
【０１７６】
及びｘとｙを入れ替えた同種の数式を意味する。このアルゴリズムの過程２ａは、Ω及びＵに対するこの連立１次方程式を最小二乗法で解き、次に、これらの解と式(8)からＺ-1nを復元する。計算量はΟ(Ｎp)である。なお、過程２ａは、その計算を∇Ｉに基づいている。測定された基準画像Ｉ0から計算した∇Ｉの値を使用すれば、Ｕ，Ω，Ｚ-1nの推定値は真のマルチフレーム推定値ではないであろう。よりよいマルチフレーム推定値を得るために、まずＩ0及び∇Ｉを全ての画像データに基づいて再計算し、次のように、式(9)によってＵ，Ωを計算するために、その結果を使用することができる。
【０１７７】
Δ(3)CHを、ΔCH≡Ｃ-1/2ΔＨTに対する最良のランク３近似であるとする。未知の回転フローによって、Δ(3)CHは、平行移動による強度シフトを与える。
【０１７８】
【数７３】
0003651590.tif 000063
【０１７９】
次に式(10)をＷに対して最小二乗法で解く、これはΔに対してΔe、Ｉ0に対してＩ0eという改良された推定値を与える。これらは次に、式(9)で使用できる。しかし、Ｉ0eのこの計算は、Ο(Ｎ-1Iω2,Ｎ-2Iτ2Ｚ-2)の誤差を与える。Ｉ0のそもそもの雑音がこれに満たなければ、過程２ａにおいてＩ0を直接使用するべきである。
【０１８０】
過程２ｂ：
【０１８１】
【数７４】
0003651590.tif 000064
【０１８２】
を得る。ここでＵTも３×３行列である。このアルゴリズムは、すでに計算したΦを使用して、ＵTに関して連立一次方程式
【０１８３】
【数７５】
0003651590.tif 000065
【０１８４】
を解き、次に、
【０１８５】
【外３２】
0003651590.tif 000066
【０１８６】
を復元するために結果を代入する(plug)ことにより、ＵT及び
【０１８７】
【外３３】
0003651590.tif 000067
【０１８８】
を復元する。
【０１８９】
［一般的注釈］
実験的には、本発明者等による画像列においては、強度パターンが画像ごとに大幅に変動することを発見した。実際には、本発明者等は、元の画像列の修正されたバージョンに、次により得られた本発明者による方法を適用する。
【０１９０】
１）エッジを強調するためにガウシアンのラプラシアンにより、画像列にフィルタ処理を施す。
【０１９１】
２）高周波成分の特徴をさらに強調するために、非線形変換（シグマ関数）を適用する。新しい画像列は、導関数を計算するために使用可能な、関心特徴(interest feature)の連続した表現を与える。この手順は、離散的な、より容易に追跡可能な特徴を強調するが、他の画像領域からの情報を除外しない。この方法は、最良のフローを与える可能性が高い導関数を決定するためにさまざまな関心特徴の中から選択するよう、そのようなさまざまな関心特徴を扱うように拡張することができる。ＢＣＥが近似的に有効であることを保証するために、空間導関数が数列全体にわたって一致していることと、現に取扱っている画像の尺度において画像強度が平面により良く近似されていることが必要である。この具体化例において、アルゴリズムは単一の尺度で適用される。
【０１９２】
本発明による他の実施の形態において、上述のアルゴリズムは、他の場合には較正が固定され既知である変化し未知である焦点距離を有するカメラや、未知の形態で画像ごとに任意に変化するカメラ較正（これは射影の場合として知られている）を取り扱うことも可能である。
【０１９３】
変化する未知の焦点距離に対して、J. Oliensis,“A Multi-frame structure from Motion Algorithm under Perspective Projection（透視投影の下での運動アルゴリズムからのマルチフレーム構造）”, IJCV, 34:2/3, 163-192, 1999及びWorkshop on Visual Scenes, 77-84, 1999に説明された方法に沿って、このアルゴリズムを修正することが可能である。修正は次のように説明される。すなわち、過程２ａにおいて、Ｈを、Ψx，Ψy，Ψzに加えてΨF≡[ｐn・∇Ｉ(ｐn)]を零化すると定義する。Ｚ-1nを復元する場合に、式(8)のΨを［Ψ,ΨF］に置換しなければならず、Ωは４×３行列となる。修正された過程２ａは、焦点距離が未知であり変化する場合には正確に復元することが本質的に難しい定数（Ｚn＝１）成分を除いて、Ｚ-1を復元することができる。
【０１９４】
射影の場合に対しては、J. Oliensis and Y. Genc,“Fast Algorithms for Projective Multi-Frame Structure from Motion（運動からの射影マルチフレーム構造に対する高速アルゴリズム）”, ICCV, 536-543, 1999に説明されている方法に沿って、このアルゴリズムは修正することができる。この場合には、８個の長さＮpのベクトルを零化にするように、Ｈを定義する。ここで、これらのベクトルのｎ番目のエントリは、ａ，ｂ，ｃ，ｄ∈[ｘ,ｙ]に対して、８つの量(∇Ｉn)a，Ｐnb(∇Ｉn)c，Ｐnd(∇Ｉn)・Ｐnにより与えられる。過程２ａは、式(8)のΨを、これらの８個のベクトルからなる行列に置換することにより、修正されなければならない。
【０１９５】
［射影のアルゴリズム］
較正された画像列に対するここで説明したアルゴリズムは、未較正の画像列を扱うために直接的な方法で一般化される。この射影の場合において、回転計算の予備的段階の代わりに、基準画像と後続の各々の画像との間の平坦な相同を計算する。射影再構成法(projective-reconstruction method)のマルチフレーム一般化により、これらの相同が容易に正確に計算できることは、指摘する価値がある。追跡点の例について、これを簡潔に説明する。シーンが平坦であると仮定すると、
【０１９６】
【数７６】
0003651590.tif 000068
【０１９７】
であり、ここで、Ｍiは３×３行列（相同）であり、Ｓmは（ある任意の座標系において）平面上のｍ番目の点の位置を与える構造３−ベクトルであり、λimは射影深度である。相同の復元の手順は次のようになる。
【０１９８】
過程１）  初期値としてλim＝１とする。
【０１９９】
過程２）  λimの現在の推定値に対して、単一の３ＮI×Ｎp行列Γの中に
【０２００】
【外３４】
0003651590.tif 000069
【０２０１】
を収集(collect)する。２つのランク３因子の積にこれを分解するために、ＳＶＤを使用する。
【０２０２】
【数７７】
0003651590.tif 000070
【０２０３】
過程３）  現在のＭ(3)，Ｓ(3)に対して、｜Γ−Ｍ(3)Ｓ(3)｜を最小にするλimを計算する。過程２に戻る。
【０２０４】
収束した後、Ｍ(3)はＭi推定値を含んでいる。基準画像を画像ｉとする相同の推定値は、Ｍ0＼Ｍiである。この手順は｜Γ(4)｜2／｜Γ｜2の局所最小値に収束することが保証される。ここで、Γ(4)はΓとその最良のランク３近似の間の差である。さらに、真の動きがあまり大きくなく、シーンが真に平坦でれば、この手順はＭiのほぼ最適な推定値を与える。
【０２０５】
線に基づくアルゴリズムのうち未較正画像列用のものは、次のように説明される。Ｋiは、零下方対角線エントリ(zero lower-diagonal entries)を有する標準の３×３行列であって、ｉ番目の画像に対する較正パラメタを与えるものであるとする。未較正の場合に対しては、Ｋ及びＫ'を第１及び第２の画像に対する較正行列として、Ｒ→(Ｋ')-TＲＫ≡ＭH及びＴ→ＫＴと置換すれば、式(3)は依然として成立する。ＭHは一般的な相同であり、点に対する類似した相同の逆転置である。
【０２０６】
【数７８】
0003651590.tif 000071
【０２０７】
のこの新しいアルゴリズムでの１次近似は、
【０２０８】
【数７９】
0003651590.tif 000072
【０２０９】
である。Ｔ→ＫＴと置換され、ω×ＡはδＭHＡで置換されている、ここでＭH≡１＋δＭHである。
【０２１０】
較正された場合の回転フローの代わりに、δＭHによる１次の相同フローを零化するように、Ｈを定義する。Ｍ-THの１次近似は１−δＭTHであるから、点と線の両方に対して１次のフローを零化するするように、Ｈを容易に定義することができる。
【０２１１】
【数８０】
0003651590.tif 000073
【０２１２】
と表現する。ここでＦは２×２であり、Ｇ及びＪは２×１である。δＭHの８つのパラメタのそれぞれによる１次の点及び線の移動は、次の表で与えられる。
【０２１３】
【表１】
0003651590.tif 000074
【０２１４】
この表の列に関連した８つのベクトルを零化するように、Ｈを定義する。示された置換の後には、未較正用のアルゴリズムは較正用のアルゴリズムと同じである。ある基準画像内の点座標を固定すれば、射影変換のもとでの残りの自由度はＺ-1→ａｘ＋ｂｙ＋ｃ＋ｄＺ-1であり、ここでｘ及びｙは画像座標である。これはスケーリングに対応し、構造に任意の平面を加える。Ｂに対しても同様の関係を導くことができる。Ｐを、Ｂにより決定される線上の３次元の点であるとすれば、Ｂ・Ｐ＝−１は、Ｂ・(ｘ,ｙ,１)＝−Ｚ-1を意味する。射影変換は、Ｂ'・Ｐ'＝−１を維持しなければならない。このことと、上記により与えられるＺ-1の変換から、スケーリングによっては、点変換に対すると同じ定数ａ，ｂ，ｃを有するＢ'−Ｂ＝(ａ,ｂ,ｃ)という結果になる。この関係は、任意の線Ｂに対して同じ定数で成立する。
【０２１５】
全体としてシフトによってのみＢが復元可能であるという事実は、未較正用のアルゴリズムからも明確である。小さい相同の１次の効果を除外するようにＨを使用する場合、Ｂの３つの定数成分による平行移動の画像フローも除外される。したがってこれらの成分は復元することができない。
【０２１６】
［具体例］
本明細書において開示されている本発明による方法は、コンピュータを制御するために適切な記憶媒体に収容されているソフトウェアにより、具体化され、完全に実行が可能であることは、当業者に明白であろう。
【０２１７】
図１において、本発明を取り入れたコンピュータ・ハードウェア・システムを、ブロック図形式で例示している。図１に示すように、このシステムはビデオ・ソース１０１を有し、ビデオ・ソース１０１の出力はデジタイザ１０２によりピクセル・マップにデジタル化される。次に、デジタル化されたビデオ・フレームは、使用中のメイン・システム・メモリがアクセスするために、システムバス１０３を介して電子的な形態で記憶装置１０４に送られる。使用中は、このシステムの動作は、デジタル化されたピクセル・マップへのアクセス及び本発明による手順を制御する中央処理装置（ＣＰＵ）１０５により制御される。コンピュータ・ハードウェア・システムは、モニタ１０６及びグラフィックス・ボード１０７のような、データ及びグラフィックスにアクセスし表示するための、当業者には公知の標準的な構成要素を有する。
【０２１８】
利用者は、キーボード１０８及び／またはマウス１０９、あるいは画面上の項目を選択しあるいはシステムの機能を指示するために使用できるトラックボールのような、他の位置検知装置により、システムを対話的に操作する。
【０２１９】
本発明の方法に関連するキータスクの実行は、システムの、ＣＰＵにより制御されるメイン・メモリに格納された命令により指示される。ＣＰＵの動作を支配する格納された命令にしたがって、本発明による方法を実行するために、ＣＰＵはメイン・メモリにアクセスして必要なステップを実行することができる。具体的には、ＣＰＵは、利用者の入力にしたがって、格納されたデジタル化されたビデオにアクセスし、関連するデジタル化されたピクセル・マップから３次元構造及びカメラの動きの情報を抽出するために、本発明の方法を具体化する命令にしたがって、選択されたビデオ画像を解析する。
【０２２０】
ここで、図２のフローチャートを参照して、本発明よる方法を説明する。まず、ステップ２０１において、基準透視点(reference perspective)においてカメラにより画像列の基準画像を撮影し、カメラを平行移動及び／または回転することにより、１以上の連続した異なる透視点(perspective)において画像列の１以上の連続した画像を撮影する。次に、３次元画像内容すなわち点、線及び画像強度の分析のために、画像はデジタル化される（ステップ２０２）。次に、基準画像に対する、連続した各々の画像の画像データ・シフト(image data shift)を決定する（ステップ２０３）。シフトは、基準透視点から連続した異なる透視点へのカメラの平行移動及び／または回転から導出される。次に、各々の画像に対する画像データ・シフトを含むシフト・データ行列を構成する（ステップ２０４）。また、特異値分解（ＳＶＤ）を使用して、シフト・データ行列から、一方が３次元構造に対応（３次元構造行列）し他方がカメラの動きに対応する、２つのランク３の因子行列を計算する（ステップ２０５）。１次方程式を解くことにより、ＳＶＤを使用して、３次元構造行列から３次元構造を復元する（ステップ２０６）。復元された３次元構造を使用して、カメラの動きに対応する行列から、カメラの動きを復元することができる（ステップ２０７）。
【０２２１】
本発明による好適実施例と考えられるものを示し説明したが、本発明による技術思想から逸脱することなく、形式あるいは細部のさまざまな修正及び変更が容易に行えることは、当然理解されよう。したがって、本発明は説明し例示した正確な形式に限定されるものではなく、添付した特許請求の範囲内に含まれるすべての修正を含むように、構成されていると解釈される。
【図面の簡単な説明】
【図１】本発明の方法を実施するためのハードウェアの具体例を概略的に示すブロック図である。
【図２】本発明による方法の一例を示すフローチャートである。
【符号の説明】
１０１  ビデオ
１０２  デジタイザ
１０３  システムバス
１０４  大容量記憶装置
１０５  ＣＰＵ
１０６  画面ディスプレイ
１０７  グラフィックス・ボード
１０８  キーボード
１０９  位置検知装置（マウス）
２０１〜２０７    ステップ
【図１】
0003651590.tif 000117
【図２】
0003651590.tif 000118