楕円モデルのハフ変換による数値地形図上のクレータ認識手法の開発／原田，平田，出村，浅田

69

楕円モデルのハフ変換による数値地形図上のクレータ認識手法の開発
原田 直人 ，平田 成 ，出村 裕英 ，浅田 智朗
1111（要旨）本論文では「かぐや」月数値地形図における大容量データ解析を見込んだ，クレータ認識手法とその準備状況について述べる．この手法で採用した個々の処理手法の原理は画像処理分野で既に広く使われているものではあるが，それらを組み合わせる事で楕円モデルのハフ変換に最適化したクレータ認識手法を確立した．クレータ検出結果を利用する観点で火星数値地形図を用いた予備的検証を行った結果は，この手法と従来研究の検出精度はほぼ同等であり，クレータサイズ頻度分布の形も保存しているため，手動カウンティングを代替できる見通しをつけることができた．

1. はじめにる問題がある[2]．このため，クレータ地形認識の自動化・省力化を図る必要性が説かれてきた．

2007年9月14日に月周回衛星「かぐや」（SELENE）

クレータ認識の自動化は，世界的にも多数研究さは打ち上げられ，本稿執筆時点では既に月周回軌道にれており，写真やDTMからの認識を目標としてい入り，観測運用に向けて準備が進められている．本稿

る．写真からクレータを認識する方法には，ハフ変換

では，「かぐや」に搭載されている地形カメラ（Terrain[4]やファジーハフ変換[6]，テンプレートマッチング

Camera，TC）が取得するデータの解析のために準備[7]，ラーニングアルゴリズム[8]といったものが代表とされている研究について，その現状を紹介する．

して挙げられる．DTMを使用した方法には，ハフ変TCはその名の通り，月の地形データを取得する．

換[3,5]がある．しかし，それらの手法では，かぐやか

本 研 究 の 目 的 は， 月 の 主 要 な 地 形 要 素 の 一 つ で あらもたらされる写真やDTMを解析するのに幾つかのる，クレータを自動的に認識するアルゴリズムの検討

問題点がある．例えば，テンプレートマッチングではと，これを実装したツールの開発である．TCの主な

あらかじめ全ての場合を尽くしたテンプレートを準備

データプロダクトには，水平分解能10ｍの数値地形図しておく必要がある．また，ラーニングアルゴリズム(Digital Terrain Model，DTM)があり，本研究ではこは教師付分類の一種と考えればこれも外部パラメータれを入力とする．クレータを楕円の凹地形として認識を必要とすることになる．このことからこれら2つのし，その諸元（中心位置，長軸径，短軸径，長軸の方位，手法は大容量の解析に向いていない．さらに，過去に深さなど）を出力する．

研究されてきたハフ変換を使用した手順[4]も前処理

クレータの自動認識の需要は大きい．例えば，クレ段階などで入力画像に依存したパラメータ調整が必要

ータを網羅的にカウンティングした結果得られるサイとされていたため，自動化されたデータ解析には至っ

ズ頻度分布はクレータ年代学の基礎データである．まていなかった．また，過去に研究された手法[3-5,7,8]

た，クレータ諸元のサイズ依存性も統計的に求められは結果の評価方法が惑星科学分野における実際の利用るので，大量の処理が必要である[1]．このような大量形態に即していなかったため，その有効性が十分検証の解析を手作業で行う場合，仮に何らかのソフトウェされていなかった．沢辺ら[6]のファジーハフ変換を採アの支援があったとしても，多大な時間を要する．さ用したアルゴリズムはパラメータの調整や利用形態にらに，手作業では，結果が作業者の熟練度に左右され即した検証の問題を解決した．しかし，沢辺らの手法

1. 会津大学コンピュータ理工学研究科情報システム学専攻

は写真を対象にしたアルゴリズムなので，かぐやのデ70

日本惑星科学会誌 Vol.17.No.1,2008

図１: 処理手順 (A)入力データ Holdenクレータおよびその周辺部のMOLA DTM (150×156ピクセル)．(B)Haarウェーブレッ
ト変換後．(C)二値化後．(D)メディアンフィルタ後．(E)ダイレーション処理後．(F) 細線化後．(G)クレータ認識結果．
(H)クレータの再フィッティング後の結果．

ータから作成されるDTMに適応できるか分からない．

非常に適している．また，我々は楕円モデルのハフ変

よって，かぐやのデータから作成されるDTMを解析換に適応した新しい前処理の手順も作成した．そして，する新しいアルゴリズムが必要となっている．

この手法の認識精度を手動認識による結果と比較し，我々は，DTMからクレータを認識する新しい手法の開発を行った．クレータ認識の部分には楕円モデルサイズ頻度分布へ与える影響などの観点からも評価した．のハフ変換を採用した．楕円モデルのハフ変換は，ハフ変換を楕円でも認識できるように拡張した手法であ2. テストデータ

り，クレータリムが欠損している場合でも高精度の認識が可能であるというハフ変換の利点も備えている．

テ ス ト デ ー タ と し てMOLA (Mars Orbiter Laser

よって，我々が解析しようとしているデータの処理にAltimeter) DTMを使用した．このデータは火星のデ楕円モデルのハフ変換による数値地形図上のクレータ認識手法の開発／原田，平田，出村，浅田

71

図２: 楕円モデルのハフ変換の手順. (A) エッジ抽出画像．(B) エッジ点の周り3ピクセルを取り出す．(C) 取り出した画像．
(D) 取り出した中のエッジから1点選び，中心のエッジ点と選んだ1点を結ぶ直線の切片と傾きを計算し，ハフ空間へ投票する．(E) 取り出した画像内の全エッジで全通り行う．(F) ハフ空間上で，最も投票数が多いものを，その取り出した画像の接線とする．(G) 画像全体に対してGを行う．(H) 接線の中から2本選ぶ ．(I) 楕円の性質 (2本の接線の交点とその2点の接点を結ぶ線の中点を結ぶと，円の中心を通る．(J) 楕円の性質を利用し，中心を通る線を書く．これを全ての接線の組み合わせで行う．中心を通る線が一定以上集まっている場所が存在すれば，それを楕円の中心とする (K) 楕円の中心とその中心に投票した全てのエッジで楕円の長軸径，短軸径，長軸の傾きを計算する．計算したパラメータを新しいハフ空間に投票する．一番投票数が集まったパラメータをそのクレータのパラメータとして出力する．

ータであるが，現在手に入るものとしては，最も解像クレータのリム抽出画像と一致させるため再フィッテ

度のよいDTMの一つであり，クレータ地形がよく表ィング処理を行う．

現されている．

3. 手法

3.1

クレータのリム抽出

本論文ではクレータのリム抽出処理に以下の5段階のステップをふんだ．(1) エッジ検出，(2) 二値化，(3)

このクレータ認識手法は大きく分けて3つの段階で構ノイズ除去，(4) エッジ連結，(5) 細線化．

成されている．(1)クレータリムの抽出．このステップ
では標高値が急激に変化するクレータのリムをエッジ

3.1.1

エッジ検出として抽出する処理を行う．(2)クレータ認識．ここでエッジ検出はHaarウェーブレット変換のハイパスは楕円モデルのハフ変換を使用し楕円のエッジをクレフィルタを採用した． Haarウェーブレット変換のハータとして認識する処理を行う．(3) 再フィッティング．

イパスフィルタは，細かい地形や複雑な地形でも検出

この最後のステップは，(2) で認識されたクレータをする事ができ，クレータを認識しやすいエッジ検出画

72

日本惑星科学会誌 Vol.17.No.1,2008

像を得ることが可能である(図1B) ．

つ目は，エッジの数を減らし投票処理にかかる時間を削減することである．細線化処理にはHilditch [8] の

3.1.2

二値化

Haarウェーブレット変換後の画像を二値化することでエッジを取り出した．二値化のための閾値を決定する方法には様々なものがあるが，既存の方法では必要なエッジだけを取り出す事はできなかった．このた細線化法を使用した (図1F)．

3.2

クレータ認識クレータ認識には楕円モデルのハフ変換を使用した認識手法を採用した[9,10]．

め，以下に述べる新しい方法を開発した．(1) 画像のヒストグラムを求める，(2) ヒストグラムの各ビンの3.2.1

楕円モデルのハフ変換

値に，前後6つのビンの値を加える，(3) ビンを大きいハフ変換とは画像上にあるエッジ点から直線を計算

方から検査し，値が全画素数の1%以下となる画素値し，その直線のパラメータをハフ空間上に投票するこを求め，そこを閾値とする，(4) 閾値以下の画素値をとで，その画像上にある直線を識別する方法である．

持つ画素をエッジとして取り出す ．この手法を用い楕円モデルのハフ変換とは，ハフ変換を楕円でも認識ることで，よい二値化結果を得ることができた (図1C)．

できるようにしたものである．

ただし，(2) における加算ビン数や，(3) における閾値の基準とする画素数は，MOLA- DTMに最適化され3.2.2

楕円モデルのハフ変換の手順ている．実際にTCによる月のDTMを解析する際には，楕円モデルのハフ変換を使用しクレータを認識する

撮像条件等に応じた最適値をそれぞれ求めてから処理手順は以下のようになる．(1) クレータのリムを抽出を行う予定である．

した画像の一部を取り出す(図2B,2C)．(2) (1) で取り出した中のエッジから2点選び，2点を結ぶ直線の切片と3.1.3

ノイズ除去

二値化された画像には細かなエッジが混ざっている．

傾きを計算し，ハフ空間へ投票する(図2D)．これを取り出した画像内の全エッジ点で全通り行う(図2E)．(3)

そこでメディアンフィルタをかけることで認識に必要ハフ空間上で，最も投票数が多いものを，その取り出でないエッジを除去した．フィルタサイズは3×3であした画像の接線とする(図2F)．(1)-(3) を繰り返し，画る．(図1D)

像全体に対して行う(図2G)．(4) 接線の中から2本選び，楕円の性質(2本の接線の交点とその2点の接点を結ぶ

3.1.4

エッジ連結線の中点を結ぶと，円の中心を通る(図2H))を利用し，エッジの膨張処理を行い，エッジを連結させた．こ

中心を通る線を書く．これを全ての接線の組み合わの処理の目的は近くに存在するエッジ同士を一つの塊せで行う(図2I)．(5) 中心を通る線が一定以上集まってとして認識させることで，エッジの欠損部分を補いクいる場所が存在すれば，それを楕円の中心とする(図レータ認識の精度を上げる事である (図1E)．この処理2J)． (6) 楕円の中心とその中心に投票した全てのエッで採用したフィルタサイズは5×5である．つまり，エジで楕円の長軸径，短軸径，長軸の傾きを計算する．

ッジ間に欠損があっても，4ピクセル以内の間隙であ(7) 計算したパラメータを新しいハフ空間に投票する．

れば埋めることが可能である．

(8) 一番投票数が集まったパラメータをそのクレータのパラメータとして出力する(図1G,2K)．

3.1.5

細線化

連結処理を行ったエッジを1ピクセルの線にするた

3.2.3

段階的な楕円認識めに細線化した．この処理の目的は２つある．一つ目

3.2.2項で示した楕円認識手順の (5) の段階での楕円は，線の太さに違いがあると楕円モデルのハフ変換の中心の投票による決定において，大きな楕円ほど多数投票処理にばらつきがでてしまい，公平な投票が行わの接線を持ち得るため，中心への投票数も多くなる．

れない．そこで細線化することで，投票はエッジ数だ

このため，楕円認識の得票数の閾値を一定とすると，けに依存し，線の太さに依存しない投票がされる．２

大きな楕円だけが選択的に認識されてしまうことにな楕円モデルのハフ変換による数値地形図上のクレータ認識手法の開発／原田，平田，出村，浅田

73

図３: 楕円モデルのハフ変換による楕円認識の欠損耐性(100×100ピクセル) ．(A)入力データ．  (B)クレータリムのエッジ抽出結果．(C) 欠損率20%のエッジデータ．(D) (C) のク
レータ認識結果 (長軸径78pixel，短軸径74pixel)．(E) 欠損率50%のエッジデータ．(F) (E) のクレータ認識結果 (長軸径72pixel，短軸径50pixel)．(G) 欠損率60%のエッジデー
タ．
(H) (G) のクレータ認識結果 (長軸径68pixel，短軸径48pixel)．

る．そのような状況に陥るのを防ぐため，あらかじめ認識対象とする楕円の径を一定範囲に制限して，小さな楕円から段階的に認識させることで，認識率を向上させ，誤認識と認識もれを減らす事が可能となる．
3.2.4
3.2.4.1欠損耐性
欠損耐性の検証方法

楕円モデルのハフ変換を使用したアルゴリズムはある程度の欠損があるエッジでも楕円を認識することができる．そこで，本手法による楕円認識の欠損耐性を調べるために，本来のエッジを一部消去したデータを入力とした試験をおこなった．クレータ一つ分のエッ

図４: アージャイルクレータ周辺部(1000×500ピクセル) に対する自動認識結果（再フィッティング前: A，再フィッティン
グ後: B）．(C) 同じ地域についての手動認識による結果．

ジを含む二値化済みのDTMを100×100ピクセルのサ検証の結果，エッジの欠損率が0%〜20%のときは，イズで切り出した．そして，エッジの一部を段階的にエッジ抽出画像とほぼ一致した楕円を認識できた．し消去し，どのように認識されるか試した．

かし欠損率が20%を超えると一致した楕円を当てはめ消去の仕方は，一部分が完全に欠損している場合をることが困難になる．そして欠損率が50%以上になる再現した．エッジが飛び飛びに欠損している場合は，とほとんど楕円を当てはめる事はできなかった．この前処理の段階で補間されたり，楕円モデルのハフ変換原因として，欠損部分が大きくなると楕円の長軸の推の投票ではエッジの連続性に関係なく投票されたりす

定ができなくなるため，欠損部分の補完が不可能になるので，このアルゴリズムでは影響は少ない．しかし，っている事が分かった．

エッジの一部が完全に損失した場合は，損失した部分を補えないので．楕円の識別が困難になる．よって今

3.3

クレータの再フィッティング

回は，最悪なケースを想定しエッジが完全に欠損する

このアルゴリズムの特性として短軸が短く認識され場合を再現する事とした．図3ではこのうち，欠損率る傾向がある(図1G)．その問題を修正するためにクレ20%，50%，60%の例を示す．

ータの再フィッティングを行った．(1)最初に決められた中心と長軸径，長軸の傾きを固定したままで短軸径3.2.4.2

検証結果を1pixel長くする．(2) エッジ抽出画像と短軸径を変更

74

日本惑星科学会誌 Vol.17.No.1,2008

した楕円との一致度を計算する．(3) (1)と(2)の作業を短軸径が長軸径と同じ長さになるまで繰り返す．(4)
最も一致度が高かった短軸径をその楕円の短軸径として採用する (図1H)．

4. クレータ自動認識の精度検証
4.1

比較対照データ

このアルゴリズムによるクレータの自動認識の精度を検証するため，同じDTMデータに対して目視による手作業のクレータカウンティングを行った．自動図５: 自動認識および手動認識による検出クレータ長軸径のヒストグラム

認識による結果を図4A（再フィッティング前）および図4B（再フィッティング後）に，また，手動による結果を図3Cにそれぞれ示した．以後はこれらを基に議論する．

4.2

検出限界

図5に，自動・手動別のサイズ（長軸径）ごとのクレータ検出数を示す．これより，このアルゴリズムによるクレータ自動認識の下限は直径10ピクセル程度であることがわかる．
下限を決める要因として考えられる要素として，第一に挙げられるのは，小さいクレータの場合，エッジからの投票数が圧倒的に少なくなり，偶然による投票とクレータ本来の投票との識別が難しくなることがある．また，ダイレーション処理などにより，小さい地形は他の地形と接合されてしまう可能性もある．

4.3

検出率

このアルゴリズムによって検出されたクレータの総数は102個，うち，長軸径10ピクセル以上の個数は101個であった．一方，手動による検出個数は309個のうち，前節で行った検出限界の長軸径10ピクセル以上

図６: 自動認識および手動認識の結果から作成したクレータサ
イズ頻度分布 (CSFD) プロット (A)，およびRプロット (B)．の個数は114個であった．従って長軸径が10ピクセル以上のクレータについて検出もれは13個．検出率は89

サイズ-累積頻度分布 (CSFD) プロットおよびRプロッ％ (検出もれ率11％ )であった．また，クレータ以外のトにまとめたものを図6に示す．いずれにプロット形地形を誤検出した例が4個あったため，誤検出率は4％式においても，自動認識の結果と手動認識の結果はグである．

ラフ形状が長軸径範囲全てにわたって，よい一致を示

4.4

クレータサイズ頻度分布による検証

自動認識の結果を，クレータ年代学で一般に用いるしている． 自動認識による結果が，手動による結果に対してやや下にプロットされており，クレータ年代が若く評価される可能性がある．しかし，年代が既知

楕円モデルのハフ変換による数値地形図上のクレータ認識手法の開発／原田，平田，出村，浅田

75の領域で検出もれを確認することで，年代値の較正が算で380Gpixelとなるので，これを今回と同じ計算機できるだろう．この結果から，本アルゴリズムによる一台で解析した場合，処理時間は約27年と見積もるこ

クレータ自動認識の出力がクレータ年代学に適応可能とができる．計算機の処理能力の向上や，複数機で同であるといえる．

時に処理を進めることを考えれば，現実的な期間内に4.5

評価

本アルゴリズムによるクレータの検出限界である直解析を終えることが可能である．

5. まとめ径10pixelという値に関しては，過去の例でDTMを対象にした手法[3,5]と同等の結果である．また，検出率

本研究では，DTMを対象とした，楕円モデルのハ89％も他の手法と比較してほぼ同等かそれ以上の結果

フ変換を用いたクレータ地形の自動認識アルゴリズムを示している．ただし，本研究を含めて，各手法におを開発した．MOLA DTMによる予備的検証では，不ける検出率の算出にあたっては，対象とするクレータ完全なエッジからでもクレータを認識することが可能のサイズ範囲，基準とするクレータリストなどが統一であることが確認されたほか，手動認識結果を正解とされていない状態であり，検出率の値の単純な比較はした場合の検出率89%を示した．この手法を実装した

難しいより正確な評価をするためにはMOLAのクレツールでは，クレータの位置の他，長軸，短軸径や長ータデータベース[11]と比較する必要があり，今後か軸の向きなどのパラメータ取得も行うことができる．

ぐやデータと平行して検証する予定である．サイズ-

本手法は，対象地形およびその大きさ毎のパラメータ

累積頻度分布プロットの形状は，手動による検出結果決定は不要であるが，撮像条件などに依存した画質なと，本アルゴリズムによる検出結果とでほとんど一致

らびに影領域に強く影響されるパラメータは求めておしており，十分クレータ年代学に使用できることが分かないといけない．具体的には，二値化アルゴリズム

かった．の箇所が該当しているが，それ以外の先見的なパラメ

本手法の前処理部分は別として，クレータ認識過程ータの決定は不要であり，「かぐや」大容量データのではテンプレートマッチングのような認識する地形に解析を一括して行うのに有利である．また，このツー依存した外部入力を必要としない．よって，テンプレ

ルはかぐやプロジェクトメンバ内での利用に限らず，ートマッチングよりは大容量データの解析に向いてい惑星科学コミュニティに向けて公開することも検討中るといえる．

である．

この手法で認識することができなかった直径10pixel以下のクレータについて，写真からの認識を謝辞行った沢辺ら[6]は，直径が2pixel以上のクレータについて認識可能であるという結果を示している．また，本稿を丁寧に査読され，有益なコメントを頂けた高写真を対象とした他の手法[7,8]でも5pixelの認識限界知大本田理恵氏に感謝申し上げます．惑星科学分野とを示している．よって，直径10pixel以下のクレータ

画像処理技術分野の境界領域成果として，綺麗にまと認識結果が必要な場合は，写真からの認識と組み合わめることができました．本手法の開発は，会津大学コせることが現実的な対応である．

ンピュータ理工学部マルチメディアシステム学講座に4.6

解析時間よる継続的な研究によって進められたもので，卒業生である松本直也氏と狩野智皓氏が試行錯誤した，楕円今 回 使 用 し た 計 算 機 は，PowerMac G5（CPU:モデルのハフ変換とエッジ検出の一部も使われておりPowerPC G5 2.5 GHz×2, RAM: DDR 8GB）である．

ます．あらためて，この二人に対してもここに感謝の入力DTMのサイズが300×250 pixelの時，平均的な意を記します．

解析時間は約170秒であった．従って，解析速度は440pixel/secとなる．かぐやによって得られる解像度10m/pixelの月全球地形モデルの全ピクセル数は，概76

日本惑星科学会誌 Vol.17.No.1,2008

引用文献[1] Pike, R. J., 1974, in Impact and Explosion
Cratering, 489.
[2] Sawabe, Y et al. 2004, 遊星人, 13, 87.
[3] Bue, B. et al., 2007, Lunar and Planetary Science
XXXVIII, Abstract #1178.
[4] Bruzzone, L. et al., 2004, Proceedings of ESAEUSC 2004, 13.
[5] Michael, G. G., 2003, Planetary and Space Science,51, 563.
[6] Sawabe, Y. et al., 2006, Advances in Space
Research, 37, 21.
[7] Bandeira, L. P. C. et al., 2006, First International
Conference on Impact Cratering in the Solar
System. 17.
[8] Honda, R. et al., 2002, Progress of Discovery
Science, Springer Verlag, 395.
[9] Hilditch C. J. et al., 1969, Machine Intelligence 4,Eginburgh Univ. Press, p.403[10] Matsumoto, N. et al., 2005, Lunar and Planetary
Science XXXVI，Abstract #1995．
[11] Salamuniccar, G. and Loncaric, S., G., 2006, 36th
COSPAR, Abstract# #450