JP3813343B2 - ３次元モデリング装置 - Google Patents
３次元モデリング装置 Download PDFInfo
- Publication number
- JP3813343B2 JP3813343B2 JP00598298A JP598298A JP3813343B2 JP 3813343 B2 JP3813343 B2 JP 3813343B2 JP 00598298 A JP00598298 A JP 00598298A JP 598298 A JP598298 A JP 598298A JP 3813343 B2 JP3813343 B2 JP 3813343B2
- Authority
- JP
- Japan
- Prior art keywords
- target object
- voxel
- standard deviation
- shape
- region
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes
Landscapes
- Physics & Mathematics (AREA)
- Geometry (AREA)
- Engineering & Computer Science (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
- Processing Or Creating Images (AREA)
- Length Measuring Devices By Optical Means (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、３次元モデリング装置およびその方法、３次元モデリングプログラムを記録した媒体、ならびに３次元形状推定プログラムを記録した媒体に関し、特に、対象物体の立体モデルを生成する３次元モデリング装置およびその方法ならびに３次元モデリングプログラムを記録した媒体に関する。
【０００２】
【従来の技術】
一般に、３次元グラフィックス画像を得るためには、３次元対象物体の形状およびテクスチャ（色彩および模様）を表現した立体モデルが作成される。
【０００３】
立体モデルを作成する従来の３次元モデリング方法としては、たとえばシルエット法やステレオ法がある。シルエット法は、対象物体を複数の視点から単眼カメラで撮影し、それにより得られた複数のシルエット画像に基づいて対象物体の立体モデルを生成する。
【０００４】
【発明が解決しようとする課題】
しかしながら、シルエット法は撮影された対象物体の輪郭に基づいて立体モデルを生成するため、対象物体の凹所を表現することができないという問題がある。
【０００５】
また、「財満久義他，『連続カラー画像を用いた双方向投票による物体形状と色の推定』，情報処理学会第５３回（平成８年後期）全国大会，２−３０１および３０２頁」には、色情報を投票して対象物体の形状を推定する方法が開示されている。この方法は、表面ボクセルにおいては同じ色が投票されるはずであるという仮定に基づいたものである。しかしながら、この方法は、すべての画像情報を一旦用いて色の投票を行なっているため、表面ボクセルに対し、裏方向の画像情報までもが投票されてしまう。したがって、通常、表面ボクセルにおいても投票された色が異なってしまう。この問題を解決するため、ボクセルに投票された候補色を絞り込む必要が生じている。この処理のため、この方法は、投票された色をすべて保持するための多大なメモリが必要となる。これに加え、候補色の中から色の同一性を抽出して表面ボクセルを決定するために、同じ部分が常に同じ色で撮影される必要がある。あるいは、同じ色であるか否かを判断するためのしきい値設定が必要になる。このしきい値設定は一般に困難な処理となる。
【０００６】
本発明は上記のような問題点を解決するためになされたもので、対象物体の凹所を表現可能な簡易かつ安価な３次元モデリング装置およびその方法ならびに３次元モデリングプログラムを記録した媒体を提供することを目的とする。
【０００７】
【課題を解決するための手段】
請求項１に係る３次元モデリング装置の発明は、対象物体の立体モデルを生成する３次元モデリング装置であって、前記対象物体を複数の視点から撮影し、複数の物体画像を得る物体画像獲得手段と、前記複数の物体画像のうち少なくとも２つの物体画像間において互いに対応する画像情報についての濃度情報及び／又は色情報の標準偏差に基づいて前記対象物体の形状を推定する形状推定手段と、前記推定された形状に基づいて前記立体モデルを生成する生成手段とを備え、前記形状推定手段は、複数の領域から構成される３次元仮想空間の領域に前記少なくとも２つの物体画像に対応する画像情報についての濃度情報及び／又は色情報を累積する画像情報累積手段と、前記領域に累積された画像情報についての濃度情報及び／又は色情報の標準偏差を算出する算出手段とを含む。
【０００８】
請求項２に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記複数の物体画像中の前記対象物体の輪郭に基づいて前記３次元仮想空間における前記対象物体の存在領域を概略的に推定する存在領域推定手段をさらに備え、前記画像情報累積手段は、前記３次元仮想空間の領域のうち前記存在領域内の領域に前記画像情報についての濃度情報及び／又は色情報を累積する。
【０００９】
請求項３に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記形状推定手段はさらに、前記対象物体の形状を推定するために、前記３次元仮想空間の領域のうち前記標準偏差の小さい領域を滑らかに結合する結合手段を含む。
【００１０】
請求項４に係る３次元モデリング装置の発明は、請求項３に記載の発明において、前記結合手段は、前記標準偏差の小さい領域を滑らかに結合するために、前記標準偏差の小ささと前記結合されるべき領域の滑らかさとを変数とする評価関数を最小化する。
【００１１】
請求項５に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記３次元仮想空間の座標空間を前記視点の１つから放射状に広がる座標軸を有する透視座標空間に変換する座標変換手段と、前記対象物体の形状を推定するために、前記透視座標空間において前記１つの視点を通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段とを含む。
【００１２】
請求項６に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記対象物体の形状を推定するために、前記３次元仮想空間において前記視点の１つを通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段を含む。
【００１３】
請求項７に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記複数の物体画像間において互いに対応する画像情報の標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含む。
【００１４】
請求項８に係る３次元モデリング装置の発明は、請求項２に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、前記累積された距離の総和が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む。
【００１５】
請求項９に係る３次元モデリング装置の発明は、請求項２に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、前記累積回数を計算する累積回数計算手段と、前記累積された距離の総和を前記計算された累積回数で除した距離平均値が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む。
【００１６】
請求項１０に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記部分表面の各々の内側にある領域に所定の第１の値を投票する投票手段と、前記投票された第１の値の総和が所定の第２の値よりも大きい領域を前記対象物体の内側と認定する認定手段とを含む。
【００１７】
請求項１１に係る３次元モデリング装置の発明は、請求項１から請求項１０のいずれかに記載の発明において、前記画像情報累積手段はさらに、前記物体画像の複数の画素を含む所定ブロックごとに前記画像情報についての濃度情報及び／又は色情報を累積する。
【００１８】
請求項１２に係る３次元モデリング装置の発明は、請求項１から請求項１１のいずれかに記載の発明において、前記領域はボクセルであり、前記３次元仮想空間はボクセル空間である。
【００５８】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して詳しく説明する。なお、図中同一または相当部分には同一符号を付してその説明は繰返さない。
【００５９】
［実施の形態１］
図１は、本発明の実施の形態１による３次元モデリング装置の全体構成を示す図である。図１を参照して、この実施の形態１による３次元モデリング装置は、回転テーブル１と、デジタルスチルカメラ２と、コンピュータ３と、ＣＤ−ＲＯＭ（コンパクトディスク−リードオンリメモリ）４とを備える。回転テーブル１は、その上に載置された対象物体（図示せず）を回転させるためのものである。ここでは回転テーブル１を用いたが、これに代えて対象物体の向きを変えることができるロボットアームなどを用いてもよい。デジタルスチルカメラ２は、立体モデルを得ようとする対象物体を複数の視点から撮影し、複数の物体画像を得るためのものである。ここでは静止画を得ることができるデジタルスチルカメラ２を用いたが、これに代えて動画を得ることができるビデオカメラなどを用いてもよい。ＣＤ−ＲＯＭ４には、コンピュータ３に、デジタルスチルカメラ２で得られた物体画像に基づいて対象物体の立体モデルを生成させるためのプログラムが記録されている。ここではプログラムを記録する媒体としてＣＤ−ＲＯＭ４を用いたが、これに代えてフロッピィディスク、ＣＤ−Ｒ、ＤＶＤ、メモリカードなどを用いてもよい。
【００６０】
なお、この３次元モデリング装置は、カメラ２を固定し、対象物体を回転させることによって対象物体を複数の視点から撮影するようにしているが、対象物体を固定し、カメラ２を移動することによって対象物体を複数の視点から撮影するようにしてもよい。ただし、カメラ２を移動する場合はカメラ２の撮影位置を求める必要がある。この３次元モデリング装置は、カメラ２を固定し、対象物体を回転させているため、各視点から対象物体までの距離は互いに等しいが、この距離は必ずしも等しくなくてもよい。距離が等しくない場合は、以下に詳述する方法で各視点から対象物体までの距離を考慮して種々の計算を行なう必要がある。
【００６１】
図２は、図１に示されたコンピュータ３の主要構成を示すブロック図である。図２を参照して、このコンピュータ３は、中央処理装置（ＣＰＵ）５と、リードオンリメモリ（ＲＯＭ）６と、ランダムアクセスメモリ（ＲＡＭ）７と、入出力インタフェース（Ｉ／Ｆ）８と、ＣＰＵ５、ＲＯＭ６、ＲＡＭ７およびＩ／Ｆ８を相互に結合するためのデータ／アドレスバス９とを備える。デジタルスチルカメラ２、およびＣＤ−ＲＯＭ４のためのドライブユニットはＩ／Ｆ８に接続され、これによりデジタルスチルカメラ２で得られた物体画像およびＣＤ−ＲＯＭ４に記録されたプログラムはＩ／Ｆ１０を介してＲＡＭ７に格納される。ＣＰＵ５は、ＲＯＭ６およびＲＡＭ７に格納されたプログラムに従って演算処理を行なう。
【００６２】
図３は、図１に示されたＣＤ−ＲＯＭ４に記録された３次元モデリングプログラムを示すフローチャートである。図４（ａ）は、図３中の対象物体および背景の撮影ステップＳ２を説明するための図である。図４（ｂ）は、図３中のシルエット画像の作成ステップＳ３を説明するための図である。図４（ｃ）は、図３中の形状推定処理ステップＳ４を説明するための図である。図４（ｄ）は、図３中のポリゴン作成ステップＳ５を説明するための図である。図４（ｅ）は、図３中のテクスチャマッピングステップＳ６を説明するための図である。
【００６３】
図１〜図４を参照して、まずステップＳ１では、回転テーブル１とカメラ２との間の相対的な位置関係を決定するためのキャリブレーションを行なう。
【００６４】
続いてステップＳ２では、対象物体１０およびその背景をカメラ２で撮影する。より具体的には、まず、対象物体１０を回転テーブル１の上に置かずに背景だけを１つの視点から撮影する。これにより１枚の背景画像が得られる。次に、対象物体１０を回転テーブル１の上に置き、回転テーブル１を回転させることによって対象物体１０を複数の視点から撮影する。たとえば回転テーブル１を２２．５°ずつ回転させると、図５に示されるように対象物体１０を１６の視点（０〜１５）から撮影することになるため、１６枚の物体画像Ａ０〜Ａ１５が得られる。回転テーブル１およびカメラ２はコンピュータ３内のＩ／Ｆ８に接続され、ＣＰＵ５によって制御される。また、得られた１枚の背景画像および１６枚の物体画像Ａ０〜Ａ１５はＩ／Ｆ８を介してＲＡＭ７に格納される。なお、ここでは対象物体１０を１６の視点から撮影しているが、視点の数は２以上であればよく、特に限定されるものではない。
【００６５】
続いてステップＳ３では、公知の画像間差分処理法に従ってシルエット画像が作成される。すなわち、物体画像Ａ０〜Ａ１５と背景画像との間で差分処理が行なわれ、図４（ｂ）に示されるように１６枚のシルエット画像Ｂ０〜Ｂ１５が作成される。差分処理においては、物体画像の色信号レベルと背景画像の色信号レベルとの差が各画素ごとに算出される。なお、ここでは画像間差分処理法を用いているが、これに代えて公知のステレオ法を用いてもよい。ステレオ法は、両眼カメラで対象物体１０を撮影し、その画像中の奥行き情報に基づいて対象物体１０を背景と区別するものである。
【００６６】
続いてステップＳ４では、形状推定処理が行なわれる。この形状推定処理のために、図６に示されるような直交座標系のボクセル空間１１が予め用意され、ＲＡＭ７に格納されている。ここでは直交座標系のボクセル空間１１を用いているが、これに代えて円筒座標系のボクセル空間を用いてもよい。ボクセル空間１１は、複数のボクセル１２から構成される。
【００６７】
図７は、図３に示された形状推定処理ステップＳ４の詳細を示すフローチャートである。まず、ステップＳ４０の概略形状推定処理では、複数のシルエット画像Ｂ１〜Ｂ１５に基づいてボクセル空間１１における対象物体１０の存在領域が概略的に推定される。次に、ステップＳ４１の詳細形状推定処理では、物体画像Ａ０〜Ａ１５間における濃度情報（輝度信号レベル）の標準偏差に基づいて対象物体１０の形状が詳細に推定される。
【００６８】
対象物体１０の存在領域を概略的に推定するステップＳ４０では公知のボーティング処理が行なわれる。すなわち、図８に示されるように、各シルエット画像Ｂ１，Ｂ２が各視点からボクセル空間１１に投影される。このときシルエット画像Ｂ１，Ｂ２中の物体部分１３が貫通するボクセル１２に対象物体１０が存在する可能性として１点が投票される。すべてのシルエット画像Ｂ１〜Ｂ１５がボクセル空間１１に投影されたとき、投票数が所定数よりも多いボクセル１２が対象物体１０の存在領域と推定される。ただし、図５に示されるような対象物体１０の凹所１４はシルエット画像Ｂ０〜Ｂ１５に現われないため、このボーティング処理によって推定された存在領域は対象物体１０の形状を概略的に表現するものにすぎない。すなわち、対象物体１０の凹所の真の表面１４は現れず、仮の表面１５が現れる。
【００６９】
通常、対象物体１０の真の表面１４はいずれの視点からも同じ色および濃度で観察されるが、仮の表面１５は視点によって異なった色または濃度で観察される。したがって、同じ部分でも視点によって異なった色または濃度で観察される部分は対象物体１０の凹所である可能性が高い。
【００７０】
そこで、対象物体１０の形状を詳細に推定するステップＳ４１では、対象物体１０の形状を概略的に推定した後、濃度情報の標準偏差に基づいて対象物体１０の形状、特に凹所の形状が詳細に推定される。より具体的には、図９に示されるように、１６枚の物体画像Ａ０〜Ａ１５のうち、中心画像として物体画像Ａ０、および周辺画像として物体画像Ａ０に隣接する物体画像Ａ１，Ａ１５を用いて形状推定を行なう場合、物体画像Ａ０，Ａ１，Ａ１５ともに、その画素の濃度情報がボクセル空間１１に投影され、各画素からの投影線１６が貫通するすべてのボクセル１２にその画素の濃度情報が累積される（色ボーティング処理ステップＳ４１０）。換言すれば、すべてのボクセル１２について、物体画像Ａ０，Ａ１，Ａ１５への投影点が求められ、その投影点の濃度情報が対応するボクセル１２に累積される。そして、各ボクセル１２に累積された濃度情報の標準偏差が算出される。標準偏差が大きいボクセル１２は、それに対応する部分が視点によって異なった濃度で観察されることを意味するので、そのようなボクセル１２は対象物体１０の凹所である可能性が高い。そのため、ボクセル空間１１内のボクセル１２のうち標準偏差が小さいボクセルを滑らかに結合することにより、図５に示された視点（０）から観察可能な対象物体１０の部分表面が推定される（部分表面推定処理ステップＳ５２）。
【００７１】
上述した色ボーティング処理ステップＳ４１０および部分表面推定処理ステップＳ５２は、物体画像Ａ０〜Ａ１５の各々を中心画像として行なわれる。したがって、推定された１６枚の部分表面が統合され、これにより対象物体１０の立体形状が得られる（部分表面統合処理ステップＳ４１２）。
【００７２】
上記ステップＳ４に続くステップＳ５では、ステップＳ４で得られた対象物体１０の立体形状に基づいて三角パッチからなるポリゴン１７が作成される。したがって、対象物体１０の立体形状は多数のポリゴン１７で表現される。ポリゴン１７で表現された立体形状はＲＡＭ７に格納される。
【００７３】
最後にステップＳ６では、ステップＳ５で作成された各ポリゴン１７に対応するテクスチャがマッピングされる。これにより、対象物体１２をモデル化した立体モデル１８が作成される。テクスチャは物体画像Ａ０〜Ａ１５から得られる色情報で、ＲＡＭ９に格納される。
【００７４】
なお、ここでは回転テーブル１を用いて対象物体１０を複数の視点から撮影して複数の物体画像Ａ０〜Ａ１５を得ているが、対象物体１０の立体形状を推定するために必要な物体画像を得ることさえできれば、上記ステップＳ１〜Ｓ３は適宜他の公知の手法に置換えてもよい。また、物体画像および推定された立体形状に基づいて立体モデル１８を生成することさえできれば、上記ステップＳ５およびＳ６は適宜他の公知の手法に置換えてもよい。また、対象物体１０の輪郭に基づいてボクセル空間１１における対象物体１０の存在領域を概略的に推定することさえできれば、上記ステップＳ４中の前半のステップは適宜他の公知の手法に置換えてもよい。
【００７５】
次に、上記のように構成された３次元モデリング装置を用いて対象物体１０の立体モデル１８を生成する方法について説明する。
【００７６】
まずステップＳ１で上述したようにキャリブレーションを行なった後、ステップＳ２で１つの視点から背景を撮影し、続いて複数の視点から対象物体１０を撮影する。これにより１枚の背景画像および複数の物体画像Ａ０〜Ａ１５が得られる。続いてステップＳ３で、物体画像Ａ０〜Ａ１５と背景画像との間で差分処理が行なわれ、これによりシルエット画像Ｂ０〜Ｂ１５が得られる。
【００７７】
続いて、ステップＳ４０でシルエット画像Ｂ０〜Ｂ１５に基づいてボーティング処理などの概略形状推定処理が行なわれた後、ステップＳ４１で物体画像Ａ０〜Ａ１５の濃度情報に基づいて詳細形状推定処理が行なわれる。図１０は、図７中の概略形状推定処理ステップＳ４０、色ボーティング処理ステップＳ４１０および部分表面推定処理ステップＳ５２の詳細を示すフローチャートである。なお、図１０中のステップＳ４２１〜ステップＳ５１が図７中の色ボーティング処理ステップＳ４１０に対応する。図１１は、図７中の部分表面統合処理ステップＳ４１２の詳細を示すフローチャートである。
【００７８】
まずステップＳ４０で、シルエット画像Ｂ０〜Ｂ１５に基づいて上述したボーティング処理などの概略形状推定処理が行なわれ、対象物体１０の初期立体形状が概略的に抽出される。より具体的には、投票数が所定数よりも多いボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域と推定され、そのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）が「１」に設定される。他方、投票数が所定数よりも少ないボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域と推定されず、そのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）が「０」に設定される。したがって、ｖ（ｘ，ｙ，ｚ）＝１のボクセルが対象物体１０の概略的な存在領域を形成する。
【００７９】
続いてステップＳ４２で、図５に示された視点（０）から対象物体１０を撮影して得られた物体画像が中心画像ｋとして選択される。
【００８０】
続いてステップＳ４２１で、変数（ｃ０，ｃ１，ｃ２）が初期化される。
続いてステップＳ４３で、中心画像ｋ（＝０）から図５上で時計方向にｎ個離れた視点（１６−ｎ（「１６」は図５上の「０」に相当する））から対象物体１０を撮影して得られた物体画像が周辺画像ｉとして選択される。たとえばｎ＝２とすると、図５に示された視点（１４）から対象物体１０を撮影して得られた物体画像（周辺画像）Ａ１４が選択される。
【００８１】
続いてステップＳ４４で、存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。
【００８２】
ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ４５で、図１２に示されるように座標（ｘ，ｙ，ｚ）に位置するボクセル１２の画像ｉ（ここではｉ＝１４）への投影点（Ｘ，Ｙ）が計算される。
【００８３】
続いてステップＳ４６で、後述するステップＳ５０で標準偏差ｓ（ｘ，ｙ，ｚ）を計算するために必要な値が計算される。より具体的には、変数ｃ０（ｘ，ｙ，ｚ）をインクリメントすることによってボクセル１２に累積される濃度情報の数がカウントされる。また、変数ｃ１（ｘ，ｙ，ｚ）に投影点（Ｘ，Ｙ）に位置する画素の濃度情報ｃ（Ｘ，Ｙ）が加算される。また、変数ｃ２（ｘ，ｙ，ｚ）にその濃度情報の自乗ｃ（Ｘ，Ｙ）2 が加算される。なお、ステップＳ４４でｖ（ｘ，ｙ，ｚ）＝０の場合、ステップＳ４５およびＳ４６の処理は行なわれない。
【００８４】
続いてステップＳ４７で、上記ステップＳ４４〜Ｓ４６の処理がすべてのボクセルについて行なわれたか否かが判定される。未だすべてのボクセルについて上記ステップＳ４４〜Ｓ４６の処理が行なわれていない場合、再び対象物体１０の存在領域内の各ボクセルの変数ｃ１（ｘ，ｙ，ｚ）およびｃ２（ｘ，ｙ，ｚ）に画像ｉ（＝１４）中の対応する投影点（Ｘ，Ｙ）の濃度情報ｃ（Ｘ，Ｙ）およびその自乗ｃ（Ｘ，Ｙ）2 がそれぞれ加算される。
【００８５】
すべてのボクセルについて上記ステップＳ４４〜Ｓ４６の処理が行なわれた場合、ステップＳ４８で、上記ステップＳ４４〜Ｓ４７の処理が中心画像ｋ（＝０）およびすべての周辺画像ｉ（＝１４，１５，１，２）について行なわれたか否かが判定される。最初は周辺画像ｉ（＝１４）について上記ステップＳ４４〜Ｓ４７の処理が行なわれたにすぎないが、ステップＳ４９で変数ｉがインクリメントされ、対象の周辺画像ｉがｉ＝１４からｉ＝１５に変更される。したがって、上記と同様に、対象物体１０の存在領域内のボクセル１２に周辺画像ｉ（＝１５）中の対応する投影点（Ｘ，Ｙ）の濃度情報ｃ（Ｘ，Ｙ）およびｃ（Ｘ，Ｙ）2 が累積される。
【００８６】
変数ｉの値が変数（ｋ＋ｎ）の値（ここでは２）を超えたとき、中心画像ｋ（＝０）および周辺画像ｉ（＝１４，１５，１，２）の濃度情報がすべて対象物体１０の存在領域内のボクセル１２に累積されたことになる。
【００８７】
ステップＳ４８で変数ｉの値が変数（ｋ＋ｎ）の値よりも大きいと判定されると、ステップＳ４９０で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５０で、ステップＳ４６で計算されたｃ０，ｃ１およびｃ２に基づいて各ボクセル１２に累積された濃度情報の標準偏差ｓ（ｘ，ｙ，ｚ）が計算される。
【００８８】
続いてステップＳ５１で、すべてのボクセルについて標準偏差ｓ（ｘ，ｙ，ｚ）が計算されたか否かが判定される。
【００８９】
すべてのボクセルについて標準偏差ｓ（ｘ，ｙ，ｚ）が計算された場合、ステップＳ５２でボクセル空間１１内において各視点（０〜１５）から観察した対象物体１０の部分表面Ｓｐ（ｋ）が推定される。上述したように、対象物体１０の凹所１４内に位置するボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）は一般に大きくなり、対象物体１０の表面に位置するボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）は小さくなる。したがって、同一のｘｙ座標をもつボクセルのうち標準偏差ｓ（ｘ，ｙ，ｚ）が最小のボクセルを結合することによって対象物体１０の表面を推定することができる。ただし、照明条件が変動したり、凹所１４の表面で同じ色が連続する場合、ノイズの影響を受け、対象物体１０の表面に位置しないボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が最小になることがある。したがって、単純に標準偏差ｓ（ｘ，ｙ，ｚ）が最小のボクセルを結合した場合、必ずしも正確な形状を推定できるとは限らない。
【００９０】
そこで、標準偏差ｓ（ｘ，ｙ，ｚ）だけでなく推定される形状の連続性も考慮して対象物体１０の形状を推定するのが望ましい。より具体的には、標準偏差ｓ（ｘ，ｙ，ｚ）および隣接するボクセル間の距離を変数とするエネルギ関数（評価関数）を規定し、そのエネルギ関数の値が最小となるボクセルを対象物体１０の表面に位置すると推定し、これらのボクセルによって構成される面を部分表面Ｓｐ（ｋ）とする。以下、この方法を「最小エネルギ法」と呼ぶ。なお、このエネルギ関数は、標準偏差ｓ（ｘ，ｙ，ｚ）、および隣接するボクセル空間の傾きと、これに対応する概略形状での表面の傾きとの差を変数とするものであってもよい。このような部分表面の推定は、たとえばダイナミックプログラミング（ＤＰ）法、緩和法などによって行なわれる。ダイナミックプログラミング法については後に図１３を用いて詳述する。
【００９１】
上記のように、中心画像ｋ（＝０）および周辺画像ｉ（＝１４，１５，１，２）に基づいて部分表面Ｓｐ（ｋ）が対象物体１０の表面候補として推定されるが、他の視点（１〜１５）からの画像を中心画像としたときも同様に部分表面Ｓｐ（ｋ）が対象物体１０の表面候補としてそれぞれ推定される（Ｓ５４１，Ｓ５４２）。したがって、計算された複数の部分表面Ｓｐ（ｋ）を統合して対象物体１０の真の表面を決定する必要がある。
【００９２】
最初にステップＳ４０で概略的に推定された対象物体１０の存在領域内に位置しかつステップＳ５２で算出された部分表面Ｓｐ（ｋ）よりも内側に位置するボクセルは対象物体１０の真の存在領域である可能性が高い。これに対し、ステップＳ４０で概略的に推定された対象物体１０の存在領域内に位置しかつステップＳ５２で計算された部分表面Ｓｐ（ｋ）よりも外側に位置するボクセルは対象物体１０の真の存在領域でない可能性がある。その可能性はステップＳ５２で推定された部分表面Ｓｐ（ｋ）から遠く離れたボクセルほど高い。そこで、部分表面Ｓｐ（ｋ）よりも外側に位置する各ボクセルに部分表面Ｓｐ（ｋ）までの距離をペナルティとして与える。そして、このペナルティの平均値が所定の値よりも大きいボクセルを対象物体１０の存在領域以外の領域とする。
【００９３】
より具体的には、図１１に示されるように、まずステップＳ５４３で変数（ｐ０，ｐ１）が初期化され、続いてステップＳ５４４で図５に示された視点（０）から対象物体１０を撮影して得られた物体画像Ａ０が中心画像ｋとして選択される。続いて、ステップＳ５４５で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５３で、ボクセル（ｘ，ｙ，ｚ）がステップＳ５２で計算された部分表面Ｓｐ（ｋ）よりも外側に位置するか否かが判定される。
【００９４】
外側に位置する場合、変数ｐ０（ｘ，ｙ，ｚ）がインクリメントされる（Ｓ５４）。これにより、そのボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの数がカウントされる。また、変数ｐ１（ｘ，ｙ，ｚ）に距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））が加算される（Ｓ５４）。距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））はボクセル（ｘ，ｙ，ｚ）からステップＳ５２で推定された対象物体１０の部分表面Ｓｐ（ｋ）までの距離である。これにより、部分表面Ｓｐ（ｋ）から遠いボクセルほど大きなペナルティが与えられる。
【００９５】
ステップＳ５５で、すべてのボクセルについて上記ステップＳ５３〜Ｓ５４の処理が行なわれたか否かが判定され、行なわれた場合は、ステップＳ５６で図５に示されたすべての視点（０〜１５）から撮影された画像を中心画像として上記ステップＳ４３〜Ｓ５５の処理が行なわれたか否か、つまり変数ｋがその最大値ｋｍａｘ（ここでは１５）に達したか否かが判定される。変数ｋがその最大値ｋｍａｘよりも小さい場合、ステップＳ５７で変数ｋがインクリメントされ、再びステップＳ５３〜Ｓ５５の処理が行なわれる。
【００９６】
他方、変数ｋがその最大値ｋｍａｘよりも小さくない場合、ステップＳ５７０で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５８〜Ｓ６０でペナルティの高いボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域ではないとされる。より具体的には、ステップＳ５８で各ボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きいか否かが判定される。ここで、ｐ０（ｘ，ｙ，ｚ）はボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの数である。
【００９７】
ペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きい場合、ステップＳ５９でそのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）は０に設定される。ペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きくない場合、ステップＳ５９の処理は行なわれない。
【００９８】
続いてステップＳ６０で、すべてのボクセルについてステップＳ５８およびＳ５９の処理が行なわれたか否かが判定される。したがって、所定の値よりも大きい平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）を有するボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域から除外される。これにより、複数の部分表面Ｓｐ（ｋ）が統合され、対象物体１０の凹所が表現され得る。以下、図１１に示された部分表面Ｓｐ（ｋ）を統合する方法を「距離ペナルティ法」と呼ぶ。
【００９９】
上記のようにして得られた物体画像Ａ０〜Ａ１５および対象物体１０の形状に基づいて、図３に示されたステップＳ５でポリゴン１７が作成され、ステップＳ６でテクスチャがマッピングされる。これにより、対象物体１０の立体モデル１８が作成される。
【０１００】
図１３は、図１０に示されたステップＳ５２において部分表面Ｓｐ（ｋ）をダイナミックプログラミングによって計算する方法を示すフローチャートである。なお、図示されていないが、図１３の処理を行なう前に、ボクセル空間を中心画像ｋに応じて、後述の図２２に示された座標変換を行なっておくのが望ましい。図１４は、図１０および図１１に示されたフローチャートに従って計算されたＹ＝ｙ０平面における標準偏差ｓを示す鳥瞰図である。図１５は、ボクセル空間１１のＹ＝ｙ０平面内のボクセルの標準偏差ｓを示す図である。図１３に示されたダイナミックプログラミング法は、図１４および図１５上で、滑らかさを保ちつつ可能な限り標準偏差ｓの小さい箇所を結ぶものである。
【０１０１】
図１３を参照して、ステップＳ５２０で座標ｙにその最小値ｙｍｉｎが設定され、続いてステップＳ５２１でｘにその最小値ｘｍｉｎが設定され、座標ｚにその最小値ｚｍｉｎが設定される。
【０１０２】
続いてステップＳ５２２〜Ｓ５２４で、Ｘ＝ｘｍｉｎの列（図１５上で最も左側の列）のボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が仮のエネルギ関数値用配列ｅ（ｚ）の初期値としてそれぞれ設定される。
【０１０３】
続いてステップＳ５２５でｘがインクリメントされ、続いてステップＳ５２６で座標ｚがその最小値ｚｍｉｎに設定される。
【０１０４】
続いてステップＳ５２７〜Ｓ５２９で、Ｘ＝ｘｍｉｎ＋１の列のエネルギ関数ｅ（ｚ）が計算される。ステップＳ５２７において、ＭＩＮ（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０），ｚ０＝ｚｍｉｎ〜ｚｍａｘ）は、現在の座標ｘ（ここではｘ＝ｘｍｉｎ＋１）において、求めようとする対象物体１０の表面が座標ｚであると仮定し、かつ、座標ｘ−１における表面が座標ｚ０（ｚ０はｚｍｉｎ〜ｚｍａｘのいずれか）であると仮定したとき、評価値（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０））のうち最小のもの（ｍｉｎ）である。ｄｉｆ（ｚ，ｚ０）は、たとえば、表面ｘの座標ｚから表面ｘ−１の座標ｚ０までの距離である。あるいは、ｄｉｆ（ｚ，ｚ０）は、表面ｘの座標ｚから表面ｘ−１の座標ｚ０までの傾きと、概略表面ｘの座標ｚ’から概略表面ｘ−１の座標ｚ０’までの傾きとの差である。このＭＩＮ（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０））にＸ＝ｘｍｉｎ＋１の列の各ボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が加算され、各ボクセルについて仮のエネルギ関数値ｅ（ｚ）が計算される。
【０１０５】
また、上記最小値ｍｉｎを与えるような座標ｚ０の値が変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）に設定される。これにより、Ｘ＝ｘｍｉｎ−１において、座標ｚの値をｚと仮定したときの、仮のエネルギ関数値ｅ（ｚ）が最小になるための座標ｘ−１における座標ｚがそれぞれ記憶される。
【０１０６】
続いてステップＳ５３０およびＳ５３１で、上記ステップＳ５２７〜Ｓ５２９の処理がｘｍｉｎ〜ｘｍａｘ−１のすべての座標ｘに対応するｚ座標列について行なわれる。したがって、エネルギ関数ｅ（ｚ）にはｘｍｉｎから現在の座標ｘまでの各ｚに対する最小エネルギ関数値が記録されている。
【０１０７】
上記の結果、計算されたエネルギ関数ｅ（ｚ）のうち最小値を与える座標ｚが算出される。この座標ｚがｘｍａｘ−１の列における対象物体１０の表面の位置である。ここで、変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）を見れば、現在よりも１つ前の座標ｘ−１における表面のｚ座標が記録されている。したがって、座標ｘを順次デクリメントし、その座標ｘを変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）に代入すれば、−Ｘ方向に向かって順次対象物体１０の表面の位置を辿っていくことができる。このようにして計算された部分表面Ｓｐ（ｋ）が図１６中に示されている。図１６から明らかなように、部分表面Ｓｐ（ｋ）と決定されたボクセルの標準偏差は各列において必ずしも最小のものではない。仮に各列において最小の標準偏差を持つボクセルが部分表面Ｓｐ（ｋ）を構成するものとして選択されたならば、その部分表面Ｓｐ（ｋ）の連続性が悪くなる。ここでは、部分表面Ｓｐ（ｋ）の連続性を考慮して表面を構成するボクセルが選択されている。
【０１０８】
続いてステップＳ５３２およびＳ５３３で、上記ステップＳ５２１〜Ｓ５３１の処理がｙｍｉｎ〜ｙｍａｘのすべての座標ｙについて行なわれる。
【０１０９】
図１７は、図１０に示されたステップＳ４０で得られた初期のＹ＝ｙ０平面におけるボクセルを示す図である。図１７中において「＊＊」で示されたボクセルは対象物体１０の概略的な存在領域に位置するものである。図１１に示されたステップＳ５４では、図１７では対象物体１０の内部であるが、図１６では外部である各ボクセルから図１６に示された対象物体１０の部分表面Ｓｐ（ｋ）までの距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））が計算される。
【０１１０】
以上のようにこの実施の形態１によれば、物体画像Ａ０〜Ａ１５間における濃度情報の標準偏差に基づいて対象物体１０の形状を推定しているため、対象物体１０の凹所を表現した立体モデル１８を作成することができる。
【０１１１】
また、標準偏差に基づいて対象物体１０の形状を推定する前に、シルエット画像Ｐ０〜Ｐ１５に基づいてボーティング処理により対象物体１０の存在領域を概略的に推定しているため、対象物体１０の立体モデル１８を迅速に生成することができる。
【０１１２】
また、対象物体１０の表面を標準偏差が小さいボクセルを滑らかに結合することによって推定しているため、照明条件などが変動しても安定して正確な立体モデル１８を生成することができる。
【０１１３】
［実施の形態２］
図１８は、本発明の実施の形態２による３次元モデリング装置における形状推定処理（部分表面統合処理を除く）を示すフローチャートである。この実施の形態２では、図１０に示されたステップＳ４５およびＳ４６に代えてステップＳ７０〜Ｓ７４の処理が行なわれる。また、図１０に示されたステップＳ５０に代えてステップＳ７５の処理が行なわれる。また、図１０に示されたステップＳ４２１では変数ｃ２が初期化されているが、ここでは変数ｃ２は用いられていないので初期化されない。その他は、図１０に示されたものと同じ処理が行なわれる。
【０１１４】
すなわち、上記実施の形態１では図１２に示されるようにボクセル１２の画像ｉへの投影点（Ｘ，Ｙ）の画素の濃度情報がボクセル１２に累積されるが、この実施の形態２では図１９に示されるようにボクセル１２の中心画像ｋへの投影点（Ｘｋ，Ｙｋ）を含む複数の画素から構成されるブロック領域Ｂｋがボクセル１２の周辺画像ｉへの投影点（Ｘｉ，Ｙｉ）を含む複数の画素から構成されるブロック領域Ｂｉと図２０に示されるように比較され、濃度情報が異なる画素の数がボクセル１２に累積される。
【０１１５】
より具体的には、ステップＳ７０でボクセル（ｘ，ｙ，ｚ）の中心画像ｋへの投影点（Ｘｋ，Ｙｋ）が計算される。続いてステップＳ７１で、投影点（Ｘｋ，Ｙｋ）を含むｍ×ｎ個の画素からなるブロック領域Ｂｋが設定される。続いてステップＳ７２で、ボクセル（ｘ，ｙ，ｚ）の周辺画像ｉへの投影点（Ｘｉ，Ｙｉ）が計算される。続いてステップＳ７３で、投影点（Ｘｉ，Ｙｉ）を含むｍ×ｎ個の画素からなるブロック領域Ｂｉが設定される。
【０１１６】
続いてステップＳ７４で、ｃ０（ｘ，ｙ，ｚ）がインクリメントされ、これにより各ボクセル（ｘ，ｙ，ｚ）に累積される濃度情報が異なる画素の数がカウントされる。また、ブロック領域Ｂｋがブロック領域Ｂｉと比較され、濃度ＬＭｋ，ＬＭｉの差の絶対値｜ＬＭｋ−ＬＭｉ｜の総和ｄｉｆ（Ｂｋ，Ｂｉ）がｃ１（ｘ，ｙ，ｚ）に加算される。たとえば図２０に示されるように、４×４個の画素からなるブロック領域ＢｋおよびＢｉを設定する場合、１６の画素のうち、３つの画素の濃度情報がブロック領域Ｂｋとブロック領域Ｂｉとの間で相違し、その差の絶対値の合計は「４」であるため、ｄｉｆ（Ｂｋ，Ｂｉ）として「４」がｃ１（ｘ，ｙ，ｚ）に加算される。
【０１１７】
上記の結果、濃度差の絶対値の総和の累積値がｃ１（ｘ，ｙ，ｚ）に設定される。
【０１１８】
その後、ステップＳ７５で、各ボクセル（ｘ，ｙ，ｚ）に累積された濃度差の絶対値の総和の平均値ｓ（ｘ，ｙ，ｚ）が計算される。この平均値ｓ（ｘ，ｙ，ｚ）が濃度情報のばらつき評価値として上記標準偏差の代わりに用いられる。
【０１１９】
以上のようにこの実施の形態２によれば、物体画像間における濃度情報のばらつきをブロックＢｋ，Ｂｉごとに評価しているため、対象物体１０の立体モデル１８をより精度良く生成することができる。
【０１２０】
なお、本実施の形態では、物体画像間における濃度情報のばらつきをブロック間で評価する際、ブロック内の画素ごとの差の絶対値｜ＬＭｋ−ＬＭｉ｜の総和を用いたが、ブロック内の画素ごとの差の自乗（ＬＭｋ−ＬＭｉ）2 の総和を用いてもよい。さらには、ブロック内の画素ごとの標準偏差の平均値を用いてもよい。
【０１２１】
［実施の形態３］
図２１は、本発明の実施の形態３による３次元モデリング装置における部分表面推定処理を示すフローチャートである。上記実施の形態では、できるだけ正確な部分表面Ｓｐ（ｋ）を得るために図１３に示されたダイナミックプログラミング法を用いているが、計算量が多くなる。したがって、以下に詳述する実施の形態３のように、単純に標準偏差が最小のボクセルを部分表面Ｓｐ（ｋ）と推定してもよい。
【０１２２】
まずステップＳ３００で、各中心画像ｋに応じてボクセル空間１１が座標変換される。ここでは、図２２に示されるように、直交座標系のボクセル空間１１が回転され、さらに並進されて再び直交座標系のボクセル空間３０に変換される。座標変換前のボクセル空間１１の各ボクセルの座標を（ｘ，ｙ，ｚ）とし、座標変換後のボクセル空間３０の各ボクセルの座標を（ｘ′，ｙ′，ｚ′）とし、回転変換のためのオイラー角をそれぞれφ、θ、ψとすると、この座標変換は次の式（１）で表わされる。
【０１２３】
【数１】
【０１２４】
なお、式（１）中の（ｘｓ，ｙｓ，ｚｓ）は、変換前のボクセル空間１１の中心から中心画像ｋの視点へのベクトルである。回転変換後のボクセル空間はこのベクトル（ｘｓ，ｙｓ，ｚｓ）に応じて並進され、変換後のボクセル空間３０となる。
【０１２５】
このような座標変換が行なわれた後、ステップＳ３０１で座標ｙ′が最小値ｙ′ｍｉｎに初期化され、ステップＳ３０２で座標ｘ′が最小値ｘ′ｍｉｎに初期化され、さらにステップＳ３０３で座標ｚ′が最小値ｚ′ｍｉｎに初期化され、標準偏差一時変数ｓｔｍｐが無限大∞に初期化され、かつ座標一時変数ｚｔｍｐが無限大∞に初期化される。
【０１２６】
続いて、ステップＳ３０３０で存否変数ｖ（ｘ′，ｙ′，ｚ′）が「０」か否か、つまりそのボクセル（ｘ′，ｙ′，ｚ′）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ′，ｙ′，ｚ′）＝０ではない場合、つまりｖ（ｘ′，ｙ′，ｚ′）＝１の場合、ステップＳ３０４で、標準偏差一時変数ｓｔｍｐの値が座標（ｘ′，ｙ′，ｚ′）のボクセルの標準偏差ｓ（ｘ′，ｙ′，ｚ′）と比較される。標準偏差一時変数ｓｔｍｐの値がそのボクセルの標準偏差ｓ（ｘ′，ｙ′，ｚ′）よりも大きい場合は、ステップＳ３０５で標準偏差一時変数ｓｔｍｐに改めてその標準偏差ｓ（ｘ′，ｙ′，ｚ′）が設定され、そのときの座標ｚ′が位置一時変数ｚｔｍｐに設定される。
【０１２７】
ステップＳ３０６およびＳ３０７により、すべての座標ｚ′（＝ｚ′ｍｉｎ〜ｚ′ｍａｘ）について上記ステップＳ３０４および３０５の処理が行なわれる。その結果、変換後のボクセル空間３０において同一のｘ′ｙ′座標をもつ複数のボクセルのうち標準偏差が最小のボクセルが特定される。このようにして特定された標準偏差が最小のボクセルのｚ′座標はステップＳ３０８で表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定される。この表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）は各中心画像ｋに対応して設定される。すなわち、この表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）は、各中心画像ｋ上に撮影されている対象物体の表面位置を示すものである。なお、この時点で位置一時変数ｚｔｍｐが∞であれば、表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に対応した物体表面は存在しないことを意味する。
【０１２８】
図２３は、変換後のボクセル空間３０におけるＹ′＝ｙ′平面を示す。ボクセル空間３０内にはボーティング処理により推定された対象物体１０の概略形状３１が示されている。この概略形状３１内の各ボクセルには標準偏差ｓ（ｘ′，ｙ′，ｚ′）が与えられている。上述したステップＳ３０４〜Ｓ３０７は、たとえば同一のｘ′ｙ′座標をもつ複数のボクセル３２の中から標準偏差ｓ（ｘ′，ｙ′，ｚ′）が最小のボクセルを特定するものである。ここでは、ｚ′＝３７のボクセルが最小の標準偏差を有している。上述したステップＳ３０８は、この最小の標準偏差をもつボクセルの座標ｚ′を表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）の対応する座標に与えるものである。
【０１２９】
続いてステップＳ３０９およびＳ３１０により、Ｙ′＝ｙ′平面上のすべての座標ｘ′（＝ｘ′ｍｉｎ〜ｘ′ｍａｘ）について上記ステップＳ３０３〜Ｓ３０８の処理が行なわれる。また、ステップＳ３１１およびＳ３１２により、すべての座標ｙ′（＝ｙ′ｍｉｎ〜ｙ′ｍａｘ）について上記ステップＳ３０２〜Ｓ３１０の処理が行なわれる。これにより、対象物体１０の部分表面Ｓｐ（ｋ）が決定され、その部分表面Ｓｐ（ｋ）のボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる。
【０１３０】
上述した方法は図１３のダイナミックプログラミング法に代えて用いられるもので、以下「座標変換による最小ばらつき探索法」と呼ぶ。
【０１３１】
この方法により各中心画像ｋごとに得られた部分表面Ｓｐ（ｋ）は上述した距離ペナルティ法または後に詳述する単純ボーティング法により１つの表面に統合され、これにより対象物体１０の形状が決定される。たとえば図２３に示されるように、●で示されるボクセルは対象物体１０の内部にあると判断され、○で表わされるボクセルは対象物体１０の外部にあると判断される。
【０１３２】
このように単純に最小の標準偏差をもつボクセルを特定することにより対象物体１０の形状を推定すれば、上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができる。
【０１３３】
［実施の形態３の変形例］
上記実施の形態においては、直交座標系のボクセル空間１１で同一のｘｙ座標をもつ複数のボクセルのうち標準偏差が最小のボクセルが必ずしも対象物体１０の真の表面１４に位置していない。これは画像ノイズに起因するものもあるが、主として計算方法に起因するものである。この理由を図９および図２４を用いて詳細に説明する。
【０１３４】
図９に示されるように、直交座標系のボクセル空間１１において、同一のｘｙ座標をもつすべてのボクセル１２から中心画像となる物体画像Ａ０への投影線１６が物体画像Ａ０上の同じ位置に投影する場合は、照明条件の著しい変動がない限り、真の表面に位置するボクセルの物体画像Ａ０への投影点の色彩が複雑であるか単一であるかに関わらず、その同一のｘｙ座標をもつ複数のボクセル１２のうち標準偏差が最小のボクセル１２が対象物体１０の真の表面上に位置すると考えてもよい。しかしながら、上述したボクセル１２から中心画像となる物体画像Ａ０への投影線１６は、透視投影に基づいているため、図２４に示されるように、同一のｘｙ座標をもつ複数のボクセル１２から物体画像Ａ０への投影線１６は物体画像Ａ０上の異なる位置に投影することが多い。この場合、真の表面１４上に位置するボクセル１２から物体画像Ａ０上への投影点付近には複雑な色彩があり、その他のボクセル１２から中心画像である物体画像Ａ０および周辺画像である物体画像Ａ１およびＡ１５上への投影点には単一の色彩があるとすると、真の表面１４上に位置しないボクセルの標準偏差が真の表面１４上に位置するボクセルの標準偏差よりも小さくなる場合がある。
【０１３５】
そこで、図２１に示されたステップＳ３００における直交座標系から直交座標系への座標変換に代えて、直交座標系から図２５に示されるような透視座標系への座標変換が行なわれてもよい。この透視座標系のボクセル空間３５は、中心画像ｋの視点を頂点とした錐体状をなしている。
【０１３６】
図２６（ａ）は図２５に示された透視座標系のボクセル空間３５をＹ方向から見た図であり、図２６（ｂ）は図６に示された直交座標系のボクセル空間１１をＹ方向から見た図である。中心画像ｋの視点からボクセル空間１１の中心までの距離をｌとすると、この直交座標系から透視座標系への座標変換は次の式（２）により表わされる。
【０１３７】
ｘ′＝ｘ・ｌ／ｚ，ｙ′＝ｙ・ｌ／ｚ，ｚ′＝ｚ …（２）
このような透視座標系のボクセル空間３５においても上記と同様に、図２７に示されるような同一のｘ′ｙ′座標をもつ複数のボクセル３６の中から標準偏差が最小のボクセルが特定され、その特定されたボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる。ただし、ここで表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる値はこの変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）が設定される中心画像からその標準偏差が最小のボクセルまでの距離である。このように、透視座標系のボクセル空間３５において各中心画像ｋの視点を通る直線状のボクセルのうち標準偏差が最小のボクセルが結合され、対象物体１０の部分表面Ｓｐ（ｋ）が決定される。
【０１３８】
なお、部分表面統合処理により、図２７中の●のボクセルは対象物体１０の内部と判定され、○のボクセルは対象物体１０の外部と判定される。
【０１３９】
このように透視座標系のボクセル空間に座標変換を行えば、最小の標準偏差をもつボクセルを特定することにより対象物体１０の形状を簡単かつ正確に推定することができる。また、このような座標変換による最小ばらつき法は上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができ、その結果、この処理に必要なメモリのワークエリアを削減するこができる。
【０１４０】
［実施の形態４］
図２８は、本発明の実施の形態４による３次元モデリング装置における部分表面推定処理を示すフローチャートである。ここでは、図２１に示されたステップＳ３０８に代えて、ステップＳ３２０〜Ｓ３２４が設けられている。この方法は上記実施の形態３の部分表面Ｓｐ（ｋ）上またはそれよりも外側に部分表面を特定するもので、以下「座標変換による準最小ばらつき探索法」と呼ぶ。
【０１４１】
上述したようにステップＳ３０４〜Ｓ３０７で同一のｘ′ｙ′座標をもつボクセルの中から標準偏差が最小のボクセルが特定される。
【０１４２】
ステップＳ３２０では、その最小の標準偏差に所定の値αを加えたものが最小標準偏差変数ｓｍｉｎに設定される。そして、ステップＳ３２０〜Ｓ３２３で、再び同一のｘ′ｙ′座標をもつ複数のボクセルの中から［最小の標準偏差＋α］以下の標準偏差をもつボクセルが特定される。このとき、ボクセルはｚ′ｍｉｎからｚ′ｍａｘに向かって順に探索される。したがって、最小の標準偏差をもつボクセルよりも外側に、［最小の標準偏差＋α］以下の標準偏差をもつボクセルが存在する場合は、そのボクセルの座標ｚ′がステップＳ３２４で表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定される。この場合、図２３または図２７に示された部分表面Ｓｐ（ｋ）よりも外側に部分表面が特定される。なお、最小の標準偏差をもつボクセルよりも外側に、［最小の標準偏差＋α］以下の標準偏差をもつボクセルが存在しない場合は、最小の標準偏差をもつボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定され、図２３または図２７に示された部分表面Ｓｐ（ｋ）上に部分表面が特定される。
【０１４３】
以上のようにこの実施の形態４によれば、最小の標準偏差よりも少しだけ大きい標準偏差をもつボクセルが特定され、これにより部分表面Ｓｐ（ｋ）が推定されているため、照明条件の変動などによりノイズが入っても、一部を欠いた立体モデルが生成されることはない。
【０１４４】
［実施の形態５］
図２９は、本発明の実施の形態５による３次元モデリング装置における部分表面推定処理を示すフローチャートである。この方法は、上記実施の形態のようにボクセル空間の座標変換を行なうことなく視点からボクセル空間への透視性を考慮して標準偏差が最小のボクセルを特定して対象物体１０の部分表面を推定するもので、以下「光線追跡による最小ばらつき探索法」と呼ぶ。
【０１４５】
まずここでは、図３０に示されるような座標（ｐｘ，ｐｙ，ｐｚ）を定義する。すなわち、中心画像ｋの視点の座標を（ｏｘ，ｏｙ，ｏｚ）と定義する。そして、この視点（ｏｘ，ｏｙ，ｏｚ）からボクセル空間１１の中心に向かう方向にＺ′座標軸を設け、中心画像ｋ上にＸ′座標軸およびＹ′座標軸を設ける。視点（ｏｘ，ｏｙ，ｏｚ）から放射状に出る光線４０が中心画像ｋを貫通するｘ′ｙ′座標を（ｐｘ，ｐｙ）とし、中心画像ｋのｚ′座標をｐｚとする。
【０１４６】
図２９に示されるように、まずステップＳ５００で座標ｐｙに最小値ｐｙｍｉｎが設定され、続いてステップＳ５０１で座標ｐｘに最小値ｐｘｍｉｎが設定される。
【０１４７】
続いてステップＳ５０２で、単位光線ベクトルｖ＝（ｖｘ，ｖｙ，ｖｚ）が計算される。ただし、ここでは次の式（３）〜（５）が成立する。
【０１４８】
ｖｚ＝１ …（３）
ｖ×（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）＝０ …（４）
ｖ・（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）≧０ …（５）
ベクトルｖとベクトル（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）との外積が０ということは、これらのベクトルが同一方向に横たわっている（向きは同じか反対か不明）ことを意味する。また、ベクトルｖとベクトル（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）との内積が０以上ということは、これらのベクトルの向きも同じであることを意味する。
【０１４９】
続いてステップＳ５０３で、整数ｎに最小値ｎｍｉｎが設定され、標準偏差一時変数ｓｔｍｐに無限大∞が設定され、かつ座標一時変数ｚｔｍｐに無限大∞が設定される。
【０１５０】
続いてステップＳ５０４〜Ｓ５０８では、図３１に示されるように、直交座標系のボクセル空間１１において視点（ｏｘ，ｏｙ，ｏｚ）を通る各光線４０上の複数のボクセルの中から標準偏差が最小のボクセルが特定される。図３１では、複数のボクセル（○，●で表わされる）のうち、●で表わされるボクセルが特定されている。
【０１５１】
より具体的には、ステップＳ５０４で視点（ｏｘ，ｏｙ，ｏｚ）を基準としてｎ倍の単位光線ベクトルｖで表わされる座標（ｘ′ｙ′ｚ′）に位置するボクセルが選択される。
【０１５２】
続いてステップＳ５０４０でボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５０５で、標準偏差一時変数ｓｔｍｐの値がその選択されたボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）と比較される。標準偏差一時変数ｓｔｍｐの値が標準偏差ｓ（ｘ，ｙ，ｚ）よりも大きい場合は、ステップＳ５０６でその標準偏差ｓ（ｘ，ｙ，ｚ）が代わりに標準偏差一時変数ｓｔｍｐに設定され、そのときの自然数ｎが位置一時変数ｚｔｍｐに設定される。
【０１５３】
続いてステップＳ５０７およびＳ５０８により、すべての自然数ｎ（＝ｎｍｉｎ〜ｎｍａｘ）について上記ステップＳ５０４〜Ｓ５０６の処理が行なわれる。そのため、光線４０上に位置する複数のボクセルの中から標準偏差が最小のボクセルが特定される。
【０１５４】
続いてステップＳ５０９で、この特定されたボクセルのｚ′座標（位置一時変数ｚｔｍｐに設定されている）が表面変数ｓｕｒｆａｃｅ（ｐｘ，ｐｙ）に設定される。
【０１５５】
続いてステップＳ５１０およびＳ５１１により、すべての座標ｐｘについて上記ステップＳ５０２〜Ｓ５０９の処理が行なわれる。
【０１５６】
続いてステップＳ５１２およびＳ５１３により、すべての座標ｐｙ（＝ｐｙｍｉｎ〜ｐｙｍａｘ）について上記ステップＳ５０１〜Ｓ５１１の処理が行なわれる。
【０１５７】
図３１では、視点（ｏｘ，ｏｙ，ｏｚ）を通るある１つの光線４０上の複数のボクセルの中から最小の標準偏差をもつ●で表わされるボクセルが特定される。この特定されたボクセルのｐｚ座標は「３７」であるので、この「３７」が中心画像ｋ上の表面変数ｓｕｒｆａｃｅ（ｘ，ｙ）の対応する座標に設定される。
【０１５８】
上記のように各光線４０上ごとに最小の標準偏差をもつボクセルのｐｚ座標が表面変数ｓｕｒｆａｃｅ（ｘ，ｙ）に設定されることにより、図３２および図３３に示されるように対象物体１０の部分表面Ｓｐ（ｋ）が特定される。
【０１５９】
このようにして特定された複数の部分表面Ｓｐ（ｋ）は、上述した距離ペナルティ法または後に詳述する単純ボーティング法により１つに統合され、これにより対象物体１０の形状が決定される。図３２に示されるように、部分表面Ｓｐ（ｋ）よりも外側に位置する○で表わされるボクセルは対象物体の外部に位置し、部分表面Ｓｐ（ｋ）よりも内側に位置する●で表わされるボクセルは対象物体１０の内部に位置している。
【０１６０】
以上のようにこの実施の形態５によれば、視点から放射上に出る光線上の複数のボクセルのうち最小の標準偏差をもつボクセルが特定されているため、対象物体１０の形状を比較的正確に推定することができる。このような光線追跡による最小ばらつき探索法は上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができ、その結果、この処理に必要なメモリのワークエリアを削減するこができる。
【０１６１】
［実施の形態６］
図３４は、本発明の実施の形態６による３次元モデリング装置における部分表面統合処理を示すフローチャートである。この方法は、上述した距離ペナルティ法に代えて用いられ、上記部分表面推定処理により得られた複数の部分表面Ｓｐ（ｋ）を１つに統合して対象物体の形状を決定するものである。以下、この方法を「単純ボーティング法」と呼ぶ。
【０１６２】
図３５（ａ）〜（ｃ）は、各観察方向からの部分表面Ｓｐ（ｋ）が決定されたボクセル空間を示す。ここでは、推定された概略形状３１よりも内側でかつ部分表面Ｓｐ（ｋ）よりも内側のボクセルに得点「１」が与えられている。このような得点がボクセル空間に投票されると、図３５（ｄ）に示されるようなボクセル空間が得られる。ここで、対象物体１０は１６個の視点から観察されているため、各ボクセルの最高得点は「１６」である。この得点数が高いボクセルほど対象物体１０の真の存在領域である可能性が高いが、ここでは表面条件の変動などによるノイズを考慮して投票数が「１４」以上のボクセルを図３５（ｅ）に示されるように対象物体１０の真の存在領域４５として推定している。
【０１６３】
より具体的には、図３４に示されるように、ステップＳ６００で変数ｃｎｔ（ｘ，ｙ，ｚ）が初期化され、ステップＳ６００１で中心画像ｋとして物体画像Ａ０が設定される。続いてステップＳ６０１〜Ｓ６０３で、座標ｙ、ｘ、ｚにそれぞれ最小値ｙｍｉｎ、ｘｍｉｎ、ｚｍｉｎが設定される。
【０１６４】
続いてステップＳ６０４で、図１０に示されたステップＳ４０の結果に従って座標（ｘ，ｙ，ｚ）のボクセルが対象物体１０の概略存在領域内か否かが判定される。概略存在領域内、つまりｖ（ｘ，ｙ，ｚ）＝１の場合は、ステップＳ６０５で中心画像ｋに応じて座標変換が行なわれる。ここでは、座標（ｘ，ｙ，ｚ）が上述した座標（ｘ′，ｙ′，ｚ′）に変換される。
【０１６５】
続いてステップＳ６０６で、このボクセルが部分表面Ｓｐ（ｋ）よりも内側にあるか否か、つまりｓｕｒｆａｃｅ（ｘ′，ｙ′）≦ｚ′か否かが判定される。このボクセルが部分表面Ｓｐ（ｋ）よりも内側にある場合は、ステップＳ６０７で変数ｃｎｔ（ｘ，ｙ，ｚ）がインクリメントされる。
【０１６６】
続いてステップＳ６０８〜Ｓ６１３により、すべての座標ｘ（＝ｘｍｉｎ〜ｘｍａｘ）、すべての座標ｙ（＝ｙｍｉｎ〜ｙｍａｘ）、すべての座標ｚ（＝ｚｍｉｎ〜ｚｍａｘ）について上記ステップＳ６０４〜Ｓ６０７の処理が行なわれる。これにより図３５（ａ）〜（ｃ）に示されるようなボクセル空間が得られる。
【０１６７】
続いてステップＳ６１４およびＳ６１５により、すべての中心画像ｋ（＝０〜ｋｍａｘ）について上記ステップＳ６０１〜Ｓ６１３の処理が行なわれる。これにより、図３５（ｄ）に示されるようなボクセル空間が得られる。
【０１６８】
続いてステップＳ６１５０でボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ６１６で、各ボクセルの変数ｃｎｔ（ｘ，ｙ，ｚ）が所定の値（ここでは「１４」）と比較される。変数ｃｎｔ（ｘ，ｙ，ｚ）が所定の値よりも小さい場合は、ステップＳ６１７でそのボクセルは対象物体１０の存在領域外、つまり変数ｖ（ｘ，ｙ，ｚ）に０が設定される。
【０１６９】
続いてステップＳ６１８により、すべてのボクセルについて上記ステップＳ６１６およびＳ６１７の処理が行なわれ、これにより図３５（ｅ）に示されるようなボクセルが得られる。これにより対象物体１０の形状４５が決定される。
【０１７０】
以上のようにこの実施の形態６によれば、推定された複数の部分表面Ｓｐ（ｋ）を単純ボーティング法により統合しているため、上述した距離ペナルティ法よりも高速で対象物体１０の形状を決定することができる。
【０１７１】
以上、本発明の実施の形態を説明したが、本発明は上記実施の形態に限定されるものではない。たとえば、上記実施の形態は画像情報として濃度を用いているが、これに代えて色を用いてもよく、さらに色および濃度の組合せを用いてもよい。また、上記実施の形態は画像情報の標準偏差を用いているが、これに代えて画像情報の分散を用いてもよく、要するに、画像情報のばらつきを評価する値を用いればよい。
【０１７２】
【発明の効果】
以上のように本発明によれば、対象物体を複数の視点から撮影して得られた複数の物体画像のうち少なくとも２つの物体画像間における画像情報のばらつきに基づいて対象物体の形状を推定しているため、対象物体の凹所を表現した立体モデルを生成することができる。
【０１７３】
また、予め対象物体の輪郭に基づいてボクセル空間における対象物体の存在領域を概略的に推定しているため、対象物体の立体モデルを迅速に生成することができる。
【０１７４】
また、ばらつきが小さいボクセルを滑らかに結合することにより対象物体の形状を推定しているため、照明条件などが変動しても正確な立体モデルを安定して生成することができる。
【０１７５】
また、３次元仮想空間の座標空間を透視座標空間に変換し、その透視座標空間において１つの視点を通る直線上の複数の領域のうち、ばらつきが最小またはその最小値よりも所定の値だけ大きい領域を結合しているため、対象物体の形状を高速で推定することができる。
【０１７６】
また、３次元仮想空間において１つの視点を通る直線上の複数の領域のうち、ばらつきが最小またはその最小値よりも所定の値だけ大きい領域を結合しているため、対象物体の形状を高速で推定することができる。
【０１７７】
また、３次元仮想空間の領域のうち存在領域の内側にありかつ部分表面の各々の外側にある領域にその領域から複数の部分表面までの距離を累積し、その累積された距離の総和が所定の値よりも大きい領域を対象物体の外側と認定しているため、対象物体の形状を比較的正確に決定することができる。
【０１７８】
また、３次元仮想空間の領域のうち部分表面の各々の内側にある領域に所定の値を投票し、その投票された値の総和が所定の値よりも大きい領域を対象物体の内側と認定しているため、対象物体の形状を高速で決定することができる。
【０１７９】
また、複数の画素を含む所定ブロックごとに画像情報を累積しているため、対象物体の立体モデルをより精度良く生成することができる。
【０１８０】
また、各物体画像内に複数の画素を含むブロックを設定し、そのブロック間において互いに対応する画素の値の差をそれぞれ算出し、その算出された差の絶対値の総和を画像情報のばらつきとして算出しているため、対象物体の形状を高速で推定することができる。
【０１８１】
また、各物体画像内に複数の画素を含むブロックを設定し、そのブロック間において互いに対応する画素の値の差をそれぞれ算出し、その算出された差の自乗の総和を画像情報のばらつきとして算出しているため、対象物体の形状を高速で推定することができる。
【図面の簡単な説明】
【図１】本発明の実施の形態１による３次元モデリング装置の全体構成を示す図である。
【図２】図１に示されたコンピュータの主要構成を示すブロック図である。
【図３】図１に示されたＣＤ−ＲＯＭに記録された３次元モデリングプログラムを示すフローチャートである。
【図４】（ａ）は図３中の対象物体および背景の撮影ステップを説明するための図であり、（ｂ）は図３中のシルエット画像作成ステップを説明するための図であり、（ｃ）は図３中の形状推定処理ステップを説明するための図であり、（ｄ）は図３中のポリゴン作成ステップを説明するための図であり、（ｅ）は図３中のテクスチャマッピングステップを説明するための図である。
【図５】図４（ａ）において対象物体を撮影する視点を説明するための図である。
【図６】図４（ｃ）の形状推定処理で用いられるボクセル空間を示す概念図である。
【図７】図３に示された形状推定処理の詳細を示すフローチャートである。
【図８】図７に示された概略形状推定処理を説明するための図である。
【図９】図７に示された色ボーティング処理を説明するための図である。
【図１０】図７に示された概略形状推定処理、色ボーティング処理、および部分表面推定処理の詳細を示すフローチャートである。
【図１１】図７に示された部分表面統合処理の詳細を示すフローチャートである。
【図１２】図１０中のステップＳ４５で用いられる投影点を説明するための図である。
【図１３】図１０中のステップＳ５２で用いられるダイナミックプログラミング法を示すフローチャートである。
【図１４】図１０に示されたフローチャートによって計算されたＹ＝ｙ０平面における標準偏差を示す鳥瞰図である。
【図１５】図１０に示されたフローチャートによって計算されたＹ＝ｙ０平面における標準偏差の数値を示す図である。
【図１６】図１３に示されたフローチャートによって決定された部分表面を示す図である。
【図１７】図１０中のステップＳ４０で得られた対象物体の存在領域の初期形状を示す図である。
【図１８】本発明の実施の形態２による３次元モデリング装置における形状推定処理を示すフローチャートである。
【図１９】図１７中のステップＳ７０およびＳ７２で用いられる投影点を説明するための図である。
【図２０】図１７中のステップＳ７４においてブロック領域間で濃度が異なる画素の数を計算する方法を説明するための図である。
【図２１】本発明の実施の形態３による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図２２】図２１中のステップＳ３００における直交座標系のボクセル空間から直交座標系のボクセル空間への座標変換を説明するための図である。
【図２３】図２１に示されたフローチャートに従って直交座標系のボクセル空間から直交座標系のボクセル空間への座標変換を行なった場合における部分表面の特定およびその統合方法を説明するための図である。
【図２４】標準偏差が最小のボクセルが必ずしも対称物体の部分表面上に位置するとは限らない理由を明らかにするために、同一のｘｙ座標をもつ複数のボクセルから中心画像への投影点が異なる場合における色ボーティング処理を説明するための図である。
【図２５】図２１中のステップＳ３００における座標変換で使用可能な透視座標系のボクセル空間を示す斜視図である。
【図２６】（ａ）は図２５に示された透視座標系のボクセル空間をＹ方向から見た上面図であり、（ｂ）は図６に示された直交座標系のボクセル空間をＹ方向から見た上面図である。
【図２７】図２１中のステップＳ３００で図２５に示されるような透視座標系のボクセル空間への座標変換が行なわれた場合における部分表面の特定およびその統合方法を説明するための図である。
【図２８】本発明の実施の形態４による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図２９】本発明の実施の形態５による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図３０】図２９に示されたフローチャートで用いられる座標系を説明するための図である。
【図３１】図２９に示された光線追跡による最小ばらつき探索法を説明するための図である。
【図３２】図２９に示された光線追跡による最小ばらつき探索法により特定された部分表面の統合方法を説明するための図である。
【図３３】図２９に示された光線追跡による最小ばらつき探索法により１つの部分表面が特定されたボクセルを示す図である。
【図３４】本発明の実施の形態６による立体モデリング装置における部分表面統合処理の詳細を示すフローチャートである。
【図３５】（ａ）〜（ｅ）は、図３４に示された単純ボーティング法を説明するための図である。
【符号の説明】
２ デジタルスチルカメラ
３ コンピュータ
４ ＣＤ−ＲＯＭ
５ 中央処理装置
７ ランダムアクセスメモリ
１０ 対象物体
１１ ボクセル空間
１２ ボクセル
１８ 立体モデル
Ａ０〜Ａ１５ 物体画像
Ｂ０〜Ｂ１５ シルエット画像
Ｂｋ，Ｂｉ ブロック領域
【発明の属する技術分野】
本発明は、３次元モデリング装置およびその方法、３次元モデリングプログラムを記録した媒体、ならびに３次元形状推定プログラムを記録した媒体に関し、特に、対象物体の立体モデルを生成する３次元モデリング装置およびその方法ならびに３次元モデリングプログラムを記録した媒体に関する。
【０００２】
【従来の技術】
一般に、３次元グラフィックス画像を得るためには、３次元対象物体の形状およびテクスチャ（色彩および模様）を表現した立体モデルが作成される。
【０００３】
立体モデルを作成する従来の３次元モデリング方法としては、たとえばシルエット法やステレオ法がある。シルエット法は、対象物体を複数の視点から単眼カメラで撮影し、それにより得られた複数のシルエット画像に基づいて対象物体の立体モデルを生成する。
【０００４】
【発明が解決しようとする課題】
しかしながら、シルエット法は撮影された対象物体の輪郭に基づいて立体モデルを生成するため、対象物体の凹所を表現することができないという問題がある。
【０００５】
また、「財満久義他，『連続カラー画像を用いた双方向投票による物体形状と色の推定』，情報処理学会第５３回（平成８年後期）全国大会，２−３０１および３０２頁」には、色情報を投票して対象物体の形状を推定する方法が開示されている。この方法は、表面ボクセルにおいては同じ色が投票されるはずであるという仮定に基づいたものである。しかしながら、この方法は、すべての画像情報を一旦用いて色の投票を行なっているため、表面ボクセルに対し、裏方向の画像情報までもが投票されてしまう。したがって、通常、表面ボクセルにおいても投票された色が異なってしまう。この問題を解決するため、ボクセルに投票された候補色を絞り込む必要が生じている。この処理のため、この方法は、投票された色をすべて保持するための多大なメモリが必要となる。これに加え、候補色の中から色の同一性を抽出して表面ボクセルを決定するために、同じ部分が常に同じ色で撮影される必要がある。あるいは、同じ色であるか否かを判断するためのしきい値設定が必要になる。このしきい値設定は一般に困難な処理となる。
【０００６】
本発明は上記のような問題点を解決するためになされたもので、対象物体の凹所を表現可能な簡易かつ安価な３次元モデリング装置およびその方法ならびに３次元モデリングプログラムを記録した媒体を提供することを目的とする。
【０００７】
【課題を解決するための手段】
請求項１に係る３次元モデリング装置の発明は、対象物体の立体モデルを生成する３次元モデリング装置であって、前記対象物体を複数の視点から撮影し、複数の物体画像を得る物体画像獲得手段と、前記複数の物体画像のうち少なくとも２つの物体画像間において互いに対応する画像情報についての濃度情報及び／又は色情報の標準偏差に基づいて前記対象物体の形状を推定する形状推定手段と、前記推定された形状に基づいて前記立体モデルを生成する生成手段とを備え、前記形状推定手段は、複数の領域から構成される３次元仮想空間の領域に前記少なくとも２つの物体画像に対応する画像情報についての濃度情報及び／又は色情報を累積する画像情報累積手段と、前記領域に累積された画像情報についての濃度情報及び／又は色情報の標準偏差を算出する算出手段とを含む。
【０００８】
請求項２に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記複数の物体画像中の前記対象物体の輪郭に基づいて前記３次元仮想空間における前記対象物体の存在領域を概略的に推定する存在領域推定手段をさらに備え、前記画像情報累積手段は、前記３次元仮想空間の領域のうち前記存在領域内の領域に前記画像情報についての濃度情報及び／又は色情報を累積する。
【０００９】
請求項３に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記形状推定手段はさらに、前記対象物体の形状を推定するために、前記３次元仮想空間の領域のうち前記標準偏差の小さい領域を滑らかに結合する結合手段を含む。
【００１０】
請求項４に係る３次元モデリング装置の発明は、請求項３に記載の発明において、前記結合手段は、前記標準偏差の小さい領域を滑らかに結合するために、前記標準偏差の小ささと前記結合されるべき領域の滑らかさとを変数とする評価関数を最小化する。
【００１１】
請求項５に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記３次元仮想空間の座標空間を前記視点の１つから放射状に広がる座標軸を有する透視座標空間に変換する座標変換手段と、前記対象物体の形状を推定するために、前記透視座標空間において前記１つの視点を通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段とを含む。
【００１２】
請求項６に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記対象物体の形状を推定するために、前記３次元仮想空間において前記視点の１つを通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段を含む。
【００１３】
請求項７に係る３次元モデリング装置の発明は、請求項１から請求項２のいずれかに記載の発明において、前記形状推定手段はさらに、前記複数の物体画像間において互いに対応する画像情報の標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含む。
【００１４】
請求項８に係る３次元モデリング装置の発明は、請求項２に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、前記累積された距離の総和が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む。
【００１５】
請求項９に係る３次元モデリング装置の発明は、請求項２に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、前記累積回数を計算する累積回数計算手段と、前記累積された距離の総和を前記計算された累積回数で除した距離平均値が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む。
【００１６】
請求項１０に係る３次元モデリング装置の発明は、請求項１に記載の発明において、前記形状推定手段はさらに、前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、前記統合手段は、前記３次元仮想空間の領域のうち前記部分表面の各々の内側にある領域に所定の第１の値を投票する投票手段と、前記投票された第１の値の総和が所定の第２の値よりも大きい領域を前記対象物体の内側と認定する認定手段とを含む。
【００１７】
請求項１１に係る３次元モデリング装置の発明は、請求項１から請求項１０のいずれかに記載の発明において、前記画像情報累積手段はさらに、前記物体画像の複数の画素を含む所定ブロックごとに前記画像情報についての濃度情報及び／又は色情報を累積する。
【００１８】
請求項１２に係る３次元モデリング装置の発明は、請求項１から請求項１１のいずれかに記載の発明において、前記領域はボクセルであり、前記３次元仮想空間はボクセル空間である。
【００５８】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して詳しく説明する。なお、図中同一または相当部分には同一符号を付してその説明は繰返さない。
【００５９】
［実施の形態１］
図１は、本発明の実施の形態１による３次元モデリング装置の全体構成を示す図である。図１を参照して、この実施の形態１による３次元モデリング装置は、回転テーブル１と、デジタルスチルカメラ２と、コンピュータ３と、ＣＤ−ＲＯＭ（コンパクトディスク−リードオンリメモリ）４とを備える。回転テーブル１は、その上に載置された対象物体（図示せず）を回転させるためのものである。ここでは回転テーブル１を用いたが、これに代えて対象物体の向きを変えることができるロボットアームなどを用いてもよい。デジタルスチルカメラ２は、立体モデルを得ようとする対象物体を複数の視点から撮影し、複数の物体画像を得るためのものである。ここでは静止画を得ることができるデジタルスチルカメラ２を用いたが、これに代えて動画を得ることができるビデオカメラなどを用いてもよい。ＣＤ−ＲＯＭ４には、コンピュータ３に、デジタルスチルカメラ２で得られた物体画像に基づいて対象物体の立体モデルを生成させるためのプログラムが記録されている。ここではプログラムを記録する媒体としてＣＤ−ＲＯＭ４を用いたが、これに代えてフロッピィディスク、ＣＤ−Ｒ、ＤＶＤ、メモリカードなどを用いてもよい。
【００６０】
なお、この３次元モデリング装置は、カメラ２を固定し、対象物体を回転させることによって対象物体を複数の視点から撮影するようにしているが、対象物体を固定し、カメラ２を移動することによって対象物体を複数の視点から撮影するようにしてもよい。ただし、カメラ２を移動する場合はカメラ２の撮影位置を求める必要がある。この３次元モデリング装置は、カメラ２を固定し、対象物体を回転させているため、各視点から対象物体までの距離は互いに等しいが、この距離は必ずしも等しくなくてもよい。距離が等しくない場合は、以下に詳述する方法で各視点から対象物体までの距離を考慮して種々の計算を行なう必要がある。
【００６１】
図２は、図１に示されたコンピュータ３の主要構成を示すブロック図である。図２を参照して、このコンピュータ３は、中央処理装置（ＣＰＵ）５と、リードオンリメモリ（ＲＯＭ）６と、ランダムアクセスメモリ（ＲＡＭ）７と、入出力インタフェース（Ｉ／Ｆ）８と、ＣＰＵ５、ＲＯＭ６、ＲＡＭ７およびＩ／Ｆ８を相互に結合するためのデータ／アドレスバス９とを備える。デジタルスチルカメラ２、およびＣＤ−ＲＯＭ４のためのドライブユニットはＩ／Ｆ８に接続され、これによりデジタルスチルカメラ２で得られた物体画像およびＣＤ−ＲＯＭ４に記録されたプログラムはＩ／Ｆ１０を介してＲＡＭ７に格納される。ＣＰＵ５は、ＲＯＭ６およびＲＡＭ７に格納されたプログラムに従って演算処理を行なう。
【００６２】
図３は、図１に示されたＣＤ−ＲＯＭ４に記録された３次元モデリングプログラムを示すフローチャートである。図４（ａ）は、図３中の対象物体および背景の撮影ステップＳ２を説明するための図である。図４（ｂ）は、図３中のシルエット画像の作成ステップＳ３を説明するための図である。図４（ｃ）は、図３中の形状推定処理ステップＳ４を説明するための図である。図４（ｄ）は、図３中のポリゴン作成ステップＳ５を説明するための図である。図４（ｅ）は、図３中のテクスチャマッピングステップＳ６を説明するための図である。
【００６３】
図１〜図４を参照して、まずステップＳ１では、回転テーブル１とカメラ２との間の相対的な位置関係を決定するためのキャリブレーションを行なう。
【００６４】
続いてステップＳ２では、対象物体１０およびその背景をカメラ２で撮影する。より具体的には、まず、対象物体１０を回転テーブル１の上に置かずに背景だけを１つの視点から撮影する。これにより１枚の背景画像が得られる。次に、対象物体１０を回転テーブル１の上に置き、回転テーブル１を回転させることによって対象物体１０を複数の視点から撮影する。たとえば回転テーブル１を２２．５°ずつ回転させると、図５に示されるように対象物体１０を１６の視点（０〜１５）から撮影することになるため、１６枚の物体画像Ａ０〜Ａ１５が得られる。回転テーブル１およびカメラ２はコンピュータ３内のＩ／Ｆ８に接続され、ＣＰＵ５によって制御される。また、得られた１枚の背景画像および１６枚の物体画像Ａ０〜Ａ１５はＩ／Ｆ８を介してＲＡＭ７に格納される。なお、ここでは対象物体１０を１６の視点から撮影しているが、視点の数は２以上であればよく、特に限定されるものではない。
【００６５】
続いてステップＳ３では、公知の画像間差分処理法に従ってシルエット画像が作成される。すなわち、物体画像Ａ０〜Ａ１５と背景画像との間で差分処理が行なわれ、図４（ｂ）に示されるように１６枚のシルエット画像Ｂ０〜Ｂ１５が作成される。差分処理においては、物体画像の色信号レベルと背景画像の色信号レベルとの差が各画素ごとに算出される。なお、ここでは画像間差分処理法を用いているが、これに代えて公知のステレオ法を用いてもよい。ステレオ法は、両眼カメラで対象物体１０を撮影し、その画像中の奥行き情報に基づいて対象物体１０を背景と区別するものである。
【００６６】
続いてステップＳ４では、形状推定処理が行なわれる。この形状推定処理のために、図６に示されるような直交座標系のボクセル空間１１が予め用意され、ＲＡＭ７に格納されている。ここでは直交座標系のボクセル空間１１を用いているが、これに代えて円筒座標系のボクセル空間を用いてもよい。ボクセル空間１１は、複数のボクセル１２から構成される。
【００６７】
図７は、図３に示された形状推定処理ステップＳ４の詳細を示すフローチャートである。まず、ステップＳ４０の概略形状推定処理では、複数のシルエット画像Ｂ１〜Ｂ１５に基づいてボクセル空間１１における対象物体１０の存在領域が概略的に推定される。次に、ステップＳ４１の詳細形状推定処理では、物体画像Ａ０〜Ａ１５間における濃度情報（輝度信号レベル）の標準偏差に基づいて対象物体１０の形状が詳細に推定される。
【００６８】
対象物体１０の存在領域を概略的に推定するステップＳ４０では公知のボーティング処理が行なわれる。すなわち、図８に示されるように、各シルエット画像Ｂ１，Ｂ２が各視点からボクセル空間１１に投影される。このときシルエット画像Ｂ１，Ｂ２中の物体部分１３が貫通するボクセル１２に対象物体１０が存在する可能性として１点が投票される。すべてのシルエット画像Ｂ１〜Ｂ１５がボクセル空間１１に投影されたとき、投票数が所定数よりも多いボクセル１２が対象物体１０の存在領域と推定される。ただし、図５に示されるような対象物体１０の凹所１４はシルエット画像Ｂ０〜Ｂ１５に現われないため、このボーティング処理によって推定された存在領域は対象物体１０の形状を概略的に表現するものにすぎない。すなわち、対象物体１０の凹所の真の表面１４は現れず、仮の表面１５が現れる。
【００６９】
通常、対象物体１０の真の表面１４はいずれの視点からも同じ色および濃度で観察されるが、仮の表面１５は視点によって異なった色または濃度で観察される。したがって、同じ部分でも視点によって異なった色または濃度で観察される部分は対象物体１０の凹所である可能性が高い。
【００７０】
そこで、対象物体１０の形状を詳細に推定するステップＳ４１では、対象物体１０の形状を概略的に推定した後、濃度情報の標準偏差に基づいて対象物体１０の形状、特に凹所の形状が詳細に推定される。より具体的には、図９に示されるように、１６枚の物体画像Ａ０〜Ａ１５のうち、中心画像として物体画像Ａ０、および周辺画像として物体画像Ａ０に隣接する物体画像Ａ１，Ａ１５を用いて形状推定を行なう場合、物体画像Ａ０，Ａ１，Ａ１５ともに、その画素の濃度情報がボクセル空間１１に投影され、各画素からの投影線１６が貫通するすべてのボクセル１２にその画素の濃度情報が累積される（色ボーティング処理ステップＳ４１０）。換言すれば、すべてのボクセル１２について、物体画像Ａ０，Ａ１，Ａ１５への投影点が求められ、その投影点の濃度情報が対応するボクセル１２に累積される。そして、各ボクセル１２に累積された濃度情報の標準偏差が算出される。標準偏差が大きいボクセル１２は、それに対応する部分が視点によって異なった濃度で観察されることを意味するので、そのようなボクセル１２は対象物体１０の凹所である可能性が高い。そのため、ボクセル空間１１内のボクセル１２のうち標準偏差が小さいボクセルを滑らかに結合することにより、図５に示された視点（０）から観察可能な対象物体１０の部分表面が推定される（部分表面推定処理ステップＳ５２）。
【００７１】
上述した色ボーティング処理ステップＳ４１０および部分表面推定処理ステップＳ５２は、物体画像Ａ０〜Ａ１５の各々を中心画像として行なわれる。したがって、推定された１６枚の部分表面が統合され、これにより対象物体１０の立体形状が得られる（部分表面統合処理ステップＳ４１２）。
【００７２】
上記ステップＳ４に続くステップＳ５では、ステップＳ４で得られた対象物体１０の立体形状に基づいて三角パッチからなるポリゴン１７が作成される。したがって、対象物体１０の立体形状は多数のポリゴン１７で表現される。ポリゴン１７で表現された立体形状はＲＡＭ７に格納される。
【００７３】
最後にステップＳ６では、ステップＳ５で作成された各ポリゴン１７に対応するテクスチャがマッピングされる。これにより、対象物体１２をモデル化した立体モデル１８が作成される。テクスチャは物体画像Ａ０〜Ａ１５から得られる色情報で、ＲＡＭ９に格納される。
【００７４】
なお、ここでは回転テーブル１を用いて対象物体１０を複数の視点から撮影して複数の物体画像Ａ０〜Ａ１５を得ているが、対象物体１０の立体形状を推定するために必要な物体画像を得ることさえできれば、上記ステップＳ１〜Ｓ３は適宜他の公知の手法に置換えてもよい。また、物体画像および推定された立体形状に基づいて立体モデル１８を生成することさえできれば、上記ステップＳ５およびＳ６は適宜他の公知の手法に置換えてもよい。また、対象物体１０の輪郭に基づいてボクセル空間１１における対象物体１０の存在領域を概略的に推定することさえできれば、上記ステップＳ４中の前半のステップは適宜他の公知の手法に置換えてもよい。
【００７５】
次に、上記のように構成された３次元モデリング装置を用いて対象物体１０の立体モデル１８を生成する方法について説明する。
【００７６】
まずステップＳ１で上述したようにキャリブレーションを行なった後、ステップＳ２で１つの視点から背景を撮影し、続いて複数の視点から対象物体１０を撮影する。これにより１枚の背景画像および複数の物体画像Ａ０〜Ａ１５が得られる。続いてステップＳ３で、物体画像Ａ０〜Ａ１５と背景画像との間で差分処理が行なわれ、これによりシルエット画像Ｂ０〜Ｂ１５が得られる。
【００７７】
続いて、ステップＳ４０でシルエット画像Ｂ０〜Ｂ１５に基づいてボーティング処理などの概略形状推定処理が行なわれた後、ステップＳ４１で物体画像Ａ０〜Ａ１５の濃度情報に基づいて詳細形状推定処理が行なわれる。図１０は、図７中の概略形状推定処理ステップＳ４０、色ボーティング処理ステップＳ４１０および部分表面推定処理ステップＳ５２の詳細を示すフローチャートである。なお、図１０中のステップＳ４２１〜ステップＳ５１が図７中の色ボーティング処理ステップＳ４１０に対応する。図１１は、図７中の部分表面統合処理ステップＳ４１２の詳細を示すフローチャートである。
【００７８】
まずステップＳ４０で、シルエット画像Ｂ０〜Ｂ１５に基づいて上述したボーティング処理などの概略形状推定処理が行なわれ、対象物体１０の初期立体形状が概略的に抽出される。より具体的には、投票数が所定数よりも多いボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域と推定され、そのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）が「１」に設定される。他方、投票数が所定数よりも少ないボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域と推定されず、そのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）が「０」に設定される。したがって、ｖ（ｘ，ｙ，ｚ）＝１のボクセルが対象物体１０の概略的な存在領域を形成する。
【００７９】
続いてステップＳ４２で、図５に示された視点（０）から対象物体１０を撮影して得られた物体画像が中心画像ｋとして選択される。
【００８０】
続いてステップＳ４２１で、変数（ｃ０，ｃ１，ｃ２）が初期化される。
続いてステップＳ４３で、中心画像ｋ（＝０）から図５上で時計方向にｎ個離れた視点（１６−ｎ（「１６」は図５上の「０」に相当する））から対象物体１０を撮影して得られた物体画像が周辺画像ｉとして選択される。たとえばｎ＝２とすると、図５に示された視点（１４）から対象物体１０を撮影して得られた物体画像（周辺画像）Ａ１４が選択される。
【００８１】
続いてステップＳ４４で、存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。
【００８２】
ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ４５で、図１２に示されるように座標（ｘ，ｙ，ｚ）に位置するボクセル１２の画像ｉ（ここではｉ＝１４）への投影点（Ｘ，Ｙ）が計算される。
【００８３】
続いてステップＳ４６で、後述するステップＳ５０で標準偏差ｓ（ｘ，ｙ，ｚ）を計算するために必要な値が計算される。より具体的には、変数ｃ０（ｘ，ｙ，ｚ）をインクリメントすることによってボクセル１２に累積される濃度情報の数がカウントされる。また、変数ｃ１（ｘ，ｙ，ｚ）に投影点（Ｘ，Ｙ）に位置する画素の濃度情報ｃ（Ｘ，Ｙ）が加算される。また、変数ｃ２（ｘ，ｙ，ｚ）にその濃度情報の自乗ｃ（Ｘ，Ｙ）2 が加算される。なお、ステップＳ４４でｖ（ｘ，ｙ，ｚ）＝０の場合、ステップＳ４５およびＳ４６の処理は行なわれない。
【００８４】
続いてステップＳ４７で、上記ステップＳ４４〜Ｓ４６の処理がすべてのボクセルについて行なわれたか否かが判定される。未だすべてのボクセルについて上記ステップＳ４４〜Ｓ４６の処理が行なわれていない場合、再び対象物体１０の存在領域内の各ボクセルの変数ｃ１（ｘ，ｙ，ｚ）およびｃ２（ｘ，ｙ，ｚ）に画像ｉ（＝１４）中の対応する投影点（Ｘ，Ｙ）の濃度情報ｃ（Ｘ，Ｙ）およびその自乗ｃ（Ｘ，Ｙ）2 がそれぞれ加算される。
【００８５】
すべてのボクセルについて上記ステップＳ４４〜Ｓ４６の処理が行なわれた場合、ステップＳ４８で、上記ステップＳ４４〜Ｓ４７の処理が中心画像ｋ（＝０）およびすべての周辺画像ｉ（＝１４，１５，１，２）について行なわれたか否かが判定される。最初は周辺画像ｉ（＝１４）について上記ステップＳ４４〜Ｓ４７の処理が行なわれたにすぎないが、ステップＳ４９で変数ｉがインクリメントされ、対象の周辺画像ｉがｉ＝１４からｉ＝１５に変更される。したがって、上記と同様に、対象物体１０の存在領域内のボクセル１２に周辺画像ｉ（＝１５）中の対応する投影点（Ｘ，Ｙ）の濃度情報ｃ（Ｘ，Ｙ）およびｃ（Ｘ，Ｙ）2 が累積される。
【００８６】
変数ｉの値が変数（ｋ＋ｎ）の値（ここでは２）を超えたとき、中心画像ｋ（＝０）および周辺画像ｉ（＝１４，１５，１，２）の濃度情報がすべて対象物体１０の存在領域内のボクセル１２に累積されたことになる。
【００８７】
ステップＳ４８で変数ｉの値が変数（ｋ＋ｎ）の値よりも大きいと判定されると、ステップＳ４９０で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５０で、ステップＳ４６で計算されたｃ０，ｃ１およびｃ２に基づいて各ボクセル１２に累積された濃度情報の標準偏差ｓ（ｘ，ｙ，ｚ）が計算される。
【００８８】
続いてステップＳ５１で、すべてのボクセルについて標準偏差ｓ（ｘ，ｙ，ｚ）が計算されたか否かが判定される。
【００８９】
すべてのボクセルについて標準偏差ｓ（ｘ，ｙ，ｚ）が計算された場合、ステップＳ５２でボクセル空間１１内において各視点（０〜１５）から観察した対象物体１０の部分表面Ｓｐ（ｋ）が推定される。上述したように、対象物体１０の凹所１４内に位置するボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）は一般に大きくなり、対象物体１０の表面に位置するボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）は小さくなる。したがって、同一のｘｙ座標をもつボクセルのうち標準偏差ｓ（ｘ，ｙ，ｚ）が最小のボクセルを結合することによって対象物体１０の表面を推定することができる。ただし、照明条件が変動したり、凹所１４の表面で同じ色が連続する場合、ノイズの影響を受け、対象物体１０の表面に位置しないボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が最小になることがある。したがって、単純に標準偏差ｓ（ｘ，ｙ，ｚ）が最小のボクセルを結合した場合、必ずしも正確な形状を推定できるとは限らない。
【００９０】
そこで、標準偏差ｓ（ｘ，ｙ，ｚ）だけでなく推定される形状の連続性も考慮して対象物体１０の形状を推定するのが望ましい。より具体的には、標準偏差ｓ（ｘ，ｙ，ｚ）および隣接するボクセル間の距離を変数とするエネルギ関数（評価関数）を規定し、そのエネルギ関数の値が最小となるボクセルを対象物体１０の表面に位置すると推定し、これらのボクセルによって構成される面を部分表面Ｓｐ（ｋ）とする。以下、この方法を「最小エネルギ法」と呼ぶ。なお、このエネルギ関数は、標準偏差ｓ（ｘ，ｙ，ｚ）、および隣接するボクセル空間の傾きと、これに対応する概略形状での表面の傾きとの差を変数とするものであってもよい。このような部分表面の推定は、たとえばダイナミックプログラミング（ＤＰ）法、緩和法などによって行なわれる。ダイナミックプログラミング法については後に図１３を用いて詳述する。
【００９１】
上記のように、中心画像ｋ（＝０）および周辺画像ｉ（＝１４，１５，１，２）に基づいて部分表面Ｓｐ（ｋ）が対象物体１０の表面候補として推定されるが、他の視点（１〜１５）からの画像を中心画像としたときも同様に部分表面Ｓｐ（ｋ）が対象物体１０の表面候補としてそれぞれ推定される（Ｓ５４１，Ｓ５４２）。したがって、計算された複数の部分表面Ｓｐ（ｋ）を統合して対象物体１０の真の表面を決定する必要がある。
【００９２】
最初にステップＳ４０で概略的に推定された対象物体１０の存在領域内に位置しかつステップＳ５２で算出された部分表面Ｓｐ（ｋ）よりも内側に位置するボクセルは対象物体１０の真の存在領域である可能性が高い。これに対し、ステップＳ４０で概略的に推定された対象物体１０の存在領域内に位置しかつステップＳ５２で計算された部分表面Ｓｐ（ｋ）よりも外側に位置するボクセルは対象物体１０の真の存在領域でない可能性がある。その可能性はステップＳ５２で推定された部分表面Ｓｐ（ｋ）から遠く離れたボクセルほど高い。そこで、部分表面Ｓｐ（ｋ）よりも外側に位置する各ボクセルに部分表面Ｓｐ（ｋ）までの距離をペナルティとして与える。そして、このペナルティの平均値が所定の値よりも大きいボクセルを対象物体１０の存在領域以外の領域とする。
【００９３】
より具体的には、図１１に示されるように、まずステップＳ５４３で変数（ｐ０，ｐ１）が初期化され、続いてステップＳ５４４で図５に示された視点（０）から対象物体１０を撮影して得られた物体画像Ａ０が中心画像ｋとして選択される。続いて、ステップＳ５４５で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５３で、ボクセル（ｘ，ｙ，ｚ）がステップＳ５２で計算された部分表面Ｓｐ（ｋ）よりも外側に位置するか否かが判定される。
【００９４】
外側に位置する場合、変数ｐ０（ｘ，ｙ，ｚ）がインクリメントされる（Ｓ５４）。これにより、そのボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの数がカウントされる。また、変数ｐ１（ｘ，ｙ，ｚ）に距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））が加算される（Ｓ５４）。距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））はボクセル（ｘ，ｙ，ｚ）からステップＳ５２で推定された対象物体１０の部分表面Ｓｐ（ｋ）までの距離である。これにより、部分表面Ｓｐ（ｋ）から遠いボクセルほど大きなペナルティが与えられる。
【００９５】
ステップＳ５５で、すべてのボクセルについて上記ステップＳ５３〜Ｓ５４の処理が行なわれたか否かが判定され、行なわれた場合は、ステップＳ５６で図５に示されたすべての視点（０〜１５）から撮影された画像を中心画像として上記ステップＳ４３〜Ｓ５５の処理が行なわれたか否か、つまり変数ｋがその最大値ｋｍａｘ（ここでは１５）に達したか否かが判定される。変数ｋがその最大値ｋｍａｘよりも小さい場合、ステップＳ５７で変数ｋがインクリメントされ、再びステップＳ５３〜Ｓ５５の処理が行なわれる。
【００９６】
他方、変数ｋがその最大値ｋｍａｘよりも小さくない場合、ステップＳ５７０で存否変数ｖ（ｘ，ｙ，ｚ）が「０」か否か、つまりそのボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝０ではない場合、つまりｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５８〜Ｓ６０でペナルティの高いボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域ではないとされる。より具体的には、ステップＳ５８で各ボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きいか否かが判定される。ここで、ｐ０（ｘ，ｙ，ｚ）はボクセル（ｘ，ｙ，ｚ）に与えられたペナルティの数である。
【００９７】
ペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きい場合、ステップＳ５９でそのボクセルの存否変数ｖ（ｘ，ｙ，ｚ）は０に設定される。ペナルティの平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）が所定の値よりも大きくない場合、ステップＳ５９の処理は行なわれない。
【００９８】
続いてステップＳ６０で、すべてのボクセルについてステップＳ５８およびＳ５９の処理が行なわれたか否かが判定される。したがって、所定の値よりも大きい平均値ｐ１（ｘ，ｙ，ｚ）／ｐ０（ｘ，ｙ，ｚ）を有するボクセル（ｘ，ｙ，ｚ）は対象物体１０の存在領域から除外される。これにより、複数の部分表面Ｓｐ（ｋ）が統合され、対象物体１０の凹所が表現され得る。以下、図１１に示された部分表面Ｓｐ（ｋ）を統合する方法を「距離ペナルティ法」と呼ぶ。
【００９９】
上記のようにして得られた物体画像Ａ０〜Ａ１５および対象物体１０の形状に基づいて、図３に示されたステップＳ５でポリゴン１７が作成され、ステップＳ６でテクスチャがマッピングされる。これにより、対象物体１０の立体モデル１８が作成される。
【０１００】
図１３は、図１０に示されたステップＳ５２において部分表面Ｓｐ（ｋ）をダイナミックプログラミングによって計算する方法を示すフローチャートである。なお、図示されていないが、図１３の処理を行なう前に、ボクセル空間を中心画像ｋに応じて、後述の図２２に示された座標変換を行なっておくのが望ましい。図１４は、図１０および図１１に示されたフローチャートに従って計算されたＹ＝ｙ０平面における標準偏差ｓを示す鳥瞰図である。図１５は、ボクセル空間１１のＹ＝ｙ０平面内のボクセルの標準偏差ｓを示す図である。図１３に示されたダイナミックプログラミング法は、図１４および図１５上で、滑らかさを保ちつつ可能な限り標準偏差ｓの小さい箇所を結ぶものである。
【０１０１】
図１３を参照して、ステップＳ５２０で座標ｙにその最小値ｙｍｉｎが設定され、続いてステップＳ５２１でｘにその最小値ｘｍｉｎが設定され、座標ｚにその最小値ｚｍｉｎが設定される。
【０１０２】
続いてステップＳ５２２〜Ｓ５２４で、Ｘ＝ｘｍｉｎの列（図１５上で最も左側の列）のボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が仮のエネルギ関数値用配列ｅ（ｚ）の初期値としてそれぞれ設定される。
【０１０３】
続いてステップＳ５２５でｘがインクリメントされ、続いてステップＳ５２６で座標ｚがその最小値ｚｍｉｎに設定される。
【０１０４】
続いてステップＳ５２７〜Ｓ５２９で、Ｘ＝ｘｍｉｎ＋１の列のエネルギ関数ｅ（ｚ）が計算される。ステップＳ５２７において、ＭＩＮ（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０），ｚ０＝ｚｍｉｎ〜ｚｍａｘ）は、現在の座標ｘ（ここではｘ＝ｘｍｉｎ＋１）において、求めようとする対象物体１０の表面が座標ｚであると仮定し、かつ、座標ｘ−１における表面が座標ｚ０（ｚ０はｚｍｉｎ〜ｚｍａｘのいずれか）であると仮定したとき、評価値（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０））のうち最小のもの（ｍｉｎ）である。ｄｉｆ（ｚ，ｚ０）は、たとえば、表面ｘの座標ｚから表面ｘ−１の座標ｚ０までの距離である。あるいは、ｄｉｆ（ｚ，ｚ０）は、表面ｘの座標ｚから表面ｘ−１の座標ｚ０までの傾きと、概略表面ｘの座標ｚ’から概略表面ｘ−１の座標ｚ０’までの傾きとの差である。このＭＩＮ（ｅ（ｚ０）＋ｄｉｆ（ｚ，ｚ０））にＸ＝ｘｍｉｎ＋１の列の各ボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）が加算され、各ボクセルについて仮のエネルギ関数値ｅ（ｚ）が計算される。
【０１０５】
また、上記最小値ｍｉｎを与えるような座標ｚ０の値が変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）に設定される。これにより、Ｘ＝ｘｍｉｎ−１において、座標ｚの値をｚと仮定したときの、仮のエネルギ関数値ｅ（ｚ）が最小になるための座標ｘ−１における座標ｚがそれぞれ記憶される。
【０１０６】
続いてステップＳ５３０およびＳ５３１で、上記ステップＳ５２７〜Ｓ５２９の処理がｘｍｉｎ〜ｘｍａｘ−１のすべての座標ｘに対応するｚ座標列について行なわれる。したがって、エネルギ関数ｅ（ｚ）にはｘｍｉｎから現在の座標ｘまでの各ｚに対する最小エネルギ関数値が記録されている。
【０１０７】
上記の結果、計算されたエネルギ関数ｅ（ｚ）のうち最小値を与える座標ｚが算出される。この座標ｚがｘｍａｘ−１の列における対象物体１０の表面の位置である。ここで、変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）を見れば、現在よりも１つ前の座標ｘ−１における表面のｚ座標が記録されている。したがって、座標ｘを順次デクリメントし、その座標ｘを変数ｐｒｅＩＤ（ｘ，ｙ，ｚ）に代入すれば、−Ｘ方向に向かって順次対象物体１０の表面の位置を辿っていくことができる。このようにして計算された部分表面Ｓｐ（ｋ）が図１６中に示されている。図１６から明らかなように、部分表面Ｓｐ（ｋ）と決定されたボクセルの標準偏差は各列において必ずしも最小のものではない。仮に各列において最小の標準偏差を持つボクセルが部分表面Ｓｐ（ｋ）を構成するものとして選択されたならば、その部分表面Ｓｐ（ｋ）の連続性が悪くなる。ここでは、部分表面Ｓｐ（ｋ）の連続性を考慮して表面を構成するボクセルが選択されている。
【０１０８】
続いてステップＳ５３２およびＳ５３３で、上記ステップＳ５２１〜Ｓ５３１の処理がｙｍｉｎ〜ｙｍａｘのすべての座標ｙについて行なわれる。
【０１０９】
図１７は、図１０に示されたステップＳ４０で得られた初期のＹ＝ｙ０平面におけるボクセルを示す図である。図１７中において「＊＊」で示されたボクセルは対象物体１０の概略的な存在領域に位置するものである。図１１に示されたステップＳ５４では、図１７では対象物体１０の内部であるが、図１６では外部である各ボクセルから図１６に示された対象物体１０の部分表面Ｓｐ（ｋ）までの距離ｄ（（ｘ，ｙ，ｚ），Ｓｐ（ｋ））が計算される。
【０１１０】
以上のようにこの実施の形態１によれば、物体画像Ａ０〜Ａ１５間における濃度情報の標準偏差に基づいて対象物体１０の形状を推定しているため、対象物体１０の凹所を表現した立体モデル１８を作成することができる。
【０１１１】
また、標準偏差に基づいて対象物体１０の形状を推定する前に、シルエット画像Ｐ０〜Ｐ１５に基づいてボーティング処理により対象物体１０の存在領域を概略的に推定しているため、対象物体１０の立体モデル１８を迅速に生成することができる。
【０１１２】
また、対象物体１０の表面を標準偏差が小さいボクセルを滑らかに結合することによって推定しているため、照明条件などが変動しても安定して正確な立体モデル１８を生成することができる。
【０１１３】
［実施の形態２］
図１８は、本発明の実施の形態２による３次元モデリング装置における形状推定処理（部分表面統合処理を除く）を示すフローチャートである。この実施の形態２では、図１０に示されたステップＳ４５およびＳ４６に代えてステップＳ７０〜Ｓ７４の処理が行なわれる。また、図１０に示されたステップＳ５０に代えてステップＳ７５の処理が行なわれる。また、図１０に示されたステップＳ４２１では変数ｃ２が初期化されているが、ここでは変数ｃ２は用いられていないので初期化されない。その他は、図１０に示されたものと同じ処理が行なわれる。
【０１１４】
すなわち、上記実施の形態１では図１２に示されるようにボクセル１２の画像ｉへの投影点（Ｘ，Ｙ）の画素の濃度情報がボクセル１２に累積されるが、この実施の形態２では図１９に示されるようにボクセル１２の中心画像ｋへの投影点（Ｘｋ，Ｙｋ）を含む複数の画素から構成されるブロック領域Ｂｋがボクセル１２の周辺画像ｉへの投影点（Ｘｉ，Ｙｉ）を含む複数の画素から構成されるブロック領域Ｂｉと図２０に示されるように比較され、濃度情報が異なる画素の数がボクセル１２に累積される。
【０１１５】
より具体的には、ステップＳ７０でボクセル（ｘ，ｙ，ｚ）の中心画像ｋへの投影点（Ｘｋ，Ｙｋ）が計算される。続いてステップＳ７１で、投影点（Ｘｋ，Ｙｋ）を含むｍ×ｎ個の画素からなるブロック領域Ｂｋが設定される。続いてステップＳ７２で、ボクセル（ｘ，ｙ，ｚ）の周辺画像ｉへの投影点（Ｘｉ，Ｙｉ）が計算される。続いてステップＳ７３で、投影点（Ｘｉ，Ｙｉ）を含むｍ×ｎ個の画素からなるブロック領域Ｂｉが設定される。
【０１１６】
続いてステップＳ７４で、ｃ０（ｘ，ｙ，ｚ）がインクリメントされ、これにより各ボクセル（ｘ，ｙ，ｚ）に累積される濃度情報が異なる画素の数がカウントされる。また、ブロック領域Ｂｋがブロック領域Ｂｉと比較され、濃度ＬＭｋ，ＬＭｉの差の絶対値｜ＬＭｋ−ＬＭｉ｜の総和ｄｉｆ（Ｂｋ，Ｂｉ）がｃ１（ｘ，ｙ，ｚ）に加算される。たとえば図２０に示されるように、４×４個の画素からなるブロック領域ＢｋおよびＢｉを設定する場合、１６の画素のうち、３つの画素の濃度情報がブロック領域Ｂｋとブロック領域Ｂｉとの間で相違し、その差の絶対値の合計は「４」であるため、ｄｉｆ（Ｂｋ，Ｂｉ）として「４」がｃ１（ｘ，ｙ，ｚ）に加算される。
【０１１７】
上記の結果、濃度差の絶対値の総和の累積値がｃ１（ｘ，ｙ，ｚ）に設定される。
【０１１８】
その後、ステップＳ７５で、各ボクセル（ｘ，ｙ，ｚ）に累積された濃度差の絶対値の総和の平均値ｓ（ｘ，ｙ，ｚ）が計算される。この平均値ｓ（ｘ，ｙ，ｚ）が濃度情報のばらつき評価値として上記標準偏差の代わりに用いられる。
【０１１９】
以上のようにこの実施の形態２によれば、物体画像間における濃度情報のばらつきをブロックＢｋ，Ｂｉごとに評価しているため、対象物体１０の立体モデル１８をより精度良く生成することができる。
【０１２０】
なお、本実施の形態では、物体画像間における濃度情報のばらつきをブロック間で評価する際、ブロック内の画素ごとの差の絶対値｜ＬＭｋ−ＬＭｉ｜の総和を用いたが、ブロック内の画素ごとの差の自乗（ＬＭｋ−ＬＭｉ）2 の総和を用いてもよい。さらには、ブロック内の画素ごとの標準偏差の平均値を用いてもよい。
【０１２１】
［実施の形態３］
図２１は、本発明の実施の形態３による３次元モデリング装置における部分表面推定処理を示すフローチャートである。上記実施の形態では、できるだけ正確な部分表面Ｓｐ（ｋ）を得るために図１３に示されたダイナミックプログラミング法を用いているが、計算量が多くなる。したがって、以下に詳述する実施の形態３のように、単純に標準偏差が最小のボクセルを部分表面Ｓｐ（ｋ）と推定してもよい。
【０１２２】
まずステップＳ３００で、各中心画像ｋに応じてボクセル空間１１が座標変換される。ここでは、図２２に示されるように、直交座標系のボクセル空間１１が回転され、さらに並進されて再び直交座標系のボクセル空間３０に変換される。座標変換前のボクセル空間１１の各ボクセルの座標を（ｘ，ｙ，ｚ）とし、座標変換後のボクセル空間３０の各ボクセルの座標を（ｘ′，ｙ′，ｚ′）とし、回転変換のためのオイラー角をそれぞれφ、θ、ψとすると、この座標変換は次の式（１）で表わされる。
【０１２３】
【数１】
なお、式（１）中の（ｘｓ，ｙｓ，ｚｓ）は、変換前のボクセル空間１１の中心から中心画像ｋの視点へのベクトルである。回転変換後のボクセル空間はこのベクトル（ｘｓ，ｙｓ，ｚｓ）に応じて並進され、変換後のボクセル空間３０となる。
【０１２５】
このような座標変換が行なわれた後、ステップＳ３０１で座標ｙ′が最小値ｙ′ｍｉｎに初期化され、ステップＳ３０２で座標ｘ′が最小値ｘ′ｍｉｎに初期化され、さらにステップＳ３０３で座標ｚ′が最小値ｚ′ｍｉｎに初期化され、標準偏差一時変数ｓｔｍｐが無限大∞に初期化され、かつ座標一時変数ｚｔｍｐが無限大∞に初期化される。
【０１２６】
続いて、ステップＳ３０３０で存否変数ｖ（ｘ′，ｙ′，ｚ′）が「０」か否か、つまりそのボクセル（ｘ′，ｙ′，ｚ′）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ′，ｙ′，ｚ′）＝０ではない場合、つまりｖ（ｘ′，ｙ′，ｚ′）＝１の場合、ステップＳ３０４で、標準偏差一時変数ｓｔｍｐの値が座標（ｘ′，ｙ′，ｚ′）のボクセルの標準偏差ｓ（ｘ′，ｙ′，ｚ′）と比較される。標準偏差一時変数ｓｔｍｐの値がそのボクセルの標準偏差ｓ（ｘ′，ｙ′，ｚ′）よりも大きい場合は、ステップＳ３０５で標準偏差一時変数ｓｔｍｐに改めてその標準偏差ｓ（ｘ′，ｙ′，ｚ′）が設定され、そのときの座標ｚ′が位置一時変数ｚｔｍｐに設定される。
【０１２７】
ステップＳ３０６およびＳ３０７により、すべての座標ｚ′（＝ｚ′ｍｉｎ〜ｚ′ｍａｘ）について上記ステップＳ３０４および３０５の処理が行なわれる。その結果、変換後のボクセル空間３０において同一のｘ′ｙ′座標をもつ複数のボクセルのうち標準偏差が最小のボクセルが特定される。このようにして特定された標準偏差が最小のボクセルのｚ′座標はステップＳ３０８で表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定される。この表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）は各中心画像ｋに対応して設定される。すなわち、この表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）は、各中心画像ｋ上に撮影されている対象物体の表面位置を示すものである。なお、この時点で位置一時変数ｚｔｍｐが∞であれば、表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に対応した物体表面は存在しないことを意味する。
【０１２８】
図２３は、変換後のボクセル空間３０におけるＹ′＝ｙ′平面を示す。ボクセル空間３０内にはボーティング処理により推定された対象物体１０の概略形状３１が示されている。この概略形状３１内の各ボクセルには標準偏差ｓ（ｘ′，ｙ′，ｚ′）が与えられている。上述したステップＳ３０４〜Ｓ３０７は、たとえば同一のｘ′ｙ′座標をもつ複数のボクセル３２の中から標準偏差ｓ（ｘ′，ｙ′，ｚ′）が最小のボクセルを特定するものである。ここでは、ｚ′＝３７のボクセルが最小の標準偏差を有している。上述したステップＳ３０８は、この最小の標準偏差をもつボクセルの座標ｚ′を表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）の対応する座標に与えるものである。
【０１２９】
続いてステップＳ３０９およびＳ３１０により、Ｙ′＝ｙ′平面上のすべての座標ｘ′（＝ｘ′ｍｉｎ〜ｘ′ｍａｘ）について上記ステップＳ３０３〜Ｓ３０８の処理が行なわれる。また、ステップＳ３１１およびＳ３１２により、すべての座標ｙ′（＝ｙ′ｍｉｎ〜ｙ′ｍａｘ）について上記ステップＳ３０２〜Ｓ３１０の処理が行なわれる。これにより、対象物体１０の部分表面Ｓｐ（ｋ）が決定され、その部分表面Ｓｐ（ｋ）のボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる。
【０１３０】
上述した方法は図１３のダイナミックプログラミング法に代えて用いられるもので、以下「座標変換による最小ばらつき探索法」と呼ぶ。
【０１３１】
この方法により各中心画像ｋごとに得られた部分表面Ｓｐ（ｋ）は上述した距離ペナルティ法または後に詳述する単純ボーティング法により１つの表面に統合され、これにより対象物体１０の形状が決定される。たとえば図２３に示されるように、●で示されるボクセルは対象物体１０の内部にあると判断され、○で表わされるボクセルは対象物体１０の外部にあると判断される。
【０１３２】
このように単純に最小の標準偏差をもつボクセルを特定することにより対象物体１０の形状を推定すれば、上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができる。
【０１３３】
［実施の形態３の変形例］
上記実施の形態においては、直交座標系のボクセル空間１１で同一のｘｙ座標をもつ複数のボクセルのうち標準偏差が最小のボクセルが必ずしも対象物体１０の真の表面１４に位置していない。これは画像ノイズに起因するものもあるが、主として計算方法に起因するものである。この理由を図９および図２４を用いて詳細に説明する。
【０１３４】
図９に示されるように、直交座標系のボクセル空間１１において、同一のｘｙ座標をもつすべてのボクセル１２から中心画像となる物体画像Ａ０への投影線１６が物体画像Ａ０上の同じ位置に投影する場合は、照明条件の著しい変動がない限り、真の表面に位置するボクセルの物体画像Ａ０への投影点の色彩が複雑であるか単一であるかに関わらず、その同一のｘｙ座標をもつ複数のボクセル１２のうち標準偏差が最小のボクセル１２が対象物体１０の真の表面上に位置すると考えてもよい。しかしながら、上述したボクセル１２から中心画像となる物体画像Ａ０への投影線１６は、透視投影に基づいているため、図２４に示されるように、同一のｘｙ座標をもつ複数のボクセル１２から物体画像Ａ０への投影線１６は物体画像Ａ０上の異なる位置に投影することが多い。この場合、真の表面１４上に位置するボクセル１２から物体画像Ａ０上への投影点付近には複雑な色彩があり、その他のボクセル１２から中心画像である物体画像Ａ０および周辺画像である物体画像Ａ１およびＡ１５上への投影点には単一の色彩があるとすると、真の表面１４上に位置しないボクセルの標準偏差が真の表面１４上に位置するボクセルの標準偏差よりも小さくなる場合がある。
【０１３５】
そこで、図２１に示されたステップＳ３００における直交座標系から直交座標系への座標変換に代えて、直交座標系から図２５に示されるような透視座標系への座標変換が行なわれてもよい。この透視座標系のボクセル空間３５は、中心画像ｋの視点を頂点とした錐体状をなしている。
【０１３６】
図２６（ａ）は図２５に示された透視座標系のボクセル空間３５をＹ方向から見た図であり、図２６（ｂ）は図６に示された直交座標系のボクセル空間１１をＹ方向から見た図である。中心画像ｋの視点からボクセル空間１１の中心までの距離をｌとすると、この直交座標系から透視座標系への座標変換は次の式（２）により表わされる。
【０１３７】
ｘ′＝ｘ・ｌ／ｚ，ｙ′＝ｙ・ｌ／ｚ，ｚ′＝ｚ …（２）
このような透視座標系のボクセル空間３５においても上記と同様に、図２７に示されるような同一のｘ′ｙ′座標をもつ複数のボクセル３６の中から標準偏差が最小のボクセルが特定され、その特定されたボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる。ただし、ここで表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に与えられる値はこの変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）が設定される中心画像からその標準偏差が最小のボクセルまでの距離である。このように、透視座標系のボクセル空間３５において各中心画像ｋの視点を通る直線状のボクセルのうち標準偏差が最小のボクセルが結合され、対象物体１０の部分表面Ｓｐ（ｋ）が決定される。
【０１３８】
なお、部分表面統合処理により、図２７中の●のボクセルは対象物体１０の内部と判定され、○のボクセルは対象物体１０の外部と判定される。
【０１３９】
このように透視座標系のボクセル空間に座標変換を行えば、最小の標準偏差をもつボクセルを特定することにより対象物体１０の形状を簡単かつ正確に推定することができる。また、このような座標変換による最小ばらつき法は上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができ、その結果、この処理に必要なメモリのワークエリアを削減するこができる。
【０１４０】
［実施の形態４］
図２８は、本発明の実施の形態４による３次元モデリング装置における部分表面推定処理を示すフローチャートである。ここでは、図２１に示されたステップＳ３０８に代えて、ステップＳ３２０〜Ｓ３２４が設けられている。この方法は上記実施の形態３の部分表面Ｓｐ（ｋ）上またはそれよりも外側に部分表面を特定するもので、以下「座標変換による準最小ばらつき探索法」と呼ぶ。
【０１４１】
上述したようにステップＳ３０４〜Ｓ３０７で同一のｘ′ｙ′座標をもつボクセルの中から標準偏差が最小のボクセルが特定される。
【０１４２】
ステップＳ３２０では、その最小の標準偏差に所定の値αを加えたものが最小標準偏差変数ｓｍｉｎに設定される。そして、ステップＳ３２０〜Ｓ３２３で、再び同一のｘ′ｙ′座標をもつ複数のボクセルの中から［最小の標準偏差＋α］以下の標準偏差をもつボクセルが特定される。このとき、ボクセルはｚ′ｍｉｎからｚ′ｍａｘに向かって順に探索される。したがって、最小の標準偏差をもつボクセルよりも外側に、［最小の標準偏差＋α］以下の標準偏差をもつボクセルが存在する場合は、そのボクセルの座標ｚ′がステップＳ３２４で表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定される。この場合、図２３または図２７に示された部分表面Ｓｐ（ｋ）よりも外側に部分表面が特定される。なお、最小の標準偏差をもつボクセルよりも外側に、［最小の標準偏差＋α］以下の標準偏差をもつボクセルが存在しない場合は、最小の標準偏差をもつボクセルの座標ｚ′が表面変数ｓｕｒｆａｃｅ（ｘ′，ｙ′）に設定され、図２３または図２７に示された部分表面Ｓｐ（ｋ）上に部分表面が特定される。
【０１４３】
以上のようにこの実施の形態４によれば、最小の標準偏差よりも少しだけ大きい標準偏差をもつボクセルが特定され、これにより部分表面Ｓｐ（ｋ）が推定されているため、照明条件の変動などによりノイズが入っても、一部を欠いた立体モデルが生成されることはない。
【０１４４】
［実施の形態５］
図２９は、本発明の実施の形態５による３次元モデリング装置における部分表面推定処理を示すフローチャートである。この方法は、上記実施の形態のようにボクセル空間の座標変換を行なうことなく視点からボクセル空間への透視性を考慮して標準偏差が最小のボクセルを特定して対象物体１０の部分表面を推定するもので、以下「光線追跡による最小ばらつき探索法」と呼ぶ。
【０１４５】
まずここでは、図３０に示されるような座標（ｐｘ，ｐｙ，ｐｚ）を定義する。すなわち、中心画像ｋの視点の座標を（ｏｘ，ｏｙ，ｏｚ）と定義する。そして、この視点（ｏｘ，ｏｙ，ｏｚ）からボクセル空間１１の中心に向かう方向にＺ′座標軸を設け、中心画像ｋ上にＸ′座標軸およびＹ′座標軸を設ける。視点（ｏｘ，ｏｙ，ｏｚ）から放射状に出る光線４０が中心画像ｋを貫通するｘ′ｙ′座標を（ｐｘ，ｐｙ）とし、中心画像ｋのｚ′座標をｐｚとする。
【０１４６】
図２９に示されるように、まずステップＳ５００で座標ｐｙに最小値ｐｙｍｉｎが設定され、続いてステップＳ５０１で座標ｐｘに最小値ｐｘｍｉｎが設定される。
【０１４７】
続いてステップＳ５０２で、単位光線ベクトルｖ＝（ｖｘ，ｖｙ，ｖｚ）が計算される。ただし、ここでは次の式（３）〜（５）が成立する。
【０１４８】
ｖｚ＝１ …（３）
ｖ×（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）＝０ …（４）
ｖ・（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）≧０ …（５）
ベクトルｖとベクトル（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）との外積が０ということは、これらのベクトルが同一方向に横たわっている（向きは同じか反対か不明）ことを意味する。また、ベクトルｖとベクトル（ｐｘ−ｏｘ，ｐｙ−ｏｙ，ｐｚ−ｏｚ）との内積が０以上ということは、これらのベクトルの向きも同じであることを意味する。
【０１４９】
続いてステップＳ５０３で、整数ｎに最小値ｎｍｉｎが設定され、標準偏差一時変数ｓｔｍｐに無限大∞が設定され、かつ座標一時変数ｚｔｍｐに無限大∞が設定される。
【０１５０】
続いてステップＳ５０４〜Ｓ５０８では、図３１に示されるように、直交座標系のボクセル空間１１において視点（ｏｘ，ｏｙ，ｏｚ）を通る各光線４０上の複数のボクセルの中から標準偏差が最小のボクセルが特定される。図３１では、複数のボクセル（○，●で表わされる）のうち、●で表わされるボクセルが特定されている。
【０１５１】
より具体的には、ステップＳ５０４で視点（ｏｘ，ｏｙ，ｏｚ）を基準としてｎ倍の単位光線ベクトルｖで表わされる座標（ｘ′ｙ′ｚ′）に位置するボクセルが選択される。
【０１５２】
続いてステップＳ５０４０でボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ５０５で、標準偏差一時変数ｓｔｍｐの値がその選択されたボクセルの標準偏差ｓ（ｘ，ｙ，ｚ）と比較される。標準偏差一時変数ｓｔｍｐの値が標準偏差ｓ（ｘ，ｙ，ｚ）よりも大きい場合は、ステップＳ５０６でその標準偏差ｓ（ｘ，ｙ，ｚ）が代わりに標準偏差一時変数ｓｔｍｐに設定され、そのときの自然数ｎが位置一時変数ｚｔｍｐに設定される。
【０１５３】
続いてステップＳ５０７およびＳ５０８により、すべての自然数ｎ（＝ｎｍｉｎ〜ｎｍａｘ）について上記ステップＳ５０４〜Ｓ５０６の処理が行なわれる。そのため、光線４０上に位置する複数のボクセルの中から標準偏差が最小のボクセルが特定される。
【０１５４】
続いてステップＳ５０９で、この特定されたボクセルのｚ′座標（位置一時変数ｚｔｍｐに設定されている）が表面変数ｓｕｒｆａｃｅ（ｐｘ，ｐｙ）に設定される。
【０１５５】
続いてステップＳ５１０およびＳ５１１により、すべての座標ｐｘについて上記ステップＳ５０２〜Ｓ５０９の処理が行なわれる。
【０１５６】
続いてステップＳ５１２およびＳ５１３により、すべての座標ｐｙ（＝ｐｙｍｉｎ〜ｐｙｍａｘ）について上記ステップＳ５０１〜Ｓ５１１の処理が行なわれる。
【０１５７】
図３１では、視点（ｏｘ，ｏｙ，ｏｚ）を通るある１つの光線４０上の複数のボクセルの中から最小の標準偏差をもつ●で表わされるボクセルが特定される。この特定されたボクセルのｐｚ座標は「３７」であるので、この「３７」が中心画像ｋ上の表面変数ｓｕｒｆａｃｅ（ｘ，ｙ）の対応する座標に設定される。
【０１５８】
上記のように各光線４０上ごとに最小の標準偏差をもつボクセルのｐｚ座標が表面変数ｓｕｒｆａｃｅ（ｘ，ｙ）に設定されることにより、図３２および図３３に示されるように対象物体１０の部分表面Ｓｐ（ｋ）が特定される。
【０１５９】
このようにして特定された複数の部分表面Ｓｐ（ｋ）は、上述した距離ペナルティ法または後に詳述する単純ボーティング法により１つに統合され、これにより対象物体１０の形状が決定される。図３２に示されるように、部分表面Ｓｐ（ｋ）よりも外側に位置する○で表わされるボクセルは対象物体の外部に位置し、部分表面Ｓｐ（ｋ）よりも内側に位置する●で表わされるボクセルは対象物体１０の内部に位置している。
【０１６０】
以上のようにこの実施の形態５によれば、視点から放射上に出る光線上の複数のボクセルのうち最小の標準偏差をもつボクセルが特定されているため、対象物体１０の形状を比較的正確に推定することができる。このような光線追跡による最小ばらつき探索法は上記ダイナミックプログラミング法よりも高速で部分表面Ｓｐ（ｋ）を決定することができ、その結果、この処理に必要なメモリのワークエリアを削減するこができる。
【０１６１】
［実施の形態６］
図３４は、本発明の実施の形態６による３次元モデリング装置における部分表面統合処理を示すフローチャートである。この方法は、上述した距離ペナルティ法に代えて用いられ、上記部分表面推定処理により得られた複数の部分表面Ｓｐ（ｋ）を１つに統合して対象物体の形状を決定するものである。以下、この方法を「単純ボーティング法」と呼ぶ。
【０１６２】
図３５（ａ）〜（ｃ）は、各観察方向からの部分表面Ｓｐ（ｋ）が決定されたボクセル空間を示す。ここでは、推定された概略形状３１よりも内側でかつ部分表面Ｓｐ（ｋ）よりも内側のボクセルに得点「１」が与えられている。このような得点がボクセル空間に投票されると、図３５（ｄ）に示されるようなボクセル空間が得られる。ここで、対象物体１０は１６個の視点から観察されているため、各ボクセルの最高得点は「１６」である。この得点数が高いボクセルほど対象物体１０の真の存在領域である可能性が高いが、ここでは表面条件の変動などによるノイズを考慮して投票数が「１４」以上のボクセルを図３５（ｅ）に示されるように対象物体１０の真の存在領域４５として推定している。
【０１６３】
より具体的には、図３４に示されるように、ステップＳ６００で変数ｃｎｔ（ｘ，ｙ，ｚ）が初期化され、ステップＳ６００１で中心画像ｋとして物体画像Ａ０が設定される。続いてステップＳ６０１〜Ｓ６０３で、座標ｙ、ｘ、ｚにそれぞれ最小値ｙｍｉｎ、ｘｍｉｎ、ｚｍｉｎが設定される。
【０１６４】
続いてステップＳ６０４で、図１０に示されたステップＳ４０の結果に従って座標（ｘ，ｙ，ｚ）のボクセルが対象物体１０の概略存在領域内か否かが判定される。概略存在領域内、つまりｖ（ｘ，ｙ，ｚ）＝１の場合は、ステップＳ６０５で中心画像ｋに応じて座標変換が行なわれる。ここでは、座標（ｘ，ｙ，ｚ）が上述した座標（ｘ′，ｙ′，ｚ′）に変換される。
【０１６５】
続いてステップＳ６０６で、このボクセルが部分表面Ｓｐ（ｋ）よりも内側にあるか否か、つまりｓｕｒｆａｃｅ（ｘ′，ｙ′）≦ｚ′か否かが判定される。このボクセルが部分表面Ｓｐ（ｋ）よりも内側にある場合は、ステップＳ６０７で変数ｃｎｔ（ｘ，ｙ，ｚ）がインクリメントされる。
【０１６６】
続いてステップＳ６０８〜Ｓ６１３により、すべての座標ｘ（＝ｘｍｉｎ〜ｘｍａｘ）、すべての座標ｙ（＝ｙｍｉｎ〜ｙｍａｘ）、すべての座標ｚ（＝ｚｍｉｎ〜ｚｍａｘ）について上記ステップＳ６０４〜Ｓ６０７の処理が行なわれる。これにより図３５（ａ）〜（ｃ）に示されるようなボクセル空間が得られる。
【０１６７】
続いてステップＳ６１４およびＳ６１５により、すべての中心画像ｋ（＝０〜ｋｍａｘ）について上記ステップＳ６０１〜Ｓ６１３の処理が行なわれる。これにより、図３５（ｄ）に示されるようなボクセル空間が得られる。
【０１６８】
続いてステップＳ６１５０でボクセル（ｘ，ｙ，ｚ）が対象物体１０の存在領域か否かが判定される。ｖ（ｘ，ｙ，ｚ）＝１の場合、ステップＳ６１６で、各ボクセルの変数ｃｎｔ（ｘ，ｙ，ｚ）が所定の値（ここでは「１４」）と比較される。変数ｃｎｔ（ｘ，ｙ，ｚ）が所定の値よりも小さい場合は、ステップＳ６１７でそのボクセルは対象物体１０の存在領域外、つまり変数ｖ（ｘ，ｙ，ｚ）に０が設定される。
【０１６９】
続いてステップＳ６１８により、すべてのボクセルについて上記ステップＳ６１６およびＳ６１７の処理が行なわれ、これにより図３５（ｅ）に示されるようなボクセルが得られる。これにより対象物体１０の形状４５が決定される。
【０１７０】
以上のようにこの実施の形態６によれば、推定された複数の部分表面Ｓｐ（ｋ）を単純ボーティング法により統合しているため、上述した距離ペナルティ法よりも高速で対象物体１０の形状を決定することができる。
【０１７１】
以上、本発明の実施の形態を説明したが、本発明は上記実施の形態に限定されるものではない。たとえば、上記実施の形態は画像情報として濃度を用いているが、これに代えて色を用いてもよく、さらに色および濃度の組合せを用いてもよい。また、上記実施の形態は画像情報の標準偏差を用いているが、これに代えて画像情報の分散を用いてもよく、要するに、画像情報のばらつきを評価する値を用いればよい。
【０１７２】
【発明の効果】
以上のように本発明によれば、対象物体を複数の視点から撮影して得られた複数の物体画像のうち少なくとも２つの物体画像間における画像情報のばらつきに基づいて対象物体の形状を推定しているため、対象物体の凹所を表現した立体モデルを生成することができる。
【０１７３】
また、予め対象物体の輪郭に基づいてボクセル空間における対象物体の存在領域を概略的に推定しているため、対象物体の立体モデルを迅速に生成することができる。
【０１７４】
また、ばらつきが小さいボクセルを滑らかに結合することにより対象物体の形状を推定しているため、照明条件などが変動しても正確な立体モデルを安定して生成することができる。
【０１７５】
また、３次元仮想空間の座標空間を透視座標空間に変換し、その透視座標空間において１つの視点を通る直線上の複数の領域のうち、ばらつきが最小またはその最小値よりも所定の値だけ大きい領域を結合しているため、対象物体の形状を高速で推定することができる。
【０１７６】
また、３次元仮想空間において１つの視点を通る直線上の複数の領域のうち、ばらつきが最小またはその最小値よりも所定の値だけ大きい領域を結合しているため、対象物体の形状を高速で推定することができる。
【０１７７】
また、３次元仮想空間の領域のうち存在領域の内側にありかつ部分表面の各々の外側にある領域にその領域から複数の部分表面までの距離を累積し、その累積された距離の総和が所定の値よりも大きい領域を対象物体の外側と認定しているため、対象物体の形状を比較的正確に決定することができる。
【０１７８】
また、３次元仮想空間の領域のうち部分表面の各々の内側にある領域に所定の値を投票し、その投票された値の総和が所定の値よりも大きい領域を対象物体の内側と認定しているため、対象物体の形状を高速で決定することができる。
【０１７９】
また、複数の画素を含む所定ブロックごとに画像情報を累積しているため、対象物体の立体モデルをより精度良く生成することができる。
【０１８０】
また、各物体画像内に複数の画素を含むブロックを設定し、そのブロック間において互いに対応する画素の値の差をそれぞれ算出し、その算出された差の絶対値の総和を画像情報のばらつきとして算出しているため、対象物体の形状を高速で推定することができる。
【０１８１】
また、各物体画像内に複数の画素を含むブロックを設定し、そのブロック間において互いに対応する画素の値の差をそれぞれ算出し、その算出された差の自乗の総和を画像情報のばらつきとして算出しているため、対象物体の形状を高速で推定することができる。
【図面の簡単な説明】
【図１】本発明の実施の形態１による３次元モデリング装置の全体構成を示す図である。
【図２】図１に示されたコンピュータの主要構成を示すブロック図である。
【図３】図１に示されたＣＤ−ＲＯＭに記録された３次元モデリングプログラムを示すフローチャートである。
【図４】（ａ）は図３中の対象物体および背景の撮影ステップを説明するための図であり、（ｂ）は図３中のシルエット画像作成ステップを説明するための図であり、（ｃ）は図３中の形状推定処理ステップを説明するための図であり、（ｄ）は図３中のポリゴン作成ステップを説明するための図であり、（ｅ）は図３中のテクスチャマッピングステップを説明するための図である。
【図５】図４（ａ）において対象物体を撮影する視点を説明するための図である。
【図６】図４（ｃ）の形状推定処理で用いられるボクセル空間を示す概念図である。
【図７】図３に示された形状推定処理の詳細を示すフローチャートである。
【図８】図７に示された概略形状推定処理を説明するための図である。
【図９】図７に示された色ボーティング処理を説明するための図である。
【図１０】図７に示された概略形状推定処理、色ボーティング処理、および部分表面推定処理の詳細を示すフローチャートである。
【図１１】図７に示された部分表面統合処理の詳細を示すフローチャートである。
【図１２】図１０中のステップＳ４５で用いられる投影点を説明するための図である。
【図１３】図１０中のステップＳ５２で用いられるダイナミックプログラミング法を示すフローチャートである。
【図１４】図１０に示されたフローチャートによって計算されたＹ＝ｙ０平面における標準偏差を示す鳥瞰図である。
【図１５】図１０に示されたフローチャートによって計算されたＹ＝ｙ０平面における標準偏差の数値を示す図である。
【図１６】図１３に示されたフローチャートによって決定された部分表面を示す図である。
【図１７】図１０中のステップＳ４０で得られた対象物体の存在領域の初期形状を示す図である。
【図１８】本発明の実施の形態２による３次元モデリング装置における形状推定処理を示すフローチャートである。
【図１９】図１７中のステップＳ７０およびＳ７２で用いられる投影点を説明するための図である。
【図２０】図１７中のステップＳ７４においてブロック領域間で濃度が異なる画素の数を計算する方法を説明するための図である。
【図２１】本発明の実施の形態３による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図２２】図２１中のステップＳ３００における直交座標系のボクセル空間から直交座標系のボクセル空間への座標変換を説明するための図である。
【図２３】図２１に示されたフローチャートに従って直交座標系のボクセル空間から直交座標系のボクセル空間への座標変換を行なった場合における部分表面の特定およびその統合方法を説明するための図である。
【図２４】標準偏差が最小のボクセルが必ずしも対称物体の部分表面上に位置するとは限らない理由を明らかにするために、同一のｘｙ座標をもつ複数のボクセルから中心画像への投影点が異なる場合における色ボーティング処理を説明するための図である。
【図２５】図２１中のステップＳ３００における座標変換で使用可能な透視座標系のボクセル空間を示す斜視図である。
【図２６】（ａ）は図２５に示された透視座標系のボクセル空間をＹ方向から見た上面図であり、（ｂ）は図６に示された直交座標系のボクセル空間をＹ方向から見た上面図である。
【図２７】図２１中のステップＳ３００で図２５に示されるような透視座標系のボクセル空間への座標変換が行なわれた場合における部分表面の特定およびその統合方法を説明するための図である。
【図２８】本発明の実施の形態４による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図２９】本発明の実施の形態５による３次元モデリング装置における部分表面推定処理の詳細を示すフローチャートである。
【図３０】図２９に示されたフローチャートで用いられる座標系を説明するための図である。
【図３１】図２９に示された光線追跡による最小ばらつき探索法を説明するための図である。
【図３２】図２９に示された光線追跡による最小ばらつき探索法により特定された部分表面の統合方法を説明するための図である。
【図３３】図２９に示された光線追跡による最小ばらつき探索法により１つの部分表面が特定されたボクセルを示す図である。
【図３４】本発明の実施の形態６による立体モデリング装置における部分表面統合処理の詳細を示すフローチャートである。
【図３５】（ａ）〜（ｅ）は、図３４に示された単純ボーティング法を説明するための図である。
【符号の説明】
２ デジタルスチルカメラ
３ コンピュータ
４ ＣＤ−ＲＯＭ
５ 中央処理装置
７ ランダムアクセスメモリ
１０ 対象物体
１１ ボクセル空間
１２ ボクセル
１８ 立体モデル
Ａ０〜Ａ１５ 物体画像
Ｂ０〜Ｂ１５ シルエット画像
Ｂｋ，Ｂｉ ブロック領域
Claims (12)
- 対象物体の立体モデルを生成する３次元モデリング装置であって、
前記対象物体を複数の視点から撮影し、複数の物体画像を得る物体画像獲得手段と、
前記複数の物体画像のうち少なくとも２つの物体画像間において互いに対応する画像情報についての濃度情報及び／又は色情報の標準偏差に基づいて前記対象物体の形状を推定する形状推定手段と、
前記推定された形状に基づいて前記立体モデルを生成する生成手段とを備え、
前記形状推定手段は、
複数の領域から構成される３次元仮想空間の領域に前記少なくとも２つの物体画像に対応する画像情報についての濃度情報及び／又は色情報を累積する画像情報累積手段と、
前記領域に累積された画像情報についての濃度情報及び／又は色情報の標準偏差を算出する算出手段とを含む、３次元モデリング装置。 - 前記複数の物体画像中の前記対象物体の輪郭に基づいて前記３次元仮想空間における前記対象物体の存在領域を概略的に推定する存在領域推定手段をさらに備え、
前記画像情報累積手段は、前記３次元仮想空間の領域のうち前記存在領域内の領域に前記画像情報についての濃度情報及び／又は色情報を累積する、請求項１に記載の３次元モデリング装置。 - 前記形状推定手段はさらに、前記対象物体の形状を推定するために、前記３次元仮想空間の領域のうち前記標準偏差の小さい領域を滑らかに結合する結合手段を含む、請求項１に記載の３次元モデリング装置。
- 前記結合手段は、前記標準偏差の小さい領域を滑らかに結合するために、前記標準偏差の小ささと前記結合されるべき領域の滑らかさとを変数とする評価関数を最小化する、請求項３に記載の３次元モデリング装置。
- 前記形状推定手段はさらに、
前記３次元仮想空間の座標空間を前記視点の１つから放射状に広がる座標軸を有する透視座標空間に変換する座標変換手段と、
前記対象物体の形状を推定するために、前記透視座標空間において前記１つの視点を通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段とを含む、請求項１から請求項２のいずれかに記載の３次元モデリング装置。 - 前記形状推定手段はさらに、
前記対象物体の形状を推定するために、前記３次元仮想空間において前記視点の１つを通る直線上の領域のうち前記標準偏差が最小またはその最小値よりも所定の値だけ大きい領域を結合する結合手段を含む、請求項１から請求項２のいずれかに記載の３次元モデリング装置。 - 前記形状推定手段はさらに、
前記複数の物体画像間において互いに対応する画像情報の標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、
前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含む、請求項１から請求項２のいずれかに記載の３次元モデリング装置。 - 前記形状推定手段はさらに、
前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、
前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、
前記統合手段は、
前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、
前記累積された距離の総和が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む、請求項２に記載の３次元モデリング装置。 - 前記形状推定手段はさらに、
前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、
前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、
前記統合手段は、
前記３次元仮想空間の領域のうち前記存在領域推定手段で推定された存在領域内にありかつ前記部分表面の各々の外側にある領域に、その領域から前記複数の部分表面までの距離を累積する距離累積手段と、
前記累積回数を計算する累積回数計算手段と、
前記累積された距離の総和を前記計算された累積回数で除した距離平均値が所定の値よりも大きい領域を前記対象物体の外側と認定する認定手段とを含む、請求項２に記載の３次元モデリング装置。 - 前記形状推定手段はさらに、
前記標準偏差に基づいて前記対象物体の複数の部分表面を推定する部分表面推定手段と、
前記複数の部分表面を統合して前記対象物体の形状を決定する統合手段とを含み、
前記統合手段は、
前記３次元仮想空間の領域のうち前記部分表面の各々の内側にある領域に所定の第１の値を投票する投票手段と、
前記投票された第１の値の総和が所定の第２の値よりも大きい領域を前記対象物体の内側と認定する認定手段とを含む、請求項１に記載の３次元モデリング装置。 - 前記画像情報累積手段はさらに、前記物体画像の複数の画素を含む所定ブロックごとに前記画像情報についての濃度情報及び／又は色情報を累積する、請求項１から請求項１０のいずれかに記載の３次元モデリング装置。
- 前記領域はボクセルであり、前記３次元仮想空間はボクセル空間である、請求項１から請求項１１のいずれかに記載の３次元モデリング装置。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP00598298A JP3813343B2 (ja) | 1997-09-09 | 1998-01-14 | ３次元モデリング装置 |
| US09/149,046 US6256036B1 (en) | 1997-09-09 | 1998-09-08 | Three-dimensional modeling apparatus for producing three-dimensional model of object of interest and method therefor |
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP24365097 | 1997-09-09 | ||
| JP9-243650 | 1997-09-09 | ||
| JP00598298A JP3813343B2 (ja) | 1997-09-09 | 1998-01-14 | ３次元モデリング装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH11149576A JPH11149576A (ja) | 1999-06-02 |
| JP3813343B2 true JP3813343B2 (ja) | 2006-08-23 |
Family
ID=26340032
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP00598298A Expired - Fee Related JP3813343B2 (ja) | 1997-09-09 | 1998-01-14 | ３次元モデリング装置 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US6256036B1 (ja) |
| JP (1) | JP3813343B2 (ja) |
Families Citing this family (33)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3954211B2 (ja) * | 1998-08-20 | 2007-08-08 | 富士通株式会社 | 三次元シーンにおける形状及び模様の復元方法及び装置 |
| US6819318B1 (en) * | 1999-07-23 | 2004-11-16 | Z. Jason Geng | Method and apparatus for modeling via a three-dimensional image mosaic system |
| JP4474743B2 (ja) * | 2000-07-03 | 2010-06-09 | ソニー株式会社 | 三次元画像生成装置および三次元画像生成方法、並びにプログラム記録媒体 |
| US7706600B2 (en) * | 2000-10-02 | 2010-04-27 | The Research Foundation Of State University Of New York | Enhanced virtual navigation and examination |
| US7423666B2 (en) * | 2001-05-25 | 2008-09-09 | Minolta Co., Ltd. | Image pickup system employing a three-dimensional reference object |
| US6914600B2 (en) * | 2001-06-26 | 2005-07-05 | Hewlett-Packard Development Company, L.P. | Volumetric warping for voxel coloring on an infinite domain |
| JP4035978B2 (ja) * | 2001-10-05 | 2008-01-23 | コニカミノルタホールディングス株式会社 | ３次元形状モデルの評価方法および生成方法並びに装置 |
| JP2003141565A (ja) * | 2001-11-06 | 2003-05-16 | Kubota Corp | ３次元形状復元システムおよび方法 |
| US7289662B2 (en) * | 2002-12-07 | 2007-10-30 | Hrl Laboratories, Llc | Method and apparatus for apparatus for generating three-dimensional models from uncalibrated views |
| FR2887058B1 (fr) * | 2005-06-10 | 2007-08-31 | Daniel Ackerman | Procede et dispositif de reconstruction 3d d'un objet a partir de plusieurs images 2d |
| WO2007069724A1 (ja) | 2005-12-16 | 2007-06-21 | Ihi Corporation | 三次元形状データの位置合わせ方法と装置 |
| CN101331379B (zh) | 2005-12-16 | 2012-04-11 | 株式会社Ihi | 自身位置辨认方法和装置以及三维形状的计测方法和装置 |
| DE112006003361T5 (de) * | 2005-12-16 | 2008-10-16 | Ihi Corporation | Verfahren und Vorrichtung für das Aufzeichnen/Anzeigen von dreidimensionalen Formdaten und Verfahren und Vorrichtung für das Messen einer dreidimensionalen Form |
| JP4839237B2 (ja) * | 2007-02-07 | 2011-12-21 | 日本電信電話株式会社 | ３次元形状復元方法，３次元形状復元装置，その方法を実装した３次元形状復元プログラム及びそのプログラムを記録した記録媒体 |
| CA2717485A1 (en) * | 2007-03-02 | 2008-09-12 | Organic Motion | System and method for tracking three dimensional objects |
| US8315477B2 (en) * | 2007-11-14 | 2012-11-20 | Integraph Software Technologies Company | Method and apparatus of taking aerial surveys |
| JP5218034B2 (ja) * | 2008-12-26 | 2013-06-26 | Ｋｄｄｉ株式会社 | マスク画像を抽出する方法及びプログラム並びにボクセルデータを構築する方法及びプログラム |
| JP5396971B2 (ja) * | 2009-03-31 | 2014-01-22 | カシオ計算機株式会社 | 位置検索システム及び位置検索方法 |
| TWI409717B (zh) * | 2009-06-22 | 2013-09-21 | Chunghwa Picture Tubes Ltd | 適用於電腦產品與影像顯示裝置的影像轉換方法 |
| JP5645010B2 (ja) * | 2010-10-14 | 2014-12-24 | 国立大学法人 千葉大学 | 立体画像を作成する方法 |
| US8587583B2 (en) * | 2011-01-31 | 2013-11-19 | Microsoft Corporation | Three-dimensional environment reconstruction |
| US8942917B2 (en) | 2011-02-14 | 2015-01-27 | Microsoft Corporation | Change invariant scene recognition by an agent |
| EP2680594A4 (en) * | 2011-02-24 | 2014-07-09 | Kyocera Corp | ELECTRONIC APPARATUS, IMAGE DISPLAY METHOD, AND IMAGE DISPLAY PROGRAM |
| JP5561786B2 (ja) * | 2011-03-30 | 2014-07-30 | Ｋｄｄｉ株式会社 | ３次元形状モデル高精度化方法およびプログラム |
| JP5777507B2 (ja) * | 2011-12-27 | 2015-09-09 | キヤノン株式会社 | 情報処理装置、情報処理方法、及びそのプログラム |
| JP6181363B2 (ja) * | 2012-11-06 | 2017-08-16 | 敏 伊木 | 画像処理装置、画像処理方法、動画作成方法、動画作成システム及びプログラム |
| US9836879B2 (en) * | 2013-04-16 | 2017-12-05 | Autodesk, Inc. | Mesh skinning technique |
| US9269018B2 (en) * | 2014-01-14 | 2016-02-23 | Microsoft Technology Licensing, Llc | Stereo image processing using contours |
| JP6816442B2 (ja) * | 2016-10-19 | 2021-01-20 | 富士ゼロックス株式会社 | データ処理装置、立体造形システム及びプログラム |
| US11100677B2 (en) | 2017-09-22 | 2021-08-24 | Nec Corporation | Information processing device, information processing method and recording medium |
| WO2019116942A1 (ja) * | 2017-12-14 | 2019-06-20 | キヤノン株式会社 | ３次元モデルの生成装置、生成方法、及びプログラム |
| JP6513169B1 (ja) | 2017-12-14 | 2019-05-15 | キヤノン株式会社 | 仮想視点画像を生成するシステム、方法及びプログラム |
| CA3061699C (en) * | 2017-12-14 | 2023-03-07 | Canon Kabushiki Kaisha | Generation device, generation method and non-transitory computer-readable meadium storing a program for generating a three-dimensional model from an image |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0896165A (ja) * | 1994-09-29 | 1996-04-12 | Fujitsu Ltd | 画像生成方法及び画像生成装置 |
| JPH08186844A (ja) * | 1994-12-28 | 1996-07-16 | Sanyo Electric Co Ltd | 立体映像発生装置および立体映像発生方法 |
| US5982375A (en) * | 1997-06-20 | 1999-11-09 | Sun Microsystems, Inc. | Floating point processor for a three-dimensional graphics accelerator which includes single-pass stereo capability |
-
1998
- 1998-01-14 JP JP00598298A patent/JP3813343B2/ja not_active Expired - Fee Related
- 1998-09-08 US US09/149,046 patent/US6256036B1/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| US6256036B1 (en) | 2001-07-03 |
| JPH11149576A (ja) | 1999-06-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3813343B2 (ja) | ３次元モデリング装置 | |
| JP6321106B2 (ja) | 現実環境内にバーチャルオブジェクトを描写する方法および装置 | |
| US6463176B1 (en) | Image recognition/reproduction method and apparatus | |
| JP2919284B2 (ja) | 物体認識方法 | |
| Yamany et al. | Free-form surface registration using surface signatures | |
| JP3735344B2 (ja) | キャリブレーション装置、キャリブレーション方法、及びキャリブレーション用プログラム | |
| JP4785880B2 (ja) | 三次元オブジェクト認識のためのシステムおよび方法 | |
| US6072903A (en) | Image processing apparatus and image processing method | |
| US8781161B2 (en) | Image processing method and apparatus for generating a 3D model of a target object | |
| WO2000036564A2 (en) | Creating a three-dimensional model from two-dimensional images | |
| JP2008145431A (ja) | ３次元表面形状再構築のための装置および方法 | |
| JP2008537190A (ja) | 赤外線パターンを照射することによる対象物の三次元像の生成 | |
| JP2010541053A (ja) | 現実環境内にバーチャルオブジェクトを描写する方法および装置 | |
| JP2001243468A (ja) | 三次元モデリング装置、方法及び媒体並びに三次形状データ記録装置、方法及び媒体 | |
| US5864779A (en) | Strict recognizing apparatus using observation point | |
| CN110546687B (zh) | 图像处理装置及二维图像生成用程序 | |
| JP2004280776A (ja) | 画像内のオブジェクトの形状を決定する方法 | |
| JP2002133413A (ja) | 画像処理を用いる３次元物体の識別方法および装置 | |
| JP2000304499A (ja) | 飛翔体の画像処理装置及び飛翔体の画像処理方法 | |
| JP2003123057A (ja) | ３次元形状モデルの生成方法および装置 | |
| JPH11328445A (ja) | 衝突判定装置および方法、および衝突判定方法を記録した媒体 | |
| JPH1040385A (ja) | 物体形状復元装置及びその方法 | |
| GB2358307A (en) | Method of determining camera projections in 3D imaging having minimal error | |
| Zhang et al. | Range image registration via probability field | |
| Soucy | View correspondence using curvature and motion consistency |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040330 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040330 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060221 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060424 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060523 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060531 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090609 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100609 Year of fee payment: 4 |
|
| S111 | Request for change of ownership or part of ownership |
Free format text: JAPANESE INTERMEDIATE CODE: R313113 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100609 Year of fee payment: 4 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100609 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110609 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110609 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120609 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120609 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130609 Year of fee payment: 7 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| LAPS | Cancellation because of no payment of annual fees |