JP3618649B2 - An extended image matching method between images using an indefinite window - Google Patents

Description

また、Ｐが投影される写真座標（投影中心の像（主点）を原点とする２次元座標）ｐ（Ｐ，Ｅ）＝^ｔ［ｘ ｙ］は次の式で表わされる。
【００４０】
【数２】

ここでｃは航空写真の焦点距離である。式（２４）は、カメラの外部標定要素と地上点の３次元座標がわかれば、それに対応する写真座標がわかることを意味する。写真座標ｐは、物理的にフィルム面やＣＣＤ素子の上の実寸の座標であるが、フィルムを読みとる際のスキャナのひずみや非線形やＣＣＤ素子の配置誤差が無視できるとすれば、画像座標ｐｉｍｇとの関係は２次元の射影変換Ｈｉｎで表すことができる。以上より、ＥおよびＨｉｎを知ることができれば、地上点の３次元座標から画像座標に変換することが可能となり、４つの座標変換関数ＦＸ１（ｘ，ｙ，ｚ）、ＦＹ１（ｘ，ｙ，ｚ）、ＦＸ２（ｘ，ｙ，ｚ）、ＦＹ２（ｘ，ｙ，ｚ）が求められたこととなる。
【００４１】
また、写真座標ｐ＝（ｘ，ｙ）であるとすると、これはカメラ座標ではｐ_ｒｐ＝（ｘ，ｙ，−ｃ）である。これにより、ｐ＝（ｘ，ｙ）に写っている点の標高Ｚ_ｐがあるとすると、この点の地上座標系における３次元座標Ｐ_ｐ＝（Ｘ，Ｙ，Ｚ_ｐ）は次の式で与えられる。
【００４２】
【数３】

画像座標ｐｉｍｇに逆射影変換Ｈｉｎ^−１を行った後、式（２５）を適用すれば、Ｚ＝Ｚ_ｐの拘束の元で画像座標から３次元座標への変換を行うことができる。つまり、座標変換関数ＧＸ１（ｘ，ｙ，ｚ）、ＧＹ１（ｘ，ｙ，ｚ）、ＧＸ２（ｘ，ｙ，ｚ）、ＧＹ２（ｘ，ｙ，ｚ）を求めることができる。
【００４３】
一方、数値地図入力部４は、既存の数値地図を入力する。この数値地図は地物（たとえば建物データ）をベクター形式（座標点列データ）で保持しているものとする。数値地図が画像データの場合は、ラスター／ベクター変換によりベクター形式に直す。
【００４４】
次に、処理対象地物データ抽出部５は、オペレータによる指示に従って数値地図データより、標高付けの対象となる地物データを抽出する。数値地図では一般に対象の種類（建物・道路等）の違いがコード付けされているので、抽出の際は対象となるコードの図形のみを抽出すればよい。
【００４５】
そして、抽出した地物Ｂｕをポリゴン点列Ｐｌ_ｉ（閉多角形、一般に最初と最後が同じ点となる点列（ｐ１，ｐ２，…，ｐｎ，ｐ１）で表される）に変換する。 このとき、対象とする地物Ｂｕが建物データの場合は普通ポリゴン点列になっているものが多いので問題はない。ポリゴン点列になっていない地物の場合は、ＣＡＤソフト等で編集してあらかじめポリゴン点列化しておく。
【００４６】
また、画像上地物計測部６は、拡張ステレオモデルＤｉの画像上に建物等の地物Ｂｕがあるが、数値地図データに対象とする地物形状が記述されていない場合には、拡張ステレオ画像Ｄｉ（輝度画像Ｉ１、Ｉ２）のどちらかを画面に表示し、その上でマウス操作により地物形状を入力する。この場合、地物形状はポリゴン点列Ｐｌ_ｉとし、画像座標を列挙した閉図形（不定形窓）として入力する。例えば、図２においては、輝度画像Ｉ１からＢ１を選択し、このＢ１の不定形窓を得る。
【００４７】
次に、標高探索部７は、標高の探索では、画像上地物計測部６で得られた各地物ポリゴン点列Ｐｌ_ｉ（ｉ＝１，…，ｎ_ｉ）について、以下の処理を繰り返す。
【００４８】
標高探索部７における標高探索範囲および探索間隔設定部７ａは、標高の探索範囲（Ｈｍｉｎ，Ｈｍａｘ）および探索間隔ｄＨを設定する。
【００４９】
このＨｍｉｎには、その地域の地表面の標高もしくはそれより小さい値を、Ｈｍａｘにはその地図の建物階数等から予想される地物の最大標高もしくはそれ以上の値を与える。
【００５０】
地形モデル（国土地理院発行の数値地図５０ｍメッシュ（標高）など）が利用できる場合は各地物位置での地表面の標高を推定することができる。
【００５１】
また、利用できない場合は、すべての地物について適当に同じＨｍｉｎとＨｍａｘを設定すればよい。またｄｚには、拡張ステレオモデルから測定される高さ方向の解像度を与えればよい。
【００５２】
たとえば、各画像間の撮影位置間隔（基線長）をＢ、拡張ステレオモデル撮影時の対地高度の平均をＨ、２つの画像の１ピクセルに対応する地上解像度平均をＤとすれば、ｄＨは次の式で与えられる。
【００５３】
【数４】

またこの時、探索標高群Ｈ_ｋを次のように定義する。
【００５４】
Ｈ_ｋ＝Ｈｍｉｎ＋ｄＨ×ｋ（ｋ＝０，１，２，…ｎ_ｋ） 式（１０）
ただし、ｎ_ｋはＨ_ｋ≦Ｈｍａｘとなる最大のｋである。
【００５５】
そして、不定形窓定義部７ｂは、地物データのそれぞれのポリゴン点列の内部をステレオマッチングの際の不定形窓として定義する。具体的には、図８に示すように、地物ポリゴン点列内に等間隔のメッシュ（ｄｘ、ｄｙ）を発生させ、これを不定形窓Ｐｌ_ｉとする。
【００５６】
この処理は、ポリゴン点列Ｐｌ_ｉが数値地図（座標が地上座標）から得られたか、それとも画像上の計測結果（座標が画像座標）から得られたかによって処理が異なる。
【００５７】
・数値地図から得られたポリゴン点列の時
地物データが数値地図起源の時は、図８に示すように、地物データＢｕのそれぞれのポリゴン点列の周辺に等間隔のメッシュ点群ＭＰｉを発生させる。メッシュは任意の地上座標に原点（ｏｘ，ｏｙ）を設定し、ｘ方向およびｙ方向にｄｘおよびｄｙの間隔で発生させ、Ｍ（ｏｘ，ｏｙ，ｄｘ，ｄｙ）とする。
【００５８】
なお、ｏｘおよびｏｙは任意に取ってかまわない。またｄｘおよびｄｙは、デジタル航空写真の１ピクセル〜数ピクセルの地上解像度に設定する。
【００５９】
ここでｉ番目の地物Ｆ_ｉを表すポリゴン点列Ｐｌ_ｉで囲まれる領域Ｐｌｉｎ_ｉ内部に存在するメッシュの交点群ＭＰ_ｉを次のように定義する。
【００６０】
ＭＰｉ＝Ｐｌｉｎ_ｉ∩Ｍ（ｏｘ，ｏｙ，ｄｘ，ｄｙ） 式（１１）
このようにして得た不定形窓Ｐｌｉを用いて、例えば図１では、不定形窓Ｐｌの高さを変化させ輝度画像Ｉ１、輝度画像Ｉ２の投影位置を求め、この投影位置における両方の画像の輝度値を比較する。
【００６１】
そして、輝度値の評価関数が最も一致（近似）したときの標高値Ｈを地物Ｂ１、Ｂ２の標高値ｚと決定する。
【００６２】
・画像上の計測結果から得られたポリゴン点列の時は、不定形窓定義部７ｃ（第２の不定形窓定義部ともいう）は、Ｐｌ_ｉが画像上の計測結果の時は、評価関数計算ループ部７ｄの中できめられる標高Ｈ_ｋにあわせてメッシュを作成する。
【００６３】
今、Ｐｌ_ｉが画像１上で計測されたと仮定し、ｉ番目の画像座標を（ｐｘ_ｊ，ｐｙ_ｊ）とすると、この点が標高Ｈ_ｋにあったときの地上
座標（Ｘ_ｊ，Ｙ_ｊ）は、式（５）および式（６）を用いて、
Ｘ_ｊ＝ＧＸ１（ｐｘ_ｊ，ｐｙ_ｊ，Ｈ_ｋ） 式（１２）
Ｙ_ｊ＝ＧＹ１（ｐｘ_ｊ，ｐｙ_ｊ，Ｈ_ｋ） 式（１３）
で表される。こうしてすべてのポリゴン点列上の点を地上座標に戻した後、数値地図の場合と同様に変換して、メッシュ交点群ＭＰ_ｉｋを得る。Ｐｌ_ｉが画像上で計測された場合も式（１２）および式（１３）によって同様にメッシュ点群ＭＰ_ｉｋを得る。
【００６４】
このような不定形窓を用いることによって、例えば図２に示すように、航空機によって得た輝度画像Ｉ１を画面に表示し、この画像上でマウス操作等で座標が既知の地物Ｂ１の輪郭に相当する閉領域（不定形窓Ｐｌｈ：ポリゴン）を与える。そして、この不定形窓Ｐｌｈを地物Ｂｕにおいて、標高Ｈ（ｚ）を変化させる。
【００６５】
前述の評価関数計算ループ部７ｄは、すべてのＨ_ｋ（ｋ＝０，１，２，…，ｎ_ｋ）について、評価値Ｅ（Ｈ_ｋ）を計算する。
【００６６】
【数５】

ここで、Ｅｋｊは、メッシュ交点群ＭＰｉのｊ番目の点（ｘｊ，ｙｊ）が標高Ｈｋの高さにあると仮定した時の評価関数値である。評価値Ｅ（Ｈｋ）は、拡張ステレオモデル上の画像の画素値の一致の程度を評価し、一致が高くなる時に小さくなるものとする。
【００６７】
すなわち、標高Ｈｋの高さにあるメッシュ交点群ＭＰｉのｊ番目の点（ｘｊ，ｙｊ，Ｈｋ）が二つの画像Ｉ１およびＩ２の上へ、各々式（１）〜式（４）により投影されるとき、その位置にある画素値Ｖ１ｋｊおよびＶ２ｋｊを

ここで、Ｖ１_ｋｊ（ｘ_ｊ，ｙ_ｊ）、Ｖ２_ｋｊ（ｘ_ｊ，ｙ_ｊ）は投影したときの画素値であり、（ｘ_ｊ，ｙ_ｊ，Ｈ_ｋ）はＨにあるときの不定形窓のポリゴン点群を示す。
【００６８】
上記（１５）式のように表せば、画素値の差の自乗Ｅｄｉｆの総和や、画素値の相関係数に（−１）をかけた値Ｅｃｏｒなどを用いることができる。
【００６９】
例えば、Ｅｄｉｆ_ｋは次の式で表される。
【００７０】
【数６】

また、Ｅｃｏｒ_ｋは次の式で表される。
【００７１】
【数７】

そして、地物標高決定部７ｅは、ｉ番目の地物Ｆ_ｉの標高Ｈ_Ｆｉを、評価関数Ｅ（Ｈ_ｋ）を最小にする高さＨ_ｋｍｉｎとする。すなわち、

次に、地物高さ計算部７ｆは、地形モデルがある場合は、２次元のポリゴン点列Ｐｌ_ｉの位置の地表面高さを標高値Ｈ_Ｆｉから引くことにより、地物比高Ｈ_Ｒｉとして出力する。また、後で信頼性の解析に用いるため、標高値Ｈ_Ｆｉでの評価関数値Ｅ_ｉも地物の属性情報として出力する。
【００７２】
そして、３次元地物データ出力部７ｅは、２次元のポリゴン点列Ｐｌ_ｉのｚ座標に標高値Ｈ_Ｆｉを付け加え、更に地表面比高Ｈ_Ｒｉを属性データとして付け加えた３次元地物データＰｌ３Ｄ_ｉとして出力する。
【００７３】
次に、３次元地図出力部８は、すべてのＰｌ３Ｄ_ｉをあわせて、標高付きの３次元地図として出力する。
【００７４】
次に、標高推定の信頼性の評価について説明を補充する。
【００７５】
次のような場合は、標高推定の信頼性の推定方法として、次の様な例が考えられる。
【００７６】
・評価関数値Ｅ_ｉが大きい時。つまり、あるしきい値をＴ２として、
Ｅ_ｉ≧Ｔ２ 式（１９）
Ｔ２の値の決定方法としては、評価関数値Ｅ_ｉの平均値より、１〜２σ大きい値を設定することが合理的である。
【００７７】
・求めた高さＨ_Ｆｉが探索範囲の両端になる場合。つまり、
Ｈ_Ｆｉ＝Ｈｍｉｎ ｏｒ Ｈ_Ｆｉ＝Ｈｍａｘ 式（２０）
この場合は、探索範囲が正しく設定されなかった可能性がある。
【００７８】
また、地物の滅失や変更の推定する場合は、次のような場合、地物の滅失や変更の可能性が高いと推定できる。
【００７９】
・地表面比高が小さい時。つまり、０に近いしきい値をＴ１として、
｜Ｈ_Ｒｉ｜≦Ｔ１ 式（２１）
また、信頼性が低くなる原因として、地物の高さ、大きさ等に変更が起こっている場合がある。
【００８０】
例えば、図９に示すように、Ａ年度の撮影時においては地物の標高がＨｏｒｇであっても、数年後（経年変化）に撮影した場合は、標高がＨｏｒｇからＨｎｅｗに変化している場合がある。
【００８１】
また、建物の大きさ（平面形状）も時間の経過に伴って変化している場合もある（Ｃ年度撮影）。
【００８２】
標高の既知の地物の異動の推定の場合は、すでに標高が既知の地物の異動は次のように推定する。
【００８３】
地物の標高をＨｏｒｇとし、標高の探索範囲を
Ｈｏｒｇ−ｄＨ≦Ｈ_ｋ≦Ｈｏｒｇ＋ｄＨにとる。
【００８４】
経年変化がある場合は、推定した標高がＨｏｒｇと大きく異なるか、評価関数値から標高推定の信頼性が低くなる。
【００８５】
また、旧・新の撮影時期の異なる画像を拡張ステレオモデルに構成し、標高推定を行い、標高推定の信頼性が低くなる場合に経年変化があると判定することも可能である。
【００８６】
また、例えば図１０に示すように、地物が複雑な形状であったり、実際には一定の高さではない場合にも対応可能である。建物について実際に航空写真測量やレーザスキャナで計測した結果と比べてみると、建物の最頻値標高に近い値を示すことが実験で確認されている。
【００８７】
また地物高さが０に近いときは、自動的に地物の滅失箇所候補として抽出し、修正の効率化を図ることができる。
【００８８】
＜実施の形態４＞
次に地図の更新について図１１のフローチャートを用いて説明する。本説明では予め３次元地図（旧）と、２次元地図（旧）と、新規撮影の拡張ステレオモデルとがディスクに記憶されているとする。
【００８９】
このような状態において、二次元地図（旧）を読み込むと共に、拡張ステレオモデルとを読み込み、上記説明のように、既知の地物の不定形窓に高さを与え、これを変化させて拡張ステレオモデルの両方の画像上の投影位置を求め、この投影位置における不定形窓で区切られた領域の輝度値をそれぞれ求めて輝度値を比較し、比較結果が最も一致するときの不定形窓の高さを拡張ステレオモデルにおける地物の高さとする拡張ステレオマッチングによる標高付与処理を行い（Ｓ１）、この３次元地図データを途中成果として記憶する（Ｓ２）。
【００９０】
そして、３次元地図（旧）と、この途中成果の３次元地図データと新規の拡張ステレオモデルとから変化を求める拡張ステレオマッチングによる変化抽出処理を行い（Ｓ３）、これを写真上（拡張ステレオモデル）に割り付ける前の途中成果の３次元地図として記憶する（Ｓ４）。
【００９１】
次に、写真上での編集処理を行う（Ｓ５）。これは、デジタル写真測量ワークステーションを用いて行うのが望ましい。
【００９２】
この編集処理は、ステップＳ４で得た途中成果の地物の３次元データを新規の拡張ステレオモデルの上にオーバレイしてモニターに表示し、３次元データと新規に撮影した拡張ステレオモデルとを比較することによってオペレータがマニュアルで編集を行う。ここでは拡張ステレオモデルとしては新規に撮影した一般的なステレオ画像を用いてもよいが、旧２次元地図を作成したときの画像と新規に撮影した画像とから構成されるステレオモデルを使用してもよい。
【００９３】
そして、この編集処理で得られた３次元データに旧の３次元データを更新する処理を行う（Ｓ６）。
【００９４】
つまり、この一連の処理によって、標高を持たない数値地図データを容易に３次元データ化がなされ、かつ既存の地図の更新ができることになる。
【００９５】
【発明の効果】
以上のように本発明によれば、２次元地図を効率的に３次元地図化することができる。特に航空写真を新しくとる必要がない場合は、コンピュータグラフィックスや都市３次元シミュレーションなどの用途に用いての都市３次元モデルを既存地図から廉価に作成することが可能となる。
【００９６】
また、数値地図を３次元化することにより、２次元地図の経年変化修正業務において、新規の撮影した地図をステレオ航空写真上に重ねることができ、修正が容易になる。
【図面の簡単な説明】
【図１】本実施の形態１の不定形窓を用いた拡張イメージング方法の概略構成図である。
【図２】実施の形態２の不定形窓を用いた拡張イメージング方法の概略構成図である。
【図３】不定形窓の高さ方向の変更を説明する説明図である。
【図４】実施の形態３の比高を説明する説明図である。
【図５】各実施の形態の具体的な構成図である。
【図６】標高探索部の構成図である。
【図７】エリアセンサの幾何学的特性を説明する説明図である。
【図８】地物データのメッシュ化の説明図である。
【図９】撮影時期の異なるときの地物データの変化を説明する説明図である。
【図１０】複雑な形状の建物の標高の求め方を説明する説明図である。
【図１１】実施の形態３の地図の更新を説明するフローチャートである。
【符号の説明】
１ デジタル画像入力部
２ ハードディスク
３ 標定要素算出部
４ 数値地図入力部
５ 処理対象地物データ抽出部
６ 画像上地物計測部
７ 標高探索部[0001]
BACKGROUND OF THE INVENTION
Using two or more types of images taken from the same area, such as aerial photographs and satellite images, two-dimensional numerical map data without elevation is converted into three-dimensional data (altitude values are assigned), or existing numerical values The present invention relates to a method for confirming the loss of a feature described in map data.
[0002]
[Prior art]
In conventional photogrammetry, the object to be measured is taken from two different positions (or more) using the same type of camera, the corresponding points are obtained on each image, and the three-dimensional space in the three-dimensional space is obtained based on the principle of triangulation Work to determine the position.
[0003]
A method for automatically performing this work is to pay attention to one point in one image constituting a stereo image or a region divided on a regular grid, or a line structure (edge), and respond to it. A stereo matching method is generally used in which a point or region or a line structure is searched for and matched with another image.
[0004]
[Problems to be solved by the invention]
However, with this method, even if there is a region whose shape is known in advance, the shape cannot be used as a search condition because, for example, the photographing altitude and the photographing date are different.
[0005]
For this reason, matching using this shape as a constraint condition cannot be performed, and the elevation cannot be determined efficiently. As a result, the three-dimensional data of the feature cannot be obtained.
[0006]
This is because, for example, when an altitude is given to a feature described in the existing numerical map data, such as a building, and data is to be converted into three-dimensional data, the points and lines measured by the stereo image are edited again. Otherwise, there is a problem that 3D numerical map data cannot be generated.
[0007]
Recently, high-resolution commercial satellites have been put to practical use, and attempts have been made to acquire digital image data by installing line scanners on airplanes. The situation of individual acquisition is beginning to spread.
[0008]
That is, unlike the conventional photogrammetry, stereo images are not intended in advance, and images obtained by photographing the same region from different positions are being acquired individually.
[0009]
In the information-oriented society toward the 21st century, with the spread of geographic information systems (GIS), it is necessary to maintain and manage data by updating numerical map data that is the information base in a timely manner.
[0010]
However, even when trying to obtain the latest numerical map in a timely manner, the timely update is difficult and the cost cannot be reduced by the conventional data update by taking a stereo aerial photograph once every several years.
[0011]
Therefore, it becomes useful to perform measurement equivalent to conventional photogrammetry or update a map using various images acquired by individually capturing the same area in a timely manner. However, in reality there is no such system.
[0012]
The present invention has been made in order to solve the above-described problems. By easily converting numerical map data having no altitude into three-dimensional data using two or more kinds of images obtained by photographing the same area. An object is to obtain an extended image matching method between images that can update an existing map.
[0013]
[Means for Solving the Problems]
The present invention relates to an imaging matching method using a combination of two types of images having different resolutions in the same region, and a step of performing an orientation operation on the two types of images and generating an extended stereo model related by a function; A polygon point sequence of known features existing in common in both images of the extended stereo model is used as a standard irregular window, and a height (z) is given to each irregular window. A step of changing the height (z), a step of obtaining a projection position on both images of the extended stereo model corresponding to the height of the irregular window, and a step of the extended stereo A step of comparing luminance values of regions divided by the irregular window at the projection positions of both images of the model, and the height (z) of the irregular window when the comparison result is the closest approximation And outputting as the height of the feature of the images of Dell, and summarized in that comprising the step of determining the three-dimensional coordinates of the feature by using the output height (z).
[0014]
DETAILED DESCRIPTION OF THE INVENTION
<Embodiment 1>
In the first embodiment, as shown in FIG. 1, any combination of an aerial photograph image, satellite photograph image, and line scanner image of the same area acquired by an aircraft (including a helicopter) and a satellite [satellite image (photograph (Including a combination of a line scanner image) and an aerial image (photo or line scanner), a combination of aerial photographs of different altitudes or a combination of different altitude scanner images, etc.] The two resolutions are combined (based on the photographing altitude and scale) and the orientation work is performed on these images.
[0015]
Then, a function for projecting the three-dimensional coordinates on the ground onto the two-dimensional coordinates on each image is determined for each image. A combination of images in which two images are related by a function for projecting the three-dimensional coordinates of the ground onto the respective image coordinates is referred to as an extended stereo model Di (luminance image I1, luminance image I2).
[0016]
When such an extended stereo model Di is used, for example, an altitude can be given to a feature on an existing numerical map described two-dimensionally without having an altitude value, and can be converted into three-dimensional data.
[0017]
That is, the feature Bu (known) having an arbitrary shape for which the altitude is to be obtained is defined as the reference indefinite window Pl (polygon array) of the closed region, and has an arbitrary altitude ground (height direction).
[0018]
Then, the positions of the luminance images I1 and I2 are obtained by using the above-described functions for the luminance image I1 and the luminance image I2 while changing the height of the irregular window Pl, and the irregular window Pl at this position is projected. The luminance values in the luminance images I1 and I2 are compared, and the height (elevation value H) of the amorphous P1 when it is determined that the evaluation function matches most closely (approximate) the height of the feature in the luminance image. decide. FIG. 1 shows an example in which the highest value is obtained when the altitude value is H3.
[0019]
Further, in FIG. 1, for easy understanding, the shape of the irregular window Pl when the height of the irregular window Pl is changed and projected onto the luminance images I1 and I2 is changed.
[0020]
<Embodiment 2>
In addition, a new building that is not on the map may be added as three-dimensional data in the map update. In this case, in the conventional method, the operator must convert the data into three-dimensional data while observing the stereo, but the operator cannot perform the stereo observation with a combination of different types of images, and there is no such system.
[0021]
However, when the extended stereo model Di is used, the luminance image I1 and the luminance image I2 are related by a function in the extended stereo model Di, and therefore, as shown in FIG. A closed region B1 corresponding to the contour of the feature Bu is given on the screen by operating the mouse on the image. Then, an irregular window Plh is generated at a position where the projection coincides with the closed region, and the altitude H (z) is changed.
[0022]
For example, when the elevation is changed in the order of 30 m, 80 m, and 150 m, the projection positions of the luminance images I1 and I2 of the irregular window Plh in the luminance images I1 and I2 when changing to 30 m, 80 m, and 150 m are obtained each time. The luminance values in the luminance images I1 and I2 when the irregular window Plh at the position is projected are compared, and the height of the irregular P1h (elevation value H) when it is determined to be the closest (approximate) by the evaluation function Is determined as the height of the feature in the luminance image. Therefore, the three-dimensional coordinates of the feature Bu can be estimated.
[0023]
Further, in FIG. 2, for easy understanding, the shape of the irregular window Plh when the height of the irregular window Plh is changed and projected onto the luminance images I1 and I2 is changed.
[0024]
That is, in the first and second embodiments, as shown in FIG. 3, in the luminance images I1 and I2 of the extended stereo model in which the above-described combination images are related by functions, By changing the height of the irregular window Plh obtained from the object Bu, the projection position in the luminance images I1 and I2 is obtained, and the evaluation function of the luminance value divided by the irregular window at this projection position is matched (approximate). The height of the feature Bu is estimated.
[0025]
<Embodiment 3>
When a numerical terrain model can be used, as shown in FIG. 4, the specific height of the feature is also obtained by subtracting the elevation at the feature position from the feature elevation z. Thus, an accurate building height can be obtained.
[0026]
Further, for example, let us focus on a building described in an existing numerical map. If this specific height is determined to be a value close to 0, it can be estimated that the building has been lost. For buildings that have already been given three-dimensional coordinates, it is also possible to estimate the loss or change of the feature by comparing the image of the position where the feature is projected between the newly captured image and the old image. .
[0027]
Next, specific examples of the first, second, and third embodiments will be described below. FIG. 5 is a configuration diagram of the method for estimating the altitude and specific height of a feature according to the present invention. FIG. 6 is a schematic configuration diagram of the elevation search unit of FIG.
[0028]
For example, the digital image input unit 1 stores two images in the form of a digital image in a hard disk 2 of a computer that is an image processing apparatus. At this time, in the case of an image supplied as a film or photographic paper such as an aerial photograph, it is read by a scanner and converted into a digital image by an A / D converter. The pixel coordinates (x _{1 pix} , Y _{1 pix} ) And (x _{2 pix} , Y _{2 pix} ) For image brightness ₁ (X _{1 pix} , Y _{1 pix} ) And I ₂ (X _{2 pix} , Y _{2 pix} ).
[0029]
Also, the hard disk 2 described above includes a ground altitude H at the time of extended stereo model shooting, a shooting position interval (baseline length B) between images, a ground resolution average D corresponding to one pixel of two images, and a known Stores polygons (indefinite windows) of buildings.
[0030]
Next, the orientation element calculation unit 3 uses, for the two luminance images I1 and I2, a coordinate conversion function FX and an inverse coordinate conversion function function GX (generic name) for calculating image coordinates of each image on which the three-dimensional coordinates on the ground are projected. (Also referred to as orientation information). This corresponds to what is called orientation work in aerial photogrammetry. The orientation work method and required parameters differ depending on the image sensor (for example, the area sensor and the linear array sensor), but finally the pixel coordinates of the left and right images on which the point (x, y, z) in the three-dimensional space is projected (X _{1 pix} , Y _{1 pix} ) And (x _{2 pix} , Y _{2 pix} ) Can be determined.
[0031]
That is, the following four coordinate transformation functions FX1 (x, y, z), FY1 (x, y, z), FX2 (x, y, z), FY2 (x, y, z) are analyzed analytically or procedurally. Can be determined.
[0032]
x _{1 pix} = FX1 (x, y, z) Formula (1)
y _{1 pix} = FY1 (x, y, z) Equation (2)
x _{2 pix} = FX2 (x, y, z) Equation (3)
y _{2 pix} = FY2 (x, y, z) Equation (4)
Further, when equations (1) and (2) are solved analytically or procedurally with z = constant, the x-coordinate and y-coordinate on the ground when the point on each image is assumed to be at a certain altitude z are obtained. Can be sought.
[0033]
x = GX1 (x _{1 pix} , Y _{1 pix} , Z) Equation (5)
y = GY1 (x _{1 pix} , Y _{1 pix} , Z) Equation (6)
Similarly, when solving for Equation (3) and Equation (4),
x = GX2 (x _{2 pix} , Y _{2 pix} , Z) Equation (7)
y = GY2 (x _{2 pix} , Y _{2 pix} , Z) Equation (8)
It is.
[0034]
For area sensor images such as aerial photographs and area CCD sensors, the functions FX1 (x, y, z), FY1 (x, y, z), FX2 (x, y, z), FY2 (x, y , Z) are described in detail in FIG.
[0035]
It should be noted that the existing method is used for the orientation work method. For example, the area CCD sensor is described in detail in Document [1], and the linear center projection sensor is described in detail in Document [2].
[0036]
[1] Japan Photogrammetry Society: Analytical Photogrammetry Revised Edition, 1989. “2” Osamu Uchida: Research on automation of three-dimensional measurement using stereo satellite images, 1989.
[0037]
Further, for example, the coordinate conversion function of the area sensor image in FIG.
The coordinates of the projection center of the central projection image are P ₀ = ^t [X ₀ Y ₀ Z ₀ ], The rotation angle with respect to the ground coordinates of the camera is T = ^t [Ω φ κ], a combination of these (external orientation element) E = ^t [TT tP ₀ ] (Left shoulder t indicates a transposed matrix or transposed vector). At this time, the ground coordinate P = ^t [X y z] is the camera coordinate system (a three-dimensional coordinate system relative to the projection center) p _r Is expressed as follows.
[0038]
p _r = ^t [X _r y _r z _r ] = R (T) (RP ₀ Formula (22)
Here, R (T) is a rotation matrix defined by the following equation.
[0039]
[Expression 1]

In addition, the photographic coordinates (two-dimensional coordinates with the origin of the projection center image (principal point)) p (P, E) = ^t [Xy] is represented by the following equation.
[0040]
[Expression 2]

Here, c is the focal length of the aerial photograph. Equation (24) means that if the external orientation elements of the camera and the three-dimensional coordinates of the ground point are known, the corresponding photographic coordinates can be known. The photographic coordinate p is physically the actual coordinate on the film surface or the CCD element, but if the distortion of the scanner when reading the film, non-linearity or the arrangement error of the CCD element can be ignored, the image coordinate pimg Can be expressed by a two-dimensional projective transformation Hin. From the above, if E and Hin can be known, it is possible to convert the three-dimensional coordinates of the ground point into image coordinates, and four coordinate conversion functions FX1 (x, y, z), FY1 (x, y, z) ), FX2 (x, y, z), FY2 (x, y, z) are obtained.
[0041]
If the photo coordinates p = (x, y), this is p in camera coordinates. _rp = (X, y, -c). As a result, the altitude Z of the point shown in p = (x, y) _p If there is, the three-dimensional coordinates P of this point in the ground coordinate system _p = (X, Y, Z _p ) Is given by:
[0042]
[Equation 3]

Reverse projection transformation Hin to image coordinates img ^-1 And then applying equation (25), Z = Z _p Conversion from image coordinates to three-dimensional coordinates can be performed under the constraints. That is, the coordinate conversion functions GX1 (x, y, z), GY1 (x, y, z), GX2 (x, y, z), and GY2 (x, y, z) can be obtained.
[0043]
On the other hand, the numerical map input unit 4 inputs an existing numerical map. It is assumed that this numerical map holds features (for example, building data) in a vector format (coordinate point sequence data). When the numerical map is image data, it is converted into a vector format by raster / vector conversion.
[0044]
Next, the processing target feature data extracting unit 5 extracts feature data to be altitude-added from the numerical map data in accordance with an instruction from the operator. Since numerical maps are generally coded with different types of objects (buildings, roads, etc.), only the figure of the target code need be extracted.
[0045]
Then, the extracted feature Bu is converted to a polygon point sequence Pl _i (Closed polygon, generally represented by a sequence of points (p1, p2,..., Pn, p1) where the first and last points are the same). At this time, if the target feature Bu is building data, there are no problems because there are many normal polygon point sequences. In the case of a feature that is not in a polygon point sequence, it is edited with CAD software or the like to create a polygon point sequence in advance.
[0046]
Further, the image top feature measurement unit 6 has a feature Bu such as a building on the image of the extended stereo model Di, but when the target feature shape is not described in the numerical map data, the extended stereo model Di is used. One of the images Di (luminance images I1 and I2) is displayed on the screen, and a feature shape is input by operating the mouse on the screen. In this case, the feature shape is a polygon point sequence Pl _i And entered as a closed figure (indefinite window) enumerating the image coordinates. For example, in FIG. 2, B1 is selected from the luminance image I1, and an irregular window of this B1 is obtained.
[0047]
Next, the elevation search unit 7 performs the feature polygon point sequence Pl obtained by the on-image feature measurement unit 6 in the elevation search. _i (I = 1, ..., n _i ), The following processing is repeated.
[0048]
The altitude search range and search interval setting unit 7a in the altitude search unit 7 sets the altitude search range (Hmin, Hmax) and the search interval dH.
[0049]
This Hmin is given an altitude on the ground surface of the area or a value smaller than that, and Hmax is given a value higher than or equal to the maximum altitude of the feature expected from the number of building floors of the map.
[0050]
If a terrain model (such as a numerical map 50m mesh (elevation) published by the Geospatial Information Authority of Japan) is available, the altitude of the ground surface at each feature position can be estimated.
[0051]
In addition, if it cannot be used, the same Hmin and Hmax may be appropriately set for all the features. Further, dz may be given a resolution in the height direction measured from the extended stereo model.
[0052]
For example, if the shooting position interval (baseline length) between each image is B, the average of the ground altitude at the time of extended stereo model shooting is H, and the average of the ground resolution corresponding to one pixel of two images is D, then dH is Is given by
[0053]
[Expression 4]

At this time, search altitude group H _k Is defined as follows.
[0054]
H _k = Hmin + dH × k (k = 0, 1, 2,... N _k Formula (10)
Where n _k Is H _k The maximum k that satisfies ≦ Hmax.
[0055]
Then, the irregular window defining unit 7b defines the inside of each polygon point sequence of the feature data as an irregular window at the time of stereo matching. Specifically, as shown in FIG. 8, meshes (dx, dy) of equal intervals are generated in the feature polygon point sequence, and this is generated as an indefinite window Pl. _i And
[0056]
This processing is performed by the polygon point sequence Pl _i The processing differs depending on whether is obtained from a numerical map (coordinates are ground coordinates) or measurement results on images (coordinates are image coordinates).
[0057]
-For polygon point sequences obtained from numerical maps
When the feature data originates from a numerical map, as shown in FIG. 8, a mesh point group MPi with equal intervals is generated around each polygon point sequence of the feature data Bu. The mesh has an origin (ox, oy) set at an arbitrary ground coordinate, and is generated at intervals of dx and dy in the x and y directions, and is defined as M (ox, oy, dx, dy).
[0058]
Note that ox and oy may be arbitrarily selected. Dx and dy are set to a ground resolution of one pixel to several pixels in the digital aerial photograph.
[0059]
Where i-th feature F _i Polygon point sequence Pl _i Area Plin surrounded by _i Intersection group MP of mesh existing inside _i Is defined as follows.
[0060]
MPi = Plin _i ∩M (ox, oy, dx, dy) Equation (11)
For example, in FIG. 1, the height of the irregular window Pl is changed to obtain the projection positions of the luminance image I1 and the luminance image I2, and the projection positions of both images at the projection position are obtained. Compare brightness values.
[0061]
Then, the altitude value H when the evaluation functions of the luminance values are most matched (approximated) is determined as the altitude value z of the features B1 and B2.
[0062]
In the case of a polygon point sequence obtained from the measurement result on the image, the irregular window definition unit 7c (also referred to as a second irregular window definition unit) _i Is the measurement result on the image, the altitude H created in the evaluation function calculation loop section 7d _k Create a mesh to match.
[0063]
Pl now _i Is measured on image 1, and the i-th image coordinate is (px _j , Py _j ), This point is at altitude H _k When we were on the ground
Coordinates (X _j , Y _j ) Is obtained by using formula (5) and formula (6),
X _j = GX1 (px _j , Py _j , H _k Formula (12)
Y _j = GY1 (px _j , Py _j , H _k Formula (13)
It is represented by In this way, all the points on the polygon point sequence are returned to the ground coordinates, and then converted in the same manner as in the case of the numerical map, and the mesh intersection group MP _ik Get. Pl _i Is also measured on the image, the mesh point group MP is similarly expressed by the equations (12) and (13). _ik Get.
[0064]
By using such an irregular window, for example, as shown in FIG. 2, a luminance image I1 obtained by an aircraft is displayed on the screen, and on this image, the contour of the feature B1 whose coordinates are known by a mouse operation or the like is displayed. A corresponding closed region (indefinite window Plh: polygon) is given. Then, the altitude H (z) of the irregular window Plh is changed in the feature Bu.
[0065]
The above-described evaluation function calculation loop unit 7d _k (K = 0, 1, 2, ..., n _k ) For evaluation value E (H _k ).
[0066]
[Equation 5]

Here, Ekj is an evaluation function value when it is assumed that the j-th point (xj, yj) of the mesh intersection group MPi is at the height of the altitude Hk. The evaluation value E (Hk) is evaluated when the degree of matching of the pixel values of the image on the extended stereo model is evaluated, and becomes small when the matching becomes high.
[0067]
That is, the j-th point (xj, yj, Hk) of the mesh intersection point group MPi at the height of the altitude Hk is projected onto the two images I1 and I2 by the expressions (1) to (4), respectively. When the pixel values V1kj and V2kj at that position are

Where V1 _kj (X _j , Y _j ), V2 _kj (X _j , Y _j ) Is the pixel value when projected, and (x _j , Y _j , H _k ) Indicates the polygon point group of the irregular window when in H.
[0068]
If expressed as the above equation (15), the sum of the squares of pixel value differences Edif, the value Ecor obtained by multiplying the correlation coefficient of pixel values by (−1), or the like can be used.
[0069]
For example, Edif _k Is represented by the following equation.
[0070]
[Formula 6]

Ecor _k Is represented by the following equation.
[0071]
[Expression 7]

Then, the feature altitude determining unit 7e reads the i-th feature F _i Altitude H _Fi Is evaluated function E (H _k ) To minimize height H _kmin And That is,

Next, the feature height calculation unit 7f, when there is a topographic model, the two-dimensional polygon point sequence Pl _i The height of the ground surface at the altitude value H _Fi Feature height H by subtracting from _Ri Output as. Also, the altitude value H is used for later analysis of reliability. _Fi Evaluation function value E _i Are also output as feature attribute information.
[0072]
Then, the three-dimensional feature data output unit 7e outputs a two-dimensional polygon point sequence Pl _i Elevation value H in z coordinate _Fi In addition, the ground surface specific height H _Ri 3D feature data Pl3D with added as attribute data _i Output as.
[0073]
Next, the 3D map output unit 8 sends all Pl3Ds. _i Are output as a 3D map with elevation.
[0074]
Next, the explanation for the evaluation of the reliability of the altitude estimation will be supplemented.
[0075]
In the following cases, the following example can be considered as a method of estimating the reliability of altitude estimation.
[0076]
・ Evaluation function value E _i When is big. In other words, let T2 be a certain threshold value.
E _i ≧ T2 Formula (19)
As a method of determining the value of T2, the evaluation function value E _i It is reasonable to set a value that is 1 to 2σ larger than the average value.
[0077]
・ Determined height H _Fi Is at both ends of the search range. That means
H _Fi = Hmin or H _Fi = Hmax Formula (20)
In this case, the search range may not be set correctly.
[0078]
In addition, when estimating the loss or change of a feature, it can be estimated that the possibility of the loss or change of the feature is high in the following cases.
[0079]
・ When the ground surface specific height is small. That is, let T1 be a threshold value close to 0.
｜ H _Ri | ≦ T1 Formula (21)
In addition, there is a case where a change in the height, size, etc. of the feature occurs as a cause of the low reliability.
[0080]
For example, as shown in FIG. 9, even when the altitude of the feature is Horg at the time of shooting in year A, the altitude changes from Horg to Hnew when shot after several years (aging). There is a case.
[0081]
In addition, the size (planar shape) of the building may also change with time (Fiscal year C).
[0082]
In the case of estimation of a change in a feature whose altitude is already known, a change in a feature whose altitude is already known is estimated as follows.
[0083]
The altitude of the feature is Horg, and the search range of the altitude is
Horg−dH ≦ H _k ≦ Horg + dH.
[0084]
When there is a secular change, the estimated altitude is significantly different from Horg, or the reliability of altitude estimation from the evaluation function value is low.
[0085]
It is also possible to construct old stereo images with different shooting times as an extended stereo model, perform altitude estimation, and determine that there is a secular change when the altitude estimation reliability decreases.
[0086]
Further, for example, as shown in FIG. 10, it is possible to deal with a case where the feature has a complicated shape or is not actually a certain height. Experiments have confirmed that the building shows a value close to the mode elevation of the building when compared with the results of actual measurements with aerial photogrammetry and laser scanners.
[0087]
When the feature height is close to 0, the feature can be automatically extracted as a candidate for a lost portion of the feature, and correction can be made more efficient.
[0088]
<Embodiment 4>
Next, map update will be described with reference to the flowchart of FIG. In this description, it is assumed that a three-dimensional map (old), a two-dimensional map (old), and a newly captured extended stereo model are stored in advance on the disc.
[0089]
In such a state, the two-dimensional map (old) is read and the extended stereo model is read. As described above, the height of the amorphous window of the known feature is given, and this is changed to change the extended stereo model. Find the projection position on both images of the model, find the brightness value of the area divided by the irregular window at this projection position, compare the brightness value, and compare the height of the irregular window when the comparison results are the best match Elevation processing is performed by extended stereo matching using the height as the height of the feature in the extended stereo model (S1), and the 3D map data is stored as an intermediate result (S2).
[0090]
Then, a change extraction process is performed by extended stereo matching to obtain a change from the 3D map (old), the 3D map data obtained in the middle of the process, and the new extended stereo model (S3). ) Is stored as a three-dimensional map of intermediate results before being assigned to ().
[0091]
Next, editing processing on the photograph is performed (S5). This is preferably done using a digital photogrammetry workstation.
[0092]
In this editing process, the 3D data of the intermediate result obtained in step S4 is overlaid on the new extended stereo model and displayed on the monitor, and the 3D data is compared with the newly captured extended stereo model. By doing so, the operator manually edits. Here, a newly taken general stereo image may be used as the extended stereo model, but a stereo model composed of an image obtained when an old 2D map was created and a newly taken image is used. Also good.
[0093]
And the process which updates old three-dimensional data to the three-dimensional data obtained by this edit process is performed (S6).
[0094]
That is, by this series of processing, the numerical map data having no elevation is easily converted into three-dimensional data, and the existing map can be updated.
[0095]
【The invention's effect】
As described above, according to the present invention, a two-dimensional map can be efficiently converted into a three-dimensional map. In particular, when it is not necessary to take a new aerial photograph, it is possible to inexpensively create a three-dimensional city model for use in computer graphics and three-dimensional city simulation.
[0096]
In addition, by converting the numerical map into a three-dimensional map, it is possible to superimpose a new photographed map on a stereo aerial photograph in a two-dimensional map aging correction work, and the correction becomes easy.
[Brief description of the drawings]
FIG. 1 is a schematic configuration diagram of an extended imaging method using an indefinite window according to the first embodiment.
FIG. 2 is a schematic configuration diagram of an extended imaging method using an indefinite window according to the second embodiment.
FIG. 3 is an explanatory diagram for explaining a change in the height direction of an indefinite window.
FIG. 4 is an explanatory diagram for explaining the specific height of the third embodiment.
FIG. 5 is a specific configuration diagram of each embodiment.
FIG. 6 is a configuration diagram of an elevation search unit.
FIG. 7 is an explanatory diagram illustrating geometric characteristics of an area sensor.
FIG. 8 is an explanatory diagram of meshing of feature data.
FIG. 9 is an explanatory diagram for explaining changes in feature data when shooting times are different;
FIG. 10 is an explanatory diagram for explaining how to obtain the altitude of a building having a complicated shape.
FIG. 11 is a flowchart illustrating map update according to the third embodiment.
[Explanation of symbols]
1 Digital image input section
2 Hard disk
3 Orientation element calculation part
4 Numerical map input part
5 Feature data extraction unit for processing
6 Image top feature measurement unit
7 Elevation search section

Claims (3)

In an imaging matching method using a combination of two types of images with different resolutions in the same region,
Performing an orientation operation on the two types of images and generating an extended stereo model related by function;
A polygon point sequence of a known feature that exists in common in both images of the extended stereo model is used as a standard amorphous window;
Providing each of the irregular windows with a height (z) and changing the height (z);
Each time the height of the irregular window is changed, obtaining a projection position on both images of the extended stereo model according to the height; and
Comparing luminance values of regions separated by the irregular window at the projection positions of both images of the extended stereo model;
Outputting the height (z) of the irregular window when the comparison result is most approximated as the height of the feature of both images of the extended stereo model;
And a step of determining a three-dimensional coordinate of the feature using the output height (z), and an extended image matching method between images using an indefinite window. 2. The method includes the step of subtracting the height from the geographical data to obtain a specific height of the feature every time the height of the feature is obtained, and having geographical data of the surface of the area. An inter-image extended imaging matching method using the described irregular window. The non-standard window according to claim 1 or 2, further comprising a step of notifying that the feature has disappeared when the image based on the non-standard window does not exist in a new image. The inter-image extended imaging matching method used.

Images