３次元サーフェス生成方法
【課題】対象物体を撮影した複数の画像から物体表面の３次元サーフェスを直接かつ高速に生成できる３次元サーフェス生成方法を提供する。
【解決手段】 対象物体を異なる視点位置から撮影して得られた複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、基準画像上に２次元三角メッシュを生成し、メッシュの全頂点奥行きパラメータとカメラパラメータとによって定まる平面射影変換によって参照画像を変形して得られた画像の画素値を要素とするベクトルと、基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、複数の画像とカメラパラメータと全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、全頂点奥行きパラメータの微小変化量の計算と、全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、コスト関数の値が最小となる、全頂点奥行きパラメータを計算する。
【解決手段】 対象物体を異なる視点位置から撮影して得られた複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、基準画像上に２次元三角メッシュを生成し、メッシュの全頂点奥行きパラメータとカメラパラメータとによって定まる平面射影変換によって参照画像を変形して得られた画像の画素値を要素とするベクトルと、基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、複数の画像とカメラパラメータと全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、全頂点奥行きパラメータの微小変化量の計算と、全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、コスト関数の値が最小となる、全頂点奥行きパラメータを計算する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象物体を撮影した複数のステレオ画像から、対象物体の３次元サーフェスを生成する３次元サーフェス生成方法に関する。
【背景技術】
【０００２】
伝統的なステレオ法では、基準画像上の特定画素に注目し、その画素の対応点を参照画像上に探索することにより奥行き推定を行う。このような方法によって得られる対象物体の３次元情報は、３次元空間中の点群であり、対象物体の密な表面形状（以下、単に、「３次元サーフェス」又は「サーフェス」とも呼ぶ。）を直接的に表現したものではない。
【０００３】
つまり、伝統的なステレオ法では、画像上における隣接する画素に対応する３次元空間中の点の関係づけがなされていないという問題が存在している。
【０００４】
また、ステレオ画像を用いて対象物体の３次元サーフェスを生成するための１つのアプローチは、ステレオ法で得られた奥行きデータを利用して、即ち、対象物体表面を３次元空間中の点の集合として求めた後に、その３次元空間中の点を、サーフェスモデルに当てはめる方法である（例えば、非特許文献１、非特許文献２、非特許文献３、非特許文献４を参照）。このような手法は、ステレオ法に限らず、レーザレンジファインダなどによって得られた３次元情報を活用する技術として、旧来から盛んに研究されてきた。
【０００５】
しかし、伝統的なステレオ法によって得られる３次元空間中の点の情報は、画像ノイズ等の影響を大きく受けて誤差を含んでいるため、このような誤差の大きい３次元情報から対象物体のサーフェスを生成する方法には、適切なサーフェスを生成できなかったり、複雑な計算を必要としたりするなどの、さまざまな問題がある。また、これらの手法のほとんどは、もとの３次元情報を補間および変更することに基づいているため、ステレオ画像の輝度情報を有効に活用していないとの問題もある。
【非特許文献１】ディー. アール. マリー（D.R. Murray）著,「パッチレツ： ア メソッド オフ インタープレッティング コリレイション ステレオ ３Ｄ データ（Patchlets: A method of interpreting correlation stereo 3D data）」,ＰｈＤ スィーザス, デパートメント オフ コンピュータ サイエンス, ザ ユニバーシティー オフ ブリティシュ コロンビア, ２００３年（PhD thesis, Department of Computer Science, The University of British Columbia, 2003）
【非特許文献２】アイ. コヘン（I. Cohen）、エル. ディー. コヘン（L.D. Cohen）、エヌ. アヤチェ（N. Ayache）共著,「イントロジューイング ニュー デフォーマブル サーフェス トゥー セグメント ３Ｄ イメージ（Introducing new deformable surfaces to segment 3d images）」,ＩＥＥＥ コンピュータ ソサイアティ カンファレンス オン コンピュータ ビジョン アンド パターン レコグニション（IEEE Computer Society Conference on Computer Vision and Pattern Recognition）,p.738-739,1991年
【非特許文献３】エー. ペントランド（A. Pentland）、エス. スクラロフ（S. Sclaroff）共著,「クローズドフォーム ソリューションズ フォー フィジカルリ ベースド シェイプ モデリング アンド レコグニション（Closed-form solutions for physically based shape modeling and recognition）」,ＩＥＥＥ トランスアクションズ オン パターン アナリシス アンド マシン インテリジェンス（IEEE Transactions on Pattern Analysis and Machine Intelligence）,第13巻,p.715-729,1991年
【非特許文献４】アール. スゼリスキ（R. Szeliski）、ディー. トンネセン（D. Tonnesen）共著,「サーフェス モデリング ウィズ オリエンティド パーティクル システムズ（Surface modeling with oriented particle systems）」,コンピュータ グラフィックス（ＳＩＧＧＲＡＰＨ）（Computer Graphics(SIGGRAPH)）,p.185-194,1992年
【非特許文献５】ピー. ジェー. ナラヤナン（P.J. Narayanan）、ピー. ダブリュ. ランデル（P.W. Rander）、ティー. カナデ（T. Kanade）共著,「コンストラクティング バーチャル ワールドズ ユジング デンス ステレオ（Constructing Virtual Worlds using Dense Stereo）」,ＩＥＥＥ インターナショナル カンファレンス オン コンピュータ ビジョン（IEEE International Conference on Computer Vision）,p.3-10,1998年
【非特許文献６】ピー. ファ（P.Fua）著,「フロム マルチプル ステレオ ビューズ トゥー マルチプル ３Ｄ サーフェス（From multiple stereo views to multiple 3D surfaces）」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International journal of Computer Vision）,第24巻,第1号,p.19-35,1997年
【非特許文献７】オー. ディー. ファゲラス（O.D. Faugeras）、アール. ケリヴェン（R. Keriven）共著,「コンプリート デンス ステレオビジョン ユジング レベル セット メソッド（Complete Dense Stereovision Using Level Set Methods）」,ヨーロピアン カンファレンス オン コンピュータ ビジョン（European Conference on Computer Vision）,第I巻,p.379-393,1998年
【非特許文献８】エス. ベイカ（S. Baker）、アイ. マシューズ（I. Matthews）共著,「ルーカス・カナデ ２０ イヤーズ オン： ア ユニファイング フレームワーク（Lucas-kanade 20 years on: A unifying framework）」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International journal of Computer Vision）,第56巻,第3号,p.221-255,2004年
【非特許文献９】ビー. ディー. ルーカス（B.D. Lucas）、ティー. カナデ（T. Kanade）共著,「アン イテラティブ イメージ レジストレイション テクニック ウィズ アン アプローチ トゥー ステレオ ビジョン（An iterative image registration technique with an approach to stereo vision）」,イメージ アンダスタンディング ワークショップ（Image Understanding Workshop）,p.121-130,1981年
【非特許文献１０】オー. ファゲラス（O. Faugeras）、エフ. ラストマン（F. Lustman）共著,「モーション アンド ストラクチャ フロム モーション イン ア ピースワイズ プレーナ エンバイアロンメント（Motion and structure from motion in a piecewise planar environment）」,レポート ディ ルシェルシュ ディ ｌ‘ＩＮＲＩＡ（Report de Recherche de l’INRIA）,1988年
【発明の開示】
【発明が解決しようとする課題】
【０００６】
かかる問題を解決するために、ステレオ画像の輝度情報を利用して、対象物体の表面形状を直接的に生成するアプローチが提案されている。このような方法として、ステレオ画像からポリゴンメッシュを生成する手法（非特許文献５を参照）や、oriented particleを最適化する手法（非特許文献６を参照）などが挙げられる。
【０００７】
しかし、これらの手法は、彫像や人形などの静止物体の正確な３次元モデルを得ることを目的にしているため、一般に、数分程度の計算時間を必要とする。これは、ロボットや自動車などのナビゲーションを考えたときには大きな問題となる。例えば、移動するロボットでは、前方の坂の勾配や段差などの形状を速やかに得る必要があり、また、ロケットやヘリコプターの着陸時には、地上のサーフェスを速やかに取得する必要がある。
【０００８】
また、ステレオ画像から直接的に物体の３次元サーフェスを生成する方法として、レベルセット法と呼ばれる手法を用い、対象物体の３次元表面形状をステレオ画像から直接生成する手法がある（例えば、非特許文献７を参照）。この手法は、物体表面形状を関数によって表現して微分方程式を繰り返し計算によって解くアプローチである。しかし、この繰り返し計算は大きな計算量を必要とするため、短時間にて対象物体の３次元サーフェスを生成できないという問題がある。
【０００９】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、このような高速な処理を目的とした効率的な３次元サーフェス生成方法、即ち、対象物体を撮影した複数のステレオ画像から、対象物体の３次元形状を点の情報として求めるのではなく、物体表面の３次元サーフェスを直接かつ高速に生成する、３次元サーフェス生成方法を提供することにある。
【課題を解決するための手段】
【００１０】
本発明は、既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法に関し、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する第１のステップと、前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータと、カメラ視点位置を含むカメラパラメータとによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルと、前記基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、前記複数の画像と、前記カメラパラメータと、前記全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、前記全頂点奥行きパラメータの微小変化量の計算と、前記全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、前記コスト関数の値が最小となる、前記全頂点奥行きパラメータを計算する第２のステップとを有することにより、或いは、前記第２のステップをより短時間内で行うために、前記第２のステップでは、前記全頂点奥行きパラメータの現在値が定める平面射影変換と、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換との合成によって真の平面射影変換を表現し、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値を要素とするベクトルと、前記全頂点奥行きパラメータの現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルとによって、前記距離と同等の距離を表現することによって、繰り返しごとに必要な情報を１度だけ計算することによって効果的に達成される。
【００１１】
また、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する、２次元メッシュ生成ステップと、前記２次元メッシュ生成ステップで生成された２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差が小さくなるように、前記全頂点奥行きパラメータを最適化手法によって推定する、全頂点奥行き推定ステップとを有することにより、或いは、前記全頂点奥行き推定ステップは、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換と、前記全頂点奥行きパラメータの現在値が定める平面射影変換とに基づいて、真の平面射影変換が合成され、前記微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値と、前記現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値の差が小さくなるように、前記微小変化量を計算することによって、繰り返しごとに必要な情報を１度だけ計算することによってより効果的に達成される。
【００１２】
更に、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角パッチによって構成された２次元メッシュを生成するステップと、前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータの初期値を定めるステップと、定めた全頂点奥行きパラメータの初期値と、カメラ視点位置を含むカメラパラメータと、前記複数の画像とに基づいて、前記２次元メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、前記全頂点奥行きパラメータを求める最適化処理を行うステップとを有することにより、或いは、前記最適化処理では、前記基準画像、及び前記カメラパラメータを用いて、前記全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し処理に必要な情報を一度だけ計算する、繰り返し処理用情報計算ステップと、前記全頂点奥行きパラメータの現在値と、前記カメラパラメータとによって定まる各三角パッチの平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差分を計算する、差分値計算処理ステップと、前記繰り返し処理用情報計算ステップで得られた必要な情報と、前記差分値計算処理ステップで得られた差分値と、前記カメラパラメータと、前記全頂点奥行きパラメータの現在値を入力として、前記全頂点奥行きパラメータの微小変化量を計算する、微小変化量計算処理ステップと、前記微小変化量計算処理ステップで得られた微小変化量と、前記全頂点奥行きパラメータの現在値とに基づいて、前記全頂点奥行きパラメータの現在値を更新する、現在値更新処理ステップとを有することにより、或いは、前記最適化処理では、前記差分値計算処理ステップ、前記微小変化量計算処理ステップ及び前記現在値更新処理ステップを、所定の収束条件を満たすまで繰り返し、繰り返し処理により、前記収束条件を満たしたときに得られた全頂点奥行きパラメータの現在値を全頂点奥行きパラメータとすることにより、前記３次元サーフェスを生成することによってより効果的に達成される。
【発明の効果】
【００１３】
本発明に係る３次元サーフェス生成方法は、対象物体を撮影した複数のステレオ画像から、対象物体の３次元形状を点の情報として求めるのではなく、物体表面の３次元サーフェスを直接かつ高速に生成する方法である。
【００１４】
本発明の３次元サーフェス生成方法では、既知の視点位置を持つ校正済みのステレオカメラで対象物体を撮影して得られたステレオ画像から、対象物体の３次元サーフェスを直接に生成するようにしているので、ステレオ画像上に隣接する画素の３次元情報は、同一の三角パッチ中にあれば同一平面上にあり、そうでなければ、それぞれが稜線（もしくは頂点）を共有する２つの平面上にあるものとして関連付けられるため、明白となる。
【００１５】
本発明によれば、複数のステレオ画像から対象物体の３次元サーフェスを短時間内に効率的に生成できるという優れた効果を奏する。
【００１６】
更に、本発明を適用してステレオ画像から対象物体の３次元サーフェスを直接に生成する際に、所定のパラメータの値を固定することにより高速化処理を行うことで、３次元サーフェス生成における更なる高速化を実現することができる。
【００１７】
従って、例えばロボットや自動車などのナビゲーションにおいても、地上のサーフェスを速やかに取得する必要があるロケットやヘリコプターの着陸時に、前方の坂の勾配や段差などの形状を速やかに得る必要がある移動するロボットに、問題なく本発明を適用することができる。
【発明を実施するための最良の形態】
【００１８】
以下、本発明を実施するための最良の形態を、図面を参照して詳細に説明する。
【００１９】
本発明は、ステレオカメラで撮影された画像（ステレオ画像）を利用して、対象物体の表面形状を高速に推定できる、３次元サーフェス生成方法に関する。
【００２０】
つまり、本発明は、既知の視点位置を持つ校正済みのカメラを用いて、それぞれ異なる視点位置から対象物体を撮影した複数のステレオ画像から、物体表面の３次元サーフェスを直接かつ高速に生成できる、３次元サーフェス生成方法である。
＜１＞本発明の着眼点及び概説
図１に示すように、２台のステレオカメラで撮影された対象物体（以下、単に、「物体」とも言う。）の３次元表面形状は、複数の三角パッチで構成されるメッシュによって表現できる。
【００２１】
３次元空間中の各三角パッチは、２つの画像に投影されたとき、画像上に２次元の三角形を描き、３次元メッシュは画像上で２次元メッシュを描く。２枚のうちの一方の画像に投影されたメッシュの全頂点は、その画像を撮影したカメラの光学中心からの距離、即ち、「奥行き」を持つ。
【００２２】
したがって、全頂点の奥行きを推定することにより、物体の３次元サーフェスを生成することができる。
【００２３】
ここで、全頂点の奥行きを定める画像を「基準画像」と呼び、もう一方の画像を「参照画像」と呼ぶ。基準画像上の１つの三角パッチの領域内の座標は、その三角パッチの３頂点の奥行きと、既知のカメラ視点位置等の情報を含むカメラパラメータを用いて、参照画像の座標に変換することができる。この座標変換を「平面射影変換（homography）」と呼ぶ。
【００２４】
すなわち、その三角パッチの３頂点の奥行きが既知であれば、基準画像上のその三角パッチ内の画素値は、その画素の座標を平面射影変換した座標における参照画像上の画素値とほぼ一致すると想定できる。
【００２５】
よって、基準画像上のメッシュの全頂点の奥行きが既知であれば、全頂点の奥行きと、カメラパラメータとを用いて、基準画像上の全ての三角パッチの平面射影変換を計算することができる。
【００２６】
このとき、各三角パッチの平面射影変換によって基準画像上の全ての座標を変換した座標における参照画像の画素値を抽出し、基準画像の座標と同じ位置に挿入することによって参照画像を変形した画像を生成すると、その画像は基準画像とほぼ一致する。
【００２７】
そこで、全頂点の奥行きと、カメラパラメータとが定める各三角パッチの平面射影変換を用いて、参照画像を変形した画像の画素値と、基準画像の画素値との差が小さくなるように、全頂点の奥行きを計算する。このようにして、全頂点の奥行きを求めることにより、物体の３次元サーフェスが生成できる。
【００２８】
要するに、本発明では、基準画像を三角パッチによるメッシュ（ポリゴンメッシュ）に分割し、各三角パッチは３次元空間中で平面をなし、物体表面の形状をよく近似できていると仮定する。このとき、それぞれのパッチ領域では、ステレオ画像間において平面射影変換（homography）による座標変換が成立する。
【００２９】
校正済みのステレオカメラを利用することにより、その座標変換は、自由度３の平面パラメータを用いて表すことができ、さらに、その自由度は三角パッチの３頂点の奥行きに対応する。つまり、全パッチの頂点の奥行きを求めることにより、物体表面の３次元形状を得ることができる。
【００３０】
そこで、本発明では、頂点の奥行きを求める際に、基準画像上のパッチと、対応する参照画像上のパッチとのＳＳＤ（Sum of Squared Differences）を考え、全パッチについて加算したＳＳＳＤ（Sum ofＳＳＤ）を最小にする奥行きを、Gauss-Newton法を用いて最適化する。
【００３１】
更に、本発明では、パッチに関するステレオ画像間のＳＳＤを、ＩＣＩＡ（Inverse Compositional Image Alignment）を利用して表現することにより、繰り返しごとの計算コストを低く抑える。本発明では、パッチ頂点の画像上の位置を、基準画像上に固定するため、１つの頂点が持つ自由度は、１（奥行き）となる。
＜２＞平面パラメータの推定
＜１＞では、本発明の着眼点及び概説を簡潔に述べたが、ここからは、図面及び数式を参照しながら、本発明に係る３次元サーフェス生成方法を詳細に説明する。
【００３２】
ここで、まず、本発明において重要な意味をもつ「平面姿勢パラメータ（以下、単に、平面パラメータと言う）」の推定方法について述べる。
【００３３】
即ち、ここでは、３次元空間中にある１つの平面に注目し、その平面を一意に表す３自由度のパラメータ（平面パラメータ）を、ＩＣＩＡ(Inverse Compositional Image Alignment)（非特許文献８参照）の考え方を利用して高速に推定する方法について、詳細に述べる。
【００３４】
ここで述べる「平面パラメータ推定方法」は、後述する「＜３＞曲面の形状推定」にて述べる「頂点の奥行き推定」において重要な役割を果たす。
＜２−１＞ＩＣＩＡ（Inverse Compositional Image Alignment）
まず、ＩＣＩＡのアルゴリズムについて概説する。
【００３５】
基準画像座標を
、参照画像座標を
とし、２つの座標間の変換（warp）を
と表す。
【００３６】
ただし、
は座標変換を表す。また、
は「座標変換パラメータベクトル」（以下、単に、「変換パラメータ」とも呼ぶ。）を表し、座標間の変換が平行移動の場合に２次元のベクトルとなり、また、座標間の変換がアフィン変換および平面射影変換の場合に、それぞれ６次元および８次元のベクトルとなる。
【００３７】
基準画像をＴ、参照画像をＩとすると、変換パラメータ
は、下記数１で示す画像間のＳＳＤ(Sum of Squared Differences)を、最小にする値として得られる。
【００３８】
【数１】
ただし、
は、基準画像における注目領域内の画素について加算することを表す。
【００３９】
数１を最小化することにより、変換パラメータ
を得る方法は、Gauss-Newton法による最適化を利用したLucas-Kanade法（非特許文献９を参照）が広く知られている。
【００４０】
ＩＣＩＡ（非特許文献８を参照）とは、Lucas-Kanade法の改良版であり、Lucas-Kanade法の精度と安定性を損なうことなく、高速化を実現したアルゴリズムである。
【００４１】
いま、座標変換
は、結合について閉じており、単位変換と逆変換を持つ（すなわち、変換群をなす）ものとする。
【００４２】
すると、座標変換
は、下記数２のように結合変換で表すことができる。
【００４３】
【数２】
ただし、
は、
の現在値（初期値）であり、
は、微小要素を持つベクトルである。また、
は、単位写像を表すものとする。
【００４４】
ＩＣＩＡでは、上記数２を利用し、下記数３を最小化する
を求める。
【００４５】
【数３】
そして、
の規則により、座標変換パラメータを更新する。この操作を、
が収束するまで繰り返し、最終的な値を
の推定値とする。
【００４６】
ＩＣＩＡの高速性は、次のように説明することができる。
【００４７】
Lucas-Kanade法（非特許文献９を参照）では、数１に対し、
付近でTaylor展開する。この場合、
が繰り返し計算のたびに更新され、それに伴い
が変化することから、
における
の値を繰り返しのたびに計算する必要がある。
【００４８】
これに対し、ＩＣＩＡでは、数３において、
付近でTaylor展開する。この場合、微分対象は基準画像Ｔであるため変化せず、また、
の値は常に固定値
において評価されることから、
を再計算する必要がない。これにより、ヘッセ行列を固定化でき、Lucas-Kanade法よりも、はるかに高速な演算を実現できる。
【００４９】
このように、ＩＣＩＡは、Lucas-Kanade法と比較して、繰り返し計算ごとの計算コストを大幅に削減できるため、近年では、画像レジストレーションのリアルタイム処理を実現するための強力な手段となっている。
＜２−２＞ＩＣＩＡを利用した平面射影変換パラメータ推定
ここで、ＩＣＩＡを利用して、８自由度の平面射影変換パラメータを推定する方法について述べる。
【００５０】
３次元空間中の平面を撮影した２枚の画像間では、下記数４のように、３×３の平面射影変換行列
による座標変換が成立する。
【００５１】
【数４】
ただし、
および
は、それぞれ
および
の同次座標を表す。また、’〜’は同値関係を表している。
【００５２】
数４は、定数倍の不定性を持つため、３×３の平面射影変換行列
の自由度は８であり、ここでは、一番右下の要素（即ち、第３行第３列要素）を固定したものとする。
【００５３】
また、数４は、
の座標変換に対応している。
【００５４】
ただし、ベクトル
は、
の８個の未知数をラスタ順に並べたものとなる。
【００５５】
ＩＣＩＡを用いた平面射影変換行列
の推定（非特許文献８を参照）では、現在値（初期値）
と微小要素を持つ行列
を用いて、
を下記数５のように表す。
【００５６】
【数５】
すなわち、
および
の２つの数式は、それぞれ、
および
による座標変換に対応する。
【００５７】
数３を最小化する
は、
をTaylor展開し、数３を
で微分しゼロとおくことにより、下記数６のように書ける。
【００５８】
【数６】
ただし、下記数７が成立する。
【００５９】
【数７】
ここで、
は、基準画像Ｔの勾配を示す１×２の列ベクトルであり、下記数８と表される。
【００６０】
【数８】
また、
は、下記数９で示す２×８の行列である。
【００６１】
【数９】
ここで、ｅは、
と表され、座標
における基準画像Ｔの画素値と、
に対応する参照画像Ｉ上の画素値との画素値の差を意味する。ｅは、全ての座標
を用いて、参照画像Ｉを変形（warp）した画像
を生成しておき、
として求める方法が効率がよいため、一般的に利用される。
および
は、注目領域内の各画素ごとに異なる値を持つ。しかし、既に述べたように、その値は一度だけ計算すればよく、繰り返しのたびに再計算する必要はない。これに伴い、ヘッセ行列の逆行列
も一度だけ計算すればよい。
＜２−３＞ＩＣＩＡを利用した平面パラメータ推定
ここでは、前述したＩＣＩＡの考え方を利用して、３自由度の平面パラメータを高速に推定する方法について説明する。
【００６２】
本発明の実施形態を説明する際にしては、煩雑さを避けるため、以後、画像座標
は内部パラメータが正規化された座標（正規化画像座標）を示すものとする。
【００６３】
図２に示すように、３次元空間中の平面
は、平面パラメータ
を用いて、
で表される。
【００６４】
ただし、
は平面
上の点の３次元座標を表し、基準カメラ座標系にて定義されるものとする。ここでは、
を求めたいパラメータベクトルとする。
【００６５】
平面射影変換行列
は、
を用いて、
と表すことができる（非特許文献１０を参照）。ただし、
は、カメラ間の回転行列および平行移動ベクトルを表し、それぞれ既知とする。
【００６６】
このとき、画像の座標変換は、
をパラメータとして、
と書ける。
【００６７】
ＩＣＩＡを利用して
を求めるには、現在値
と微小変化量
を用いて、座標変換を結合写像
で表現し、下記数１０を最小にする
を求めればよい。
【００６８】
【数１０】
結合写像を表現するために、
を用いて、平面射影変換行列
を下記数１１のようにおく。
【００６９】
【数１１】
数１１を数５の形式に変形すると、
および
は、それぞれ下記数１２、数１３で表される。
【００７０】
【数１２】
【００７１】
【数１３】
このとき、数１２および数１３は、それぞれ
および
による座標変換に対応する。
【００７２】
数１３は、
の関数として記述したことを意味するため、数１０を最小にする
は、微分の連鎖法則を利用して、下記数１４のように得られる。
【００７３】
【数１４】
ただし、下記数１５、数１６及び数１７が成立する。
【００７４】
【数１５】
【００７５】
【数１６】
【００７６】
【数１７】
ここで、
である。また、
は、数６の場合と同様に、基準画像Ｔの勾配を示す１×２の列ベクトルである。そして、
は、下記数１８で示す２×９の行列である。
【００７７】
【数１８】
数９と数１８における
のサイズが異なるのは、
で表現された
には不定性が存在しないためである。
【００７８】
また、
にて評価した結果であり、行列
は、下記数１９で表される。
【００７９】
【数１９】
ただし、下記数２０、数２１、及び数２２が成立する。
【００８０】
【数２０】
【００８１】
【数２１】
【００８２】
【数２２】
ここで、{ｒ１，ｒ２，…，ｒ９}は、
の要素をラスタ順に並べたものであり、また、
である。
【００８３】
κは画像座標
に依存しないため、数１４を改めて書き直すと、下記数２３のようになる。
【００８４】
【数２３】
ただし、下記数２４、数２５及び数２６が成立する。
【００８５】
【数２４】
【００８６】
【数２５】
【００８７】
【数２６】
＜２−２＞で述べた平面射影変換パラメータの推定と、ここにて述べた平面パラメータ推定との大きな相違点は、前者では、
が変換群をなすのに対し、後者では、
が変換群をなさないことにある。
【００８８】
そのため、数１３が
に依存してしまい、
の値は繰り返し計算ごと変化する。しかし、その変化は、数２３におけるスカラーκに集約されており、
および
は変化しないことを考慮すると、変換群でないことによって生じる繰り返しごとの計算コストの増加は、無視できるほど小さい。
【００８９】
これにより、ＩＣＩＡを利用した平面パラメータ推定は、高速な演算にて実現できることが分かる。
＜３＞曲面の形状推定
＜２−３＞では、３次元空間中の１つの平面に注目し、その平面パラメータを高速に推定する方法について述べた。ここでは、＜２−３＞にて述べた方法を利用して、三角パッチで構成されたポリゴンメッシュ（三角メッシュ）の全頂点の奥行きを求める方法について述べる。
＜３−１＞三角メッシュによる表面形状モデル
まず、三角メッシュを用いた物体表面形状モデルについて述べる。
【００９０】
本発明では、基準画像を小さな三角形で分割してメッシュを生成する必要がある。
【００９１】
メッシュを生成する方法としては、例えば、画像全体を方形で分割した後に方形の対角線を結び三角形に再分割する方法や、画像上から特徴点を抽出した後に特徴点を母点としてドロネー分割する方法など、様々な方法が考えられる。
【００９２】
本発明では、画像領域全体（基準画像領域全体）が微小な三角形で分割されていれば、上記のメッシュ生成方法に限らず、どのような方法を用いてもよい。ここで、正三角形によって構成されたメッシュの例を図３に示す。
【００９３】
メッシュの頂点を
と表し、
は対象物体表面の３次元形状に対応した奥行きｄｍを持つものとする。また、頂点によって構成される三角形をパッチと呼び、Ｃｎ（ｎ＝１，２，…，Ｎ）と表す。各パッチは３次元空間中において平面を構成するものとする。
【００９４】
すなわち、パッチによって構成されるメッシュはポリゴンメッシュであり、メッシュが十分細かければ、任意の３次元形状を精度よく近似できると考えられる。
【００９５】
本発明では、メッシュの細かさに関する検討は特に行わない。ただし、画像レジストレーション分野では広く知られているように、粗いメッシュ（即ち、大きなメッシュ）を利用した場合には、より推定が安定になると考えられるため、＜３−３＞にて述べる全頂点の奥行き推定では、最初は粗いメッシュを用いて奥行きを推定し、徐々にメッシュを細かくする粗密探索のアプローチが有効である。
【００９６】
また、画像のＲＭＳＥの値を利用して細かい起伏を持つ領域を推定し、その領域に対応するパッチのみを細かくするなど、再分割（subdivision）の考え方を利用したパッチの分割方法も有効であると考えられる。
＜３−２＞三角パッチの３頂点に対する奥行き推定
次に、＜２−３＞にて述べた平面パラメータ推定方法を利用して、１つの三角パッチの３頂点の奥行きを求める方法について述べる。
【００９７】
ある三角パッチの３頂点の画像上の座標を
とし、それら頂点の奥行きをｄｍとする。このとき、当該三角パッチの３頂点の３次元空間中の座標
は、図４に示すように、
と表すことができる。
【００９８】
３次元空間中の三角パッチの平面パラメータを
とすると、ある三角パッチについて、平面パラメータ
と３頂点の奥行きｄｍ（ｍ＝ｉ,ｊ,ｋ）の関係は、下記数２７で表される。
【００９９】
【数２７】
ただし、下記数２８及び数２９が成立する。
【０１００】
【数２８】
【０１０１】
【数２９】
これより、平面パラメータ
と、奥行きパラメータγは比例関係にあり、その比例関係を決定する行列
は、画像上の三角パッチの位置によって一意に定まる。
【０１０２】
なお、「奥行き」と「奥行きパラメータ」の関係は、上記数２９に示す通りである。つまり、γは、ある三角パッチの３頂点の奥行きパラメータを表し、ｄｍ（ｍ＝ｉ,ｊ,ｋ）はそれら頂点の奥行きを表す。
【０１０３】
よって、数２３に微分の連鎖法則を適用することにより、平面パラメータ
の微小変化量
の代わりに、奥行きパラメータγの微小変化量を求めることができる。
【０１０４】
すなわち、γ０およびγΔを、それぞれ奥行きパラメータγの現在値および微小変化量とし、γΔ＝γ０＋γΔとすると、下記数３０で表される式を最小化するγΔは、数２３により、下記数３１で得られる。
【０１０５】
【数３０】
【０１０６】
【数３１】
ただし、下記数３２及び数３３が成立する。
【０１０７】
【数３２】
【０１０８】
【数３３】
ここで、
である。また、κは、数２７により平面パラメータの現在値
を求め、求めた
を数２４に代入することにより、得られる。そして、
および
は、数２３と同じものである。
【０１０９】
数２３と同様に、数３１においても、
は繰り返し計算のたびに変化しない。
【０１１０】
よって、平面パラメータ推定の場合とほぼ同等の計算コストで、奥行きパラメータγを推定することが可能となる。
＜３−３＞三角メッシュの全頂点の奥行き推定
ここで、三角メッシュの全頂点の奥行きを最適化する方法について説明する。
【０１１１】
＜３−２＞では、１つの三角パッチに注目したときの基準画像と参照画像間のＳＳＤを数３０によって表し、そのＳＳＤを最小にする３頂点の奥行きパラメータγを求めた。
【０１１２】
ここでは、３自由度の平面パラメータと３頂点の奥行きの関係を利用して、＜３−１＞にて述べた三角メッシュにおける全パッチについてＳＳＤを加算することによって得られるＳＳＳＤ(Sum of ＳＳＤ)を最小にすることにより、メッシュの全頂点の奥行きを求める。
【０１１３】
＜３−１＞にて述べたように、三角メッシュにおけるパッチ数をＮ、頂点の数をＭとする。また、全頂点の奥行きｄｍ，（ｍ＝１，２，…，Ｍ）の逆数を並べ、
とする。
【０１１４】
ここでは、
を「三角メッシュにおける全頂点の奥行きパラメータ」とし、以下、単に、「全頂点奥行きパラメータ」とも称する。
【０１１５】
また、
および
を、それぞれ全頂点奥行きパラメータ
の現在値および微小変化量とし、
とおく。
【０１１６】
そして、下記数３４に示すように、数３０を三角メッシュにおける全パッチについて加算したＳＳＳＤ(Sum ofＳＳＤ)を最小化する問題を考える。
【０１１７】
【数３４】
をTaylor展開し、数３４を
で微分してゼロとおくことにより、数３４を最小にする
は、下記数３５となる。
【０１１８】
【数３５】
ただし、下記数３６及び数３７が成立する。
【０１１９】
【数３６】
【０１２０】
【数３７】
ここで、κｎは数２４に準じたものであるが、平面パラメータの現在値
がパッチごと異なり、それに伴いκの値がパッチごと異なるため、添え字ｎを付与して表記している。
【０１２１】
また、
は数２８に準じているが、その要素はパッチの画像上の座標
に依存しており、パッチごと異なるため、κと同様に添え字ｎを付与した。
【０１２２】
そして、表１に、数３５の計算に必要な行列やベクトルを示している。
【０１２３】
【表１】
数３５により、全頂点奥行きパラメータの微小変化量
を求める場合は、＜２−３＞で述べた単一平面のパラメータ推定や、＜３−２＞で述べた三角パッチの３頂点の奥行き推定とは異なり、数３６の
に関する計算を繰り返しのたびに行う必要がある。
【０１２４】
ただし、
の計算に注目すると、行列
は、各パッチごと異なるものの、繰り返しごとに変化せず、さらに、この行列は実質的に３×３＝９個の要素しか持っていないため、全パッチに関して加算しても計算量は比較的小さい。
【０１２５】
また、
の計算は、実質的に画素数に依存しており、単一平面のパラメータ推定や３頂点の奥行き推定を、同じ画素数に関して計算した場合とほぼ同じ計算量となる。
【０１２６】
すなわち、数３５の計算量全体を評価すると、同じ画素数を用いて単一平面のパラメータ推定をした場合と比較して、
の逆行列計算の分だけ計算コストの増加と言える。
【０１２７】
ただし、
の行列サイズは、Ｍ×Ｍであり、三角メッシュにおいて頂点の数が多くなるほど、その計算量は大きくなる。本発明では、計算速度に関して、評価の高いCLapackライブラリを利用して逆行列計算を実装し、計算時間の短縮を実現している。
【０１２８】
ここで、上述したように、三角メッシュにおける全頂点の奥行き推定を行うことにより、対象物体の表面形状を推定する、「曲面形状推定」、すなわち、「微小平面モデルを用いた曲面形状推定」のアルゴリズムの一例を図５に示す。
【０１２９】
図５に例示されるアルゴリズムは、初期化処理と推定処理に分けられる。初期化処理では、推定処理にて利用する固定値（繰り返しごとに変化しない行列やベクトル）を計算している。また、推定処理では、数３５の計算に必要な値を各パッチについて求め、
の計算と
の更新を収束するまで繰り返す。
【０１３０】
なお、数３４は、基準画像１枚と参照画像１枚を利用した場合、即ち２枚の画像を利用した場合のＳＳＳＤを示している。参照画像を２枚以上利用する場合は、即ち基準画像を含めて３枚以上の画像を利用する場合は、数３４は、基準画像と参照画像のペアの数だけ存在する。そこで、数３４を基準画像と参照画像のペアの数だけ加算することによって得られる新たなＳＳＳＤを利用して、全頂点の奥行きパラメータを推定する。
＜３−４＞高速化処理
＜３−３＞で述べた「三角メッシュの全頂点の奥行き推定方法」において、三角メッシュにおける頂点数が非常に大きい場合は、数３５に必要な
の逆行列を得るために、繰り返しごとに必要な計算時間が長くなり、つまり、計算コストが高くなることがある。
【０１３１】
この逆行列の計算が必要となるのは、κｎが繰り返し計算ごとに異なることが原因である。そこで、計算時間をより短くする方法として、κｎを固定化する方法が考えられる。
【０１３２】
数２４から明らかなように、κｎの変化は、ｎ番目のパッチの平面パラメータ
が繰り返し計算ごとに変化することによって生じる。しかし、
の繰り返しごとの変化は、非常に小さいため、κｎを固定しても、計算結果に大きな違いは生じない。
【０１３３】
よって、特定の繰り返し回数の間は、全パッチにおけるκｎを固定化し、
の再計算を行わずに計算時間を短くすることが可能である。例えば、１〜５回、６〜１０回の繰り返しの間では、
を固定し、１回目と６回目だけ
とその逆行列を再計算する方法などが考えられる。
【０１３４】
さらに、もし、ｎ番目のパッチにおける平面パラメータ
の現在値が、
とほぼ直交している場合、もしくは、
で表される平面と基準カメラ座標原点との距離が大きければ、数２４により
となることが分かる。
【０１３５】
そこで、はじめからκｎ＝−１,（ｎ＝１,２,…,Ｎ）として、対象物体の表面モデルを生成することも可能である。この場合、１度だけ
を求めればよく、計算時間をさらに短縮することができる。
＜３−５＞本発明による３次元サーフェスを生成するための処理の流れ
図６に、本発明を適用して、対象物体を撮影した複数のステレオ画像から、対象物体の３次元サーフェスを生成するための基本的な処理の流れを示す。
【０１３６】
なお、複数のステレオ画像のうち、一枚のみを基準画像とし、残りの画像を参照画像とする。また、ステレオ画像を撮影した多眼ステレオカメラについて、カメラのレンズ歪や焦点距離などの内部パラメータと、ステレオカメラ間の配置を示す外部パラメータとをあわせて、「カメラパラメータ」としている。
【０１３７】
説明を簡単にするために、以下、２眼ステレオカメラを用いた場合について説明する。つまり、撮影したステレオ画像が２枚で、１枚は基準画像とし、もう１枚は参照画像とする。
【０１３８】
図６に示されるように、この処理における入力は、複数のステレオ画像と、ステレオカメラのカメラパラメータである。
【０１３９】
即ち、既知の視点位置にある校正済みの多眼ステレオカメラを用いて、対象物体を撮影して、ステレオ画像を取得すれば、本発明を適用して対象物体の３次元サーフェスを直接かつ高速に生成することができる。
【０１４０】
処理の手順として、図６に示されるように、まず、所定の方法（＜３−１＞の記述を参照）により、基準画像上に２次元三角メッシュ（三角メッシュ）を生成する。
【０１４１】
次に、生成された三角メッシュの全頂点に対して、奥行きパラメータの初期値（全頂点奥行きパラメータの初期値）を定める。
【０１４２】
そして、定めた全頂点奥行きパラメータの初期値と、カメラパラメータと、基準画像及び参照画像で構成される複数のステレオ画像（２眼ステレオカメアを用いた場合に、一枚の基準画像と一枚の参照画像で構成する２枚のステレオ画像になる。）とを入力とし、基準画像上のパッチと対応する参照画像上のパッチとのＳＳＤを考え、三角メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、全頂点奥行きパラメータを求める最適化処理（全頂点奥行きパラメータ最適化処理）を行う。
【０１４３】
即ち、全パッチについて加算したＳＳＳＤを最小にする、最適な全頂点奥行きパラメータは、全頂点奥行きパラメータの初期値とカメラパラメータとを入力とする最適化手法を用いて、求めることができる。最適化手法では、全頂点奥行きパラメータの微小変化量の計算と、全頂点奥行きパラメータの現在値の更新とを繰り返し行う必要がある。
【０１４４】
最適化計算における微小変化量の計算は、コスト関数を微小変化量で表現した後に、コスト関数を微小変化量で微分してその結果をゼロとおくことにより、コスト関数の極値（極小値）となるパラメータを求める方程式をたてる。このとき、その方程式を解くために必要な情報は、繰り返しごとに求める必要があるため、一般に短時間で計算を行うことは困難である。
【０１４５】
しかし、＜３−４＞高速化処理でも述べたように、全頂点奥行きパラメータの微小変化量の計算において、必要な情報を固定化することができれば、これらの情報は繰り返し毎に計算する必要はなくなり、一度だけ計算すればよいことになり、極めて短時間で計算することが可能となる。
【０１４６】
全頂点の真の奥行きパラメータが定める平面射影変換は、全頂点奥行きパラメータの現在値が定める平面射影変換と、全頂点奥行きパラメータの微小変化量が定める平面射影変換の合成によって表すことができる。そして、全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって基準画像を変形して得られた画像の画素値と、全頂点奥行きパラメータの現在値が定める平面射影変換によって参照画像を変形して得られた画像の画素値との差を用いて、コスト関数を表現すると、繰り返しごとに必要な多くの情報は固定化でき、一度だけ計算すればよいことになる。これにより、より短時間での計算を実現できる。
【０１４７】
図７は、図６に示した最適化処理（全頂点奥行きパラメータ最適化処理）の基本フローである。
【０１４８】
図７に示されるように、最適化処理（全頂点奥行きパラメータ最適化処理）では、まず、基準画像及び基準画像を撮影した基準カメラのカメラパラメータを用いて、全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し計算に必要な情報を一度だけ計算する。
【０１４９】
次に、全頂点奥行きパラメータの現在値と、カメラ視点位置を含むカメラパラメータとによって定まる各三角形（各三角パッチ）の平面射影変換によって参照画像を変形して得られた画像の画素値と、基準画像の画素値との差分を計算する。以下、この処理を「差分値計算処理」とも称する。
【０１５０】
差分値が得られた後に、全頂点奥行きパラメータの微小変化量を計算する。つまり、図７に示されるように、繰り返し計算に必要な情報と、差分値と、カメラパラメータと、全頂点奥行きパラメータの現在値を入力として、全頂点奥行きパラメータの微小変化量を計算する。以下、この処理を「全頂点奥行きパラメータの微小変化量計算処理」とも称する。
【０１５１】
そして、得られた全頂点奥行きパラメータの微小変化量と、全頂点奥行きパラメータの現在値とに基づいて、全頂点奥行きパラメータの現在値を更新する。以下、この処理を「全頂点奥行きパラメータの現在値更新処理」とも称する。
【０１５２】
「最適化処理」では、「差分値計算処理」、「全頂点奥行きパラメータの微小変化量計算処理」及び「全頂点奥行きパラメータの現在値更新処理」を、所定の収束条件を満たすまで繰り返す。繰り返し処理（繰り返し計算）により、所定の収束条件を満たしたときに最終的に得られた全頂点奥行きパラメータの現在値を「全頂点奥行きパラメータ」とし、よって、三角メッシュにおける全頂点の奥行きが定められることになり、対象物体の３次元サーフェスを生成できる訳である。
＜４＞実験結果
以下では、本発明の有効性を確認するために行った実験結果を示す。まず、合成画像を利用して本発明を適用した様々なシミュレーション実験の結果について示し、次に、実画像を用いて本発明を適用した実画像実験の結果を示す。
＜４−１＞シミュレーション実験
シミュレーション実験（合成画像を用いた実験）では、５１２［pix］×５１２［pix］の画像を利用し、また、ステレオカメラの設定パラメータは、以下に示す値を利用した。
【０１５３】
カメラの内部パラメータは、基準カメラ及び参照カメラとも同じであり、下記数３８の
で示す。
【０１５４】
【数３８】
また、カメラの外部パラメータは、
（ここで、
は単位行列である）、
とした。ここで、
は、カメラ間の回転行列と平行移動ベクトルをそれぞれ示している。
【０１５５】
図８には、シミュレーション実験用に作成した合成画像を示す。本発明では、推定時に使用するメッシュを生成するのに、即ち、基準画像上に三角メッシュを生成するのに、前述したように、例えば、固定のメッシュを利用する方法や、抽出された画像特徴点を母点としたドロネー分割などの方法を利用することができるが、シミュレーション実験では、図８に示すように、５０［pix］の辺を持つ正三角形を用いてメッシュを作成し、テクスチャを張った平面を１０ｍ先に配置した。
【０１５６】
よって、図８の合成画像を用いた場合、頂点の奥行きは全て１０ｍであり、全パッチの平面パラメータは、
である。
【０１５７】
図８（Ａ）及び図８（Ｂ）には、作成された参照画像と基準画像をそれぞれ示している。図８（Ａ）及び図８（Ｂ）の画像は、５１２［pix］×５１２［pix］の画像における中央部分の４２０［pix］×４２０［pix］領域を表示したものである。このようにして作成した参照画像及び基準画像に対し、図８（Ｃ）に示したメッシュの頂点に関して、奥行き推定を行った。
【０１５８】
図８（Ｃ）におけるメッシュの頂点数は６１で、パッチ数は９６である。奥行き推定に利用した初期値は、全頂点において、１５ｍとした。
【０１５９】
本実験では、収束条件を特に設定せず、３０回まで繰り返し計算を行った。プログラムはＣ言語で実装し、逆行列計算にはCLapackライブラリを利用した。
【０１６０】
図８（Ｃ）で示したメッシュを利用した場合、Pentium（登録商標）ＩＶ(２.８ＧＨｚ）でLinux（登録商標）ＯＳを備えたパソコン上において、３０回の繰り返し計算に要する時間は、約０.９秒であった。
【０１６１】
メッシュが収束してゆく様子を図９に示す。図９では、基準画像上のメッシュを赤線で示しており、また、真値から計算される参照画像上の頂点を、推定された奥行きによって得られる平面射影変換を用いて基準画像上に投影し、その点によって構成されるメッシュを緑線で示している。
【０１６２】
よって、正しい奥行きが推定されれば、赤線のメッシュと緑線のメッシュは、完全に一致する。図９から、奥行き推定を繰り返すことにより、赤線のメッシュと緑線のメッシュが一致してゆくことがよく分かる。
【０１６３】
また、図１０には、画像のＲＭＳＥ（Root Mean Squared Error）と奥行きについて、繰り返し計算ごとの値の変化を示している。即ち、図１０の横軸は、最適化計算（最適化処理）における繰り返し回数を表している。
【０１６４】
図１０（Ａ）に、いくつかのパッチに関する画像ＲＭＳＥの変化と、全パッチについての画像ＲＭＳＥの変化（＊付きの実線）を示す。図１０（Ａ）から、３０回の繰り返し計算までに、適切に収束していることが分かる。
【０１６５】
そして、図１０（Ｂ）は、いくつかの頂点の奥行きの変化を示している。図１０（Ｂ）から、最終的に、真値（１０ｍ）に収束していることが分かる。
【０１６６】
次に、対象を球体として生成した合成画像を図１１に示す。図１１では、光軸上の１５ｍ先に中心を持ち、半径約７ｍの球体にテクスチャを貼り付けた。図１１（Ａ）及び図１１（Ｂ）は、作成された参照画像と基準画像を示している。図１１の画像は、図８の場合と同様に、５１２［pix］×５１２［pix］の画像における中央部分の４２０［pix］×４２０［pix］領域を表示したものである。このようにして作成した参照画像及び基準画像に対し、図１１（Ｃ）に示したメッシュの頂点に関して、奥行き推定を行った。
【０１６７】
図８の場合と同様に、図１１（Ｃ）におけるメッシュの頂点数は６１で、パッチ数は９６である。奥行き推定に利用した初期値は、全頂点において、１０ｍとした。
【０１６８】
図１２に、図１１の物体の実形状と、推定された形状を示す。図１２により、奥行き推定が適切に行われていることがよく分かる。なお、図１２（Ａ）、図１２（Ｂ）、図１２（Ｃ）及び図１２（Ｄ）は、全て同一視点から見たものである。
【０１６９】
また、繰り返し計算における画像ＲＭＳＥの変化、および、各頂点の真の奥行きと推定された奥行きとの絶対誤差（以下、単に、奥行き推定誤差の絶対値と呼ぶ。）を、図１３に示す。図１３（Ａ）では、いくつかのパッチにおける画像ＲＭＳＥと、全パッチに関する画像ＲＭＳＥを示し、図１３（Ｂ）では、いくつかの頂点における奥行き推定誤差の絶対値と、全頂点に関する奥行きのＲＭＳＥを示している。図１３（Ａ）及び図１３（Ｂ）から、ほぼ２０回の繰り返し計算にて、物体の形状が正しく推定できていることが分かる。
【０１７０】
更に、球体を用いた実験において、κｎ＝−１,（ｎ＝１,２,…,Ｎ）と固定した場合における繰り返しごとの画像ＲＭＳＥの変化と奥行きの絶対誤差の変化を図１４に示す。
【０１７１】
図１４（Ａ）及び図１４（Ｂ）は、図１３（Ａ）及び図１３（Ｂ）と、同様である。繰り返しごとの画像ＲＭＳＥ、および奥行き推定誤差の絶対値は、図１３（Ａ）及び図１３（Ｂ）のものとはやや異なるものの、ほぼ２０回の繰り返し計算にて、物体の形状が正しく推定できていることが分かる。
【０１７２】
すなわち、κｎ＝−１と固定して計算を行っても、κｎの値を固定しない場合との違いは、僅かである。よって、κｎの値を固定することにより、計算時間を一層短縮できることが分かる。
＜４−２＞実画像実験
実画像実験（実画像を用いた実験）に利用された実画像を図１５に示す。図１５の画像は、月面を模擬した１２８０［pix］×９６０［pix］のステレオ画像である。
【０１７３】
本実験では、粗いメッシュを徐々に細かくしてゆく粗密探索的な方法を用いて、表面モデルを生成した。
【０１７４】
ここでは、一辺の長さが４６４［pix］の三角パッチから始め、次の推定では、パッチのサイズを半分にするとともに、前回の頂点の奥行きを初期値とした推定を行った。
【０１７５】
図１５（Ａ）の基準画像では、最後にパッチのサイズとして使用した２９［pix］のメッシュを示している。このとき、メッシュの頂点数は８１７で、パッチ数は１５３６であった。
【０１７６】
図１５（Ａ）の基準画像及び、図１５（Ｂ）の参照画像に基づいて、本発明を適用して生成された表面モデルを図１６に示す。
【０１７７】
図１６には、図１５の実画像を用いた実験において、別々の視点から見た、本発明を適用して推定された３次元形状を示している。図１６（Ａ）、図１６（Ｃ）及び図１６（Ｅ）は、メッシュとテクスチャを同時に表示したもので、図１６（Ｂ）、図１６（Ｄ）及び図１６（Ｆ）は、テクスチャのみを表示したものである。
【０１７８】
図１６から、メッシュの端に位置する頂点は、その外側にメッシュが存在しないため、他のメッシュによる十分な拘束を受けることができず、推定が不安定になる傾向が見られる。しかし、メッシュ内部の領域の推定は安定であり、メッシュ中央部に見られるエッジを持った勾配や、左上部にみられる小さなクレータの穴が良好に復元されている。
以上のように、合成画像を用いた実験、及び実画像を用いた実験を通じて、本発明の有効性、即ち、本発明を適用することにより、効率的な表面形状推定が可能であることを確認した。
【０１７９】
また、上記の実験では、２台のステレオカメラを利用して得られた２枚のステレオ画像に基づいて、本発明を適用し、対象物体の３次元サーフェスを生成したが、本発明はそれに限定されず、３台以上のステレオカメラを利用して得られた３枚以上のステレオ画像に基づいても、本発明を適用し、対象物体の３次元サーフェスを生成することができる。
【図面の簡単な説明】
【０１８０】
【図１】対象物体とその３次元形状を説明するための模式図である。
【図２】３次元空間中の平面を説明するための模式図である。
【図３】基準画像上に生成した三角メッシュの例を示す図である。
【図４】三角パッチの頂点の３次元座標を説明するための模式図である。
【図５】本発明において、曲面形状推定のアルゴリズムを示すプログラムフローである。
【図６】本発明を適用して３次元サーフェスを生成するための基本的な処理の流れを示す図である。
【図７】図６に示した最適化処理の基本フローである。
【図８】シミュレーション実験用に作成した合成画像を示す図である。図８（Ａ）に参照画像を、図８（Ｂ）に基準画像を、図８（Ｃ）に基準画像上に生成されたメッシュ（推定時使用メッシュ）をそれぞれ示している。
【図９】図８の合成画像を用いた実験におけるメッシュが収束してゆく様子を説明するための図である。
【図１０】図８の合成画像を用いた実験において、画像のＲＭＳＥと奥行きについて、繰り返し計算ごとの値の変化を示すグラフである。
【図１１】シミュレーション実験用に対象を球体として生成された合成画像を示す図である。図１１（Ａ）に参照画像を、図１１（Ｂ）に基準画像を、図１１（Ｃ）に基準画像上に生成されたメッシュ（推定時使用メッシュ）をそれぞれ示している。
【図１２】図１１の合成画像を用いた実験における実形状と、本発明を適用して推定された形状を示す図である。
【図１３】図１１の合成画像を用いた実験において、画像のＲＭＳＥと奥行き推定誤差の絶対値について、繰り返し計算ごとの値の変化を示すグラフである。
【図１４】図１１の合成画像を用いた実験において、κｎの値を「−１」と固定した場合に、画像のＲＭＳＥと奥行き推定誤差の絶対値について、繰り返し計算ごとの値の変化を示すグラフである。
【図１５】実画像実験に利用された実画像（ステレオ画像）を示す図である。図１５（Ａ）に基準画像を、図１５（Ｂ）に参照画像をそれぞれ示している。
【図１６】図１５の実画像を用いた実験において、本発明を適用して推定された３次元形状を示す図である。
【技術分野】
【０００１】
本発明は、対象物体を撮影した複数のステレオ画像から、対象物体の３次元サーフェスを生成する３次元サーフェス生成方法に関する。
【背景技術】
【０００２】
伝統的なステレオ法では、基準画像上の特定画素に注目し、その画素の対応点を参照画像上に探索することにより奥行き推定を行う。このような方法によって得られる対象物体の３次元情報は、３次元空間中の点群であり、対象物体の密な表面形状（以下、単に、「３次元サーフェス」又は「サーフェス」とも呼ぶ。）を直接的に表現したものではない。
【０００３】
つまり、伝統的なステレオ法では、画像上における隣接する画素に対応する３次元空間中の点の関係づけがなされていないという問題が存在している。
【０００４】
また、ステレオ画像を用いて対象物体の３次元サーフェスを生成するための１つのアプローチは、ステレオ法で得られた奥行きデータを利用して、即ち、対象物体表面を３次元空間中の点の集合として求めた後に、その３次元空間中の点を、サーフェスモデルに当てはめる方法である（例えば、非特許文献１、非特許文献２、非特許文献３、非特許文献４を参照）。このような手法は、ステレオ法に限らず、レーザレンジファインダなどによって得られた３次元情報を活用する技術として、旧来から盛んに研究されてきた。
【０００５】
しかし、伝統的なステレオ法によって得られる３次元空間中の点の情報は、画像ノイズ等の影響を大きく受けて誤差を含んでいるため、このような誤差の大きい３次元情報から対象物体のサーフェスを生成する方法には、適切なサーフェスを生成できなかったり、複雑な計算を必要としたりするなどの、さまざまな問題がある。また、これらの手法のほとんどは、もとの３次元情報を補間および変更することに基づいているため、ステレオ画像の輝度情報を有効に活用していないとの問題もある。
【非特許文献１】ディー. アール. マリー（D.R. Murray）著,「パッチレツ： ア メソッド オフ インタープレッティング コリレイション ステレオ ３Ｄ データ（Patchlets: A method of interpreting correlation stereo 3D data）」,ＰｈＤ スィーザス, デパートメント オフ コンピュータ サイエンス, ザ ユニバーシティー オフ ブリティシュ コロンビア, ２００３年（PhD thesis, Department of Computer Science, The University of British Columbia, 2003）
【非特許文献２】アイ. コヘン（I. Cohen）、エル. ディー. コヘン（L.D. Cohen）、エヌ. アヤチェ（N. Ayache）共著,「イントロジューイング ニュー デフォーマブル サーフェス トゥー セグメント ３Ｄ イメージ（Introducing new deformable surfaces to segment 3d images）」,ＩＥＥＥ コンピュータ ソサイアティ カンファレンス オン コンピュータ ビジョン アンド パターン レコグニション（IEEE Computer Society Conference on Computer Vision and Pattern Recognition）,p.738-739,1991年
【非特許文献３】エー. ペントランド（A. Pentland）、エス. スクラロフ（S. Sclaroff）共著,「クローズドフォーム ソリューションズ フォー フィジカルリ ベースド シェイプ モデリング アンド レコグニション（Closed-form solutions for physically based shape modeling and recognition）」,ＩＥＥＥ トランスアクションズ オン パターン アナリシス アンド マシン インテリジェンス（IEEE Transactions on Pattern Analysis and Machine Intelligence）,第13巻,p.715-729,1991年
【非特許文献４】アール. スゼリスキ（R. Szeliski）、ディー. トンネセン（D. Tonnesen）共著,「サーフェス モデリング ウィズ オリエンティド パーティクル システムズ（Surface modeling with oriented particle systems）」,コンピュータ グラフィックス（ＳＩＧＧＲＡＰＨ）（Computer Graphics(SIGGRAPH)）,p.185-194,1992年
【非特許文献５】ピー. ジェー. ナラヤナン（P.J. Narayanan）、ピー. ダブリュ. ランデル（P.W. Rander）、ティー. カナデ（T. Kanade）共著,「コンストラクティング バーチャル ワールドズ ユジング デンス ステレオ（Constructing Virtual Worlds using Dense Stereo）」,ＩＥＥＥ インターナショナル カンファレンス オン コンピュータ ビジョン（IEEE International Conference on Computer Vision）,p.3-10,1998年
【非特許文献６】ピー. ファ（P.Fua）著,「フロム マルチプル ステレオ ビューズ トゥー マルチプル ３Ｄ サーフェス（From multiple stereo views to multiple 3D surfaces）」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International journal of Computer Vision）,第24巻,第1号,p.19-35,1997年
【非特許文献７】オー. ディー. ファゲラス（O.D. Faugeras）、アール. ケリヴェン（R. Keriven）共著,「コンプリート デンス ステレオビジョン ユジング レベル セット メソッド（Complete Dense Stereovision Using Level Set Methods）」,ヨーロピアン カンファレンス オン コンピュータ ビジョン（European Conference on Computer Vision）,第I巻,p.379-393,1998年
【非特許文献８】エス. ベイカ（S. Baker）、アイ. マシューズ（I. Matthews）共著,「ルーカス・カナデ ２０ イヤーズ オン： ア ユニファイング フレームワーク（Lucas-kanade 20 years on: A unifying framework）」,インターナショナル ジャーナル オフ コンピュータ ビジョン（International journal of Computer Vision）,第56巻,第3号,p.221-255,2004年
【非特許文献９】ビー. ディー. ルーカス（B.D. Lucas）、ティー. カナデ（T. Kanade）共著,「アン イテラティブ イメージ レジストレイション テクニック ウィズ アン アプローチ トゥー ステレオ ビジョン（An iterative image registration technique with an approach to stereo vision）」,イメージ アンダスタンディング ワークショップ（Image Understanding Workshop）,p.121-130,1981年
【非特許文献１０】オー. ファゲラス（O. Faugeras）、エフ. ラストマン（F. Lustman）共著,「モーション アンド ストラクチャ フロム モーション イン ア ピースワイズ プレーナ エンバイアロンメント（Motion and structure from motion in a piecewise planar environment）」,レポート ディ ルシェルシュ ディ ｌ‘ＩＮＲＩＡ（Report de Recherche de l’INRIA）,1988年
【発明の開示】
【発明が解決しようとする課題】
【０００６】
かかる問題を解決するために、ステレオ画像の輝度情報を利用して、対象物体の表面形状を直接的に生成するアプローチが提案されている。このような方法として、ステレオ画像からポリゴンメッシュを生成する手法（非特許文献５を参照）や、oriented particleを最適化する手法（非特許文献６を参照）などが挙げられる。
【０００７】
しかし、これらの手法は、彫像や人形などの静止物体の正確な３次元モデルを得ることを目的にしているため、一般に、数分程度の計算時間を必要とする。これは、ロボットや自動車などのナビゲーションを考えたときには大きな問題となる。例えば、移動するロボットでは、前方の坂の勾配や段差などの形状を速やかに得る必要があり、また、ロケットやヘリコプターの着陸時には、地上のサーフェスを速やかに取得する必要がある。
【０００８】
また、ステレオ画像から直接的に物体の３次元サーフェスを生成する方法として、レベルセット法と呼ばれる手法を用い、対象物体の３次元表面形状をステレオ画像から直接生成する手法がある（例えば、非特許文献７を参照）。この手法は、物体表面形状を関数によって表現して微分方程式を繰り返し計算によって解くアプローチである。しかし、この繰り返し計算は大きな計算量を必要とするため、短時間にて対象物体の３次元サーフェスを生成できないという問題がある。
【０００９】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、このような高速な処理を目的とした効率的な３次元サーフェス生成方法、即ち、対象物体を撮影した複数のステレオ画像から、対象物体の３次元形状を点の情報として求めるのではなく、物体表面の３次元サーフェスを直接かつ高速に生成する、３次元サーフェス生成方法を提供することにある。
【課題を解決するための手段】
【００１０】
本発明は、既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法に関し、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する第１のステップと、前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータと、カメラ視点位置を含むカメラパラメータとによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルと、前記基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、前記複数の画像と、前記カメラパラメータと、前記全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、前記全頂点奥行きパラメータの微小変化量の計算と、前記全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、前記コスト関数の値が最小となる、前記全頂点奥行きパラメータを計算する第２のステップとを有することにより、或いは、前記第２のステップをより短時間内で行うために、前記第２のステップでは、前記全頂点奥行きパラメータの現在値が定める平面射影変換と、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換との合成によって真の平面射影変換を表現し、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値を要素とするベクトルと、前記全頂点奥行きパラメータの現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルとによって、前記距離と同等の距離を表現することによって、繰り返しごとに必要な情報を１度だけ計算することによって効果的に達成される。
【００１１】
また、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する、２次元メッシュ生成ステップと、前記２次元メッシュ生成ステップで生成された２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差が小さくなるように、前記全頂点奥行きパラメータを最適化手法によって推定する、全頂点奥行き推定ステップとを有することにより、或いは、前記全頂点奥行き推定ステップは、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換と、前記全頂点奥行きパラメータの現在値が定める平面射影変換とに基づいて、真の平面射影変換が合成され、前記微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値と、前記現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値の差が小さくなるように、前記微小変化量を計算することによって、繰り返しごとに必要な情報を１度だけ計算することによってより効果的に達成される。
【００１２】
更に、本発明の上記目的は、前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角パッチによって構成された２次元メッシュを生成するステップと、前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータの初期値を定めるステップと、定めた全頂点奥行きパラメータの初期値と、カメラ視点位置を含むカメラパラメータと、前記複数の画像とに基づいて、前記２次元メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、前記全頂点奥行きパラメータを求める最適化処理を行うステップとを有することにより、或いは、前記最適化処理では、前記基準画像、及び前記カメラパラメータを用いて、前記全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し処理に必要な情報を一度だけ計算する、繰り返し処理用情報計算ステップと、前記全頂点奥行きパラメータの現在値と、前記カメラパラメータとによって定まる各三角パッチの平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差分を計算する、差分値計算処理ステップと、前記繰り返し処理用情報計算ステップで得られた必要な情報と、前記差分値計算処理ステップで得られた差分値と、前記カメラパラメータと、前記全頂点奥行きパラメータの現在値を入力として、前記全頂点奥行きパラメータの微小変化量を計算する、微小変化量計算処理ステップと、前記微小変化量計算処理ステップで得られた微小変化量と、前記全頂点奥行きパラメータの現在値とに基づいて、前記全頂点奥行きパラメータの現在値を更新する、現在値更新処理ステップとを有することにより、或いは、前記最適化処理では、前記差分値計算処理ステップ、前記微小変化量計算処理ステップ及び前記現在値更新処理ステップを、所定の収束条件を満たすまで繰り返し、繰り返し処理により、前記収束条件を満たしたときに得られた全頂点奥行きパラメータの現在値を全頂点奥行きパラメータとすることにより、前記３次元サーフェスを生成することによってより効果的に達成される。
【発明の効果】
【００１３】
本発明に係る３次元サーフェス生成方法は、対象物体を撮影した複数のステレオ画像から、対象物体の３次元形状を点の情報として求めるのではなく、物体表面の３次元サーフェスを直接かつ高速に生成する方法である。
【００１４】
本発明の３次元サーフェス生成方法では、既知の視点位置を持つ校正済みのステレオカメラで対象物体を撮影して得られたステレオ画像から、対象物体の３次元サーフェスを直接に生成するようにしているので、ステレオ画像上に隣接する画素の３次元情報は、同一の三角パッチ中にあれば同一平面上にあり、そうでなければ、それぞれが稜線（もしくは頂点）を共有する２つの平面上にあるものとして関連付けられるため、明白となる。
【００１５】
本発明によれば、複数のステレオ画像から対象物体の３次元サーフェスを短時間内に効率的に生成できるという優れた効果を奏する。
【００１６】
更に、本発明を適用してステレオ画像から対象物体の３次元サーフェスを直接に生成する際に、所定のパラメータの値を固定することにより高速化処理を行うことで、３次元サーフェス生成における更なる高速化を実現することができる。
【００１７】
従って、例えばロボットや自動車などのナビゲーションにおいても、地上のサーフェスを速やかに取得する必要があるロケットやヘリコプターの着陸時に、前方の坂の勾配や段差などの形状を速やかに得る必要がある移動するロボットに、問題なく本発明を適用することができる。
【発明を実施するための最良の形態】
【００１８】
以下、本発明を実施するための最良の形態を、図面を参照して詳細に説明する。
【００１９】
本発明は、ステレオカメラで撮影された画像（ステレオ画像）を利用して、対象物体の表面形状を高速に推定できる、３次元サーフェス生成方法に関する。
【００２０】
つまり、本発明は、既知の視点位置を持つ校正済みのカメラを用いて、それぞれ異なる視点位置から対象物体を撮影した複数のステレオ画像から、物体表面の３次元サーフェスを直接かつ高速に生成できる、３次元サーフェス生成方法である。
＜１＞本発明の着眼点及び概説
図１に示すように、２台のステレオカメラで撮影された対象物体（以下、単に、「物体」とも言う。）の３次元表面形状は、複数の三角パッチで構成されるメッシュによって表現できる。
【００２１】
３次元空間中の各三角パッチは、２つの画像に投影されたとき、画像上に２次元の三角形を描き、３次元メッシュは画像上で２次元メッシュを描く。２枚のうちの一方の画像に投影されたメッシュの全頂点は、その画像を撮影したカメラの光学中心からの距離、即ち、「奥行き」を持つ。
【００２２】
したがって、全頂点の奥行きを推定することにより、物体の３次元サーフェスを生成することができる。
【００２３】
ここで、全頂点の奥行きを定める画像を「基準画像」と呼び、もう一方の画像を「参照画像」と呼ぶ。基準画像上の１つの三角パッチの領域内の座標は、その三角パッチの３頂点の奥行きと、既知のカメラ視点位置等の情報を含むカメラパラメータを用いて、参照画像の座標に変換することができる。この座標変換を「平面射影変換（homography）」と呼ぶ。
【００２４】
すなわち、その三角パッチの３頂点の奥行きが既知であれば、基準画像上のその三角パッチ内の画素値は、その画素の座標を平面射影変換した座標における参照画像上の画素値とほぼ一致すると想定できる。
【００２５】
よって、基準画像上のメッシュの全頂点の奥行きが既知であれば、全頂点の奥行きと、カメラパラメータとを用いて、基準画像上の全ての三角パッチの平面射影変換を計算することができる。
【００２６】
このとき、各三角パッチの平面射影変換によって基準画像上の全ての座標を変換した座標における参照画像の画素値を抽出し、基準画像の座標と同じ位置に挿入することによって参照画像を変形した画像を生成すると、その画像は基準画像とほぼ一致する。
【００２７】
そこで、全頂点の奥行きと、カメラパラメータとが定める各三角パッチの平面射影変換を用いて、参照画像を変形した画像の画素値と、基準画像の画素値との差が小さくなるように、全頂点の奥行きを計算する。このようにして、全頂点の奥行きを求めることにより、物体の３次元サーフェスが生成できる。
【００２８】
要するに、本発明では、基準画像を三角パッチによるメッシュ（ポリゴンメッシュ）に分割し、各三角パッチは３次元空間中で平面をなし、物体表面の形状をよく近似できていると仮定する。このとき、それぞれのパッチ領域では、ステレオ画像間において平面射影変換（homography）による座標変換が成立する。
【００２９】
校正済みのステレオカメラを利用することにより、その座標変換は、自由度３の平面パラメータを用いて表すことができ、さらに、その自由度は三角パッチの３頂点の奥行きに対応する。つまり、全パッチの頂点の奥行きを求めることにより、物体表面の３次元形状を得ることができる。
【００３０】
そこで、本発明では、頂点の奥行きを求める際に、基準画像上のパッチと、対応する参照画像上のパッチとのＳＳＤ（Sum of Squared Differences）を考え、全パッチについて加算したＳＳＳＤ（Sum ofＳＳＤ）を最小にする奥行きを、Gauss-Newton法を用いて最適化する。
【００３１】
更に、本発明では、パッチに関するステレオ画像間のＳＳＤを、ＩＣＩＡ（Inverse Compositional Image Alignment）を利用して表現することにより、繰り返しごとの計算コストを低く抑える。本発明では、パッチ頂点の画像上の位置を、基準画像上に固定するため、１つの頂点が持つ自由度は、１（奥行き）となる。
＜２＞平面パラメータの推定
＜１＞では、本発明の着眼点及び概説を簡潔に述べたが、ここからは、図面及び数式を参照しながら、本発明に係る３次元サーフェス生成方法を詳細に説明する。
【００３２】
ここで、まず、本発明において重要な意味をもつ「平面姿勢パラメータ（以下、単に、平面パラメータと言う）」の推定方法について述べる。
【００３３】
即ち、ここでは、３次元空間中にある１つの平面に注目し、その平面を一意に表す３自由度のパラメータ（平面パラメータ）を、ＩＣＩＡ(Inverse Compositional Image Alignment)（非特許文献８参照）の考え方を利用して高速に推定する方法について、詳細に述べる。
【００３４】
ここで述べる「平面パラメータ推定方法」は、後述する「＜３＞曲面の形状推定」にて述べる「頂点の奥行き推定」において重要な役割を果たす。
＜２−１＞ＩＣＩＡ（Inverse Compositional Image Alignment）
まず、ＩＣＩＡのアルゴリズムについて概説する。
【００３５】
基準画像座標を
、参照画像座標を
とし、２つの座標間の変換（warp）を
と表す。
【００３６】
ただし、
は座標変換を表す。また、
は「座標変換パラメータベクトル」（以下、単に、「変換パラメータ」とも呼ぶ。）を表し、座標間の変換が平行移動の場合に２次元のベクトルとなり、また、座標間の変換がアフィン変換および平面射影変換の場合に、それぞれ６次元および８次元のベクトルとなる。
【００３７】
基準画像をＴ、参照画像をＩとすると、変換パラメータ
は、下記数１で示す画像間のＳＳＤ(Sum of Squared Differences)を、最小にする値として得られる。
【００３８】
【数１】
ただし、
は、基準画像における注目領域内の画素について加算することを表す。
【００３９】
数１を最小化することにより、変換パラメータ
を得る方法は、Gauss-Newton法による最適化を利用したLucas-Kanade法（非特許文献９を参照）が広く知られている。
【００４０】
ＩＣＩＡ（非特許文献８を参照）とは、Lucas-Kanade法の改良版であり、Lucas-Kanade法の精度と安定性を損なうことなく、高速化を実現したアルゴリズムである。
【００４１】
いま、座標変換
は、結合について閉じており、単位変換と逆変換を持つ（すなわち、変換群をなす）ものとする。
【００４２】
すると、座標変換
は、下記数２のように結合変換で表すことができる。
【００４３】
【数２】
ただし、
は、
の現在値（初期値）であり、
は、微小要素を持つベクトルである。また、
は、単位写像を表すものとする。
【００４４】
ＩＣＩＡでは、上記数２を利用し、下記数３を最小化する
を求める。
【００４５】
【数３】
そして、
の規則により、座標変換パラメータを更新する。この操作を、
が収束するまで繰り返し、最終的な値を
の推定値とする。
【００４６】
ＩＣＩＡの高速性は、次のように説明することができる。
【００４７】
Lucas-Kanade法（非特許文献９を参照）では、数１に対し、
付近でTaylor展開する。この場合、
が繰り返し計算のたびに更新され、それに伴い
が変化することから、
における
の値を繰り返しのたびに計算する必要がある。
【００４８】
これに対し、ＩＣＩＡでは、数３において、
付近でTaylor展開する。この場合、微分対象は基準画像Ｔであるため変化せず、また、
の値は常に固定値
において評価されることから、
を再計算する必要がない。これにより、ヘッセ行列を固定化でき、Lucas-Kanade法よりも、はるかに高速な演算を実現できる。
【００４９】
このように、ＩＣＩＡは、Lucas-Kanade法と比較して、繰り返し計算ごとの計算コストを大幅に削減できるため、近年では、画像レジストレーションのリアルタイム処理を実現するための強力な手段となっている。
＜２−２＞ＩＣＩＡを利用した平面射影変換パラメータ推定
ここで、ＩＣＩＡを利用して、８自由度の平面射影変換パラメータを推定する方法について述べる。
【００５０】
３次元空間中の平面を撮影した２枚の画像間では、下記数４のように、３×３の平面射影変換行列
による座標変換が成立する。
【００５１】
【数４】
ただし、
および
は、それぞれ
および
の同次座標を表す。また、’〜’は同値関係を表している。
【００５２】
数４は、定数倍の不定性を持つため、３×３の平面射影変換行列
の自由度は８であり、ここでは、一番右下の要素（即ち、第３行第３列要素）を固定したものとする。
【００５３】
また、数４は、
の座標変換に対応している。
【００５４】
ただし、ベクトル
は、
の８個の未知数をラスタ順に並べたものとなる。
【００５５】
ＩＣＩＡを用いた平面射影変換行列
の推定（非特許文献８を参照）では、現在値（初期値）
と微小要素を持つ行列
を用いて、
を下記数５のように表す。
【００５６】
【数５】
すなわち、
および
の２つの数式は、それぞれ、
および
による座標変換に対応する。
【００５７】
数３を最小化する
は、
をTaylor展開し、数３を
で微分しゼロとおくことにより、下記数６のように書ける。
【００５８】
【数６】
ただし、下記数７が成立する。
【００５９】
【数７】
ここで、
は、基準画像Ｔの勾配を示す１×２の列ベクトルであり、下記数８と表される。
【００６０】
【数８】
また、
は、下記数９で示す２×８の行列である。
【００６１】
【数９】
ここで、ｅは、
と表され、座標
における基準画像Ｔの画素値と、
に対応する参照画像Ｉ上の画素値との画素値の差を意味する。ｅは、全ての座標
を用いて、参照画像Ｉを変形（warp）した画像
を生成しておき、
として求める方法が効率がよいため、一般的に利用される。
および
は、注目領域内の各画素ごとに異なる値を持つ。しかし、既に述べたように、その値は一度だけ計算すればよく、繰り返しのたびに再計算する必要はない。これに伴い、ヘッセ行列の逆行列
も一度だけ計算すればよい。
＜２−３＞ＩＣＩＡを利用した平面パラメータ推定
ここでは、前述したＩＣＩＡの考え方を利用して、３自由度の平面パラメータを高速に推定する方法について説明する。
【００６２】
本発明の実施形態を説明する際にしては、煩雑さを避けるため、以後、画像座標
は内部パラメータが正規化された座標（正規化画像座標）を示すものとする。
【００６３】
図２に示すように、３次元空間中の平面
は、平面パラメータ
を用いて、
で表される。
【００６４】
ただし、
は平面
上の点の３次元座標を表し、基準カメラ座標系にて定義されるものとする。ここでは、
を求めたいパラメータベクトルとする。
【００６５】
平面射影変換行列
は、
を用いて、
と表すことができる（非特許文献１０を参照）。ただし、
は、カメラ間の回転行列および平行移動ベクトルを表し、それぞれ既知とする。
【００６６】
このとき、画像の座標変換は、
をパラメータとして、
と書ける。
【００６７】
ＩＣＩＡを利用して
を求めるには、現在値
と微小変化量
を用いて、座標変換を結合写像
で表現し、下記数１０を最小にする
を求めればよい。
【００６８】
【数１０】
結合写像を表現するために、
を用いて、平面射影変換行列
を下記数１１のようにおく。
【００６９】
【数１１】
数１１を数５の形式に変形すると、
および
は、それぞれ下記数１２、数１３で表される。
【００７０】
【数１２】
【００７１】
【数１３】
このとき、数１２および数１３は、それぞれ
および
による座標変換に対応する。
【００７２】
数１３は、
の関数として記述したことを意味するため、数１０を最小にする
は、微分の連鎖法則を利用して、下記数１４のように得られる。
【００７３】
【数１４】
ただし、下記数１５、数１６及び数１７が成立する。
【００７４】
【数１５】
【００７５】
【数１６】
【００７６】
【数１７】
ここで、
である。また、
は、数６の場合と同様に、基準画像Ｔの勾配を示す１×２の列ベクトルである。そして、
は、下記数１８で示す２×９の行列である。
【００７７】
【数１８】
数９と数１８における
のサイズが異なるのは、
で表現された
には不定性が存在しないためである。
【００７８】
また、
にて評価した結果であり、行列
は、下記数１９で表される。
【００７９】
【数１９】
ただし、下記数２０、数２１、及び数２２が成立する。
【００８０】
【数２０】
【００８１】
【数２１】
【００８２】
【数２２】
ここで、{ｒ１，ｒ２，…，ｒ９}は、
の要素をラスタ順に並べたものであり、また、
である。
【００８３】
κは画像座標
に依存しないため、数１４を改めて書き直すと、下記数２３のようになる。
【００８４】
【数２３】
ただし、下記数２４、数２５及び数２６が成立する。
【００８５】
【数２４】
【００８６】
【数２５】
【００８７】
【数２６】
＜２−２＞で述べた平面射影変換パラメータの推定と、ここにて述べた平面パラメータ推定との大きな相違点は、前者では、
が変換群をなすのに対し、後者では、
が変換群をなさないことにある。
【００８８】
そのため、数１３が
に依存してしまい、
の値は繰り返し計算ごと変化する。しかし、その変化は、数２３におけるスカラーκに集約されており、
および
は変化しないことを考慮すると、変換群でないことによって生じる繰り返しごとの計算コストの増加は、無視できるほど小さい。
【００８９】
これにより、ＩＣＩＡを利用した平面パラメータ推定は、高速な演算にて実現できることが分かる。
＜３＞曲面の形状推定
＜２−３＞では、３次元空間中の１つの平面に注目し、その平面パラメータを高速に推定する方法について述べた。ここでは、＜２−３＞にて述べた方法を利用して、三角パッチで構成されたポリゴンメッシュ（三角メッシュ）の全頂点の奥行きを求める方法について述べる。
＜３−１＞三角メッシュによる表面形状モデル
まず、三角メッシュを用いた物体表面形状モデルについて述べる。
【００９０】
本発明では、基準画像を小さな三角形で分割してメッシュを生成する必要がある。
【００９１】
メッシュを生成する方法としては、例えば、画像全体を方形で分割した後に方形の対角線を結び三角形に再分割する方法や、画像上から特徴点を抽出した後に特徴点を母点としてドロネー分割する方法など、様々な方法が考えられる。
【００９２】
本発明では、画像領域全体（基準画像領域全体）が微小な三角形で分割されていれば、上記のメッシュ生成方法に限らず、どのような方法を用いてもよい。ここで、正三角形によって構成されたメッシュの例を図３に示す。
【００９３】
メッシュの頂点を
と表し、
は対象物体表面の３次元形状に対応した奥行きｄｍを持つものとする。また、頂点によって構成される三角形をパッチと呼び、Ｃｎ（ｎ＝１，２，…，Ｎ）と表す。各パッチは３次元空間中において平面を構成するものとする。
【００９４】
すなわち、パッチによって構成されるメッシュはポリゴンメッシュであり、メッシュが十分細かければ、任意の３次元形状を精度よく近似できると考えられる。
【００９５】
本発明では、メッシュの細かさに関する検討は特に行わない。ただし、画像レジストレーション分野では広く知られているように、粗いメッシュ（即ち、大きなメッシュ）を利用した場合には、より推定が安定になると考えられるため、＜３−３＞にて述べる全頂点の奥行き推定では、最初は粗いメッシュを用いて奥行きを推定し、徐々にメッシュを細かくする粗密探索のアプローチが有効である。
【００９６】
また、画像のＲＭＳＥの値を利用して細かい起伏を持つ領域を推定し、その領域に対応するパッチのみを細かくするなど、再分割（subdivision）の考え方を利用したパッチの分割方法も有効であると考えられる。
＜３−２＞三角パッチの３頂点に対する奥行き推定
次に、＜２−３＞にて述べた平面パラメータ推定方法を利用して、１つの三角パッチの３頂点の奥行きを求める方法について述べる。
【００９７】
ある三角パッチの３頂点の画像上の座標を
とし、それら頂点の奥行きをｄｍとする。このとき、当該三角パッチの３頂点の３次元空間中の座標
は、図４に示すように、
と表すことができる。
【００９８】
３次元空間中の三角パッチの平面パラメータを
とすると、ある三角パッチについて、平面パラメータ
と３頂点の奥行きｄｍ（ｍ＝ｉ,ｊ,ｋ）の関係は、下記数２７で表される。
【００９９】
【数２７】
ただし、下記数２８及び数２９が成立する。
【０１００】
【数２８】
【０１０１】
【数２９】
これより、平面パラメータ
と、奥行きパラメータγは比例関係にあり、その比例関係を決定する行列
は、画像上の三角パッチの位置によって一意に定まる。
【０１０２】
なお、「奥行き」と「奥行きパラメータ」の関係は、上記数２９に示す通りである。つまり、γは、ある三角パッチの３頂点の奥行きパラメータを表し、ｄｍ（ｍ＝ｉ,ｊ,ｋ）はそれら頂点の奥行きを表す。
【０１０３】
よって、数２３に微分の連鎖法則を適用することにより、平面パラメータ
の微小変化量
の代わりに、奥行きパラメータγの微小変化量を求めることができる。
【０１０４】
すなわち、γ０およびγΔを、それぞれ奥行きパラメータγの現在値および微小変化量とし、γΔ＝γ０＋γΔとすると、下記数３０で表される式を最小化するγΔは、数２３により、下記数３１で得られる。
【０１０５】
【数３０】
【０１０６】
【数３１】
ただし、下記数３２及び数３３が成立する。
【０１０７】
【数３２】
【０１０８】
【数３３】
ここで、
である。また、κは、数２７により平面パラメータの現在値
を求め、求めた
を数２４に代入することにより、得られる。そして、
および
は、数２３と同じものである。
【０１０９】
数２３と同様に、数３１においても、
は繰り返し計算のたびに変化しない。
【０１１０】
よって、平面パラメータ推定の場合とほぼ同等の計算コストで、奥行きパラメータγを推定することが可能となる。
＜３−３＞三角メッシュの全頂点の奥行き推定
ここで、三角メッシュの全頂点の奥行きを最適化する方法について説明する。
【０１１１】
＜３−２＞では、１つの三角パッチに注目したときの基準画像と参照画像間のＳＳＤを数３０によって表し、そのＳＳＤを最小にする３頂点の奥行きパラメータγを求めた。
【０１１２】
ここでは、３自由度の平面パラメータと３頂点の奥行きの関係を利用して、＜３−１＞にて述べた三角メッシュにおける全パッチについてＳＳＤを加算することによって得られるＳＳＳＤ(Sum of ＳＳＤ)を最小にすることにより、メッシュの全頂点の奥行きを求める。
【０１１３】
＜３−１＞にて述べたように、三角メッシュにおけるパッチ数をＮ、頂点の数をＭとする。また、全頂点の奥行きｄｍ，（ｍ＝１，２，…，Ｍ）の逆数を並べ、
とする。
【０１１４】
ここでは、
を「三角メッシュにおける全頂点の奥行きパラメータ」とし、以下、単に、「全頂点奥行きパラメータ」とも称する。
【０１１５】
また、
および
を、それぞれ全頂点奥行きパラメータ
の現在値および微小変化量とし、
とおく。
【０１１６】
そして、下記数３４に示すように、数３０を三角メッシュにおける全パッチについて加算したＳＳＳＤ(Sum ofＳＳＤ)を最小化する問題を考える。
【０１１７】
【数３４】
をTaylor展開し、数３４を
で微分してゼロとおくことにより、数３４を最小にする
は、下記数３５となる。
【０１１８】
【数３５】
ただし、下記数３６及び数３７が成立する。
【０１１９】
【数３６】
【０１２０】
【数３７】
ここで、κｎは数２４に準じたものであるが、平面パラメータの現在値
がパッチごと異なり、それに伴いκの値がパッチごと異なるため、添え字ｎを付与して表記している。
【０１２１】
また、
は数２８に準じているが、その要素はパッチの画像上の座標
に依存しており、パッチごと異なるため、κと同様に添え字ｎを付与した。
【０１２２】
そして、表１に、数３５の計算に必要な行列やベクトルを示している。
【０１２３】
【表１】
数３５により、全頂点奥行きパラメータの微小変化量
を求める場合は、＜２−３＞で述べた単一平面のパラメータ推定や、＜３−２＞で述べた三角パッチの３頂点の奥行き推定とは異なり、数３６の
に関する計算を繰り返しのたびに行う必要がある。
【０１２４】
ただし、
の計算に注目すると、行列
は、各パッチごと異なるものの、繰り返しごとに変化せず、さらに、この行列は実質的に３×３＝９個の要素しか持っていないため、全パッチに関して加算しても計算量は比較的小さい。
【０１２５】
また、
の計算は、実質的に画素数に依存しており、単一平面のパラメータ推定や３頂点の奥行き推定を、同じ画素数に関して計算した場合とほぼ同じ計算量となる。
【０１２６】
すなわち、数３５の計算量全体を評価すると、同じ画素数を用いて単一平面のパラメータ推定をした場合と比較して、
の逆行列計算の分だけ計算コストの増加と言える。
【０１２７】
ただし、
の行列サイズは、Ｍ×Ｍであり、三角メッシュにおいて頂点の数が多くなるほど、その計算量は大きくなる。本発明では、計算速度に関して、評価の高いCLapackライブラリを利用して逆行列計算を実装し、計算時間の短縮を実現している。
【０１２８】
ここで、上述したように、三角メッシュにおける全頂点の奥行き推定を行うことにより、対象物体の表面形状を推定する、「曲面形状推定」、すなわち、「微小平面モデルを用いた曲面形状推定」のアルゴリズムの一例を図５に示す。
【０１２９】
図５に例示されるアルゴリズムは、初期化処理と推定処理に分けられる。初期化処理では、推定処理にて利用する固定値（繰り返しごとに変化しない行列やベクトル）を計算している。また、推定処理では、数３５の計算に必要な値を各パッチについて求め、
の計算と
の更新を収束するまで繰り返す。
【０１３０】
なお、数３４は、基準画像１枚と参照画像１枚を利用した場合、即ち２枚の画像を利用した場合のＳＳＳＤを示している。参照画像を２枚以上利用する場合は、即ち基準画像を含めて３枚以上の画像を利用する場合は、数３４は、基準画像と参照画像のペアの数だけ存在する。そこで、数３４を基準画像と参照画像のペアの数だけ加算することによって得られる新たなＳＳＳＤを利用して、全頂点の奥行きパラメータを推定する。
＜３−４＞高速化処理
＜３−３＞で述べた「三角メッシュの全頂点の奥行き推定方法」において、三角メッシュにおける頂点数が非常に大きい場合は、数３５に必要な
の逆行列を得るために、繰り返しごとに必要な計算時間が長くなり、つまり、計算コストが高くなることがある。
【０１３１】
この逆行列の計算が必要となるのは、κｎが繰り返し計算ごとに異なることが原因である。そこで、計算時間をより短くする方法として、κｎを固定化する方法が考えられる。
【０１３２】
数２４から明らかなように、κｎの変化は、ｎ番目のパッチの平面パラメータ
が繰り返し計算ごとに変化することによって生じる。しかし、
の繰り返しごとの変化は、非常に小さいため、κｎを固定しても、計算結果に大きな違いは生じない。
【０１３３】
よって、特定の繰り返し回数の間は、全パッチにおけるκｎを固定化し、
の再計算を行わずに計算時間を短くすることが可能である。例えば、１〜５回、６〜１０回の繰り返しの間では、
を固定し、１回目と６回目だけ
とその逆行列を再計算する方法などが考えられる。
【０１３４】
さらに、もし、ｎ番目のパッチにおける平面パラメータ
の現在値が、
とほぼ直交している場合、もしくは、
で表される平面と基準カメラ座標原点との距離が大きければ、数２４により
となることが分かる。
【０１３５】
そこで、はじめからκｎ＝−１,（ｎ＝１,２,…,Ｎ）として、対象物体の表面モデルを生成することも可能である。この場合、１度だけ
を求めればよく、計算時間をさらに短縮することができる。
＜３−５＞本発明による３次元サーフェスを生成するための処理の流れ
図６に、本発明を適用して、対象物体を撮影した複数のステレオ画像から、対象物体の３次元サーフェスを生成するための基本的な処理の流れを示す。
【０１３６】
なお、複数のステレオ画像のうち、一枚のみを基準画像とし、残りの画像を参照画像とする。また、ステレオ画像を撮影した多眼ステレオカメラについて、カメラのレンズ歪や焦点距離などの内部パラメータと、ステレオカメラ間の配置を示す外部パラメータとをあわせて、「カメラパラメータ」としている。
【０１３７】
説明を簡単にするために、以下、２眼ステレオカメラを用いた場合について説明する。つまり、撮影したステレオ画像が２枚で、１枚は基準画像とし、もう１枚は参照画像とする。
【０１３８】
図６に示されるように、この処理における入力は、複数のステレオ画像と、ステレオカメラのカメラパラメータである。
【０１３９】
即ち、既知の視点位置にある校正済みの多眼ステレオカメラを用いて、対象物体を撮影して、ステレオ画像を取得すれば、本発明を適用して対象物体の３次元サーフェスを直接かつ高速に生成することができる。
【０１４０】
処理の手順として、図６に示されるように、まず、所定の方法（＜３−１＞の記述を参照）により、基準画像上に２次元三角メッシュ（三角メッシュ）を生成する。
【０１４１】
次に、生成された三角メッシュの全頂点に対して、奥行きパラメータの初期値（全頂点奥行きパラメータの初期値）を定める。
【０１４２】
そして、定めた全頂点奥行きパラメータの初期値と、カメラパラメータと、基準画像及び参照画像で構成される複数のステレオ画像（２眼ステレオカメアを用いた場合に、一枚の基準画像と一枚の参照画像で構成する２枚のステレオ画像になる。）とを入力とし、基準画像上のパッチと対応する参照画像上のパッチとのＳＳＤを考え、三角メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、全頂点奥行きパラメータを求める最適化処理（全頂点奥行きパラメータ最適化処理）を行う。
【０１４３】
即ち、全パッチについて加算したＳＳＳＤを最小にする、最適な全頂点奥行きパラメータは、全頂点奥行きパラメータの初期値とカメラパラメータとを入力とする最適化手法を用いて、求めることができる。最適化手法では、全頂点奥行きパラメータの微小変化量の計算と、全頂点奥行きパラメータの現在値の更新とを繰り返し行う必要がある。
【０１４４】
最適化計算における微小変化量の計算は、コスト関数を微小変化量で表現した後に、コスト関数を微小変化量で微分してその結果をゼロとおくことにより、コスト関数の極値（極小値）となるパラメータを求める方程式をたてる。このとき、その方程式を解くために必要な情報は、繰り返しごとに求める必要があるため、一般に短時間で計算を行うことは困難である。
【０１４５】
しかし、＜３−４＞高速化処理でも述べたように、全頂点奥行きパラメータの微小変化量の計算において、必要な情報を固定化することができれば、これらの情報は繰り返し毎に計算する必要はなくなり、一度だけ計算すればよいことになり、極めて短時間で計算することが可能となる。
【０１４６】
全頂点の真の奥行きパラメータが定める平面射影変換は、全頂点奥行きパラメータの現在値が定める平面射影変換と、全頂点奥行きパラメータの微小変化量が定める平面射影変換の合成によって表すことができる。そして、全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって基準画像を変形して得られた画像の画素値と、全頂点奥行きパラメータの現在値が定める平面射影変換によって参照画像を変形して得られた画像の画素値との差を用いて、コスト関数を表現すると、繰り返しごとに必要な多くの情報は固定化でき、一度だけ計算すればよいことになる。これにより、より短時間での計算を実現できる。
【０１４７】
図７は、図６に示した最適化処理（全頂点奥行きパラメータ最適化処理）の基本フローである。
【０１４８】
図７に示されるように、最適化処理（全頂点奥行きパラメータ最適化処理）では、まず、基準画像及び基準画像を撮影した基準カメラのカメラパラメータを用いて、全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し計算に必要な情報を一度だけ計算する。
【０１４９】
次に、全頂点奥行きパラメータの現在値と、カメラ視点位置を含むカメラパラメータとによって定まる各三角形（各三角パッチ）の平面射影変換によって参照画像を変形して得られた画像の画素値と、基準画像の画素値との差分を計算する。以下、この処理を「差分値計算処理」とも称する。
【０１５０】
差分値が得られた後に、全頂点奥行きパラメータの微小変化量を計算する。つまり、図７に示されるように、繰り返し計算に必要な情報と、差分値と、カメラパラメータと、全頂点奥行きパラメータの現在値を入力として、全頂点奥行きパラメータの微小変化量を計算する。以下、この処理を「全頂点奥行きパラメータの微小変化量計算処理」とも称する。
【０１５１】
そして、得られた全頂点奥行きパラメータの微小変化量と、全頂点奥行きパラメータの現在値とに基づいて、全頂点奥行きパラメータの現在値を更新する。以下、この処理を「全頂点奥行きパラメータの現在値更新処理」とも称する。
【０１５２】
「最適化処理」では、「差分値計算処理」、「全頂点奥行きパラメータの微小変化量計算処理」及び「全頂点奥行きパラメータの現在値更新処理」を、所定の収束条件を満たすまで繰り返す。繰り返し処理（繰り返し計算）により、所定の収束条件を満たしたときに最終的に得られた全頂点奥行きパラメータの現在値を「全頂点奥行きパラメータ」とし、よって、三角メッシュにおける全頂点の奥行きが定められることになり、対象物体の３次元サーフェスを生成できる訳である。
＜４＞実験結果
以下では、本発明の有効性を確認するために行った実験結果を示す。まず、合成画像を利用して本発明を適用した様々なシミュレーション実験の結果について示し、次に、実画像を用いて本発明を適用した実画像実験の結果を示す。
＜４−１＞シミュレーション実験
シミュレーション実験（合成画像を用いた実験）では、５１２［pix］×５１２［pix］の画像を利用し、また、ステレオカメラの設定パラメータは、以下に示す値を利用した。
【０１５３】
カメラの内部パラメータは、基準カメラ及び参照カメラとも同じであり、下記数３８の
で示す。
【０１５４】
【数３８】
また、カメラの外部パラメータは、
（ここで、
は単位行列である）、
とした。ここで、
は、カメラ間の回転行列と平行移動ベクトルをそれぞれ示している。
【０１５５】
図８には、シミュレーション実験用に作成した合成画像を示す。本発明では、推定時に使用するメッシュを生成するのに、即ち、基準画像上に三角メッシュを生成するのに、前述したように、例えば、固定のメッシュを利用する方法や、抽出された画像特徴点を母点としたドロネー分割などの方法を利用することができるが、シミュレーション実験では、図８に示すように、５０［pix］の辺を持つ正三角形を用いてメッシュを作成し、テクスチャを張った平面を１０ｍ先に配置した。
【０１５６】
よって、図８の合成画像を用いた場合、頂点の奥行きは全て１０ｍであり、全パッチの平面パラメータは、
である。
【０１５７】
図８（Ａ）及び図８（Ｂ）には、作成された参照画像と基準画像をそれぞれ示している。図８（Ａ）及び図８（Ｂ）の画像は、５１２［pix］×５１２［pix］の画像における中央部分の４２０［pix］×４２０［pix］領域を表示したものである。このようにして作成した参照画像及び基準画像に対し、図８（Ｃ）に示したメッシュの頂点に関して、奥行き推定を行った。
【０１５８】
図８（Ｃ）におけるメッシュの頂点数は６１で、パッチ数は９６である。奥行き推定に利用した初期値は、全頂点において、１５ｍとした。
【０１５９】
本実験では、収束条件を特に設定せず、３０回まで繰り返し計算を行った。プログラムはＣ言語で実装し、逆行列計算にはCLapackライブラリを利用した。
【０１６０】
図８（Ｃ）で示したメッシュを利用した場合、Pentium（登録商標）ＩＶ(２.８ＧＨｚ）でLinux（登録商標）ＯＳを備えたパソコン上において、３０回の繰り返し計算に要する時間は、約０.９秒であった。
【０１６１】
メッシュが収束してゆく様子を図９に示す。図９では、基準画像上のメッシュを赤線で示しており、また、真値から計算される参照画像上の頂点を、推定された奥行きによって得られる平面射影変換を用いて基準画像上に投影し、その点によって構成されるメッシュを緑線で示している。
【０１６２】
よって、正しい奥行きが推定されれば、赤線のメッシュと緑線のメッシュは、完全に一致する。図９から、奥行き推定を繰り返すことにより、赤線のメッシュと緑線のメッシュが一致してゆくことがよく分かる。
【０１６３】
また、図１０には、画像のＲＭＳＥ（Root Mean Squared Error）と奥行きについて、繰り返し計算ごとの値の変化を示している。即ち、図１０の横軸は、最適化計算（最適化処理）における繰り返し回数を表している。
【０１６４】
図１０（Ａ）に、いくつかのパッチに関する画像ＲＭＳＥの変化と、全パッチについての画像ＲＭＳＥの変化（＊付きの実線）を示す。図１０（Ａ）から、３０回の繰り返し計算までに、適切に収束していることが分かる。
【０１６５】
そして、図１０（Ｂ）は、いくつかの頂点の奥行きの変化を示している。図１０（Ｂ）から、最終的に、真値（１０ｍ）に収束していることが分かる。
【０１６６】
次に、対象を球体として生成した合成画像を図１１に示す。図１１では、光軸上の１５ｍ先に中心を持ち、半径約７ｍの球体にテクスチャを貼り付けた。図１１（Ａ）及び図１１（Ｂ）は、作成された参照画像と基準画像を示している。図１１の画像は、図８の場合と同様に、５１２［pix］×５１２［pix］の画像における中央部分の４２０［pix］×４２０［pix］領域を表示したものである。このようにして作成した参照画像及び基準画像に対し、図１１（Ｃ）に示したメッシュの頂点に関して、奥行き推定を行った。
【０１６７】
図８の場合と同様に、図１１（Ｃ）におけるメッシュの頂点数は６１で、パッチ数は９６である。奥行き推定に利用した初期値は、全頂点において、１０ｍとした。
【０１６８】
図１２に、図１１の物体の実形状と、推定された形状を示す。図１２により、奥行き推定が適切に行われていることがよく分かる。なお、図１２（Ａ）、図１２（Ｂ）、図１２（Ｃ）及び図１２（Ｄ）は、全て同一視点から見たものである。
【０１６９】
また、繰り返し計算における画像ＲＭＳＥの変化、および、各頂点の真の奥行きと推定された奥行きとの絶対誤差（以下、単に、奥行き推定誤差の絶対値と呼ぶ。）を、図１３に示す。図１３（Ａ）では、いくつかのパッチにおける画像ＲＭＳＥと、全パッチに関する画像ＲＭＳＥを示し、図１３（Ｂ）では、いくつかの頂点における奥行き推定誤差の絶対値と、全頂点に関する奥行きのＲＭＳＥを示している。図１３（Ａ）及び図１３（Ｂ）から、ほぼ２０回の繰り返し計算にて、物体の形状が正しく推定できていることが分かる。
【０１７０】
更に、球体を用いた実験において、κｎ＝−１,（ｎ＝１,２,…,Ｎ）と固定した場合における繰り返しごとの画像ＲＭＳＥの変化と奥行きの絶対誤差の変化を図１４に示す。
【０１７１】
図１４（Ａ）及び図１４（Ｂ）は、図１３（Ａ）及び図１３（Ｂ）と、同様である。繰り返しごとの画像ＲＭＳＥ、および奥行き推定誤差の絶対値は、図１３（Ａ）及び図１３（Ｂ）のものとはやや異なるものの、ほぼ２０回の繰り返し計算にて、物体の形状が正しく推定できていることが分かる。
【０１７２】
すなわち、κｎ＝−１と固定して計算を行っても、κｎの値を固定しない場合との違いは、僅かである。よって、κｎの値を固定することにより、計算時間を一層短縮できることが分かる。
＜４−２＞実画像実験
実画像実験（実画像を用いた実験）に利用された実画像を図１５に示す。図１５の画像は、月面を模擬した１２８０［pix］×９６０［pix］のステレオ画像である。
【０１７３】
本実験では、粗いメッシュを徐々に細かくしてゆく粗密探索的な方法を用いて、表面モデルを生成した。
【０１７４】
ここでは、一辺の長さが４６４［pix］の三角パッチから始め、次の推定では、パッチのサイズを半分にするとともに、前回の頂点の奥行きを初期値とした推定を行った。
【０１７５】
図１５（Ａ）の基準画像では、最後にパッチのサイズとして使用した２９［pix］のメッシュを示している。このとき、メッシュの頂点数は８１７で、パッチ数は１５３６であった。
【０１７６】
図１５（Ａ）の基準画像及び、図１５（Ｂ）の参照画像に基づいて、本発明を適用して生成された表面モデルを図１６に示す。
【０１７７】
図１６には、図１５の実画像を用いた実験において、別々の視点から見た、本発明を適用して推定された３次元形状を示している。図１６（Ａ）、図１６（Ｃ）及び図１６（Ｅ）は、メッシュとテクスチャを同時に表示したもので、図１６（Ｂ）、図１６（Ｄ）及び図１６（Ｆ）は、テクスチャのみを表示したものである。
【０１７８】
図１６から、メッシュの端に位置する頂点は、その外側にメッシュが存在しないため、他のメッシュによる十分な拘束を受けることができず、推定が不安定になる傾向が見られる。しかし、メッシュ内部の領域の推定は安定であり、メッシュ中央部に見られるエッジを持った勾配や、左上部にみられる小さなクレータの穴が良好に復元されている。
以上のように、合成画像を用いた実験、及び実画像を用いた実験を通じて、本発明の有効性、即ち、本発明を適用することにより、効率的な表面形状推定が可能であることを確認した。
【０１７９】
また、上記の実験では、２台のステレオカメラを利用して得られた２枚のステレオ画像に基づいて、本発明を適用し、対象物体の３次元サーフェスを生成したが、本発明はそれに限定されず、３台以上のステレオカメラを利用して得られた３枚以上のステレオ画像に基づいても、本発明を適用し、対象物体の３次元サーフェスを生成することができる。
【図面の簡単な説明】
【０１８０】
【図１】対象物体とその３次元形状を説明するための模式図である。
【図２】３次元空間中の平面を説明するための模式図である。
【図３】基準画像上に生成した三角メッシュの例を示す図である。
【図４】三角パッチの頂点の３次元座標を説明するための模式図である。
【図５】本発明において、曲面形状推定のアルゴリズムを示すプログラムフローである。
【図６】本発明を適用して３次元サーフェスを生成するための基本的な処理の流れを示す図である。
【図７】図６に示した最適化処理の基本フローである。
【図８】シミュレーション実験用に作成した合成画像を示す図である。図８（Ａ）に参照画像を、図８（Ｂ）に基準画像を、図８（Ｃ）に基準画像上に生成されたメッシュ（推定時使用メッシュ）をそれぞれ示している。
【図９】図８の合成画像を用いた実験におけるメッシュが収束してゆく様子を説明するための図である。
【図１０】図８の合成画像を用いた実験において、画像のＲＭＳＥと奥行きについて、繰り返し計算ごとの値の変化を示すグラフである。
【図１１】シミュレーション実験用に対象を球体として生成された合成画像を示す図である。図１１（Ａ）に参照画像を、図１１（Ｂ）に基準画像を、図１１（Ｃ）に基準画像上に生成されたメッシュ（推定時使用メッシュ）をそれぞれ示している。
【図１２】図１１の合成画像を用いた実験における実形状と、本発明を適用して推定された形状を示す図である。
【図１３】図１１の合成画像を用いた実験において、画像のＲＭＳＥと奥行き推定誤差の絶対値について、繰り返し計算ごとの値の変化を示すグラフである。
【図１４】図１１の合成画像を用いた実験において、κｎの値を「−１」と固定した場合に、画像のＲＭＳＥと奥行き推定誤差の絶対値について、繰り返し計算ごとの値の変化を示すグラフである。
【図１５】実画像実験に利用された実画像（ステレオ画像）を示す図である。図１５（Ａ）に基準画像を、図１５（Ｂ）に参照画像をそれぞれ示している。
【図１６】図１５の実画像を用いた実験において、本発明を適用して推定された３次元形状を示す図である。
【特許請求の範囲】
【請求項１】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する第１のステップと、
前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータと、カメラ視点位置を含むカメラパラメータとによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルと、前記基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、
前記複数の画像と、前記カメラパラメータと、前記全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、前記全頂点奥行きパラメータの微小変化量の計算と、前記全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、前記コスト関数の値が最小となる、前記全頂点奥行きパラメータを計算する第２のステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項２】
前記第２のステップをより短時間内で行うために、
前記第２のステップでは、
前記全頂点奥行きパラメータの現在値が定める平面射影変換と、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換との合成によって真の平面射影変換を表現し、
前記全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値を要素とするベクトルと、前記全頂点奥行きパラメータの現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルとによって、前記距離と同等の距離を表現することによって、繰り返しごとに必要な情報を１度だけ計算する請求項１に記載の３次元サーフェス生成方法。
【請求項３】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する、２次元メッシュ生成ステップと、
前記２次元メッシュ生成ステップで生成された２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差が小さくなるように、前記全頂点奥行きパラメータを最適化手法によって推定する、全頂点奥行き推定ステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項４】
前記全頂点奥行き推定ステップは、
前記全頂点奥行きパラメータの微小変化量が定める平面射影変換と、前記全頂点奥行きパラメータの現在値が定める平面射影変換とに基づいて、真の平面射影変換が合成され、
前記微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値と、前記現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値の差が小さくなるように、前記微小変化量を計算することによって、繰り返しごとに必要な情報を１度だけ計算する請求項３に記載の３次元サーフェス生成方法。
【請求項５】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角パッチによって構成された２次元メッシュを生成するステップと、
前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータの初期値を定めるステップと、
定めた全頂点奥行きパラメータの初期値と、カメラ視点位置を含むカメラパラメータと、前記複数の画像とに基づいて、前記２次元メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、前記全頂点奥行きパラメータを求める最適化処理を行うステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項６】
前記最適化処理では、
前記基準画像、及び前記カメラパラメータを用いて、前記全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し処理に必要な情報を一度だけ計算する、繰り返し処理用情報計算ステップと、
前記全頂点奥行きパラメータの現在値と、前記カメラパラメータとによって定まる各三角パッチの平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差分を計算する、差分値計算処理ステップと、
前記繰り返し処理用情報計算ステップで得られた必要な情報と、前記差分値計算処理ステップで得られた差分値と、前記カメラパラメータと、前記全頂点奥行きパラメータの現在値を入力として、前記全頂点奥行きパラメータの微小変化量を計算する、微小変化量計算処理ステップと、
前記微小変化量計算処理ステップで得られた微小変化量と、前記全頂点奥行きパラメータの現在値とに基づいて、前記全頂点奥行きパラメータの現在値を更新する、現在値更新処理ステップと、
を有する請求項５に記載の３次元サーフェス生成方法。
【請求項７】
前記最適化処理では、前記差分値計算処理ステップ、前記微小変化量計算処理ステップ及び前記現在値更新処理ステップを、所定の収束条件を満たすまで繰り返し、繰り返し処理により、前記収束条件を満たしたときに得られた全頂点奥行きパラメータの現在値を全頂点奥行きパラメータとすることにより、前記３次元サーフェスを生成する請求項６に記載の３次元サーフェス生成方法。
【請求項１】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する第１のステップと、
前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータと、カメラ視点位置を含むカメラパラメータとによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルと、前記基準画像の画素値を要素とするベクトルとの距離をコスト関数の１つの項とし、
前記複数の画像と、前記カメラパラメータと、前記全頂点奥行きパラメータの初期値とを入力とする最適化手法を用いて、前記全頂点奥行きパラメータの微小変化量の計算と、前記全頂点奥行きパラメータの現在値の更新を、所定の条件まで繰り返し行うことによって、前記コスト関数の値が最小となる、前記全頂点奥行きパラメータを計算する第２のステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項２】
前記第２のステップをより短時間内で行うために、
前記第２のステップでは、
前記全頂点奥行きパラメータの現在値が定める平面射影変換と、前記全頂点奥行きパラメータの微小変化量が定める平面射影変換との合成によって真の平面射影変換を表現し、
前記全頂点奥行きパラメータの微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値を要素とするベクトルと、前記全頂点奥行きパラメータの現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値を要素とするベクトルとによって、前記距離と同等の距離を表現することによって、繰り返しごとに必要な情報を１度だけ計算する請求項１に記載の３次元サーフェス生成方法。
【請求項３】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角形によって構成された２次元メッシュを生成する、２次元メッシュ生成ステップと、
前記２次元メッシュ生成ステップで生成された２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータによって定まる各三角形の平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差が小さくなるように、前記全頂点奥行きパラメータを最適化手法によって推定する、全頂点奥行き推定ステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項４】
前記全頂点奥行き推定ステップは、
前記全頂点奥行きパラメータの微小変化量が定める平面射影変換と、前記全頂点奥行きパラメータの現在値が定める平面射影変換とに基づいて、真の平面射影変換が合成され、
前記微小変化量が定める平面射影変換の逆変換によって前記基準画像を変形して得られた画像の画素値と、前記現在値が定める平面射影変換によって前記参照画像を変形して得られた画像の画素値の差が小さくなるように、前記微小変化量を計算することによって、繰り返しごとに必要な情報を１度だけ計算する請求項３に記載の３次元サーフェス生成方法。
【請求項５】
既知の視点位置にある校正済みのカメラを用い、対象物体を異なる視点位置から撮影して得られた複数の画像から、前記対象物体の３次元サーフェスを直接に生成する、３次元サーフェス生成方法であって、
前記複数の画像のうち、１枚の画像を基準画像とし、残りの画像を参照画像とした上で、前記基準画像上に複数の三角パッチによって構成された２次元メッシュを生成するステップと、
前記２次元メッシュの全頂点の奥行きを表す全頂点奥行きパラメータの初期値を定めるステップと、
定めた全頂点奥行きパラメータの初期値と、カメラ視点位置を含むカメラパラメータと、前記複数の画像とに基づいて、前記２次元メッシュにおける全パッチについて加算したＳＳＳＤを最小にする、前記全頂点奥行きパラメータを求める最適化処理を行うステップと、
を有することを特徴とする３次元サーフェス生成方法。
【請求項６】
前記最適化処理では、
前記基準画像、及び前記カメラパラメータを用いて、前記全頂点奥行きパラメータの微小変化量を求めるために利用される情報、即ち、繰り返し処理に必要な情報を一度だけ計算する、繰り返し処理用情報計算ステップと、
前記全頂点奥行きパラメータの現在値と、前記カメラパラメータとによって定まる各三角パッチの平面射影変換によって前記参照画像を変形して得られた画像の画素値と、前記基準画像の画素値との差分を計算する、差分値計算処理ステップと、
前記繰り返し処理用情報計算ステップで得られた必要な情報と、前記差分値計算処理ステップで得られた差分値と、前記カメラパラメータと、前記全頂点奥行きパラメータの現在値を入力として、前記全頂点奥行きパラメータの微小変化量を計算する、微小変化量計算処理ステップと、
前記微小変化量計算処理ステップで得られた微小変化量と、前記全頂点奥行きパラメータの現在値とに基づいて、前記全頂点奥行きパラメータの現在値を更新する、現在値更新処理ステップと、
を有する請求項５に記載の３次元サーフェス生成方法。
【請求項７】
前記最適化処理では、前記差分値計算処理ステップ、前記微小変化量計算処理ステップ及び前記現在値更新処理ステップを、所定の収束条件を満たすまで繰り返し、繰り返し処理により、前記収束条件を満たしたときに得られた全頂点奥行きパラメータの現在値を全頂点奥行きパラメータとすることにより、前記３次元サーフェスを生成する請求項６に記載の３次元サーフェス生成方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【公開番号】特開２００８−１２３０１９（Ｐ２００８−１２３０１９Ａ）
【公開日】平成２０年５月２９日（２００８．５．２９）
【国際特許分類】
【出願番号】特願２００６−３０２７１７（Ｐ２００６−３０２７１７）
【出願日】平成１８年１１月８日（２００６．１１．８）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
【公開日】平成２０年５月２９日（２００８．５．２９）
【国際特許分類】
【出願日】平成１８年１１月８日（２００６．１１．８）
【出願人】（３０４０２１４１７）国立大学法人東京工業大学 (1,821)
【Ｆターム（参考）】
[ Back to top ]