(57)【要約】【課題】カメラの内部パラメータが所与の場合に、４つの特徴点の座標情報に基づいてカメラの外部パラメータの推定を行うカメラキャリブレーション手法を提供する。【解決手段】ワールド座標系とカメラ座標系が同一の場合に、ワールド座標系における４つの特徴点が作る２つの角度を基にして、カメラ座標系での各特徴点の座標位置を推定する。視線ベクトルＡを延ばした直線で形成される距離Ｔを８分割して、８分割した各地点のうち特徴点Ａが存在する座標位置に最も近い地点がいずれかを推定する。この推定は、８分割した各地点を基点にして視線ベクトルＢ、Ｃ、Ｄ上の点のうち辺ＡＢ、辺ＢＣ、辺ＤＡの長さと等しくなる点Ｂ´、Ｃ´、Ｄ´をそれぞれ求め、該基点、点Ｂ´、Ｃ´、Ｄ´で作られる角度Ｂ´、角度Ｄ´の角度Ｂ、角度Ｄとの差異を各地点で求め、特徴点Ａ、Ｂ、Ｃ、Ｄに最も近似する８分割した地点のうちの一つの点、およびその点に対応するＢ´、Ｃ´、Ｄ´を求める。【選択図】図４
2010239515.tif 000002
【特許請求の範囲】
【請求項１】
ワールド座標系（world coordinate system）とカメラ座標系（camera coordinate system）の原点、Ｘ軸、Ｙ軸およびＺ軸が同一の場合に、ワールド座標系における４つの特徴点が作る２つの角度を基にして、カメラ座標系での各特徴点の座標位置を推定することによってカメラのキャリブレーションを行う方法であって、
コンピュータが実行する前記方法は、
（ａ）ワールド座標系における４つの特徴点Ａ、Ｂ、Ｃ、Ｄを選択するステップ（４０１）と、
（ｂ）前記選択した特徴点を結んだ辺ＡＢ、辺ＢＣ、辺ＣＤ、辺ＤＡの長さを計算するステップ（４０２）と、
（ｃ）前記辺ＡＢおよび辺ＢＣが作る角度Ｂ、並びに前記辺ＣＤおよび辺ＤＡが作る角度Ｄを計算するステップ（４０３）であって、角度Ｂおよび角度Ｄは０度より大きく９０度以下であるステップと、
（ｄ）ワールド座標系での辺ＡＢをＸ軸とする場合の回転行列Ｒを計算するステップ（４０４）と、
（ｅ）前記原点と、前記４つの特徴点Ａ、Ｂ、Ｃ、Ｄをスクリーン座標系に変換した変換後の各点とを結ぶ４つの視線ベクトルＡ、Ｂ、Ｃ、Ｄを計算するステップ（４０５）と、
（ｆ）前記視線ベクトルＡを延ばした直線で形成される距離Ｔを８分割して、該８分割した各地点のうち前記特徴点Ａが存在する座標位置に最も近い地点がいずれかを推定するステップ（４０６）であって、該推定は、前記８分割した各地点のうち原点に最も近い地点から始め、各地点を基点にして視線ベクトルＢ、Ｄ上の点のうち前記辺ＡＢ、辺ＤＡの長さと等しくなる点Ｂ´、Ｄ´をそれぞれ求め、点Ｂ´と視線ベクトルＣ上の点を結ぶ点のうち前記辺ＢＣの長さと等しくなる点Ｃ´を求め、該基点、点Ｂ´、Ｃ´、Ｄ´で作られる角度Ｂ´、角度Ｄ´の前記角度Ｂ、角度Ｄとの差異を前記８分割した各地点で求め、特徴点Ａ、Ｂ、Ｃ、Ｄに最も近似する８分割した地点のうちの一つの基点、および該一つの基点に対応するＢ´、Ｃ´、Ｄ´を求めるステップと、
（ｇ）前記距離Ｔを２分の１に収束した距離を求めるステップであって、前記収束された新たな距離Ｔは、前記推定された地点から距離Ｔの４分の１の距離だけ原点方向および視線ベクトルＡを延ばした方向に離れた地点で形成されるステップ（４０７）と、
（ｈ）前記新たな距離Ｔが予め定められた値を下回る（４０８）まで、ステップ（ｆ）および（ｇ）を繰返すステップと
を含み、
前記新たな距離Ｔが予め定められた値を下回った後に、
（ｉ）ワールド座標系からカメラ座標系への変換行列Ｎを求めるステップ（４０９）であって、該変換行列Ｎは、前記回転行列Ｒおよびカメラの姿勢と位置の情報に基づいて計算されるステップと、
（ｊ）前記変換行列Ｎに基づいて、カメラの外部パラメータであるカメラの回転行列およびカメラの移動ベクトルを計算するステップ（４１０）と、
（ｋ）前記カメラの外部パラメータの誤差補正を行うステップ（４１１）と
をさらに含むことを特徴とする方法。
【請求項２】
前記ステップ（ｆ）は、
（ｆ１）前記８分割された視線ベクトルＡ上の地点の一つの基点を点ｖＡとするステップ（８０１）と、
（ｆ２）前記点ｖＡからの長さが前記辺ＡＢの長さと等しくなる視線ベクトルＢ上の１つまたは複数の点Ｂ´を求めるステップ（８０３）と、
（ｆ３）前記点ｖＡからの長さが前記辺ＤＡの長さと等しくなる視線ベクトルＤ上の１つまたは複数の点Ｄ´を求めるステップ（８０４）と、
（ｆ４）前記点Ｂ´からの長さが前記辺ＢＣの長さと等しくなる視線ベクトルＣ上の１つまたは複数の点Ｃ´を求めるステップ（８０５）と、
（ｆ５）前記点ｖＡ、点Ｂ´、点Ｃ´によって作られる１つまたは複数の角度Ｂ´を求め、前記角度Ｂに最も近似する前記点ｖＡ、点Ｂ´、点Ｃ´の組み合わせを決定するステップであって、該決定された組み合わせでの角度Ｂ´と前記角度Ｂの差異を第１の差異値として記憶しておくステップ（８０６）と、
（ｆ６）前記点ｖＡ、点Ｄ´、前記一つに決定された点Ｃ´によって作られる１つまたは複数の角度Ｄ´を求め、前記角度Ｄに最も近似する前記点ｖＡ、該点Ｃ´、点Ｄ´の組み合わせを決定するステップであって、該決定された組み合わせでの角度Ｄ´と前記角度Ｄの差異である第２の差異値を前記第１の差異値に加算して合計値を算出するステップ（８０７）と、
（ｆ７）前記一つに決定された点ｖＡ、点Ｂ´、点Ｃ´、点Ｄ´の組み合わせ、並びに前記合計値を格納するステップと
をさらに含むことを特徴とする請求項１に記載の方法。
【請求項３】
前記合計値が値ゼロに等しい場合、前記一つに決定された点ｖＡ、点Ｂ´、点Ｃ´、点Ｄ´が４つの特徴点Ａ、Ｂ、Ｃ、Ｄにそれぞれ対応することを特徴とする請求項２に記載の方法。
【請求項４】
ワールド座標系（world coordinate system）とカメラ座標系（camera coordinate system）の原点、Ｘ軸、Ｙ軸およびＺ軸が同一の場合に、ワールド座標系における４つの特徴点が作る２つの角度を基にして、カメラ座標系での各特徴点の座標位置を推定する方法であって、
コンピュータが実行する前記方法は、
（ａ）ワールド座標系における４つの特徴点Ａ、Ｂ、Ｃ、Ｄを選択するステップ（４０１）と、
（ｂ）前記選択した特徴点を結んだ辺ＡＢ、辺ＢＣ、辺ＣＤ、辺ＤＡの長さを計算するステップ（４０２）と、
（ｃ）前記辺ＡＢおよび辺ＢＣが作る角度Ｂ、並びに前記辺ＣＤおよび辺ＤＡが作る角度Ｄを計算するステップ（４０３）であって、角度Ｂおよび角度Ｄは０度より大きく９０度以下であるステップと、
（ｄ）前記原点と、前記４つの特徴点Ａ、Ｂ、Ｃ、Ｄをスクリーン座標系に変換した変換後の各点とを結ぶ４つの視線ベクトルＡ、Ｂ、Ｃ、Ｄを計算するステップ（４０５）と、
（ｅ）前記視線ベクトルＡを延ばした直線で形成される距離Ｔを８分割して、該８分割した各地点のうち前記特徴点Ａが存在する座標位置に最も近い地点がいずれかを推定するステップ（４０６）であって、該推定は、前記８分割した各地点のうち原点に最も近い地点から始め、各地点を基点にして視線ベクトルＢ、Ｄ上の点のうち前記辺ＡＢ、辺ＤＡの長さと等しくなる点Ｂ´、Ｄ´をそれぞれ求め、点Ｂ´と視線ベクトルＣ上の点を結ぶ点のうち前記辺ＢＣの長さと等しくなる点Ｃ´を求め、該基点、点Ｂ´、Ｃ´、Ｄ´で作られる角度Ｂ´、角度Ｄ´の前記角度Ｂ、角度Ｄとの差異を前記８分割した各地点で求め、特徴点Ａ、Ｂ、Ｃ、Ｄに最も近似する８分割した地点のうちの一つの基点、および該一つの基点に対応するＢ´、Ｃ´、Ｄ´を求めるステップと、
（ｆ）前記距離Ｔを２分の１に収束した距離を求めるステップであって、前記収束された新たな距離Ｔは、前記推定された地点から距離Ｔの４分の１の距離だけ原点方向および視線ベクトルＡを延ばした方向に離れた地点で形成されるステップ（４０７）と、
（ｇ）前記新たな距離Ｔが予め定められた値を下回る（４０８）まで、ステップ（ｅ）および（ｆ）を繰返すステップと
を含むことを特徴とする方法。
【請求項５】
前記ステップ（ｅ）は、
（ｅ１）前記８分割された視線ベクトルＡ上の地点の一つの基点を点ｖＡとするステップ（８０１）と、
（ｅ２）前記点ｖＡからの長さが前記辺ＡＢの長さと等しくなる視線ベクトルＢ上の１つまたは複数の点Ｂ´を求めるステップ（８０３）と、
（ｅ３）前記点ｖＡからの長さが前記辺ＤＡの長さと等しくなる視線ベクトルＤ上の１つまたは複数の点Ｄ´を求めるステップ（８０４）と、
（ｅ４）前記点Ｂ´からの長さが前記辺ＢＣの長さと等しくなる視線ベクトルＣ上の１つまたは複数の点Ｃ´を求めるステップ（８０５）と、
（ｅ５）前記点ｖＡ、点Ｂ´、点Ｃ´によって作られる１つまたは複数の角度Ｂ´を求め、前記角度Ｂに最も近似する前記点ｖＡ、点Ｂ´、点Ｃ´の組み合わせを決定するステップであって、該決定された組み合わせでの角度Ｂ´と前記角度Ｂの差異を第１の差異値として記憶しておくステップ（８０６）と、
（ｅ６）前記点ｖＡ、点Ｄ´、前記一つに決定された点Ｃ´によって作られる１つまたは複数の角度Ｄ´を求め、前記角度Ｄに最も近似する前記点ｖＡ、該点Ｃ´、点Ｄ´の組み合わせを決定するステップであって、該決定された組み合わせでの角度Ｄ´と前記角度Ｄの差異である第２の差異値を前記第１の差異値に加算して合計値を算出するステップ（８０７）と、
（ｅ７）前記一つに決定された点ｖＡ、点Ｂ´、点Ｃ´、点Ｄ´の組み合わせ、並びに前記合計値を格納するステップと
をさらに含むことを特徴とする請求項４に記載の方法。
【請求項６】
前記合計値が値ゼロに等しい場合、前記一つに決定された点ｖＡ、点Ｂ´、点Ｃ´、点Ｄ´が４つの特徴点Ａ、Ｂ、Ｃ、Ｄにそれぞれ対応することを特徴とする請求項５に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、カメラキャリブレーションの算出方法に関する。より詳細には、本発明は、内部パラメータが既知の場合に、少なくとも４つの特徴点が存在すれば、同一平面上、非同一平面上の区別なくカメラの外部パラメータの推定を行うことができる方法に関する。
【背景技術】
【０００２】
まず、カメラキャリブレーションの概要および既存のカメラキャリブレーションについて説明する。
【０００３】
まず、３次元空間を「ワールド座標」とし、そのワールド座標上に絶対的な空間上の位置ＸＹＺが複数存在し、その位置をスクリーンに投影したスクリーン座標での位置ＵＶが存在するとする。ワールド座標上の位置ＸＹＺとスクリーン座標での位置ＵＶはそれぞれ１対１で対応しており、これらが複数（少なくとも５または６個）存在すれば、カメラの姿勢（３×３行列）と移動ベクトル（ワールド座標上でのカメラ位置）を推定することができる。ここで、カメラの姿勢と移動ベクトルは「外部パラメータ」と呼ばれる。また、ワールド座標上の位置ＸＹＺとしてサンプリングしたそれぞれの点を「特徴点」と呼ぶことにする。特徴点は、カメラ情報を推定するために使用されるが、特徴点は、ワールド座標において同一平面上にある場合もあれば、非同一平面上にある場合もある。
【０００４】
図１は、従来の、カメラから見たときのスクリーン座標での特徴点位置を示す図であり、図１の座標は右手系を示し、スクリーン手前に向かってＺ軸が伸びている。
【０００５】
図１では、カメラから見たときのスクリーン座標上での特徴点が８つ与えられている。この例において、８つの特徴点は、それぞれワールド座標上の位置ＸＹＺとスクリーン上に投影された位置ＵＶが確定しているものとする。カメラの手前に投影面があり、投影された座標値を（U1, V1）、（U2, V2）、〜（U8, V8）とし、対応するワールド座標を（Wx1, Wy1, Wz1）、（Wx2, Wy2, Wz2）、〜（Wx8, Wy8, Wz8）とする。
【０００６】
また、ワールド座標からスクリーン座標への変換は、
ワールド座標  →  ビュー座標  →  透視座標  →  スクリーン座標
の順番にて行われる。座標間のそれぞれの変換は、４×４行列にて表すことができる。ここで、ビュー座標は、視点を中心とした座標系のことであり、透視座標は、３次元空間を２次元空間に投影する座標系のことである。本明細書の説明において、ビュー座標から透視座標、スクリーン座標に変換する行列はあらかじめ算出できているものとする。
【０００７】
特に、ビュー座標から透視座標変換に必要な要素としては、スクリーンサイズ、カメラの視野角度、カメラから投影面までの距離、ゆがみ係数、が存在する。これらはカメラの「内部パラメータ」と呼ばれる。
【０００８】
このうち、ワールド座標からビュー座標に変換する行列を求める手法をカメラキャリブレーションと呼ぶ（厳密に言うと、内部パラメータを推定する部分もカメラキャリブレーションに含まれるが、今回の発明ではこの部分は既に計算できているものとして考えるため、本明細書の説明においてはその部分を省略する）。
【０００９】
ワールド座標上のカメラの外部パラメータとしては３×３の姿勢行列（回転行列）と移動ベクトル（ワールド座標上の位置）が存在する。カメラでの姿勢をあらわすＸＹＺ軸を、CX(CX.x, CX.y, CX.z)、CY(CYx, CY.y, CY.z)、CZ(CZ.x, CZ.y, CZ.z)の正規直交基底として表した場合、以下の３×３行列にて回転行列を表すことができる。
【００１０】
【数１】
2010239515.tif 000003
【００１１】
また、カメラの移動ベクトルをV(Vx, Vy, Vz)とすると、カメラの姿勢と移動ベクトルより、ワールド座標からビュー座標の変換行列は以下の４×４行列にてあらわすことができる。
【００１２】
【数２】
2010239515.tif 000004
【００１３】
図２は、上記で説明したワールド座標上の８つの特徴点とカメラの関係を示す図である。図２では、カメラの前に投影面（スクリーン座標）があり、８つの特徴点を上述した変換によってスクリーン座標に投影した図が示されている。
【００１４】
カメラキャリブレーションを求める手法としては、従来から様々な手法（例えば、ＤＬＴ、Ｔｓａｉ、Ｚｈａｎｇ、マーカーを使用したカメラキャリブレーション）が知られている。
【００１５】
ＤＬＴ手法（非特許文献１参照）は、三次元空間上の特徴点とスクリーンに投影された点を６つ以上与えることにより、内部パラメータ・外部パラメータの算出を行い、カメラキャリブレーションを行う手法である。ＤＬＴ手法では、同一平面上に特徴点が存在する場合は正しい推定を行うことができず、また、計算後に視線方向（光軸方向）に誤差が出る可能性が高いため、補正を行う必要がある。
【００１６】
Ｔｓａｉ手法（非特許文献２参照）は、三次元空間上の特徴点とスクリーンに投影された点を６つ以上与えることにより、内部パラメータ・外部パラメータの算出を行い、カメラキャリブレーションを行う手法である。Ｔｓａｉ手法では、同一平面上に特徴点が存在する場合に使用できる。また、非同一平面上に特徴点がある場合でも推定できるが、その場合には同一平面上にしか特徴点がない場合に比べて大きな誤差が出てしまう。このため、同一平面上か非同一平面上のいずれに特徴点が存在するか、をあらかじめ判定しておく必要があり、また、計算結果の精度はそれほど高くはなく、計算後に視線方向（光軸方向）に誤差が出る可能性が高いため、補正を行う必要がある。
【００１７】
Ｚｈａｎｇ手法は、同一平面上の格子状のマーカーよりカメラキャリブレーションを行い、内部パラメータ・外部パラメータの算出を行う手法である。Ｚｈａｎｇ手法は、カメラキャリブレーションの精度が高いが、同一平面上に特徴点が格子状に存在している必要がある。また、カメラの内部パラメータの推定にて利用される。
【００１８】
また、マーカーを使用したカメラキャリブレーションの手法は、特定の条件を満たす左右非対称の模様の書かれた正方形のマーカーを登録してパターンマッチングによりマーカーを識別し、これよりカメラキャリブレーションを行う手法である。この手法は、同一平面限定で適用可能であるという制限があり、また、あらかじめマーカーを登録しておく必要がある。
【先行技術文献】
【非特許文献】
【００１９】
【非特許文献１】DLT Method（Young-Hoo Kwon 著）、http://www.kwon3d.com/theory/dlt/dlt.html
【非特許文献２】Tsai's Camera Model（Steven Shafer, Takeo Kanade 著）、http://www-2.cs.cmu.edu/~rgw/TsaiDesc.html
【非特許文献３】マーカー追跡に基づく拡張現実感システムとそのキャリブレーション（加藤博一，Mark Billinghurst，浅野浩一，橘啓八郎 著）、http://intron.kz.tsukuba.ac.jp/tvrsj/4.4/kato/p-99_VRSJ4_4.pdf
【非特許文献４】任意視点画像に適応したカメラキャリブレーションの研究（久保智裕 著）、http://www433.elec.ryukoku.ac.jp/Houbun/YOUSHI/2003/2003_master_thesis_kubo.pdf
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかしながら、従来技術は、特徴点が同一平面上に存在する場合のみ適用可能なもの、非同一平面上に存在する場合のみ適用可能なもの、あるいは両方の場合に適用できるが上述したような不都合も同時に内包しているという課題がある。
【課題を解決するための手段】
【００２１】
このような課題を解決するため、本発明の一実施形態に係る方法は、ワールド座標系（world coordinate system）とカメラ座標系（camera coordinate system）の原点、Ｘ軸、Ｙ軸およびＺ軸が同一の場合に、ワールド座標系における４つの特徴点が作る２つの角度を基にして、カメラ座標系での各特徴点の座標位置を推定することによってカメラのキャリブレーションを行う方法であることを特徴とする。この方法はコンピュータによって実行される。また、この方法は、（ａ）ワールド座標系における４つの特徴点Ａ、Ｂ、Ｃ、Ｄを選択するステップ（４０１）と、（ｂ）選択した特徴点を結んだ辺ＡＢ、辺ＢＣ、辺ＣＤ、辺ＤＡの長さを計算するステップ（４０２）と、（ｃ）辺ＡＢおよび辺ＢＣが作る角度Ｂ、並びに辺ＣＤおよび辺ＤＡが作る角度Ｄを計算するステップ（４０３）であって、角度Ｂおよび角度Ｄは０度より大きく９０度以下であるステップと、（ｄ）ワールド座標系での辺ＡＢをＸ軸とする場合の回転行列Ｒを計算するステップ（４０４）と、（ｅ）原点と、４つの特徴点Ａ、Ｂ、Ｃ、Ｄをスクリーン座標系に変換した変換後の各点とを結ぶ４つの視線ベクトルＡ、Ｂ、Ｃ、Ｄを計算するステップ（４０５）と、（ｆ）視線ベクトルＡを延ばした直線で形成される距離Ｔを８分割して、該８分割した各地点のうち特徴点Ａが存在する座標位置に最も近い地点がいずれかを推定するステップ（４０６）であって、該推定は、８分割した各地点のうち原点に最も近い地点から始め、各地点を基点にして視線ベクトルＢ、Ｄ上の点のうち辺ＡＢ、辺ＤＡの長さと等しくなる点Ｂ´、Ｄ´をそれぞれ求め、点Ｂ´と視線ベクトルＣ上の点を結ぶ点のうち辺ＢＣの長さと等しくなる点Ｃ´を求め、該基点、点Ｂ´、Ｃ´、Ｄ´で作られる角度Ｂ´、角度Ｄ´の角度Ｂ、角度Ｄとの差異を８分割した各地点で求め、特徴点Ａ、Ｂ、Ｃ、Ｄに最も近似する８分割した地点のうちの一つの基点、および該一つの基点に対応するＢ´、Ｃ´、Ｄ´を求めるステップと、（ｇ）距離Ｔを２分の１に収束した距離を求めるステップであって、収束された新たな距離Ｔは、推定された地点から距離Ｔの４分の１の距離だけ原点方向および視線ベクトルＡを延ばした方向に離れた地点で形成されるステップ（４０７）と、（ｈ）新たな距離Ｔが予め定められた値を下回る（４０８）まで、ステップ（ｆ）および（ｇ）を繰返すステップと、を含むことを特徴とする。さらに、本発明に係る方法は、新たな距離Ｔが予め定められた値を下回った後に、（ｉ）ワールド座標系からカメラ座標系への変換行列Ｎを求めるステップ（４０９）であって、該変換行列Ｎは、回転行列Ｒおよびカメラの姿勢と位置の情報に基づいて計算されるステップと、（ｊ）変換行列Ｎに基づいて、カメラの外部パラメータであるカメラの回転行列およびカメラの移動ベクトルを計算するステップ（４１０）と、（ｋ）カメラの外部パラメータの誤差補正を行うステップ（４１１）と、を含むことを特徴とする。
【００２２】
また、他の実施形態に係る方法は、ワールド座標系とカメラ座標系の原点、Ｘ軸、Ｙ軸およびＺ軸が同一の場合に、ワールド座標系における４つの特徴点が作る２つの角度を基にして、カメラ座標系での各特徴点の座標位置を推定する方法であることを特徴とする。この方法はコンピュータによって実行される。また、この方法は、（ａ）ワールド座標系における４つの特徴点Ａ、Ｂ、Ｃ、Ｄを選択するステップ（４０１）と、（ｂ）選択した特徴点を結んだ辺ＡＢ、辺ＢＣ、辺ＣＤ、辺ＤＡの長さを計算するステップ（４０２）と、（ｃ）辺ＡＢおよび辺ＢＣが作る角度Ｂ、並びに辺ＣＤおよび辺ＤＡが作る角度Ｄを計算するステップ（４０３）であって、角度Ｂおよび角度Ｄは０度より大きく９０度以下であるステップと、（ｄ）原点と、４つの特徴点Ａ、Ｂ、Ｃ、Ｄをスクリーン座標系に変換した変換後の各点とを結ぶ４つの視線ベクトルＡ、Ｂ、Ｃ、Ｄを計算するステップ（４０５）と、（ｅ）視線ベクトルＡを延ばした直線で形成される距離Ｔを８分割して、該８分割した各地点のうち特徴点Ａが存在する座標位置に最も近い地点がいずれかを推定するステップ（４０６）であって、該推定は、８分割した各地点のうち原点に最も近い地点から始め、各地点を基点にして視線ベクトルＢ、Ｄ上の点のうち辺ＡＢ、辺ＤＡの長さと等しくなる点Ｂ´、Ｄ´をそれぞれ求め、点Ｂ´と視線ベクトルＣ上の点を結ぶ点のうち辺ＢＣの長さと等しくなる点Ｃ´を求め、該基点、点Ｂ´、Ｃ´、Ｄ´で作られる角度Ｂ´、角度Ｄ´の角度Ｂ、角度Ｄとの差異を前記８分割した各地点で求め、特徴点Ａ、Ｂ、Ｃ、Ｄに最も近似する８分割した地点のうちの一つの基点、および該一つの基点に対応するＢ´、Ｃ´、Ｄ´を求めるステップと、（ｆ）距離Ｔを２分の１に収束した距離を求めるステップであって、収束された新たな距離Ｔは、推定された地点から距離Ｔの４分の１の距離だけ原点方向および視線ベクトルＡを延ばした方向に離れた地点で形成されるステップ（４０７）と、（ｇ）新たな距離Ｔが予め定められた値を下回る（４０８）まで、ステップ（ｅ）および（ｆ）を繰返すステップと、を含むことを特徴とする。
【発明の効果】
【００２３】
本発明に係るカメラキャリブレーション手法によれば、内部パラメータは既に計算済みであるという前提条件に立つと、少なくとも４つの特徴点が存在すれば、同一平面上、非同一平面上の区別なくカメラの外部パラメータの推定を行うことができる。また、４つの特徴点のみでカメラの外部パラメータを推定するため、従来の計算負荷を大幅に軽減できる。さらに、既存の特徴点情報にて補正を行うため、計算結果の精度も確保することができる。
【図面の簡単な説明】
【００２４】
【図１】従来技術の、カメラから見たときのスクリーン座標での特徴点位置を示す図である。
【図２】従来技術の、ワールド座標上の８つの特徴点とカメラの関係を示す図である。
【図３】本発明のカメラキャリブレーション算出方法を行うためのシステムのハードウェア構成の一例を示すブロック図である。
【図４】本発明のカメラキャリブレーション算出方法の全体的な流れを示す図である。
【図５】ワールド座標上の４つの特徴点をランダムに選択した例を示す図である。
【図６】視点位置から４つのスクリーン座標上の特徴点を結ぶ視線ベクトルの概念を示す図である。
【図７】Ａの視線ベクトルを伸ばして特徴点Ａのワールド座標位置を推定する概念を示す図である。
【図８】任意の特徴点のワールド座標上の位置を推定する処理の詳細を流れ図である。
【図９】点Ｂ´から直線Ｃに下ろした直線のうちワールド座標でのＢＣの長さと同一のものを計算するための概念図である。
【図１０】minT、maxTを更新するための流れ図である。
【図１１】カメラキャリブレーションにおいて解が２つ以上存在するパターンを示す概念図である。
【図１２】カメラの外部パラメータの誤差を補正する処理を加えたカメラキャリブレーションを示す流れ図である。
【図１３】カメラの外部パラメータの誤差を補正する処理を加えたカメラキャリブレーションを示す流れ図である。
【発明を実施するための形態】
【００２５】
図３は、本発明のカメラキャリブレーション算出方法を行うためのシステムのハードウェア構成の一例を示すブロック図である。このシステムは、装置３００およびカメラ３０１を含んでおり、両者は有線／無線媒体３０２を介して接続されている。
【００２６】
カメラ３０１は、ＣＣＤ（Charge Coupled Devices）やＣＭＯＳ（Complementary Metal Oxide Semiconductor）などを搭載した一般的なカメラで構わない。また、カメラ３０１は、撮像した映像データを一時的に記憶しておくためのメモリ（例えば、フラッシュメモリ）を内蔵し、その映像データを装置３００に送信する機能を有する。
【００２７】
装置３００は、制御部３１１、主記憶部３１２、補助記憶部３１３、入出力インターフェース（Ｉ／Ｆ）部３１４、及び出力表示部３１５を含み、これらの構成要素はバス３１６によって相互に接続されており、制御部３１１によって各構成要素３１２〜３１５は統括して制御される。なお、この装置３００は、一般的なパーソナルコンピュータ（ＰＣ）などのコンピュータであって構わない。
【００２８】
制御部３１１は、中央処理装置（ＣＰＵ）を含んでおり、上述したように各構成要素３１２〜３１５の制御やデータの演算を行い、また、補助記憶部３１３に格納されている各種プログラムを主記憶部３１２に読み出して実行する。なお、制御部３１１は、３Ｄグラフィックスの表示を行うために必要な演算処理を行うＧＰＵ（Graphics Processing Unit）を有していても良い。
【００２９】
主記憶部３１２は、メインメモリとも呼ばれ、装置３００に入力されたデータ、コンピュータ実行可能な命令及び当該命令による演算処理後のデータなどを記憶する。補助記憶部３１３は、ハードディスク（ＨＤＤ）などに代表される記憶装置であり、データやプログラムを長期的に保存する際に使用される。つまり、制御部３１１がプログラムを実行してデータの演算を行う場合には、補助記憶部３１３から必要なデータだけ主記憶部３１２に読み出し、演算結果のデータを長期的に記憶・保存するには制御部３１１が補助記憶部３１３に演算結果のデータを書き込む。
【００３０】
入出力インターフェース（Ｉ／Ｆ）部３１４は、カメラ３０１との間でデータ送受信を行うインターフェースであり、有線／無線媒体３０２に接続する。出力表示部３１５は、カメラ３０１によって撮像した映像データを表示するディスプレイなどの表示装置であり、本発明に係るカメラキャリブレーションの詳細処理も同時に表示可能である。
【００３１】
図４は、本発明に係るカメラキャリブレーション手法の全体的な流れを示す図である。以後、図４を参照しながら、本発明の実施形態に係るカメラキャリブレーション手法を使用して、カメラの姿勢（３×３の回転行列）と移動ベクトル（ワールド座標上の位置）を計算する手順を詳細に説明する。以下の説明では、ｎ個のサンプリングポイントが存在する場合、ワールド座標上の特徴点位置を（Wx1，Wy1，Wz1）、（Wx2，Wy2，Wz2）、〜（Wxn，Wyn，Wzn）とし、対応するスクリーン上の位置を（U1，V1）、（U2，V2）、〜（Un，Vn）とする。
【００３２】
ステップＳ４０１では、ワールド座標における４つの特徴点Ａ、Ｂ、Ｃ、Ｄを、以下の条件、
０°＜  ∠Ｂ  ≦９０°
０°＜  ∠Ｄ  ≦９０°
を満たすようにランダムに選択する。図５ａ、ｂは、ワールド座標上の４つの特徴点Ａ、Ｂ、Ｃ、Ｄをランダムに選択した例を示すが、図５ｂのように、４つの特徴点Ａ、Ｂ、Ｃ、Ｄを結んだそれぞれの辺が交差していてもかまわない。
【００３３】
この場合の４組分のワールド座標値とスクリーン座標値が、以下に説明するカメラキャリブレーション算出手法にて使用される。４つの特徴点Ａ、Ｂ、Ｃ、Ｄのワールド座標上の空間の位置ＸＹＺは、それぞれ、w[0]、w[1]、w[2]、w[3]で表され、対応するスクリーン座標上の位置ＸＹは、それぞれ、s[0]、s[1]、s[2]、s[3]で表されるものとする。
【００３４】
ステップＳ４０２では、ワールド座標上での各辺の長さを計算する。ＡＢ、ＢＣ、ＣＤ、ＤＡの長さは、それぞれ以下の計算式のようになる。
【００３５】
【数３】
2010239515.tif 000005
【００３６】
【数４】
2010239515.tif 000006
【００３７】
【数５】
2010239515.tif 000007
【００３８】
【数６】
2010239515.tif 000008
【００３９】
ステップＳ４０３では、ワールド座標上の特徴点Ｂ、Ｄにおける角度を計算する。∠Ｂを求めるには、まず、辺ＡＢ、ＢＣを正規化した上で、内積を計算する。そうすることで∠Ｂをｃｏｓ値で求めることができる。
【００４０】
【数７】
2010239515.tif 000009
【００４１】
【数８】
2010239515.tif 000010
【００４２】
【数９】
2010239515.tif 000011
【００４３】
同様にして、辺ＡＤ、ＣＤを正規化した上で、内積を計算する。そして、∠Ｄのｃｏｓ値を求める。
【００４４】
【数１０】
2010239515.tif 000012
【００４５】
【数１１】
2010239515.tif 000013
【００４６】
【数１２】
2010239515.tif 000014
【００４７】
ステップＳ４０４では、ワールド座標上の辺ＡＢをＸ軸方向と仮定したときの回転行列Ｒを計算する。以下の式で「cross」はベクトルの外積を求める関数である。  
【００４８】
【数１３】
2010239515.tif 000015
【００４９】
その上で、ｖＸ，ｖＹ，ｖＺをベクトルを正規化する関数「normalize」を使用して、
【００５０】
【数１４】
2010239515.tif 000016
【００５１】
とし、３×３の回転行列Ｒを以下のように計算する。
【００５２】
【数１５】
2010239515.tif 000017
【００５３】
ステップＳ４０５では、ビュー座標とワールド座標が等しいものとして、かつビュー座標を姿勢ベクトル＝単位行列、視点位置＝（0, 0, 0）として、視点位置（0, 0, 0）から４つのスクリーン座標上の特徴点を結ぶ視線ベクトルを計算する。図６はステップＳ４０５の処理内容の概念図を表している。本ステップでは、既知のカメラ内部パラメータより、透視座標からビュー座標への変換行列（ビュー座標から透視座標に変換する行列の逆行列）をあらかじめ求めておき、このときの方向ベクトルを正規化することで視線方向のベクトルを求めることができる。それぞれの正規化されたベクトルをvSDir[0]、vSDir[1]、vSDir[2]、vSDir[3]とする。それぞれのベクトルを伸ばした位置に、対象となる特徴点Ａ、Ｂ、Ｃ、Ｄのワールド座標値が存在することになる。
【００５４】
ステップＳ４０６では、４つの特徴点のうちひとつを選択し、その特徴点に対応する視線ベクトルを延長した先のどの位置にワールド座標上の特徴点が存在するかを推定する。ここでは、まず特徴点Ａのワールド座標上での位置を求めることにする。スクリーン座標でのＡの視線ベクトルvSDir[0]を伸ばした先には必ず特徴点Ａのワールド座標位置が存在する。本発明に係る方法では、以下に説明するように、Ａの視線ベクトルを十分な距離伸ばしてその距離を８分割することを繰り返しながら、最適な位置を収束させる。また、走査する視点からの距離をminTからmaxTまで変化させるとし、例えば、初期値は、
minT = 5.0
maxT = 99999.0
などの十分小さい値と十分大きい値を設定する。そして、
dT = (maxT - minT) / 8.0
とし、視点位置（0, 0, 0）からの距離ｔは、minTから始まり、dTずつ増加させてmaxTまで走査していくものとする。
【００５５】
次に、図８を参照しながら、ステップＳ４０６で行われる処理の詳細を説明する。
【００５６】
ステップＳ８０１では、視線ベクトルＡを伸ばした先のワールド座標上の仮の位置を計算する。スクリーン座標上の点Ａを伸ばした先の仮のワールド座標上の点は、仮のワールド座標上の点の位置をｖＡとし、視点からの距離をｔとした場合に、以下の式によって計算される。
vA = vSDir[0] * t
上述したように、距離ｔは、minTから始まるので、t = minTとする。
【００５７】
ステップＳ８０２では、ｔがmaxTであるか否かを判定する。ｔがmaxTに到達した場合には、本ステップの処理を終了する。
【００５８】
ステップＳ８０３では、点ｖＡから視線ベクトルＢ上に、距離ＡＢの拘束条件で下ろした位置を算出する。図７に示したように、視線ベクトルＢ（vSDir[1]）で表される直線に対して点ｖＡから下ろした垂線の足をｐとする。また、点ｖＡから垂線の足ｐまでの距離をｈとし、ステップＳ４０２において既に求めているワールド座標上のＡＢの距離と距離ｈとを比較する。
【００５９】
ｈ＞ＡＢの関係が成立する場合は、視点位置（0, 0, 0）から見た仮のワールド座標上の点の位置ｖＡが実際のワールド座標上の点Ａの位置よりも遠くにあることを意味するので、この段階で計算を中断し、ステップＳ４０７での収束計算にてmaxTの値をさらに小さいものにしてからステップＳ８０１の処理からやり直す。
【００６０】
ｈ＝ＡＢの場合は、ｖＡがワールド座標上の特徴点Ａの位置であり、かつ垂線の足ｐがワールド座標上の特徴点Ｂの位置と仮定できる。
【００６１】
また、ｈ＜ＡＢの場合は、直線Ｂ上に仮の点ｖＡからの距離がＡＢとなる位置を計算して当てはめる。
【００６２】
【数１６】
2010239515.tif 000018
【００６３】
計算された直線Ｂ上の点Ｂ´１とｖＡの距離、およびＢ´２とｖＡの距離は、共にＡＢの長さと等しくなる。
【００６４】
ステップＳ８０４では、点ｖＡから視線ベクトルＤ上に、距離ＡＤの拘束条件で下ろした位置を算出する。ステップＳ８０３と同様に、視線ベクトルＤ（vSDir[3]）で表される直線に対して点ｖＡから下ろした垂線の足をｐ１とする。また、点ｖＡから垂線の足ｐ１までの距離をｈ１として、ＡＤの距離と距離ｈ１とを比較する。ｈ１＜ＡＤの場合は、直線Ｄ上にてｖＡとの距離がＡＤとなる位置を計算し、以下の式のように、Ｄ´１、Ｄ´２とする。
【００６５】
【数１７】
2010239515.tif 000019
【００６６】
ステップＳ８０５では、点Ｂ´から視線ベクトルＣ上に、距離ＢＣの拘束条件で下ろした位置を算出する。点Ｂ´は、ステップＳ８０３において求めたＢ´１またはＢ´２である。
【００６７】
図９に示したように、点Ｂ´から視線ベクトルＣ（vSDir[2]）で表される直線に下ろした垂線の足をｐ２とする。点Ｂ´とｐ２間の距離ｈ２と、ステップＳ４０２において既に求めているワールド座標上のＢＣの距離とを対比する。なお、本ステップでは、ｈ２がＢＣの距離より小さい値であることが条件となる（大きい値の場合は、ステップＳ８０３で既に処理が中断している）。
【００６８】
ｈ２＝ＢＣとなる場合は、垂線の足ｐ２がワールド座標上の点Ｃの位置と仮定できる。一方、ｈ２＜ＢＣの場合は、直線Ｃ上で点Ｂ´からの距離がＢＣとなる位置を計算して当てはめる。
【００６９】
【数１８】
2010239515.tif 000020
【００７０】
計算された直線Ｃ上の点Ｃ´１とＢ´の距離、およびＣ´２とＢ´の距離は、共にＢＣの長さと等しくなる。
【００７１】
ステップＳ８０６では、上記ステップまでで求めた仮の点ＡＢＣ（ｖＡ、Ｂ´、Ｃ´）が作る角度を計算する。上記のステップまでで以下の点が計算できている。
【００７２】
スクリーン上の点Ａを伸ばしたベクトルvSDir[0]上の点ｖＡ、
スクリーン上の点Ｂを伸ばしたベクトルvSDir[1]上の点Ｂ´１、Ｂ´２、
スクリーン上の点Ｃを伸ばしたベクトルvSDir[2]上の点Ｃ´１、Ｃ´２。
【００７３】
このため、最大で以下の４つの組み合わせが存在するので、その組み合わせより、仮の点Ｂでの角度を計算する。  
（ｖＡ，Ｂ´１，Ｃ´１）、
（ｖＡ，Ｂ´１，Ｃ´２）、
（ｖＡ，Ｂ´２，Ｃ´１）、
（ｖＡ，Ｂ´２，Ｃ´２）
【００７４】
４つの組み合わせの内積値より、仮の点Ｂをはさむ辺の角度（ｃｏｓ値）を計算し、それぞれ、angle[0]、angle[1]、angle[2]、angle[3]とする。ステップＳ４０３での計算より、ワールド座標上のオリジナルの点Ｂでの角度はangleBであるため、
【００７５】
【数１９】
2010239515.tif 000021
【００７６】
で求めた４つの絶対値のうちゼロ値（0.0）に最も近いものを採用し、仮の点Ｂ、Ｃの位置とする（これをＢ´、Ｃ´とする）。また、このときのv[i]の最小値を誤差sValとして格納する。なお、誤差sValの初期値は9999.0などの大きな値を予め代入しておく。
【００７７】
ステップＳ８０７では、上記で求めた仮の点ＡＤＣが作る角度を計算する。ステップＳ８０６において、点Ｂ´、Ｃ´は固定されたので、ステップＳ８０４で計算された直線Ｄ上の点Ｄ´１、Ｄ´２を使って、仮の点ＡＤＣで作られる点Ｄでの角度を計算する。ここで、以下の２つの組み合わせに分かれる。  
（ｖＡ，Ｄ´１，Ｃ´）、
（ｖＡ，Ｄ´２，Ｃ´）
２つの組み合わせの内積値より、仮の点Ｄをはさむ辺の角度（ｃｏｓ値）を計算し、angle[4]、angle[5]とする。ステップＳ４０３での計算より、ワールド座標上のオリジナルの点Ｄでの角度はangleDであるため、
【００７８】
【数２０】
2010239515.tif 000022
【００７９】
で求めた２つの絶対値のうちゼロ値（0.0）に最も近いものを採用し、仮のＤの位置とする（これをＤ´とする）。また、このときのv[i]の最小値をステップＳ８０６で格納した誤差sValに加算する（sVal = sVal + v[i]の最小値）。
【００８０】
このときの誤差値であるｓＶａｌが完全に0.0になったとき、そのときのｖＡ、Ｂ´、Ｃ´、Ｄ´が求めるグローバル座標上での座標位置になる。だが、ピンポイントで誤差が収束することは少ないため、誤差が収束するように上記ステップを反復させることになる。
【００８１】
ステップＳ８０８では、視点からの距離ｔと誤差、求まった座標位置Ａ、Ｂ´、Ｃ´、Ｄ´を配列に格納する。すなわち、tA[i]に距離ｔを、sValA[i]に誤差の値sValを、Apos[i]、Bpos[i]、Cpos[i]、Dpos[i]にワールド座標値を格納していく。ここで、ｉは配列の添え字でminTからmaxTの距離を８分割したため０〜８の値が入る。
【００８２】
最後に、ステップＳ８０９では、ｔにｄＴを加算してステップＳ８０２に戻る。ｔがmaxTになるまで図８の処理を繰り返す。
【００８３】
再び、図４のステップＳ４０７に戻る。ステップＳ４０７では、収束計算を行う。図８で説明した各ステップの処理の結果、sValA[i]に誤差の値が格納されているので、この配列のうちで、一番値が小さくなるｉの位置を求める。以後、この誤差の配列にて谷となる部分を狭めていく処理を繰り返すことにより、誤差が0.0に近づくように収束させることになる。
【００８４】
ここで、図１０に示したように、谷の部分での視点位置からの距離tA[i]を中心にminT、maxTを更新する。すなわち、一番値が小さくなるｉの位置に相当する距離tA[i]と、従来のminT、maxTの間の距離の４分の１の距離とを各々減算、加算して、新たにminT、maxTとする。なお、新たなminTが初期値より小さくなる場合は初期値にし、新たなmaxTが初期値より大きくなる場合は初期値にする。新たなminT、maxTの間の距離を８等分した値をdTとし、距離ｔを新たなminTとする。その後、図８のステップＳ８０１に処理に戻る。
【００８５】
なお、maxTとminTの差が1.0より小さくなれば（ステップＳ４０８）、そのときのsValA[i]が最終的な誤差値であり、このときの特徴点位置がApos[i]、Bpos[i]、Cpos[i]、Dpos[i]と確定される。本発明によれば、距離は１回の処理のループにて２分割されていくため、例えば、初期値minT = 5.0、maxT = 99999.0の場合は１４回の繰り返し処理にて最小の位置で誤差は収束する。
【００８６】
ステップＳ４０９では、ワールド座標からビュー座標に変換する行列を求める。ステップＳ４０６〜４０７およびステップＳ８０１〜Ｓ８０９にて求められた３次元上の座標値Apos[i]、Bpos[i]、Cpos[i]、Dpos[i]を改めてw2[0]、w2[1]、w2[2]、w2[3]とおく。このw2[i]は、カメラを姿勢ベクトル＝単位行列、視点位置＝（0, 0, 0）としたときの座標値である。これらの座標値w2[0]、w2[1]、w2[2]、w2[3]を実際のワールド座標位置w[0]、w[1]、w[2]、w[3]の座標系に変換するときの変換行列からカメラの姿勢と移動ベクトルが算出することができる。回転行列は以下の式にて計算することができる。
【００８７】
【数２１】
2010239515.tif 000023
【００８８】
ここで、４×４行列にてオリジナルの姿勢と位置を含めたカメラ情報を以下のようにする。
【００８９】
【数２２】
2010239515.tif 000024
【００９０】
同様に、推定された回転行列Ｒと位置w2[0]より、推定後の４×４行列は以下のようになる。
【００９１】
【数２３】
2010239515.tif 000025
【００９２】
ここで、Ｎ＝inv（ｏｒｇＭ）＊Ｍを計算する（invは、４×４行列の逆行列を計算する関数）。計算された４×４行列Ｎが、ワールド座標からビュー座標への変換行列に相当する。
【００９３】
ステップＳ４１０では、ワールド座標からビュー座標に変換する行列Ｎからカメラの回転行列と移動ベクトルを逆算する。
【００９４】
【数２４】
2010239515.tif 000026
【００９５】
ここで求まったrMatがカメラの回転行列、tVecがカメラの移動ベクトルとなる。
ステップＳ４１１では、ステップＳ４１０にて推定されたカメラの外部パラメータの誤差を補正する。上述した各ステップの処理にて、カメラの外部パラメータ（回転行列、移動ベクトル）が求まるが、特徴点４点のみの計算であるため、誤差が発生している可能性がある。また、図１１に示すように、カメラキャリブレーションでは解が２つ以上存在するパターンもあり、その場合は軸が大きく反転することになる。カメラ位置と姿勢が原点を中心に反対に位置し、かつ、カメラ姿勢の軸が反転している場合は、ワールド座標上の特徴点とカメラの投影面に投影した座標とは、共に正しく対応している関係になる。本発明では、この現象の回避策として、１つ前のカメラ姿勢でのＸ軸とＹ軸と、新しいカメラ姿勢でのＸ軸とＹ軸との傾きを調べ、傾きが小さい方を採用する手法をとることができる。
【００９６】
一方、誤差の補正処理は、上記で求めたカメラの外部パラメータ（回転行列、移動ベクトル）よりワールド座標からビュー座標への変換行列を計算し、その行列を使って特徴点のワールド座標位置（Wxn, Wyn, Wzn）をスクリーン座標値に変換する。このときのスクリーン上の点の位置を(U'n, V'n)とし、オリジナルのスクリーン上の点（Un, Vn）との距離の二乗和を評価のための誤差とする。
【００９７】
【数２５】
2010239515.tif 000027
【００９８】
この二乗和の値Ｆが0.0に近いほど推定されたカメラの外部パラメータの精度が高いことを示す。
【００９９】
また、本発明では、この誤差値が大きい場合は、スクリーン上の位置（Un, Vn）に微小のノイズを加えて揺らすことにより収束を早める、という手段をとることができる。スクリーン上の特徴点位置にノイズを加える理由は、スクリーン上の特徴点は特徴点検出などですでにノイズが加わった状態で与えられる場合があるためである。
【０１００】
このような誤差の補正処理を加えたカメラキャリブレーションを図１２及び図１３のフローチャートで示す。
【産業上の利用可能性】
【０１０１】
本発明は、一般的なカメラキャリブレーションの算出方法として利用することができる。より詳細に言えば、本発明は、拡張現実（ＡＲ：Augmented Reality）の分野に利用することができる。
【０１０２】
ＡＲとは、カメラで撮影した動画に対してコンピュータで生成した三次元空間をオーバーラップさせることにより、現実空間に人工的な付加価値をつけることができる表現手法である。ＡＲを利用することにより、カメラを通して見た市街地の個々の建物に対して、その建物に入居している会社や店舗情報をオーバーラップして表示したり（タグ付け）、待ち合わせ場所までの道順を仮想的なラインを引いて導いたり、あるいは、室内のインテリアの配置や壁の装飾を仮想的に生成・配置して、実画像を通しながら仕上がりを確認する、などの応用を実現することができる。
【０１０３】
本発明に係るカメラキャリブレーションは、ＡＲでのカメラの位置と姿勢を推定する部分にて使用することが可能である。従来技術のようにカメラキャリブレーションの精度が悪いと、３次元空間を構築するときにずれやカメラのぶれが生じてしまうという不都合が発生するが、本発明によりそのような従来技術の不都合を解消することができる。
【符号の説明】
【０１０４】
３００  装置
３０１  カメラ
３０２  有線／無線媒体
３１１  制御部
３１２  主記憶部
３１３  補助記憶部
３１４  入出力Ｉ／Ｆ
３１５  出力表示部
【図１】
2010239515.tif 000028
【図２】
2010239515.tif 000029
【図３】
2010239515.tif 000030
【図４】
2010239515.tif 000031
【図５】
2010239515.tif 000032
【図６】
2010239515.tif 000033
【図７】
2010239515.tif 000034
【図８】
2010239515.tif 000035
【図９】
2010239515.tif 000036
【図１０】
2010239515.tif 000037
【図１１】
2010239515.tif 000038
【図１２】
2010239515.tif 000039
【図１３】
2010239515.tif 000040