(57)【要約】【課題】この発明は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を自動的に算出することができるチルト角算出方法を提供することを目的とする。【解決手段】カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を算出する方法であって、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求める第１ステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求める第２ステップを備えている。【選択図】　図９
2004046573.tif 000002
【特許請求の範囲】
【請求項１】
カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を算出する方法であって、
撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求める第１ステップ、および
上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求める第２ステップ、
を備えていることを特徴とするチルト角算出方法。
【請求項２】
第１ステップは、
撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および
上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップ、
を備えていることを特徴とする請求項１に記載のチルト角算出方法。
【請求項３】
カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成する円筒パノラマ画像合成方法であって、
撮像された複数の画像からチルト角を算出する第１ステップ、得られたチルト角を考慮して各画像を円筒面投影画像に変換する第２ステップ、得られた各円筒面投影画像を接合する第３ステップを備えており、
第１ステップは、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求めるステップを備えていることを特徴とする円筒パノラマ画像合成方法。
【請求項４】
撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップは、
撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、　一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および
上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップ、
を備えていることを特徴とする請求項３に記載の円筒パノラマ画像合成方法。
【請求項５】
カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成するための円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
撮像された複数の画像からチルト角を算出する第１ステップ、得られたチルト角を考慮して各画像を円筒面投影画像に変換する第２ステップ、得られた各円筒面投影画像を接合する第３ステップをコンピュータに実行させるための円筒パノラマ画像合成処理プログラムを記録しており、
第１ステップは、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求めるステップを備えていることを特徴とする、円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項６】
撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップは、
撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、　一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および
上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップ、
を備えていることを特徴とする請求項５に記載の円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
この発明は、チルト角算出方法、円筒パノラマ画像合成方法、ならびに円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
【０００２】
【従来の技術】
画像間で重なりを持つ複数の画像を貼り合わせて、広視野かつ高解像度の画像を作成する技術が活発に研究されている。古典的な手法では、衛星写真や航空写真の貼り合わせに使用されているが、最近では、実シーンを撮像した複数の静止画像を円筒面に投影し、貼り合わせることで１枚のパノラマ画像を作成する手法が注目されている（　文献［１］　参照）　。
【０００３】
文献［１］　：Ｓ．Ｅ．Ｃｈｅｎ，”Ｑｕｉｃｋ　Ｔｉｍｅ　ＶＲ　−　Ａｎ　Ｉｍａｇｅ−Ｂａｓｅｄ　Ａｐｐｒｏａｃｈ　ｔｏ　ＶｉｒｔｕａｌＥｎｖｉｒｏｎｍｅｎｔ　Ｎａｖｉｇａｔｉｏｎ，”　ＳＩＧＧＲＡＰＨ’９５，　ｐｐ．　２９−３８，　１９９５．
【０００４】
これにより、三脚等を用いてカメラ運動をパンのみに制限して撮影した画像から、水平方向３６０度のパノラマ画像を合成することができる。この画像を用いて、ユーザは希望する視線方向の画像を対話的に表示し、臨場感のある仮想環境を体験することができる。
【０００５】
従来の全方位を撮像可能なパノラマ合成手法として、カメラを自由に回転させて撮影した画像を球面に投影して貼り付ける手法がある（　文献［２］　参照）　。しかし、幾何変換パラメータが８つあり、自由度が高いため、単一色の壁や空など濃淡の変化の少ない画像を精度よく合成できないという問題がある。
【０００６】
文献［２］　：Ｒ．　Ｓｚｅｌｉｓｋｉ　ａｎｄ　Ｈ．　Ｙ．　Ｓｈｕｍ，　”Ｃｒｅａｔｉｎｇ　Ｆｕｌｌ　Ｖｉｅｗ　Ｐａｎｏｒａｍｉｃ　Ｉｍａｇｅ　Ｍｏｓａｉｃｓ　ａｎｄ　Ｅｎｖｉｒｏｎｍｅｎｔ　Ｍａｐｓ，”ＳＩＧＧＲＡＰＨ’９７，　ｐｐ．　２５１−２５８，　１９９７．
【０００７】
一方、円筒パノラマ合成は、カメラ運動をパンのみに制限するため自由度は１つと少ない。そのため、濃淡の変化の少ない画像でも合成精度は高い。しかし、垂直方向の撮影範囲が制限されるという問題がある。ただし、円筒パノラマ合成を一定のチルト角を含んだままパンする場合に拡張するのは容易である。この際には、チルト角を求める必要がある。
【０００８】
なお、一定のチルト角を含む撮影条件は、高い対象物を近くから見上げて撮影する場合や、高い場所からパノラマ風景を見下ろしながら撮影する場合に発生する。
【０００９】
ところで、チルト角を求める従来の方法には、手動で求める方法と、自動で求める方法とがある。
【００１０】
手動でチルト角を求める場合には、特殊な治具を使用する必要があり、試行錯誤を繰り返す必要がある。そのため、手間がかかるという問題がある。
【００１１】
Ｓｚｅｌｉｓｋｉ（　上記文献［２］　参照）　やＳｅｏ（文献［３］　参照）　は自動でチルト角を算出する手法を提案している。
【００１２】
文献［３］　：Ｙ．　Ｓｅｏ，”Ａ　Ｍｕｌｔｉｐｌｅ　Ｖｉｅｗ　Ａｐｐｒｏａｃｈ　ｆｏｒ　Ａｕｔｏ−ｃａｌｉｂｒａｔｉｏｎ　ｏｆ　ａ　Ｒｏｔａｔｉｎｇ　ａｎｄ　Ｚｏｏｍｉｎｇ　Ｃａｍｅｒａ，”ＩＥＣＩＥ　Ｔｒａｎｓ．　Ｉｎｆ．　＆　Ｓｙｓｔ．，　Ｖｏｌ．　Ｅ８３−Ｄ，　Ｎｏ．７，ｐｐ．１３７５−１３８５，２０００．
【００１３】
これは、キャリブレーション問題と呼ばれ、チルト角だけでなく、カメラ間の相対的な３つの回転角や焦点距離などのカメラの内部パラメータを求める手法である。しかし、これらの手法はカメラ間の相対的なチルト角しか求めていない。
そのため、カメラをチルト角一定のままパンしたときのチルト角は求まらない。
【００１４】
【発明が解決しようとする課題】
この発明は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を自動的に算出することができるチルト角算出方法を提供することを目的とする。
【００１５】
また、この発明は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成できる円筒パノラマ画像合成方法および円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することを目的とする。
【００１６】
【課題を解決するための手段】
この発明によるチルト角算出方法は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を算出する方法であって、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求める第１ステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求める第２ステップを備えていることを特徴とする。
【００１７】
第１ステップとしては、たとえば、撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップを備えているものが用いられる。
【００１８】
この発明による円筒パノラマ画像合成方法は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成する円筒パノラマ画像合成方法であって、撮像された複数の画像からチルト角を算出する第１ステップ、得られたチルト角を考慮して各画像を円筒面投影画像に変換する第２ステップ、得られた各円筒面投影画像を接合する第３ステップを備えており、第１ステップは、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求めるステップを備えていることを特徴とする。
【００１９】
撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップとしては、たとえば、撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップを備えていることを特徴とする。
【００２０】
この発明による記録媒体は、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成するための円筒パノラマ画像合成処理プログラムを記録したコンピュータ読み取り可能な記録媒体であって、撮像された複数の画像からチルト角を算出する第１ステップ、得られたチルト角を考慮して各画像を円筒面投影画像に変換する第２ステップ、得られた各円筒面投影画像を接合する第３ステップをコンピュータに実行させるための円筒パノラマ画像合成処理プログラムを記録しており、第１ステップは、撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップ、および上記２つの画像間の対応点が射影的に等しいことから導出される、パン角度とチルト角との関係式と、上記２つの画像間の２組以上の対応点の座標値とに基づいて、チルト角を求めるステップを備えていることを特徴とする。
【００２１】
撮影された重なり部を有する２つの画像間の２組以上の対応点の座標値を求めるステップは、撮影された重なり部を有する２つの画像間の重なり部を抽出するステップ、一方の画像における他方の画像との重なり部分から、両画像間のオプティカルフローによる追跡に有効な複数の部分画像を特徴点として抽出するステップ、および上記一方の画像上の各特徴点に対応する上記他方の画像上の点を、両画像間のオプティカルフローに基づいて追跡するステップを備えている。
【００２２】
【発明の実施の形態】
この発明の実施の形態について説明する前に、この発明の実施の形態で用いられるオプティカルフローの算出方法について説明しておく。
【００２３】
まず、従来のオプティカルフローの算出方法について説明する。
【００２４】
〔１〕従来のオプティカルフローの算出方法の説明
２枚の画像からオプティカルフローを計算し、得られたオプティカルフローに基づいて、２枚の画像間での位置合わせを行う技術が知られている。従来のオプテカルフローの算出方法について説明する。
【００２５】
（１）Ｌｕｃａｓ−Ｋａｎａｄｅ法
従来から、動画像における運動物体の見かけの速度場（オプティカルフロー）を計算する手法が数多く提案されている。中でも局所勾配法であるＬｕｃａｓ−Ｋａｎａｄｅ法は、最も良い手法の一つである。その理由は、処理が高速、実装が容易、結果が信頼度を持つことである。
【００２６】
Ｌｕｃａｓ−Ｋａｎａｄｅ法の詳細については、文献［４］　を参照のこと。
文献［４］　：　Ｂ．Ｌｕｃａｓ　ａｎｄ　Ｔ．Ｋａｎａｄｅ，”Ａｎ　Ｉｔｅｒａｔｉｖｅ　Ｉｍａｇｅ　Ｒｅｇｉｓｔｒａｔｉｏｎ　Ｔｅｃｈｎｉｑｕｅ　ｗｉｔｈ　ａｎ　Ａｐｐｌｉｃａｔｉｏｎ　ｔｏ　Ｓｔｅｒｅｏ　Ｖｉｓｉｏｎ，”Ｉｎ　Ｓｅｖｅｎｔｈ　Ｉｎｔｅｒｎａｔｉｏｎａｌ　Ｊｏｉｎｔ　Ｃｏｎｆｅｒｅｎｃｅ　ｏｎ　Ａｒｔｉｆｉｃｉａｌ　Ｉｎｔｅｌｌｉｇｅｎｃｅ（ＩＪＣＡＩ−８１），　ｐｐ．　６７４−９７９，　１９８１．
【００２７】
以下に、Ｌｕｃａｓ−Ｋａｎａｄｅ法の概要を述べる。ある時刻ｔの画像座標ｐ＝（ｘ，ｙ）の濃淡パターンＩ（ｘ，ｙ，ｔ）が、ある微小時間後（δｔ）に座標（ｘ＋δｘ，ｙ＋δｙ）に、その濃淡分布を一定に保ったまま移動した時、次のオプティカルフロー拘束式（１）が成り立つ。
【００２８】
【数１】
2004046573.tif 000003
【００２９】
２次元画像でオプティカルフロー｛ｖ＝（δｘ／δｔ，δｙ／δｔ）＝（ｕ，ｖ）｝を計算するには、未知パラメータ数が２個であるので、もう一個拘束式が必要である。Ｌｕｃａｓ　とＫａｎａｄｅ（　金出）　は、同一物体の局所領域では、同一のオプティカルフローを持つと仮定した。
【００３０】
例えば、画像上の局所領域ω内で、オプティカルフローが一定であるとすると、最小化したい濃淡パターンの二乗誤差Ｅは、
Ｉ０　（ｐ）＝Ｉ（ｘ，ｙ，ｔ），
Ｉ１　（ｐ＋ｖ）＝Ｉ（ｘ＋ｕ，ｙ＋ｖ，ｔ＋δｔ）
と書き改めると、次式（２）で定義できる。
【００３１】
【数２】
2004046573.tif 000004
【００３２】
ここで、ｖが微少な場合には、テーラー展開の２次以上の項を無視できるので、次式（３）の関係が成り立つ。
【００３３】
【数３】
2004046573.tif 000005
【００３４】
ここで、ｇ（ｐ）は、Ｉ１　（ｐ）の一次微分である。
【００３５】
誤差Ｅが最小になるのは、Ｅのｖに対する微分値が０の時であるので、次式（４）の関係が成り立つ。
【００３６】
【数４】
2004046573.tif 000006
【００３７】
故にオプティカルフローｖは次式（５）で求められる。
【００３８】
【数５】
2004046573.tif 000007
【００３９】
更に、次式（６）に示すように、ニュートン・ラフソン的な反復演算によって精度良く求めることができる。
【００４０】
【数６】
2004046573.tif 000008
【００４１】
（２）階層的推定法
Ｌｕｃａｓ−Ｋａｎａｄｅ法を含む勾配法の最も大きな問題点は、良好な初期値が必要なために、大きな動きに対しては適用できないことである。そこで、従来からピラミッド階層構造型に数段回の解像度の異なる画像を作成して解決する方法が提案されている。
【００４２】
これは、まず、２枚の連続した画像から、予めそれぞれの画像の数段階の解像度の異なる画像を作成する。次に、最も解像度の低い画像間において、おおまかなオプティカルフローを計算する。そして、この結果を参考にして、一段解像度の高い画像間においてより精密なオプティカルフローを計算する。この処理を最も解像度の高い画像間まで順次繰り返す。
【００４３】
図４は原画像を、図３は図４の原画像より解像度の低い画像を、図２は図３の低解像度画像より解像度の低い画像を、図１は図２の低解像度画像より解像度の低い画像を、それぞれ示している。図１〜図４において、Ｓは、１つのパッチを示している。
【００４４】
図１の画像（階層１の画像）、図２の画像（階層２の画像）、図３の画像（階層３の画像）および図４の画像（階層４の画像）の順番で段階的にオプティカルフローが求められる。図１〜図４において矢印は、パッチ毎に求められたオプティカルフローベクトルを示している。
【００４５】
しかしながら、ここでの問題点は、実画像では、十分な模様（テクスチャ）を含む領域が少なく、信頼性のあるオプティカルフローが得られないことにある。
【００４６】
次に、この実施の形態で採用されるオプティカルフロー推定方法について説明する。
【００４７】
〔２〕この実施の形態で採用されるオプティカルフロー推定方法についての説明この実施の形態で採用されるオプティカルフロー推定方法は、ピラミッド階層型に数段回の解像度の異なる画像を作成して、オプティカルフローを段階的に計算する階層的推定を前提としている。オプティカルフローの計算方法は、Ｌｕｃａｓ−Ｋａｎａｄｅ法等の勾配法に従う。つまり、階層構造化した勾配法によるオプティカルフロー推定法を前提としている。ここでは、勾配法としてＬｕｃａｓ−Ｋａｎａｄｅ法が用いられている。
【００４８】
この実施例で採用されるオプティカルフロー推定方法の特徴は、階層構造化したＬｕｃａｓ−Ｋａｎａｄｅ法によるオプティカルフロー推定法の各段階において得られたオプティカルフローを、膨張処理によって補完することにある。以下、これについて詳しく説明する。
【００４９】
Ｌｕｃａｓ−Ｋａｎａｄｅ法の長所の一つは、追跡結果が信頼性を持つことである。Ｔｏｍａｓｉと　Ｋａｎａｄｅ　とは、ある領域の追跡可能性が、以下のように微分画像から算出できることを示した（Ｃ．Ｔｏｍａｓｉ　ａｎｄ　Ｔ．Ｋａｎａｄｅ，”Ｓｈａｐｅ　ａｎｄ　Ｍｏｔｉｏｎ　ｆｒｏｍ　Ｉｍａｇｅ　Ｓｔｒｅａｍｓ：　ａ　Ｆａｃｔｏｒｉｚａｔｉｏｎ　ｍｅｔｈｏｄ−Ｐａｒｔ　３　Ｄｅｔｅｃｔｉｏｎ　ａｎｄ　Ｔｒａｃｋｉｎｇ　ｏｆ　Ｐｏｉｎｔ　Ｆｅａｔｕｒｅｓ　，”ＣＭＵ−ＣＳ−９１−１３２，　Ｃａｒｎｅｇｉｅ　Ｍｅｌｌｏｎ　Ｕｎｉｖｅｒｓｉｔｙ，　１９９１．）　。
【００５０】
ある領域画像ωの垂直・水平方向の微分の２乗を要素に持つ次式（７）の２×２の係数行列Ｇから、その固有値を計算することで、その領域の追跡可能性を決定することができる。
【００５１】
【数７】
2004046573.tif 000009
【００５２】
この行列Ｇの固有値が両方とも大きい場合には、その領域は直交方向に変化を持ち、一意の位置決めが可能である。従って、小さい方の固有値λｍｉｎ　と、追跡後の領域間の濃淡残差Ｅから、追跡結果の信頼度γを次式（８）によって得ることができる。
【００５３】
【数８】
2004046573.tif 000010
【００５４】
本発明者らは、オプティカルフローの同一階層内で信頼度の高い結果を用いて、信頼度の低い領域を補間する方法を開発した。これは、一段階粗い階層での結果を、追跡の初期値だけに用いて、着目している現段階の階層の結果には何も利用しない。代わりに、テクスチャの少ない領域のオプティカルフローはその周囲のオプティカルフローに近い値を持つと仮定し、モルフォロジー処理によりフロー場を補完するものである。
【００５５】
図５にフローベクトルの膨張処理の様子を示す。
【００５６】
左図は、フローベクトルの信頼度のマップを濃淡で表したものである。ここで、黒ければ黒い程信頼度が高いとする。
【００５７】
まず、得られたフローをしきい値処理する。白い部分は、結果の信頼度が低いために、しきい値処理されたものである。
【００５８】
次に、２値画像でのモルフォロジー演算による穴埋め処理を模して、フロー場において結果の膨張処理を次のように行う。ある領域ｉ，ｊのフローベクトルｕ（ｉ，ｊ）は、その４近傍のフローベクトルから信頼度γに応じて重み付けを行って次式（９）のように計算できる。
【００５９】
【数９】
2004046573.tif 000011
【００６０】
この処理を、しきい値処理されたすべての信頼度の低い領域が埋まるまで、繰り返す。この補完処理を、各階層において行う。なお、ある領域ｉ，ｊのフローベクトルｕ（ｉ，ｊ）を、その８近傍のフローベクトルから信頼度γに応じて重み付けを行って算出するようにしてもよい。
【００６１】
図６（ａ）は、ある階層の画像に対してしきい値処理されたオプティカルフローを示し、図６（ｂ）は補完後のオプティカルフローを示している。図６（ａ）において、矢印はしきい値処理によって信頼度が高いと判定されたオプティカルフローベクトルであり、×印は信頼度が低いとされた部分を示している。
【００６２】
〔３〕円筒パノラマ画像合成装置についての説明
【００６３】
この実施の形態では、図７に示すように、カメラ１を垂直方向に所定角度回転（チルト）させたまま水平回転（パン）させながら撮影された画像から円筒パノラマ画像を合成する円筒パノラマ画像合成装置について説明する。チルト角Ψは一定である。図７において、２は撮像対象物を示している。また、Ｏ−ＸＹＺは世界座標系を示している。ただし、Ｏは視点である。
【００６４】
図８は、円筒パノラマ画像合成を実現するための装置を示している。
【００６５】
円筒パノラマ画像合成は、パーソナルコンピュータ（以下ＰＣという）によって実現される。
【００６６】
ＰＣ１０には、ディスプレイ２１、マウス２２およびキーボード２３が接続されている。ＰＣ１０は、ＣＰＵ１１、メモリ１２、ハードディスク１３、ＣＤ−ＲＯＭのようなリムーバブルディスクのドライブ（ディスクドライブ）１４を備えている。
【００６７】
ハードディスク１３には、ＯＳ（オペレーティングシステム）等の他、円筒パノラマ画像合成を行なうための画像合成処理プログラムが格納されている。画像合成処理プログラムは、例えば、それが格納されたＣＤ−ＲＯＭ２０を用いて、ハードディスク１３にインストールされる。また、ハードディスク１３には、図７に示すように、カメラ１を垂直方向に所定角度回転（チルト）させたまま水平回転（パン）させながら撮影された複数の画像が予め格納されているものとする。
【００６８】
〔４〕ＣＰＵ１１によって行なわれる円筒パノラマ画像合成処理手順の説明
【００６９】
図９は、ＣＰＵ１１によって行なわれる円筒パノラマ画像合成処理の全体的な処理手順を示している。
【００７０】
円筒パノラマ画像合成処理では、画像読み込み処理（ステップ１）、重なり部の抽出処理（ステップ２）、特徴点抽出処理（ステップ３）、特徴点追跡処理（ステップ４）、チルト角の算出処理（ステップ５）、チルト角を考慮した円筒面投影変換処理（ステップ６）および画像接合処理（ステップ７）が行なわれる。
【００７１】
以下、これらの各ステップについて説明する。
【００７２】
〔４．１〕画像読み込み処理（ステップ１）の説明
ステップ１では、ユーザによって指定された接合されるべき複数の撮像画像がメモリに読み込まれる。
【００７３】
〔４．２〕重なり部の抽出処理（ステップ２）の説明
ステップ２では、重なり部を有する２枚の撮像画像（第１画像Ａ１及び第２画像Ａ２）の重なり部が抽出される。この処理は、たとえば、ＳＳＤ法（Ｓｕｍ　ｏｆ　Ｓｑｕａｒｅｄ　Ｄｉｆｆｅｒｅｎｃｅ）、正規化相互相関法に基づいて行われる。
【００７４】
（ａ）ＳＳＤ法の説明
ＳＳＤ法では、まず、重なり部を抽出すべき２枚の画像Ａ１、Ａ２それぞれについて、原画像より解像度の低い画像Ｉ１　，Ｉ２　が生成される。２枚の低解像度画像Ｉ１　，Ｉ２　の重なり部分ω（サイズ：Ｍ×Ｎ）が、次式（１０）に示すように画素あたりの２乗誤差Ｅを用いて求められる。画像間の移動量（ｄ）が可能な範囲で変化せしめられ、Ｅが最も小さい移動量（ｄ）から、重なり部が抽出される。
【００７５】
【数１０】
2004046573.tif 000012
【００７６】
（ｂ）正規化相互相関法の説明
正規化相互相関法では、まず、重なり部を抽出すべき２枚の画像Ａ１、Ａ２それぞれについて、原画像より解像度の低い画像Ｉ１　，Ｉ２　が生成される。２枚の低解像度画像Ｉ１　，Ｉ２　の重なり部分ω（サイズ：Ｍ×Ｎ）が、次式（１１）に示すように正規化相互相関係数Ｃを用いて求められる。画像間の移動量（ｄ）が可能な範囲で変化せしめられ、Ｃが最も大きな移動量（ｄ）から、重なり部が抽出される。
【００７７】
【数１１】
2004046573.tif 000013
【００７８】
式（１１）において、Ｉ１　￣、Ｉ２　￣は、第１画像を固定させ、第２画像をｄだけ移動させたときの、両画像の重なり部における各画像それぞれの濃淡値の平均である。また、σ１　、σ２　は、第１画像Ｉ１　を固定させ、第２画像Ｉ２　をｄだけ移動させたときの、両画像の重なり部における各画像それぞれの濃淡値の分散である。
【００７９】
〔４．３〕特徴点抽出処理（ステップ３）の説明
ステップ３では、まず、第１画像Ａ１における第２画像Ａ２との重なり部分から、追跡に有効な複数の部分画像（矩形領域）が特徴点として抽出される。ただし、各特徴点は互いに重ならないように抽出される。具体的には、上述した固有値λｍｉｎ　（式（８）参照）の高い部分が特徴点として抽出される。
【００８０】
〔４．４〕特徴点追跡処理（ステップ４）の説明
ステップ４では、抽出された第１画像Ａ１上の特徴点に対応する第２画像Ａ２上の位置が追跡される。
【００８１】
具体的には、まず、本出願人が開発したオプティカルフロー推定方法（上記〔２〕のこの実施の形態で採用されるオプティカルフロー推定方法）によって、適当なサイズ（例えば、１３×１３）のパッチ毎のオプティカルフローベクトルが求められる。第１画像Ａ１上の特徴点に対応する第２画像Ａ２上の位置は、第１画像Ａ１上の特徴点の４近傍のパッチのフローベクトルから線形補間により画素単位以下で求められる。これにより、第１画像Ａ１と第２画像Ａ２との重なり部において、両画像の対応点の座標が得られる。
【００８２】
ここでは、重なり部を有する各組の画像毎に両画像の対応点（特徴点）の座標が得られたものとする。
【００８３】
〔４．５〕チルト角の算出処理（ステップ５）の説明
ステップ５では、ステップ４で得られた対応点（特徴点）の座標に基づいて、チルト角が算出される。
【００８４】
まず、２つの画像間の２以上の対応点（特徴点）に基づいて、チルト角を算出する方法について説明する。
【００８５】
図１０に示すように、３次元空間上の点Ｐを、１台のカメラでチルト角Ψを一定としてかつパン角度φだけ離れた２つの位置において撮像したとき、チルト角Ψはそれらのカメラ位置の関係から求まる。
【００８６】
点Ｐを各撮像画像Ｉ１　、Ｉ２　に投影し、世界座標系Ｏ−ＸＹＺで表した点　　　　Ｐ１　＝（Ｘ１，Ｙ１，Ｚ１　）Ｔ　およびＰ２　＝（Ｘ２，Ｙ２，Ｚ２　）Ｔ　は次のように表記できる。
【００８７】
まず、点Ｐ１　は、次のように表せる。図１１に示すように、世界座標系Ｏ−ＸＹＺをΨ度チルトした座標系を第１のカメラ座標系Ｏ−Ｕ１　Ｖ１　Ｗ１　とする。このとき、世界座標系の点Ｐ１　は、この点を第１のカメラ座標系で表した点ｐ１　＝（ｕ１，ｖ１，−ｆ）Ｔ　およびＸ軸周りの回転行列ＲＸ　を用いて次式（１２）のように表記できる。
【００８８】
【数１２】
2004046573.tif 000014
【００８９】
ただし、（ｕ１，ｖ１　）は点ｐ１　を画像座標で表した点の座標であり、ｆは焦点距離を表している。
【００９０】
この式（１２）は、各要素を表示すると、次式（１３）のように表記できる。
【００９１】
【数１３】
2004046573.tif 000015
【００９２】
一方、点Ｐ２　は、次のように表せる。図１２に示すように、世界座標系Ｏ−ＸＹＺをΨ度チルトしたまま、さらにφ度パンした座標系を第２のカメラ座標系Ｏ−Ｕ２　Ｖ２　Ｗ２　とする。このとき、世界座標系の点Ｐ２　は、この点を第２のカメラ座標系で表した点ｐ２　＝（ｕ２，ｖ２，−ｆ）Ｔ　、Ｘ軸周りの回転行列ＲＸ　およびＹ軸周りの回転行列ＲＹ　を用いて次式（１４）のように表記できる。
【００９３】
【数１４】
2004046573.tif 000016
【００９４】
ただし、（ｕ２，ｖ２　）は画像座標（ｕ１，ｖ１　）に対応する画像Ｉ２　上の画像座標である。なお、対応する画像座標（ｕ１，ｖ１　）と（ｕ２，ｖ２　）としては、図９のステップ２、３、４によって得られた対応点（特徴点）が用いられる。また、焦点距離ｆはユーザによって与えられる既知の定数である。
【００９５】
この式（１４）は、各要素を表示すると、次式（１５）のように表記できる。
【００９６】
【数１５】
2004046573.tif 000017
【００９７】
点Ｐ１　と点Ｐ２　とは、射影的に等しいので、次式（１６）が成り立つ。
【００９８】
【数１６】
2004046573.tif 000018
【００９９】
ただし、ｋは０でない定数である。式（１６）は、Ｐ１　とＰ２　の各要素間の比が等しいことを表しているので、係数ｋを削除して次式（１７）のように整理できる。
【０１００】
【数１７】
2004046573.tif 000019
【０１０１】
上記式（１３）および（１５）を、上記式（１７）に代入すると、次式（１８）が得られる。
【０１０２】
【数１８】
2004046573.tif 000020
【０１０３】
上記式（１８）は、チルト角Ψに関して、次式（１９）および（２０）のように整理できる。
【０１０４】
【数１９】
2004046573.tif 000021
【０１０５】
【数２０】
2004046573.tif 000022
【０１０６】
ここで、
ａ１１＝（ｆ２　−ｖ１　ｖ２　）　ｓｉｎφ
ａ１２＝ｆ（ｖ１　＋ｖ２　）　ｓｉｎφ
ａ１３＝（ｕ１　ｖ２　−ｕ２　ｖ１　）　ｃｏｓφ
ａ１４＝ｆ（ｕ２　−ｕ１　）　ｃｏｓφ
ａ２１＝ｆ（ｖ１　＋ｖ２　）（　ｃｏｓφ−１）
ａ２２＝（ｖ１　ｖ２　−ｆ２　）（　ｃｏｓφ−１）
ａ２３＝−ｆｕ２ｓｉｎφ
ａ２４＝−ｕ２　ｖ１ｓｉｎφ
ｂ１　＝（ｆ２　＋ｕ１　ｕ２　）　ｓｉｎφ
ｂ２　＝ｆ（ｖ１ｃｏｓφ−ｖ２　）
とすると、上記式（１９）および（２０）は、次式（２１）のように書き換えることができる。
【０１０７】
【数２１】
2004046573.tif 000023
【０１０８】
ただし、係数ａ１１，ａ１２…，ａ２４およびｂ１　，ｂ２　は定数であり、上記のように、パン角φ、焦点距離ｆ、対応する画像座標（ｕ１，ｖ１　），（ｕ２，ｖ２　）から求まる。上記式（２１）からわかるように、パン角φはチルト角Ψに依存する変数であるが、チルト角の初期値をΨ＝０として反復的に求めることができる。
【０１０９】
つまり、上記式（１９）からパン角φを反復的に求めることができる。上記式（１９）をパン角φに関して整理すると、次式（２２）が求まる。
【０１１０】
【数２２】
2004046573.tif 000024
【０１１１】
対応する特徴点座標（ｕ１，ｖ１　），（ｕ２，ｖ２　）および焦点距離ｆは既知であり、未知数はチルト角Ψおよびパン角φである。このため、上記式（２２）はチルト角Ψが求まれば、パン角φが求まることを意味している。そこで、次式（２３）を用いて、チルト角の初期値をΨ０　＝０としてパン角φｉ　（ｉ＝０，１，２，…）を計算する。
【０１１２】
【数２３】
2004046573.tif 000025
【０１１３】
そして、パン角とチルト角の一方を固定して、他方を求めることを繰り返すことで、パン角を精度よく求めることができる。
【０１１４】
ところで、上記式（２１）は２画像間のすべての対応する点に対して成り立つので、次式（２４）のように表すことができる。
【０１１５】
【数２４】
2004046573.tif 000026
【０１１６】
ただし、ｎ≧２は点の対応の組数を表す。求めたいのは、チルト角Ψであり、上記式（２４）の列ベクトルＸから求まる。列ベクトルＸは最小２乗法で求まる。
【０１１７】
最後に、列ベクトルＸで求まった三角関数からチルト角Ψを計算する。ここで、チルト角Ψに関する三角関数は４つ（ｓｉｎ２Ψ、　ｓｉｎΨ　ｃｏｓΨ、　ｓｉｎΨおよび　ｃｏｓΨ）あるが、この実施の形態では　ｓｉｎΨからΨを求める。これは次の理由による。実際の撮影条件を考慮すると、チルト角はおよそ−４５度から４５度の間の値である。このような場合には、上記４つの三角関数の中で、　ｓｉｎΨが最も線形に近いからである。
【０１１８】
図１３は、チルト角の計算処理手順を示している。
【０１１９】
ここでは、チルト角およびパン角をΨｉ　、φｉ　（ｉ＝０，１，２，…）で表すことにする。
【０１２０】
まず、ｉ＝０とするとともに、チルト角の初期値Ψ０　を０とする（ステップ１１）。そして、Ψｉ　を上記式（２３）に代入して、パン角φｉ　を求める（ステップ１２）。次に、得られたパン角φｉ　を上記式（２４）に代入して、チルト角Ψｉ＋１　を算出する（ステップ１３）。
【０１２１】
そして、ｉが予め定められた所定値Ｎ（たとえば、２０）に達しているか否かを判別する（ステップ１４）。ｉがＮに達していなければ、ｉを１だけインクリメント（ｉ＝ｉ＋１）した後（ステップ１５）、ステップ１２に戻る。
【０１２２】
このようにして、ステッブ１２〜ステップ１３の処理が所定回繰り返され、ステップ１４において、ｉがＮに達すると、処理を終了する。最後に求まったチルト角Ψが、２つの画像間の対応点から得られたチルト角とする。
【０１２３】
重なり部を有する画像が複数組存在する場合には、各組毎にそれらの対応点を用いてチルト角を算出し、各組毎に得られたチルト角の平均値をチルト角とする。
【０１２４】
なお、重なり部を有する１組の画像の対応点からのみチルト角を算出するようにしてもよい。この場合には、図９のステップ１〜４の処理では、重なり部を有する１組の画像の対応点のみを抽出すればよい。
【０１２５】
〔４．６〕チルト角を考慮した円筒面投影変換処理（ステップ６）の説明
ステップ６では、各撮影画像を円筒面に投影する。カメラの焦点距離をｆ画素とすると、図１４に示すように、視点ＯからΨ度チルトして撮影した画像Ｉ上の点Ｐ（ｕ，ｖ）は、次式（２５）を用いて、Ｘ−Ｚ平面に垂直で半径ｆの円筒面上の点Ｃ＝（ｘ，ｙ）に投影される。なお、円筒面上の座標系ｘｙは、Ｚ軸と円筒面との交点を原点とし、垂直方向にｙ座標を、円筒面の周方向にｘ軸を取った座標系である。また、ｘ，ｙの各座標値の単位は画素単位である。
【０１２６】
【数２５】
2004046573.tif 000027
【０１２７】
ステップ６では、上記式（２５）を用いて、各撮影画像を円筒面画像に変換する。
【０１２８】
〔４．７〕画像接合処理（ステップ７）の説明
ステップ７では、円筒面投影変換された各画像が接合される。
【０１２９】
図１５は、画像接合処理手順を示している。
説明の便宜上、ここでは、重なり部を有する２枚の円筒面投影変換画像（第１画像Ｂ１及び第２画像Ｂ２）を接合する場合について説明する。
【０１３０】
第１画像Ｂ１と第２画像Ｂ２との重なり部の抽出処理が行われる（ステップ２１）。この重なり部の抽出処理は、図９のステップ２の処理と同様であるので、その説明を省略する。
【０１３１】
次に、特徴点抽出処理が行われる（ステップ２２）。つまり、第１画像Ｂ１における第２画像Ｂ２との重なり部分から、追跡に有効な複数の部分画像（矩形領域）が特徴点として抽出される。この特徴点抽出処理は、図９のステップ３の処理と同様であるので、その説明を省略する。
【０１３２】
次に、特徴点追跡処理が行われる（ステップ２３）。つまり、抽出された第１画像Ｂ１上の特徴点に対する第２画像Ｂ２上の位置が追跡される。この特徴点追跡処理は、図９のステップ４の処理と同様であるので、その説明を省略する。
【０１３３】
次に、幾何学変換行列の算出処理が行なわれる（ステップ２４）。第１画像Ｂ１に第２画像Ｂ２を接合する場合には、上記ステップ２３で求められた第１画像Ｂ１と第２画像Ｂ２との対応点の座標に基づいて、第２画像Ｂ２内の各画素の座標を第１画像Ｂ１上の座標に変換するための幾何変換行列（幾何変換係数）を算出する。
【０１３４】
なお、第２画像Ｂ２に第１画像Ｂ１を接合する場合には、第１画像Ｂ１内の各画素の座標を第２画像Ｂ２上の座標に変換するための幾何変換行列（幾何変換係数）を算出すればよい。
【０１３５】
幾何変換行列としては、２次元平行移動のための行列または２次元剛体変換行列が用いられる。
【０１３６】
この段階では、第２画像Ｂ２内の各画素の座標を第１画像Ｂ１上の座標に変換しないが、第２画像Ｂ２内の各画素の座標（ｘ，ｙ，１）を第１画像Ｂ１の座標（ｘ’，ｙ’，１）に変換する式は、次式（２６）または次式（２７）で表される。
【０１３７】
【数２６】
2004046573.tif 000028
【０１３８】
【数２７】
2004046573.tif 000029
【０１３９】
上記式（２６）は、幾何変換行列として２次元平行移動のための行列を用いた場合の変換式を示し、上記式（２７）は、幾何変換行列として２次元剛体変換行列を用いた場合の変換式を示している。
【０１４０】
上記式（２６）において、（ｔｘ　，ｔｙ　）は平行移動量を示し、ｋは零でない任意の定数を示している。また、上記式（２７）において、（ｔｘ　，ｔｙ　）は平行移動量を示し、θは回転角を示し、ｋは零でない任意の定数を示している。
【０１４１】
次に、画素値調合の重み（画素値調合係数）を算出する（ステップ２５）。
【０１４２】
この実施の形態のように、１台のカメラをチルト角を一定してパンさせて複数の画像を撮像した場合には、自動露出調整機能が動作することにより、画像間の明るさが異なってしまう。
【０１４３】
そこで、接合されるべき第１画像Ｂ１と第２画像Ｂ２とが重なる領域では、各画像の画素値を調合する処理を後処理として行なう。画素値調合処理は、両画像の重なり領域の画素値を画像間で滑らかに変化させることにより、明るさの差異を目立たなくさせる処理である。また、この処理は、同時に、幾何変換による合成誤差も目立たなくさせることができる。
【０１４４】
この実施の形態では、各画像の重心からの距離に基づいて両画像の重なり領域の画素値を調合する。
【０１４５】
まず、第１画像の重心位置Ｇ１　を求める。また、第２画像を第１画像に接合した合成画像上で、第２画像の重心位置Ｇ２　を、接合のために使用される変換式（上記式（２６）式または上記式（２７）式）に基づいて求める。
【０１４６】
次に、第１画像の４頂点（４隅）の位置をそれぞれ求める。また、第２画像を第１画像に接合した合成画像上で、第２画像の４頂点の位置を、接合に使用される変換式（上記式（２６）式または上記式（２７）式）に基づいて求める。これにより、合成画像上における第１画像の４頂点の位置および第２画像の４頂点の位置が分かるので、合成画像上において第１画像と第２画像との重なり領域を抽出する。
【０１４７】
合成画像上での第１画像と第２画像との重なり領域内の各画素（画素値調合対象画素）から第１画像と第２画像のそれぞれの画像の重心位置Ｇ１　、Ｇ２　までの距離に応じて、当該画素値調合対象画素の画素値の重みを決定する。この重みは、画素値調合対象画素から重心までの距離が近い方の画像の画素値を大きくするように決定する。
【０１４８】
つまり、合成画像上での第１画像と第２画像との重なり領域内の画素値調合対象画素の座標から、合成画像上での第１画像の重心位置Ｇ１　までのユークリッド距離をｄ１　、合成画像上での第２画像の重心位置Ｇ２　までのユークリッド距離をｄ２　とすると、画素値調合対象画素における第１画像の画素値に対する重みｗ１　および第２画像の画素値に対する重みｗ２　は、次式（２８）で表される。
【０１４９】
【数２８】
2004046573.tif 000030
【０１５０】
各画素値調合対象画素に対して式（２５）に基づいて重みｗ１　およびｗ２　を算出し、画素値調合対象画素に対する重みｗ１　およびｗ２　のテーブル（重みテーブル）を作成する。
【０１５１】
次に、画像を合成する（ステップ２６）。
【０１５２】
つまり、上記式（２６）または上記式（２７）を用いて第２画像を第１画像の座標系に描画する。
【０１５３】
この際、ステップ２５で得られた重みテーブルに基づいて、第１画像と第２画像の重なり領域の画素値Ｐ１　、Ｐ２　を画素調合する。すなわち、次式（２９）により、各画素値調合対象画素の画素値Ｐｉ　を算出する。
【０１５４】
【数２９】
2004046573.tif 000031
【０１５５】
〔２〕画像接合処理の変形例の説明
【０１５６】
図９のステップ７の画像接合処理の変形例について説明する。
【０１５７】
ここでは、図９のステップ１〜４によって、重なり部を有する画像の各組み合わせ毎に、それらに対応する特徴点が抽出されているものとする。
【０１５８】
図１６は、画像接合処理の変形例の処理手順を示している。
【０１５９】
上述したように、図１５で示した画像接合処理では、第１の円筒投影変換画像Ｂ１と第２の円筒投影変換画像Ｂ２との対応点を、第１の円筒投影変換画像Ｂ１と第２の円筒投影変換画像Ｂ２に対して、ステップ２１の重なり部抽出処理、ステップ２２の特徴点抽出処理およびステップ２３の特徴点追跡処理を行なうことによって求めている。そして、得られた第１の円筒投影変換画像Ｂ１と第２の円筒投影変換画像Ｂ２との対応点の座標に基づいて、幾何変換行列を算出している。
【０１６０】
これに対して、この変形例では、チルト角を算出するために図９のステップ　１〜４によって求められた対応点（特徴点）を利用して、第１の円筒投影変換画像Ｂ１と第２の円筒投影変換画像Ｂ２との対応点を求める。
【０１６１】
つまり、チルト角を算出するために求められた第１画像Ａ１と第２画像Ａ２との対応点（第１画像Ａ１上の特徴点と第２画像上の特徴点）を、上記式（２５）を用いて円筒面上の座標に変換する（ステップ３１）。そして、得られた第１画像Ａ１上の特徴点の円筒面上の座標と、それに対応する第２画像Ａ２上の特徴点の円筒面上の座標とを、第１の円筒投影変換画像Ｂ１と第２の円筒投影変換画像Ｂ２との対応点として、幾何変換行列を求める（ステップ３２）。
【０１６２】
この後、図９のステップ６で得られた円筒投影変換画像に対して、次のような処理が行なわれる。つまり、図１５のステップ２５と同様な方法により、画素値調合の重みを算出し（ステップ３３）、最後に図１５のステップ２６と同様な方法により、画像を合成する（ステップ３４）。
【０１６３】
【発明の効果】
この発明によれば、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像に基づいて、チルト角を自動的に算出できるようになる。
【０１６４】
また、この発明によれば、カメラを一定角度でチルトした状態でカメラをパンさせて撮影した複数の画像から円筒パノラマ画像を合成できるようになる。
【図面の簡単な説明】
【図１】図１は、階層的推定法を説明するための図であって、階層１の画像を示す模式図である。
【図２】図２は、階層的推定法を説明するための図であって、階層２の画像を示す模式図である。
【図３】図３は、階層的推定法を説明するための図であって、階層３の画像を示す模式図である。
【図４】図４は、階層的推定法を説明するための図であって、階層４の画像を示す模式図である。
【図５】図５は、実施の形態で採用されるオプティカルフロー推定方法において行われる膨張処理を説明するための模式図である。
【図６】図６（ａ）は、ある階層の画像に対してしきい値処理されたオプティカルフローの例を示す模式図であり、図６（ｂ）は、補完後のオプティカルフローを示す模式図である。
【図７】カメラ１を垂直方向に回転（チルト）したまま水平回転（パン）して撮影対象物を撮影する様子を示す模式図である。
【図８】円筒パノラマ画像合成を実現するためのパーソナルコンピュータの構成を示すブロック図である。
【図９】ＣＰＵ１１によって行なわれる円筒パノラマ画像合成処理の全体的な処理手順を示すフローチャートである。
【図１０】３次元空間上の点Ｐを、１台のカメラでチルト角Ψを一定としてかつパン角度φだけ離れた２つの位置において撮像したときの、画像の配置例を示す模式図である。
【図１１】第１のカメラ座標系と世界座標系との関係を示す模式図である。
【図１２】第２のカメラ座標系と世界座標系との関係を示す模式図である。
【図１３】チルト角の計算処理手順を示すフローチャートである。
【図１４】チルト角を考慮した円筒面投影変換処理を説明するための模式図である。
【図１５】画像接合処理手順を示すフローチャートである。
【図１６】画像接合処理の変形例の処理手順を示すフローチャートである。
【符号の説明】
１１　ＣＰＵ
１２　メモリ
１３　ハードディスク
【図１】
2004046573.tif 000032
【図２】
2004046573.tif 000033
【図３】
2004046573.tif 000034
【図４】
2004046573.tif 000035
【図５】
2004046573.tif 000036
【図６】
2004046573.tif 000037
【図７】
2004046573.tif 000038
【図８】
2004046573.tif 000039
【図９】
2004046573.tif 000040
【図１０】
2004046573.tif 000041
【図１１】
2004046573.tif 000042
【図１２】
2004046573.tif 000043
【図１３】
2004046573.tif 000044
【図１４】
2004046573.tif 000045
【図１５】
2004046573.tif 000046
【図１６】
2004046573.tif 000047