LEVEL SET METHODの高速化および医療画像への応用に関する研究

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

LEVEL SET METHODの高速化および医療画像への応用 に関する研究

由井, 俊太郎

https://doi.org/10.15017/1441253

出版情報：Kyushu University, 2013, 博士（工学）, 課程博士 バージョン：

権利関係：Fulltext available.

LEVEL SET METHOD

の高速化および 医療画像への応用に関する研究

2013

年

由井 俊太郎

目 次

第

章 序論

1.1

本研究の背景

1.1.1

セグメンテーション手法の応用分野

1.1.2

距離画像におけるセグメンテーション手法

1.1.3

医用画像におけるセグメンテーション手法

1.1.4

課題

1.2

本研究の目的

1.3

本論文の構成

第

章

と動的輪郭モデルの概要

2.2 Level Set Method . . . . 11

2.2.1

概要

2.2.2 Level Set Method

の実装

2.2.3

再初期化

2.2.4

成長速度場

2.3 Level Set Method

の高速化手法

2.3.1 Narrow Band Method . . . . 19

2.3.2 Fast Marching Method . . . . 20

2.3.3 Hermes Algorithm . . . . 22

2.4 Level Set Method

を用いた動的輪郭モデル

2.4.1 Geometric Active Contours (GAC) . . . . 24

2.4.2 Geodesic Active Contours . . . . 26

2.4.3 Chan-Vese

法（C-V 法）

第

章

の高速化手法

の提案

3.1.1

距離変換を用いた

の生成

3.1.2

距離マップを用いた最近傍点の探索

第

章 距離画像に対する動的輪郭モデルの高速化

距離画像に対する

を用いた高速な曲面再構成法の提案

4.1.1 FNB

を用いた高速な曲面再構成法

4.1.2

解像度制御型

4.2

実験

4.2.1 FNB

による高速化

4.2.2

解像度制御型

による高速化

第

章 サイズ適応型原発性肝がん検出法

従来手法と課題

5.1.1 CT

画像を用いた腫瘍の検出と抽出の従来手法と問題

5.1.2 Level Set Method

を用いた一般的な医療画像処理

5.1.3 Level Set Method

を用いた腫瘍検出手法とその問題

5.2

サイズ適応型原発性肝がん検出法の提案

5.2.1

白色雑音起因の精度劣化の改善

5.2.2

低コントラスト起因の精度劣化の改善

5.2.3

血管起因の精度劣化の改善

5.2.4 Level Set Method

の改良

5.3

実験

5.3.1

白色雑音起因の精度劣化改善手法の定性評価

iii

5.3.2

低コントラスト起因の精度劣化改善手法の定性評価

5.3.3

血管起因の精度劣化改善手法の定性評価

5.3.4

提案手法の定量評価

5.4

考察

第

章 結論

謝 辞

参考文献

付 録

付録

幾何モデリングの関連研究

データ計測

A.2

距離画像の位置合わせ

付録

条件

付録

モフォロジーを用いた曲面形状進化法

モフォロジー演算

C.2

モフォロジーを用いた曲面形状進化の理論

C.2.1

モフォロジーを用いた曲面形状進化

付録

形状強調フィルタ

iv

第

^{章 序論}

カメラ，CTなどを用いて得られる画像情報や，レーザレンジファインダ，ステレオカメラ などを用いて得られる距離画像情報から，ユーザが注目する領域を自動的に抽出するセグメ ンテーション手法の一つに，Snakesに代表される動的輪郭モデルがある．従来の動的輪郭モ デルは，対象領域の分離，結合など位相変化への対応が困難であることが問題であった．これ に対し，位相変化に対応可能な動的輪郭モデルとして， Osher，SethianらによりLevel Set

Methodが提案され，医療画像処理やカメラ画像内の移動体検出・追跡などで利用されてい

る．しかし，Level Set Methodには多くの計算量が必要であり，実応用には計算量の削減に よる高速化が課題であった．さらに，医療画像において腫瘍領域の自動検出・抽出にLevel Set

Methodを適用した場合，輪郭の変形を制御する従来の定式化では対象領域の検出・抽出精度

が低いため，その実応用には検出・抽出精度の向上が課題であった．

本論文は，Level Set Methodに基づくセグメンテーション手法に対し，計算の高速化，お よび医療画像における高精度化を目指した研究を纏めたものである．本章ではまず，研究の 背景を示した後，研究の目的と課題について述べ，最後に論文の構成を示す．

1.1

本研究の背景

1.1.1 セグメンテーション手法の応用分野

セグメンテーション手法は，コンピュータビジョン・医療分野・バーチャルリアリティな ど多岐にわたる分野に適用されている．

• ^{コンピュータビジョン}(Computer Vision)

特徴毎に画像を領域分割するセグメンテーションは，意味理解や形状計測に必要不可 欠な基礎技術である．例えば，コンピュータビジョンの典型的なアプリケーションとし て，カメラ画像内での特定物体の探索や，その位置姿勢を推定するモデルベース物体認 識技術（Model Based Object Recognition）がある．この手法は，予め計算機内部に3 次元物体モデルを構築するモデル構築ステップと，モデル照合ステップの2ステップか ら構成される．このうちモデル照合ステップでは，まず，カメラ画像，あるいはレーザ レンジファインダなどから得られた対象物体の距離画像を，領域毎にセグメンテーショ

1.1.

本研究の背景

ンして対象の表面形状を抽出する．次に得られた表面形状から幾何学的，あるいは見え の特徴を抽出して，モデルとの照合により物体を認識する．

• ^医用画像(Medical Image Processing)

疾病の早期発見・治療を目指し，コンピュータ断層法(computed tomography:CT)^や X線を用いて得られた画像処理技術が盛んである．例えば，胸部X線画像から癌領域 をセグメンテーションすることで，診断支援を実現している[1]．この診断法では，過 去に撮影された同一被検者の胸部画像を比較することで，病変の検出や病状の変化を 得る．

また，CTなどで得た一枚の画像は2次元画像であるが，複数枚のCT^{画像を用いて}3 次元画像をつくり出すことができる．この3次元画像は人体内部の構造に関する多くの 情報を含んでおり，3次元画像を使った臓器形状のセグメンテーションも試みられてい る．この技術を発展させて，仮想臓器を切開することで，手術をシミュレーションする 研究もなされている[2]．

• ^{バーチャルリアリティ}(Virtual Reality)

バーチャルリアリティは，人間と計算機のインターフェースとして非常に重要である．

その中でも，全てを仮想世界のみで構成するのではなく，仮想物体を現実背景画像に埋 め込む複合現実感は，仮想世界のリアリティ向上に特に効果的な手法である．複合現実 感システムは，仮想モール・ゲーム・ロボットプランニング・文化財のデジタルコンテ ンツ化といった幅広い応用が考えられる．しかし，バーチャルリアリティで用いられる 様々なモデルの作成作業は，その大部分をオペレータの人手に頼っている．このためコ ンテンツの制作に多大な時間と費用が必要となり，これがバーチャルリアリティの普及 を妨げる一因となっている．従って，ビジョンセンサなどにより取得したシーンに対し てセグメンテーション手法などを適用することで，モデル作成作業を自動化・効率化で きれば，バーチャルリアリティシステムの低価格化が可能となり，より一層の普及が期 待される．

1.1.2 距離画像におけるセグメンテーション手法

距離画像を用いた，実物体の３次元曲面を再構成する一般的な手順を以下に示す．

(1) 距離情報の計測

まず，レーザレンジファインダ等の計測装置を用いて対象物体を計測し，ある視点から 観測した対象物体までの奥行き情報を画素値とする距離画像を取得する．

1.1.

本研究の背景

(2) 複数距離画像の位置合わせ

1枚の距離画像には，対象物体表面のうち，ある視点から観測可能な領域の3次元デー タしか含まれない．そのため，複数視点から対象物体を計測し複数枚の距離画像を取得 することが必要である．ここで取得した距離画像はセンサの視線方向に依存した座標 系で表されているため，画像間の位置合わせを行い各座標系を統一した基準座標系に 変換する．

(3) 曲面再構成（データ統合，曲面形状の復元）

(2)により位置合わせされた距離画像は３次元点群のデータであり，ノイズなど対象物 以外の点群データも含まれる．この点群データからノイズを除去し，さらに点間の隣 接関係を構築することで，対象物体の表面形状を抽出する．つまり，曲面の再構成をす る．曲面再構成により得られる表面の法線情報は，物体認識・見えの解析・仮想現実・

複合現実など様々なアプリケーションで重要となる．

以下では(3)に着目し，(1)(2)は付録に掲載した．

曲面再構成に関して，これまでに多くの手法が提案されている．Turk[3]らは，2つの距離 画像の重複部を取り除き，境界をつなぎ合わせることで距離画像を統合した．しかし，複数 の距離画像を単純に縫合するだけであるため，例えば，複数の計測装置を用いた場合などは，

物体の部位によってメッシュの密度が極端に異なるモデルを獲得するという問題がある．ま た，得られるデータには必ずノイズが含まれるため，距離画像の位置合わせの際に誤差が発 生し，不自然なつなぎ目が生じる問題もある．この他に，距離画像から陰表面関数(implicit

surface function)を構成し，そのゼロ等位面を抽出することで曲面再構成する手法も提案さ

れている[4][5][6]．Hoppeら[5]の提案手法では，まず点の集合から接平面を推定し，ボクセ ルからその接平面までの符号付距離(signed distance)を計算して陰表面関数を構築する．次 に，その陰表面関数をmarching cubes algorithm[7]に適用してモデルを生成した．しかし，

ノイズに非常に弱いという問題がある．Wheeler[6]らはこのノイズの問題に対して合致表面 法(consensus surface algorithm)を提案した．この手法は，局所的に法線の連続性を考慮す るため，ノイズに頑強である．しかしこれらの陰表面関数を用いた手法は，薄い物体を両面 から計測した場合，測定誤差のため正しい符号付距離関数が得られず，滑らかで閉じた曲面 が必ずしも得られないという問題がある．

これに対し，Snakes[8]に代表される動的輪郭モデルが提案されており，3次元に拡張した 手法（変形可能曲面(deformable surfaces)[9]）も提案されている[10][11][12]．これらの手法 ではまず，対象物体の内部または対象物体を囲むように閉曲面を設置する．次に，曲面の関 数として定義されたエネルギー関数を減少するように曲面の位置・形状が繰り返し更新され，

1.1.

本研究の背景

エネルギー関数を最小(^極小)とする曲面を対象領域の境界として求められる．

また，Chenらは動的バルーンモデルを用いた手法を提案した[13]．具体的な手法を以下に 示す．まず，多視点の距離画像を統合して得られた全周型距離画像を入力データとし，その 中心に初期状態として小さな三角形メッシュモデルをおく．このメッシュの頂点に膨張力を 働かせる．これによりあたかも風船が膨張するかのように各頂点が放射状に移動する．つま り，距離画像上のデータからの外力を持たせずに，バルーンモデルを構成する三角形メッシュ の頂点にのみ膨張力と張力を与える．これにより三角形メッシュを物体表面に張り付ける(図 1.1)^．

Fig.1.1:

バルーンモデルの概略

これまでに示した動的輪郭モデルは，滑らかで閉じた曲面が必ず得られるという特徴を 有するため，広く利用されてきた．しかし，図1.2に示すように曲面の分裂や統合など位相

(topology)が変化する場合の対応が困難である．同様に，複数ある物体においても，モデリ

ングの自動化は困難である．これらの問題に対応して正しいモデルを生成するためには，対 象物体の位相に関する知識を事前に準備する必要がある．

1.1.

本研究の背景

How to connect?

initial model

inflation

boundary data

(a) (b)

Fig.1.2:

バルーンモデルでモデリングできない例

1.1.

本研究の背景

この問題を解決するため，動的輪郭モデルの概念を拡張し，曲面形状進化（波面伝播(front propagation)とも呼ばれる）を用いた曲面再構成の研究が盛んに行われている[14][15][16]．曲 面形状進化を実現する数値計算法の一つとして，Level Set Methodが提案されている[17][18]．

Level Set Methodは，予め設定した運動方程式に従う曲面の形状進化を追跡する手法であり，

曲面を一つ次元の高い補助関数のゼロ等高面とみなし，運動方程式に従って補助関数を更新 することで，ゼロ等高面に対応する曲面の形状を進化する．Level Set Methodの最大の利点 は，曲面形状進化の過程で分裂や統合といった位相の変化が可能であり，柔軟な曲面形状の 進化が実現できる点である．そのため，Level Set Method を位相適応型のセグメンテーショ ンに適用する研究が盛んになっている．さらにLevel Set Methodは曲面形状進化の過程で位 相変化が可能なことからセグメンテーションにとどまらず，移動物体追跡[19]・結晶形成シ ミュレーション[20]・煙，水，炎といった自然現象のシミュレーション[21]など様々な分野で 利用されている．

Level Set Methodを用いた曲面再構成の手法は，geometric active contoursを用いたもの と，geodesic active contours[22][23]を3次元に拡張したminimal surfaces[24]を用いたもの とに分類できる．geometric active contours^はSnakes^{に代わる手法として，}Caselles^ら[25]

やMalladiら[14]が提案したものである．これらの手法ではまず，対象物体の内部または対象

物体を囲むように閉曲面を設置し，それを初期曲面として対象物体の表面近傍まで膨張，収縮 させて面モデルを生成する．ここでは，Snakesのように頂点間距離などから定義されるエネ ルギー関数が最小になるような曲面形状進化を実現するのではなく，曲面の存在する場所が 対象物体の輪郭かどうかを判別しながら対象物体の境界に向かって曲面形状が進化する．こ の曲面形状進化をLevel Set Methodにより実現しているため，曲面形状の位相変化が可能で ある．つまり，たとえ対象物体が複数個あっても，事前知識を用いずに同時に曲面再構成が 可能になる．一方minimal surfacesでは，Snakesと同様，曲面の関数として定義されたエネ ルギー関数が減少するように曲面の位置・形状が繰り返し更新され，エネルギー関数を最小 (極小)とする曲面が対象物体の輪郭として求められる．さらに，Level Set Methodを用いて いるため，位相変化が可能な曲面形状進化を実現している．しかし，エネルギー関数の最小 化に計算コストを要するというという問題がある．本論文では，曲面再構成の高速化と精度 向上の両立を目指し，計算時間を要する曲面形状の位相変化はgeometric active contoursを，

曲面形状の補正はSnakes^を3^{次元に拡張した}Deformable Surfaceをベースに曲面再構成を 実現する．

1.1.

本研究の背景

1.1.3 医用画像におけるセグメンテーション手法

医療診断の分野で使われる画像処理技術の一つに，CTや磁気共鳴イメージング（Magnetic Resonance Imaging: MRI）があり，その有効性はよく知られている．しかしCTやMRIな どの高度化により，医用画像は２次元画像から３次元画像に広がり始め，医療の現場からは

『膨大なデータの読影に今まで以上に読影時間がかかるにも関わらず，装置の進歩により１日 の検査件数は飛躍的に増加し全体の読影時間が著しく増加している．また，１検査当たりの データ増加により，一人の医師が短時間の診療時間内に見落としなく読影することがかなり 困難になっている．このため，放射線科医の負担が著しく増加している．この問題を解決す るには短時間で正確な省力化読影法の開発が必須であり，コンピュータ支援画像診断開発な どの緊急の課題となっている[26]．』との医師の悲痛な叫び声が年々増している．実際，残念 ながら医師の画像診断（読影）において，病変の見落としや判断ミスはゼロではない．例え ば，検診マンモグラフィ（乳房X線写真）における病変の見落としは約20^％[27]^{，胸部単純} X線写真を用いた肺ノジュール検出ではその割合は約30％という報告がある[28]．Kundelら は，このような見落としの原因の多くが人間のミスと指摘しており，コンピュータ支援診断 に対する期待が高まっている[29]．

病変部の検出を支援するため，動的輪郭モデルを医用画像に適用し，病変の候補位置を自 動検出することができれば，病変の見落としや判断ミスを減少することが可能になる．特に

Level Set Methodを用いた動的輪郭モデルは，位相変化が可能であるため，CTなどの医用画

像に適用し，複数臓器や腫瘍を同時に検出・抽出する応用研究が盛んである．例えばMalladi らは，CT画像から肝臓などの臓器形状を抽出する手法を提案した[14]．

1.1.4 課題

Level Set Methodを用いたセグメンテーション手法は優れた手法であるが，計算量が膨大

であるため，莫大な時間を要するという問題がある．これに対し，これまでにNarrow Band Method (NB) [30]やFast Marching Method[31]などのLevel Set Methodの高速化手法が提 案されている．しかしNBでは，ゼロ等高面近傍の空間を格子状に分割したボクセルにおい て最近傍曲面を探索しなければならず，実応用の観点から計算量の削減が必要である．一方，

Fast Marching Methodの計算量はNBに比べ改善されているが，曲面の膨張と収縮が同時に

実現できず，曲面形状が一方向しか進化できないという制限がある．そこで，曲面形状進化 の制限がない，柔軟な曲面形状進化が可能なLevel Set Methodの高速化が必要である．

さらに，医療画像への適用を考えた場合，医用画像は人間の複雑な臓器構造を反映してい

るため，Level Set Methodを用いても病変の候補位置を自動検出する事が難しいという問題

1.2.

本研究の目的

がある．特に，３大疾患の一つである悪性腫瘍の早期診断は重要であり，肺がんや乳がんな ど一部の悪性腫瘍に対する自動検出手法は確立されているものの，５年生存率が低い肝がん の自動検出手法は確立されていない．従来，CT画像を用いた肝がん検出では，1cm以上の 肝がんを自動検出する研究が盛んであり，これらの研究では造影剤投入後の撮影タイミング が異なる2種類の画像（動脈相，遅延相）が用いられてきた．このうち，遅延相は正確な肝 がんのサイズ計測が可能であるが，コントラストが低いため1cm以下の肝がん（小結節）の 検出が難しい．一方，動脈相では小結節の検出が可能になるが，雑音と肝がんの輝度値が類 似しているため誤検出が生じやすい，という問題があった．そこで，CT^{画像を用いた，様々} なサイズの肝がん（以下，腫瘍）に対応可能な腫瘍検出法の確立が必要である．

1.2

本研究の目的

本論文の目的は，Level Set Methodに基づくセグメンテーション手法に対し，計算の高速 手法，および医療画像におけるセグメンテーションの高精度手法を開発することである．本 論文では，特に以下の２つの課題に対して検討を行う．

課題1 柔軟な曲面形状進化が可能なLevel Set Methodの高速化

課題2 Level Set Methodを用いた肝臓がんの自動検出手法の確立

課題１を解決するため，NBをさらに高速化したFast Narrow Band Method(FNB)を提案 する．本手法では，最近傍曲面を探索する際，ゼロ等高面からの距離情報を備えた距離マッ

プ(distance map)を利用し，探索が必要な曲面点を大幅に限定する．これにより，効率的な

最近傍曲面点の探索処理を実現し，曲面形状進化の高速化を実現する．また，FNBを用いた 曲面形状進化は，処理が高速であるにもかかわらず，従来手法であるNarrow Band (NB)に 対して原理上，精度の低下がない．さらに，距離画像を用いた動的輪郭モデルの高速化を目 的に，解像度を変化させながら曲面形状を進化させる手法も提案する．これは，始めに物体 の概形を抽出し，徐々に細部の輪郭を抽出することで，高速化を実現する．

課題２を解決するため，CT画像を用いた，様々なサイズの腫瘍が検出可能な腫瘍の自動検 出法を提案する．低コントラスト画像や雑音に対しても頑健な腫瘍検出を実現するため，輪 郭の変形を制御する成長速度の新たな設定法を考案した．これにより，Level Set Methodの 特徴である腫瘍個数の事前知識が不要になることに加えて，様々な形状，サイズの腫瘍が検 出可能になる．

1.3.

本論文の構成

1.3

本論文の構成

本論文は，第1章「序論」，第2章「Level Set Methodと動的輪郭モデルの概要」，第3章

「Level Set Methodの高速化手法:FNB」，第4章「距離画像に対する動的輪郭モデルの高速 化」，第5章「サイズ適応型原発性肝がん検出法」，及び第6章「結論」から構成されている．

以下に各章の概要を述べる．

第2章では，Level Set MethodとLevel Set Methodを用いた動的輪郭モデルの関連研究 について述べる．第3^章では，Level Set Method^{の高速化手法として}FNB^{を提案する．第} 4章では，距離画像に対する動的輪郭モデルにFNBを適用し，解像度を変化させながら曲面 形状進化をすることで，より高速化を実現する手法を提案する．さらに，提案手法の有効性 を実験的に示す．第5章では，Level Set Methodを医療画像に適用するため，サイズ適応型 の原発性肝がん検出手法を提案し，その有効性を実験的に示す．第6章では，本論文で得ら れた結果を統括し，今後の課題を述べる．

第

^章

と動的輪郭モデルの概要

本章では，まず動的輪郭モデルの代表的な手法であるSnakesの概要と問題点を述べる．次に，

Level Set Method[17][18]の概要を述べる．続いて，これまでに提案されたLevel Set Method の高速化手法を紹介する．最後に，従来の動的輪郭モデルの概要と問題点を述べる．

2.1 Snakes

Kass^{らが提案した}Snakes[8]^を3次元に拡張した変形可能曲面(deformable surfaces)[9]^は，

種々の対象に適用でき，かつ，計測ノイズを最小限に抑えることができるため，広く利用さ れてきた．これは，Snakesやdeformable surfacesにおける対象物体の輪郭抽出では，ある評 価関数(エネルギー)が最小になるように曲面形状を修正することで，曲面を滑らかにすると 同時に曲面が対象物体の輪郭に引きつけられるためである．本章では，2次元の濃淡画像にお ける輪郭抽出の方法と，本論文の対象の一つである距離画像データから曲面を再構成する手 法について述べる．

パラメータ化された閉曲線を[C : [0,1]→R², p7→C(p)]，対象画像をI :Z⁺×Z⁺→R⁺ とし，対象画像から物体形状の輪郭を抽出する問題を考える．Snakes^{では，次のエネルギー} 関数の値を最小化することにより，対象物体の輪郭を抽出する．

E[(C)(p)] =α

∫ ₁

0

∂C

∂p

²dp+β

∫ ₁

0

∂²C

∂p^{2 2}dp

| {z }

正則項

+ (−γ

∫ ₁

0 |∇I(C(p))|dp)

| {z }

引力項

(2.1)

始めの２項は内部エネルギー(E_int)と呼ばれる項で，この値は曲線Cの長さや形状の複雑 さが減少するにつれ小さくなる．一方，最後の項は外部エネルギー(E_ext)と呼ばれる項で，

この値はエッジ強度が増すにつれて小さくなる．つまり，内部エネルギー(Eint)の値が小さ くなるにつれ曲線は滑らかになり，外部エネルギー(Eext)の値が小さくなるにつれて曲線が 対象物体の輪郭に近づくことを意味する．

このSnakesは3次元に拡張することができ，deformable surfaces[9]と呼ばれ, 3次元物体 の輪郭抽出などに使われてきた．ここでパラメータ化された閉曲面を

v(r,s) = (x(r.s),y(r.s),z(r.s)) (r,s)∈[0,1]×[0,1]

2.2. LEVEL SET METHOD 11

と定義する．このとき，最小化するエネルギー関数は次のように与えられる．

E[v] =

∫

Ω

[ ω₁₀∂v

∂s

²+ω₀₁∂v

∂r

²+ω₂₀∂²v

∂s^{2 2} +2ω₁₁∂²v

∂s∂r

²+ω₀₂∂²v

∂r²

²+P(v(s,r)) ]

drds (2.2)

最後のP(v(s,r))はポテンシャル関数を表し，外部エネルギー(Eext) に対応し，用途に応 じて決定される．濃淡画像の輪郭抽出の場合，P(v(s,r)) =−||∇I(C(p))||^{2 となる．このエ} ネルギー関数を最小化する曲面v(s,r)^は

− ∂

∂s (

ω₁₀∂v

∂s )− ∂

∂r (

ω₀₁∂v

∂r )

+ ∂²

∂s² (

ω₂₀∂²v

∂s² ) + ∂²

∂s∂r (

2ω₁₁ ∂²v

∂s∂r )

+ ∂²

∂r² (

ω₀₂∂²v

∂r² )

= −∇P(v(s,r)) (2.3) のオイラーラグランジェ方程式(Euler-Lagrange equations)を満たすことが知られている．

距離画像からの曲面再構成にdeformable surfacesを用いた手法が提案されている[10][11][12][32]． これらの手法でのポテンシャル関数は

P(v(s,r)) =







−1 d(v(s,r)) =0

− 1

d(v(s,r)) otherwise (2.4)

のように定義される．ただし，d(v(s,r))は距離画像データと曲面v(s,r)との距離である．[10]

では，有限要素法を用いて式(2.3)を解き，曲面再構成を実現した．詳細は[33]などが詳しい．

このSnakesモデルでは，初期曲面の近傍にある適切な輪郭が得られ，雑音に強い滑らかな

曲線が得られる．しかし，曲面形状進化の過程で分裂や統合といった位相の変化に対応でき ないため，個数の事前知識がない状態では複数個の物体の輪郭を同時に自動抽出することが できない．そのため，対象物体と同じ位相を持った初期曲面を与える必要があった．

2.2 Level Set Method

2.2.1 概要

まず，C(x,0)を滑らかな初期閉曲面とし，この初期閉曲面がt秒後に外向き法線方向nの 成長速度F(図2.1)でC(x, t)になる場合を考える．ここで，nは外向き単位法線ベクトルで あり，成長速度Fは．通常，曲率κなど曲面形状に依存した関数である．以上述べた事を式 で表現すると，







∂C

∂t(x) = F(κ)n C(x,0) = C₀(x)

(2.5)

2.2. LEVEL SET METHOD 12

A

A‘

normal vector at A

F A move on A‘

Initial Curve

Deformable Curve

Fig.2.1:

成長速度

となる．上記の式に従って曲面形状進化を実現するためには，差分方程式を利用したラグラ ンジェ法を用いればよい．しかし，曲面形状進化の過程で位相の変化に対応できない，安定 解が得られないという問題があった．

この問題を解決するため，OsherとSethianによりLevel Set Methodが提案された[17]． Level Set Methodでは，曲面C(x, t)を一つ次元の高い補助関数ψのゼロ等高面であると考 え，補助関数ψ(x, t)を次のように定義する．

ψ(x, t) =±d (2.6)

dは，xから曲面C(x, t)までの最短距離である．また，xが曲面C(x, t)の外側にある時はプ ラスの符号を選び，逆に曲面の内側にある時はマイナスの符号を選ぶ．この曲面C(x, t)^は，

定義により

C(x, t) ={x|ψ(x, t) = 0} (2.7) を満足する．このため，補助関数ψ(x, t)を追跡することができれば，ゼロ等高面(ψ= 0)で ある曲面C(x, t)の運動を追跡することができる．

ここで，x^が曲面C(x, t)上の点であると仮定すると，式(2.7)^より，

ψ(x(t), t) = 0 (2.8)

となる．これを時間で偏微分すると，連鎖法則より，

ψ_t+

∑N i=1

ψ_xix_it = 0 (2.9)

2.2. LEVEL SET METHOD 13

となる．ただし，x= (x₁, x₂,· · ·, x_n), ψ_xi = ∂ψ

∂x_i, x_it = ∂xi

∂t ^{である．式}(2.9)を内積表現に すると，

ψt+∇ψ(x, t)·xt= 0 (2.10) となる．ここで，曲面C(x, t)の外向き法線の単位ベクトルnは ∇ψ(x, t)

|∇ψ(x, t)|^{であり，成長速度} はF =x_t·nであることに注意すると，曲面形状進化を表す次式のハミルトンヤコビ型偏微 分方程式が得られる．

ψt+F|∇ψ|= 0 (2.11)

この偏微分方程式は，時間依存型等高面方程式(Time-dependent Level Set Equation)と 呼ばれている[17]．Level Set Methodでは，このように初期曲面C₀(x)と成長速度を与えて，

式(2.11)に基づいて曲面の形状進化を実現する．

Level Set Methodには従来の曲面形状進化法と比べ様々な利点がある．まず，Level Set Methodでは，曲面C(x, t)を補助関数のゼロ等高面(ψ= 0)として扱っているので，曲面形 状進化の過程で分裂や統合といった位相の変化が生じても，曲面の運動の追跡が容易である．

つまり，曲面C(x, t)を直接扱うのではなく，補助関数ψ^{を更新し，曲面を}ψ= 0^{とみなすこ}

とで(式(2.7))，位相変化に対応した曲面形状進化が実現できる．例えば，図2.2は２次元平

面上の曲線の形状進化の過程で位相が変化する様子を示している．まず，式(2.6)に従い，初 期曲線の補助関数を生成する(図2.2(b)→(a))．次に，式(2.11)に従い微小時間後の補助関数 を得る(図2.2(a)→(c))．微小時間後の曲線形状は，式(2.7)に従い補助関数ψのゼロ等高面 (ψ= 0)を抽出することによって得られる(図2.2(c)→(d))．

また，曲面の幾何学的特徴がψから比較的簡単に求めることができるという利点がある．

代表的な幾何学的特徴である法線ベクトルはn=∇ψ(x, t)^{で表される．また曲率}κ^は，単位 法線ベクトルを考慮することで次式のようになる．

κ = ∇ · ∇ψ

|∇ψ|

= −ψxxψ_y²−2ψxψyψxy+ψyyψ²_x

(ψ_x²+ψ²_y)^3/2 (2.12)

2.2.2 Level Set Methodの実装

Level Set Method^{の実装法について，}3次元空間での適用を例に述べる．実装のためには，

連続空間で定義された偏微分方程式(2.11)から，離散空間で実装される差分方程式への変換方 法を検討する必要がある．この偏微分方程式から差分方程式を導く手法は[18]にていくつか紹 介されている．ここでは，Upwind Schemeと呼ばれる手法について述べる．Upwind Scheme

2.2. LEVEL SET METHOD 14

ƒˆ(t) ƒ (t)

ƒˆ(t+1) ƒ (t+1)

(a) (b)

(c) (d)

Fig.2.2:

等高面法の概略

2.2. LEVEL SET METHOD 15

とは，曲面形状進化を曲面情報の伝播と考え，その情報が伝わる方向を考慮して導いた手法 である[18]．

Level Set Methodの実装では，まず，対象空間をI×J×K個のボクセルに分割する．次 に，成長速度Fが常に凸関数であるとき，式(2.11)は次式のような差分方程式に近似できる．

ψ_i,j,kⁿ⁺¹ =ψⁿ_i,j,k−∆t(max(F_ijk,0)∇⁺+ min(F_ijk,0)∇⁻) (2.13) なお，h^{を空間グリッド幅，}∆tを時間ステップ幅とすると，ψ_i,j,kⁿ =ψ(ih, jh, kh, n∆t)^であ る．ここで，補助関数の一次近似微分を用いる場合，∇^+や∇⁻は以下のように表せる．

∇⁺ = (max(D⁻_ijk^x,−D_ijk^+x,0)² + max(D⁻_ijk^y,−D^+y_ijk,0)²

+ max(D⁻_ijk^z,−D^+z_ijk,0)²)¹² (2.14)

∇⁻ = (max(D^+x_ijk,−D_ijk^−x,0)² + max(D^+y_ijk,−D⁻_ijk^y,0)²

+ max(D^+z_ijk,−D⁻_ijk^z,0)²)¹² (2.15) ただし，

D^+x_ijk= ψ_i+1,j,kⁿ −ψ_i,j,kⁿ

h D^+y_ijk= ψⁿ_i,j+1,k−ψ_i,j,kⁿ

h D^+z_ijk= ψ_i,j,k+1ⁿ −ψ_i,j,kⁿ h D⁻_ijk^x= ψ_i,j,kⁿ −ψⁿ_i−1,j,k

h D⁻_ijk^y = ψⁿ_i,j,k−ψ_i,jⁿ_−1,k

h D⁻_ijk^z = ψ_i,j,kⁿ −ψⁿ_i,j,k−1 h

(2.16)

である．ここでは，補助関数の一次近似微分を用いているが，より滑らかな曲面を実現する ためには，補助関数の二次近似微分を用いればよい．

2.2.3 再初期化

補助関数ψを繰り返し更新する際，補助関数ψが微分不可能になる場合があり，差分演算 (式(2.13))をそのまま実行したのでは補助関数が距離関数の性質(式(2.6))を失う恐れがあ る．また，等高面方程式を差分方程式で近似しているため，補助関数の更新とともに誤差が 積算されることからも，補助関数が距離関数の性質を失う恐れがある．そこで，ある決まっ た回数だけ補助関数を更新するたびに，式(2.6)に従って補助関数ψを再構築する再初期化 (reinitialization)が通常行われる[14]．しかし，再初期化による補助関数の値を得る際，現在 の曲面点からの距離を求めるには，そのボクセルからの最近傍曲面点の探索処理を全ボクセ ルで行わなければならない．そのため，この再初期化処理の計算量が膨大になるという問題

2.2. LEVEL SET METHOD 16

がある．ここで，3次元空間における再初期化処理の計算量は，各座標方向に格子点がn^個と すると，曲面点の数がn²，全ボクセル数はn³にそれぞれ比例するので，O(n²×n³) =O(n⁵) となる．そこで，距離関数の性質を利用した高速な再初期化処理手法が提案されている．具 体的には，距離関数は|∇ψ|= 1というEikonal方程式を満たすことが知られており，Eikonal 方程式の高速な数値解法として知られているFast Marching Methodを用いることで，再初 期化の高速化を実現した[18]．しかし，あるボクセルの距離関数を決定するのに，logn² に比 例する計算時間がかかるため，再初期化処理の計算量は，O(n³×logn²) =O(n³logn²)とな り，依然として計算コストは大きい．そのため，補助関数の安定解を得るため，CFL^条件(^付 録.B)を満たした小さな時間ステップ幅∆tで補助関数の更新を行い距離関数の性質をできる だけ保つことで，再初期化処理の回数を少なくする必要があった．

2.2.4 成長速度場

Level Set Method^は式(2.5)を実現する手法であるため，成長速度F^{は全空間において定} 義されているものではなく，曲面点でのみ定義されており，{(i, j, k)|ψ_i,j,kⁿ 6= 0} ^{を満たすボ} クセルでは成長速度が存在しない．そこで全てのボクセルの成長速度を決定するために，局 所成長速度場（No extension velocity field）と拡張成長速度場（Extension velocity field）を 構築する2つの手法が提案されている．局所成長速度場とは，曲面点でないボクセルも曲面 点であると仮定して成長速度を決定する．つまり，各ボクセルにおける補助関数の成長速度 を決定する際に，当該ボクセル自身の状態，例えば対象物体輪郭かどうかや補助関数の曲率 などの情報を用いるものである．一方，拡張成長速度場とは，曲面点上のボクセルにおいて 成長速度を決定し，その他のボクセルでは最も近い曲面点のボクセルにおける成長速度をコ ピーして用いるものである．具体的には，全空間における成長速度を以下のように定義する [14]．まず，曲面点Q(ψ= 0)における成長速度を定義する．次に，ゼロ等高面にないボクセ ルP(ψ6= 0)の成長速度FをボクセルPに最も近い曲面点Qにおける成長速度と同じ値とす る(図2.3)．

2次元における両者の曲面形状進化の様子を図2.4に示す．局所成長速度（図2.4（a））は，

曲面点が対象物体輪郭に到達しその部分の成長速度が減少しても，その情報が周辺に伝わら ず，周辺の速度には影響がない．そのため，形状進化を続けるうちに補助関数場が境界外側 では密に，内側では粗になり，さらに境界周辺以外では関数場は変化し続け，停止することは ない．一方，拡張成長速度（図2.4（b））は，曲面点が対象物体境界に到達すると，その情報 が周囲に伝わって周囲の進行速度も減少させるため，形状進化を続けても局所成長速度の場 合のように補助関数場が不均一になることがない．Sethian[18]も多くの計算機実験から，拡 張成長速度場のほうが補助関数の解が安定に求まり，最終的に得られる曲面形状も高精度で

2.3. LEVEL SET METHOD

の高速化手法

ψ=0

P(x,y) ψ=C

Q(x,y)

Q is closest to P

Fig.2.3:

成長速度の拡張

あることを示している．しかし拡張成長速度場を構築するには，再初期化処理と同様に最近 傍曲面点の探索処理を全ボクセルで行わなければならないため，計算コストが莫大になると いう問題があった．

これらのことから，補助関数の更新の際に解を安定的に得るには，

(1) CFL条件を満たした小さな時間ステップ幅∆tで補助関数の更新を行う．

(2) ^{再初期化処理を行う．}

(3) 拡張成長速度場を用いる．

が必要である．しかし，これまで提案されてきたLevel Set Methodのアプリケーションの多 くは，計算コストの問題から，1.(CFL条件)のみ行い，再初期化処理は時々行われ，局所成 長速度場が用いられてきた．

2.3 Level Set Method

の高速化手法

これまでに，多くのLevel Set Methodの高速化手法が提案されており，代表的なものと して，Narrow Band Method (NB) [30]とFast Marching Method[31]が挙げられる．NBで は，全てのボクセルにおいて成長速度場を構築したり補助関数ψの更新をするのではなく，曲

面C(x, t)の近傍にあるボクセルにおいてのみ成長速度場の構築や補助関数ψの更新を行う．

2.3. LEVEL SET METHOD

の高速化手法

Zero level set

Zero level set

(a) No extension velocity field

(b) Extension velocity field

Fig.2.4:

局所成長速度と拡張成長速度での曲線形状進化の比較

2.3. LEVEL SET METHOD

の高速化手法

Fast Marching Method^は，Eikonal^方程式(|∇T|F = 1)の数値解析法として提案されたもの である．この方程式を解くには，収束計算による繰り返し処理のために膨大な計算時間が必 要になるが，Fast Marching Methodでは，成長速度の符号を固定するという条件を加えるこ とで，繰り返し処理が不要になり計算量が削減された．しかし，成長速度の符号に制約があ り，曲面の膨張と収縮を同時に行うといった自由な曲面進化ができない．これらの問題を解 決するため，様々な手法が提案されている．Hermes Algorithm[19]では，成長速度Fの絶対 値が最も大きい曲面点を選択し，その近傍でのみ曲面形状進化を行う．曲面形状進化の過程

では，式(2.11)の正しい解は得られないが，最終的には正しい形状が得られ，自由な曲面変

化が可能である．また，この計算量はFast Marching Methodとほぼ同じであることが実験 的に確かめられている．本項では，それぞれの手法とその問題点について述べる．

2.3.1 Narrow Band Method

Narrow Band Method(以下，NB)はLevel Set Methodの代表的な高速化手法であり，最 も広く利用されている手法の一つである[30]．NBでは，全てのボクセルにおいて成長速度や 補助関数ψを決定するのではなく，曲面C(x, t)の近傍にあるボクセルにおいてのみ成長速度 や補助関数ψを構築し補助関数の更新をする．具体的には，{x| |ψ(x, t)|< δ} ^{を満たすボク}

セルをNarrow Bandと定義し，この領域内でのみ成長速度Fや補助関数ψを更新する．こ

こで，δはバンド幅と呼ばれている．2次元平面上のゼロ等高線周辺のNarrow Bandを図2.5 に示す．

ψ=0 δ

ψ=-δ

ψ=δ

Fig.2.5: Narrow Band

NBはLevel Set Methodに比べると高速な手法であるが，拡張成長速度を用いたり再初期