2014 年度修士論文

(1)

2014 年度修士論文

高階エネルギー最小化による 3D _医用画像の多臓器同時セグメンテーション

指導教員石川博教授

早稲田大学大学院基幹理工学研究科情報理工学専攻 5113B013-8 岡川明日翔

2015 年 2 月 6 日提出

(2)

図目次

2.1 人体と軸の関係 . . . . 6

2.2 人体胸部付近のCT 画像 . . . . 6

2.3 胸部から腹部における臓器の 3D モデル . . . . 7

2.4 複数臓器の同時領域分割 . . . . 9

2.5 fusion-move . . . . 14

2.6 1階と高階の違い(2階の例) . . . . 15

3.1 肝臓が上手くセグメンテーションされていない例 . . . . 17

3.2 相互情報量 . . . . 18

3.3 直角三角錐状に並んだクリーク . . . . 19

3.4 離れたボクセルを用いたクリーク . . . . 20

3.5 高階平滑化項のクリークとラベル値の組み合わせの例 . . . . 21

3.6 離れた高階平滑化項のクリークとラベル値の組み合わせの例 . . . . 22

3.7 x 方向に並ぶクリークの例 . . . . 24

3.8 ラベルのパターンと頻度 . . . . 25

3.9 ラベルのパターンと頻度(v₀ のラベルのの条件付き) . . . . 26

3.10 エッジとその上下に登録されたクリーク . . . . 27

4.1 Ground Truth画像 . . . . 29

4.2 隣接した L字型の 3 ボクセル . . . . 30

4.3 L 字型の3 ボクセル間の相互情報量の計算結果 . . . . 31

4.4 2ボクセル間の相互情報量の計算結果 . . . . 32

4.5 エッジ抽出の例 . . . . 35

4.6 直角三角錐状に並んだクリークを用いた高階ポテンシャルを加えたセグメンテーション結果 . . . . 36

4.7 離れたボクセルを用いたクリークのポテンシャルを加えたセグメンテーション結果 . . . . 37

(5)

4.8 距離と提案エネルギーの重みを様々に変えてセグメンテーションした結果 38 4.9 基準となるボクセルのz 座標を固定した場合とそうでない場合の比較 . 39

4.10 E_smooth を E_HSL に変えた時とそうでない時の比較 . . . . 40

4.11 Esmooth を EHSI に変えた時とそうでない時の比較 . . . . 40

4.12 E_smooth を E_HSC に変えた時とそうでない時の比較 . . . . 41

4.13 E_smooth をボクセル間が離れたクリークの E_HSC に変えた時とそうでない時の比較 . . . . 41

4.14 肝臓を含むと考えられる領域 . . . . 42

4.15 セグメンテーション結果(L^′ =Lliver) . . . . 44

4.16 異なるスライス面のセグメンテーション結果 (L^′ =L_liver) . . . . 44

4.17 異なる重みパラメータの時のセグメンテーション結果の例(L^′ =L_liver) . 45 4.18 セグメンテーション結果(L^′ =L_liver，クリークは3 方向に登録) . . . . 46

4.19 異なるスライス面のセグメンテーション結果(L^′ =L_liver，クリークは3方向に登録) . . . . 46

4.20 セグメンテーション結果(L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}) . . . . . 48

4.21 セグメンテーション結果(L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}) . . . . . 48

4.22 セグメンテーション結果(L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}，クリークは 3方向に登録) . . . . 49

4.23 セグメンテーション結果(L^′ ={Lliver, Lstomach, Lspleen, Lpancreas}，クリークは 3方向に登録) . . . . 49

4.24 E_HA を加えた時のセグメンテーション結果 . . . . 50

4.25 E_HAF を加えた時のセグメンテーション結果 . . . . 51

4.26 E_HAC を加えた時のセグメンテーション結果 . . . . 51

4.27 EEA を加えた時のセグメンテーション結果 . . . . 52

(6)

表目次

2.1 臓器名とラベル値の対応 . . . . 8 4.1 全症例の Jaccard index の平均(L^′ =L_liver) . . . . 43 4.2 半分の症例のJaccard index の平均 (L^′ =L_liver) . . . . 43 4.3 全症例の Jaccard index の平均(L^′ =L_liver，クリークは3 方向に登録) . 45 4.4 半分の症例の Jaccard indexの平均 (L^′ =L_liver，クリークは 3方向に登録) 45 4.5 全症例の Jaccard index の平均(L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}) . 47 4.6 半分の症例のJaccard index の平均 (L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}) 47 4.7 全症例のJaccard indexの平均 (L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}，ク

リークは 3方向に登録) . . . . 48 4.8 半分の症例のJaccard indexの平均(L^′ ={L_liver, L_stomach, L_spleen, L_pancreas}，

クリークは 3方向に登録) . . . . 48

(7)

第 1 _{章はじめに}

1.1 背景

セグメンテーションとは，特徴や意味が類似している部分領域に画像を分割する処理である．医用画像の場合，人体の臓器や病変などの異常部位といった特定の領域を画像から抽出するために用いられる[1]．特定の領域の輪郭を囲ってセグメンテーション結果として示す場合と，画像内の全画素あるいはボクセルのうち，特定の領域に属する画素あるいはボクセルのみを分類することによってセグメンテーション結果とする場合に大別される．

CT ， MR ， PET ， CAD などの技術の急速な発展に伴ない，それらの技術を使用した医用画像診断装置や診断支援システムから得られるデータの高精度化が近年進んでいる．医用画像診断装置を用いて医療現場で得られた高精度なデータに，コンピュータによる医用画像処理を施すことにより，医師の診断支援を目的とした研究が数多く行われている．医用画像処理では，セグメンテーションなどが主に使われる[1]．

医用画像の高精度化に伴ない，仮想現実感・拡張現実感などの技術と組み合わせた手術・治療支援システムの応用が拡大している．医療現場では，手術中・治療中(術中と呼ぶ)に，術中画像を用いて，実空間の患者と手術前・治療前(術前と呼ぶ)の 3D臓器モデルの位置合わせを行うことで，術前と術中の情報を統一し，統一空間内において手術器具・治療装置の位置を提示することで手術支援が行われる[1]．従って，手術支援を行ったり，手術・治療計画を立てるためには，臓器の 3Dモデルや臓器の位置をあらかじめ取得しておく必要がある．臓器の 3D モデルを得るため，術前に撮影された 3D 医用画像にセグメンテーションが施される．仮想現実感・拡張現実感などの技術と組み合わせた手術・治療支援システムにおいて，正確なシミュレーションや手術支援を行うには，高精度なセグメンテーションを行う必要がある．

セグメンテーションの対象となる医用画像には， 2D 医用画像と 3D 医用画像，および 4D 医用画像が存在する．体の断層画像を 2D 医用画像，断層画像を積み重ねたものを 3D 医用画像，さらに時間軸を追加して考慮した画像を 4D 医用画像と呼ぶ[1]．

セグメンテーションの対象とされる画像の次元は研究によって異なり，[2]では 2D 医

(8)

用画像，[3][4][5][6][7][8]では 3D 医用画像，[9][10]では 4D 医用画像を扱っている．3D 医用画像は， 2D 医用画像と比べて z 軸方向の情報も用いることができるため，セグメンテーションに用いる情報の選択の幅が広がる．4D 医用画像では時間軸が追加され，[9]

では造影剤投与からの時間による輝度値の違いをセグメンテーションの情報として用いており，[10]では心臓の時間による形状の違いをセグメンテーションの情報として用いている．

セグメンテーションに用いられる医用画像は研究によって異なり， CT 画像を扱う研究[3][6][7][9][11]や MRI 画像を扱う研究[12][13][14]がその他の画像に比べて多く存在する．一般的に， CT 画像は MRI 画像に比べて空間分解能が優れており， MRI 画像は CT 画像に比べてコントラスト分解能が優れている．たとえば，急性期脳梗塞のような病変は，コントラストの高い MRI画像で撮影するほうがふさわしく，骨や肺などの内部構造の細かい部分も撮影したい場合は CT 画像がふさわしいため，撮影したい対象によって使い分けられている．

また，頭部領域の画像を対象とする研究[8][12][13][14]と腹部・胸部領域の画像を対象とする研究[3][6][7][9][11][15][16]に大別することもできる．上述した CT 画像と MRI 画像の違いのため，頭部領域を対象とした研究では CT 画像より MRI 画像のほうが多く，

逆に胸部・腹部領域を対象とした研究では MRI 画像より CT 画像のほうが多い．

セグメンテーションの対象物は研究ごとに異なり多岐にわたるが，癌・腫瘍などの病巣を対象とする研究[12][14]，脳の領域の一部を対象とする研究[8]，臓器を対象とする研究[3][4][6][7][9][11][15][17]が挙げられる．[12][14]では脳腫瘍がセグメンテーションの対象であり，[8]は海馬や脳回，[3]は肝臓，[11]は腎臓，[15]は冠状動脈と心臓周辺の血管，

[4][6][7][9][17]は複数の腹部臓器がセグメンテーションの対象である．

臓器のセグメンテーションの場合，単一臓器のみのセグメンテーション[3][11][15]または複数臓器のセグメンテーション[4][6][7][9][17]に分けられる．臓器のセグメンテーションでは，人体内で臓器が存在する位置は異なる症例間でも大体同じであるという事前知識を利用して，臓器ごとの確率アトラスというものが作られる．確率アトラスとは，ある位置に各臓器が存在する確率の統計を取ることにより得られるモデルであり，確率アトラスを用いることにより，セグメンテーションの精度は上がる．一般的に，臓器の輝度値情報やアトラスを用いてセグメンテーションは行われる．各研究ではそれらの情報に加えて，新たな情報を用いることにより，精度向上を目指している．

単一臓器のみのセグメンテーションの場合，対象物を一つの臓器に絞れるので，臓器ごとに特化した手法を考えることができる．[3]では肝臓の形状の制約情報を加えてセグメンテーションを行っている．[11]では腎臓の形状推定とセグメンテーションを繰り返し

(9)

行うことによって精度を高めている．[15]では心臓の表面に対して平行に血管が存在する事前知識を利用してセグメンテーションを行っている．

複数臓器のセグメンテーションの場合，研究ごとに対象とする臓器数と臓器の種類は異なっている．臓器によっては症例間でアトラスに大きな差が生じる．[4][7]では，あらかじめ画像領域を分割してアトラスを階層的に作成し，症例間による臓器の形状や位置の違いを可能な限り吸収することで，セグメンテーション精度の低下を防いでいる．[6]

では，複数のアトラスのうち，セグメンテーションを施す画像に似ているアトラスを選択することにより，セグメンテーション精度の低下を防いでいる．[9]では，時間軸を考慮することにより，複数の臓器全体のセグメンテーション精度の向上を目指している．[17]

では，一部の臓器が失われている場合でも，セグメンテーションができる手法を提案している．

既存研究で行われている様々なセグメンテーションの手法は，より高精度なセグメンテーション結果を得ることや，セグメンテーション実行時間の高速化が主な目的である．

セグメンテーションの手法は，次の 5つに分類される[1]．領域拡張法やウォータシェッド法に代表される二値化手法，空間フィルタなどを利用したエッジ抽出に基づく手法，統計的分類器・クラスタリングを用いた手法， Snakes 法やレベルセット法に代表される可変形状モデルを用いた手法，グラフカットなどのその他の手法である．領域拡張法やウォータシェッド法に代表される二値化手法は，他の手法に比べてアルゴリズムが単純であるという利点があるが，領域のCT 値のみに基づく手法のため，CT 値が類似する隣接領域の誤った分割や，対象領域の過分割という問題点がある．空間フィルタなどを利用したエッジ抽出に基づく手法は，画像中のノイズに対して頑強な手法を構築することが容易ではない．統計的分類器・クラスタリングを用いた手法は，画素単位で判別が行われ，アトラスと組み合わせて使うことで医学的に不自然な位置や形状の領域が抽出されるのを防ぐことが出来る．可変形状モデルを用いた手法は，現在，医用画像処理の分野で最もよく用いられているセグメンテーション手法の一つである[1]．統計的形状特徴に関する事前知識を用いて，セグメンテーション対象の輪郭に力学的特性を仮定してエネルギーを定義し，最小化することでセグメンテーションを行う．可変形状モデルを用いた手法では，エネルギー最小化の際に，局所解に陥る問題に苦しめられることに比べ，その他の手法であるグラフカットは，劣モジュラ性を満たすという条件のもとでエネルギーの最小化が保証されているという利点が挙げられる．劣モジュラ性を満たさない場合でも，QPBO [18]のようなアルゴリズムを用いることによって，近似的にエネルギーを最小化することができる．また，二値ラベルだけでなく，多値ラベルを扱うための α 拡張[19]や融合移動[20]などの近似アルゴリズムが提案されている．

(10)

エネルギーの定義に用いる画素数またはボクセル数により，エネルギーの階数が決まる．従来のグラフカットでは，高階のエネルギーを最小化することができなかったため，

これまでエネルギーは 1階のものしか使われてこなかったが，最近になって高階のエネルギーを最小化できるアルゴリズムが開発された[21]．高階になると，エネルギーに用いる画素数またはボクセル数が増えるため，隣接する画素またはボクセル間の関係性のみならず，複数の画素またはボクセルを繋げることによって離れた画素またはボクセルの複雑な関係性を再現できるようになる．その一方，扱える画素数またはボクセル数が増えることで，高階のエネルギーの設計は複雑になり難しくなる．[21]では，高階のエネルギーを用いることにより， 1 階のエネルギーを用いた時と比べ，ノイズ除去の結果が良くなったことを示している．また，[22]では，形状に関する高階エネルギーを定義することにより，血管のセグメンテーションを行っている．同様に，胸部から腹部の領域における多臓器の同時セグメンテーションの場合でも，高階エネルギーの導入により，結果が改善されることを期待できる．

1.2 目的

本研究では，胸部から腹部の領域が撮影された 3D の CT 画像を入力としてセグメンテーションを行う．複数臓器の同時セグメンテーションをグラフカットで解く際に使用するエネルギーに高階エネルギーを導入する．高階エネルギーを導入することで，導入前と比較して，より高精度なセグメンテーションになると期待される．これまで胸部から腹部領域での医用画像セグメンテーションに高階エネルギーは用いられたことがないため，医用画像セグメンテーションに有効な高階エネルギーを設計する一般的な方法は存在しない．本研究では，様々な高階エネルギーを提案し，実験によりその有用性を評価する．

また，高階エネルギーを設計する際に，良いエネルギーを見つけるための尺度として，

相互情報量の導入が有効であるかの検討を行う．

(11)

第 2 章医用画像セグメンテーション

2.1 3D CT 画像

3D CT画像内の方向は解剖学と同じく，被撮影者の頭のある方が上，足のある方が下，

被撮影者から見た左右が左右，顔の向いている方が前，背部の向いている方が後ろと表現される．図2.1のように，右から左に向かってが x軸，前から後ろに向かってがy 軸，

下から上に向かってがz 軸となる．各断面は，xy平面が水平面，yz 平面が矢状面，zx 平面が冠状面と呼ばれる．本研究で対象とする医用画像は，水平面の断層 CT 画像であり，その断面の各臓器や組織の吸収するX 線の比率を基にグレースケールで表されている．この断層画像を z 軸方向に沿って一定の間隔で撮影することによって3 次元とみなしている．CT 画像は，図2.2のように，グレースケールで表示される．

CT 画像はボクセルごとに実数値をもっており，その数値が X 線の吸収率を表している．X 線の吸収率の単位は HU（ハウンスフィールドユニット）であり，空気は −1000 HU ，水は0 HUとなるように正規化したものである．例えば，各臓器は 20 HU から70 HU ，骨は数百 HU となる[1]．本論文では， CT 画像を図2.2のように，この値が高ければ高いほど白く表示する．

2.2 セグメンテーションの対象

本研究では胸部から腹部の範囲に含まれるいくつかの主要な臓器をセグメンテーションの対象としているので，胸部から腹部領域全体が撮影された CT 画像を用いる．胸部から腹部領域に含まれるのは，図2.3で示すような臓器である．この範囲に含まれる主要な臓器は，右肺，左肺，心臓，肝臓，胆嚢，胃，脾臓，右腎臓，左腎臓，膵臓などである．

この中でも，胆嚢や脾臓，膵臓は，症例間で臓器の位置の差が大きいため，セグメンテーションが難しい．

症例ごとに撮影時の体の位置や画像の解像度，被撮影者の年齢，性別，体格が異なるので，セグメンテーションを行うにあたり，正規化が必要となる．本論文のこれ以降の話では， CT 画像の正規化がされていることを前提とする．

(12)

図 2.1: 人体と軸の関係．体の左右方向の軸をx軸，体の前後方向の軸を y 軸，体の上下方向の軸を z 軸としている．

図 2.2: 人体胸部付近のCT 画像．人体胸部付近をグレースケールで表示している．白はより高い CT 値である．

(13)

図 2.3: 胸部から腹部における臓器の 3D モデル． 3D Slicer を用いてラベルをボリュームレンダリングして作成している．セグメンテーション対象の主要な臓器が含まれている．各臓器と色の対応は，右肺 (Blue)，左肺 (Lime) ，心臓(Red) ，大動脈 (Yellow)，食道 (Cyan) ，肝臓 (Coral) ，胃 (Green) ，脾臓 (Orange) ，右腎臓 (Purple) ，左腎臓 (Sienna) ，下大静脈 (Chocolate) ，膵臓 (Brown) である．

2.3 セグメンテーションの定義

CT の撮影対象である立方体の領域を各方向ごとに格子状に分割したボクセルの集合を考えV とすると，CT 画像は各ボクセル v ∈ V に CT 値を与える関数X: V →R で表される．CT 装置によって得られる CT 画像は，通常，撮影される人間の位置やボクセルの大きさが異なっているが，ここでは適当な方法で位置合せやリサンプリングなどの処理が行われているとする．本研究では，ラベルが割り当てられていない臓器や組織，

人体外部の空気や衣類，装置の全てを，背景ラベルとして定義する．背景及びM 個の異なる臓器を表すラベル値の集合を L :={0,1, . . . , M} とし，0 を背景ラベルとする．臓器名とラベル値の対応は表2.1に示した通りである．図2.4のように，画像 X に対して，

各ボクセルに臓器ラベルを与えるラベリングL: V → Lを，X のセグメンテーションといい，このとき Lv =L(v) と表す．

(14)

表 2.1: 臓器名とラベル値の対応

ラベル値臓器名本論文内での表記

0 背景 L_background

1 右肺 L_{right lung}

2 左肺 L_{left lung}

3 心臓 L_heart

4 大動脈 L_aorta

5 食道 L_esophagus

6 食道内腔 Lesophageal lumen

8 肝臓 L_liver

9 胆嚢 Lgallbladder

10 胃+十二指腸 L_stomach 11 胃+十二指腸内腔 Lstomach lumen

12 胃+十二指腸内容物 Lstomach contents

13 脾臓 L_spleen

14 右腎臓 Lright kidney

15 左腎臓 Lleft kidney

16 下大静脈 Linferior vena cava

17 門脈+脾静脈+上腸間膜静脈 Lportal vein

18 膵臓 Lpancreas

19 膀胱 L_bladder

20 前立腺 L_prostate

21 子宮 L_womb

(15)

⼊⼒

3D CT 画像 ܺ

出⼒

ラベリング画像 ܮ

図 2.4: 複数臓器の同時領域分割．画像 X に対して，各ボクセルに臓器ラベルを与えている．

2.4 エネルギー最小化による定式化

2.4.1 MRF

各ボクセルの臓器ラベルを決める問題は確率的に解くことができる．あるラベリング Lの各ボクセル v に対するラベルの確率変数を L_v として，そのラベリングを取りうる確率 P(L) = P(L₁,· · · , L_n) が得られるモデルに Markov random ﬁeld (MRF)がある[23]．

ここで，nは総ボクセル数を示している．MRF では，各確率変数L_v をノードとし，ノード間の依存関係を表すエッジを与えたグラフを構成することができる．

データとして，確率変数 X が与えられることを考える．確率変数X が与えられた場合のラベリングLを取りうる確率P(L|X)が得られるモデルとして，Conditional Random Field (CRF)がある[23]．MRF と同様に， CRF もグラフを構成できる．各ノード i にラベル L_i を与えられる確率がノード間で互いに独立していると考えるとき，

P(L|X) =∏

i∈V

P(L_i |X_i) (2.1)

と表される．また，ノード iにラベルL_i が与えられる確率とノードi に隣接するノード j にラベル L_j が与えられる確率に何らかの依存があるとき，

P(L|X) = ∏

i∈V

P(L_i |X_i)· ∏

{i,j}∈C2

P(L_i, L_j |X_i, X_j) (2.2)

(16)

と表される．ここで，C2 は，隣接するボクセルの組の集合を表す．依存関係のあるノードが3 つ以上の場合も同様に，条件付き確率の積で表すことができる．

2.4.2 MAP-MRF

画像 X に対してラベリング L が得られる確率 P(L| X) を事後確率という．X のセグメンテーションとして，事後確率を最大にするラベリング

LMAP:= argmax

L

P(L|X) (2.3)

を与えることを X の最大事後確率推定 (MAP 推定) という．エネルギー関数

E(L;X) :=−lnP(L|X) (2.4)

を使うと，

argmax

L

P(L|X) = argmin

L

E(L;X) (2.5)

となるので，エネルギー E の最小化問題として解くことができる．

2.4.3 ポテンシャル関数

一般的に，エネルギー関数は複数のポテンシャル関数から成り，エネルギー関数は E(L) = ∑

C∈C

f_C(L_C) (2.6)

と表される．ここで， C はクリークと呼ばれるボクセル V の部分集合の族であり， f_C はクリークC 内のボクセルに割り振ったラベルL_C に依存するポテンシャル関数である．

ポテンシャル関数は，クリークに用いるボクセルの最大の数により，階数が定義される．

階数はクリークに用いるボクセルの最大の数より 1 小さい数である．最小化されるエネルギー関数は，主に次のように表すことができる[24]．第一の和はデータ項，第二の和は平滑化項と呼ばれる．

E(L) =∑

v∈V

g_v(Lv) + ∑

(u,v)∈C2

huv(Lu, Lv). (2.7)

(17)

例えば，式 (2.7) において，データ項は 0 階，平滑化項は 1 階となる．エネルギー関数の設計の仕方により，セグメンテーションの精度は左右されるので，どのようなエネルギー関数を設計するかが重要となってくる．

医用画像セグメンテーションにおけるポテンシャル関数としては，アトラス項，データ項，平滑化項が一般的によく用いられる[1]．

アトラス項

各ボクセルに各臓器が存在する確率 P(Lv)を用いて定義されるポテンシャル関数である．アトラス項は，

E_atlas(L) =∑

v∈V

f_v(L_v) = −∑

v∈V

lnP(L_v) (2.8)

で表される．

データ項

CT 値ごとの各臓器となりうる確率 P(X_v|L_v) を用いて定義されるポテンシャル関数である．データ項は，

E_data(L;X) = ∑

v∈V

g(L_v;X_v) = −∑

v∈V

lnP(X_v|L_v) (2.9)

で表される．

平滑化項

隣接するボクセル同士をなるべく同じ臓器とするポテンシャル関数である．隣接するボクセルの CT 値の差が大きければ，その 2 ボクセルは異なる臓器である可能性が高いので，与えるエネルギーを小さくすることにより，異なるラベルが割り振られやすくなっている．平滑化項は，

E_smooth(L;X) = ∑

(u,v)∈C²

h(L_u, L_v, X_u, X_v) (2.10)

h(L_u, L_v, X_u, X_v) =







0 (L_u =L_v)

1

|X_u−X_v|+ 1 (Lu ̸=Lv)

(2.11)

(18)

で表される．

実際は，これらを異なる比率wで重み付けし，総和したエネルギー関数をセグメンテーションに用いる．

E(L;X) =w_atlasE_atlas(L) +w_dataE_data(L;X) +w_smoothE_smooth(L;X). (2.12) その他の項には，臓器の形状に関する事前知識を用いてポテンシャル関数を定義する研究がある[3]．

エネルギー最小化問題を解く手法として，[24]ではグラフカットを挙げている．グラフカットは，ラベルが二値の場合，つまり，ここでは M = 1の場合において用いられるアルゴリズムである．

2.4.4 グラフカットによるエネルギー最小化手法

グラフカットは，ラベルが二値(M = 1)の場合において，劣モジュラ条件

h_v₀_v₁(0,0) +h_v₀_v₁(1,1)≤h_v₀_v₁(0,1) +h_v₀_v₁(1,0) (2.13) を満たす場合，エネルギー大域最小解が得られるアルゴリズムとして知られている．

劣モジュラ条件を満たさない場合は，QPBO [18]というアルゴリズムを用いることにより近似解を得る．また，複数の対象物をセグメンテーションする場合，つまり，ラベルが多値(M ≥2)の場合に拡張するには，α 拡張[19]，融合移動[20]といったアルゴリズムを用いる必要がある．

QPBO

QPBO は，二値問題を解く過程において，どちらか一方のラベルを選択する以外に，

どちらのラベルを与えるか不明という状態を与える．

2 値をとる変数の実数値関数 f を擬ブール関数 (PBF Pseudo-Boolean Function)という．変数の番号の集合を U ={1, . . . , n} とすると，任意の PBF は次のように多項式の形に変換できる．

f(x1, . . . , xn) = ∑

S⊂U

cs

∏

i∈S

xi. (2.14)

ここで，xi ∈ {0,1}である．

(19)

一般に，二値ポテンシャルの式は，PBFに変換することができる．QPBOでは式(2.14) において S= 2 の場合，つまり 1階に相当する場合ならエネルギーの大域最小解を得ることができるが，S ≥ 3の場合，つまり階数が 2 階以上の場合は解を得ることは出来ない．従って，QPBOでエネルギーの大域最小解を得るためには，階数を削減する必要がある．[21]で提案されたアルゴリズムにより，階数の削減が可能になり，階数が 3 階以上の場合でも， QPBOで解けるようになった．

融合移動

融合移動は Fusion-move とも呼ばれ， α 拡張アルゴリズムを一般化したものである．

各ボクセルにおいて，任意のラベル二つを固定し，そのどちらを選択するかの二値問題を繰り返すことで解く．一般の融合移動は劣モジュラでないため実用的ではなかったが，

QPBOと組み合わせて用いることで，多値のラベルの場合でも近似的にエネルギーを最小化することができる．QPBO において，状態が不明と判断された場合，現在のラベルのままにしておき，エネルギーが増加しないことを保障するというアルゴリズムである．

図2.5においては，現在のラベリングC に対して，任意のラベリング P への変更を提案する．現在のラベリング C が採用されたボクセルを c ，提案のラベリングP が採用されたボクセルを p ，どちらを採用すればよいか不明なボクセルをφとすると，結果のラベリング C^′ は提案のラベリング p が採用された部分のみ提案されたラベルに置き換わっており，また φの部分は現在のラベルのままになっている．

2.4.5 高階エネルギーと高階グラフカット

ポテンシャル関数が最大k個のボクセルに依存するとき，そのポテンシャル関数をk−1 階のポテンシャル関数と呼ぶ．エネルギーが高々k 階のポテンシャル関数を含む時，そのエネルギーの階数を k とする．この階数が 2以上の時を高階と呼び，そのときのエネルギーを高階エネルギーと呼ぶ．従来のグラフカットでは，最小化が可能なエネルギーは一般に 2 階であり，また全ての場合において解けるわけではない[21]ので，実際はほぼ 1 階のエネルギーの使用が主流であった．近年になって，任意の高階エネルギーを 1 階のエネルギーに階数削減してからエネルギーを最小化するアルゴリズムである高階グラフカットが開発された[21]．これにより，高階エネルギーが実用的になった．本研究では，この高階グラフカットを用いて，高階エネルギーの最小化を行う．

図2.6に示すように，隣接する 2 ボクセルだと臓器の境界面の情報しかセグメンテーションに用いることが出来ないが，3ボクセル以上（高階）だと入り組んだ臓器の情報を

(20)

現在のラベリング

提案のラベリング

結果のラベリング′

：不明

図 2.5: fusion-move (2D の例)．現在のラベリング C に提案のラベリング P を与えてラベリングC^′ を得る．各ボクセルにおいて C とP のどちらのラベルを採用するかの二値問題は QPBOによって解かれる．

(21)

図 2.6: 1 階と高階の違い(2階の例)．高階を用いた方がセグメンテーションに，より多くの情報を用いることができる．

用いることや， 3 つの離れた臓器の情報を用いることが可能になるということがメリットとして挙げられる．デメリットとしては，高階エネルギーを用いると，セグメンテーションの実行時間がかかることや，実行に必要なメモリが階数に応じて増える，エネルギーの設計が複雑になることが挙げられる．

(22)

第 3 章高階エネルギーの提案

3.1 高階ポテンシャルの設計

胸部から腹部領域での多臓器セグメンテーションを扱う既存の研究において，セグメンテーションに用いるポテンシャル関数は，前章で述べたデータ項，アトラス項，平滑化項が主であり，また，それらのポテンシャル関数は最大 1 階であるものばかりであった．本研究では，既存研究では用いられていない，階数が 2 以上の高階ポテンシャルを提案する．高階ポテンシャルは， 1 階のポテンシャルと比べ，ポテンシャルの設計に用いるボクセルの数が増えるため，セグメンテーションに利用できる情報量が増える．本研究では，様々な高階ポテンシャルを提案し，実験による結果と高階ポテンシャルの有用性の評価を次章にて述べる．

データ項，アトラス項，平滑化項の三つだけを用いてセグメンテーションを行った場合，一部の臓器が正しくセグメンテーションされないことがある．図3.1は，臓器の一部のセグメンテーションが出来ていない例を示している．この例では，肝臓の左葉部分が正しくセグメンテーション出来ていない．CT 画像を見ると，肝臓領域内ではCT 値がほぼ一定であり，隣接ボクセル同士でも似たような CT 値をとっているように見えるにもかかわらず，セグメンテーションは正しく行われていない．セグメンテーションが正しく行われない理由として，アトラス項やデータ項は単一ボクセルの情報しか得られないため，隣接ボクセルで似たような CT 値を取っていても同じラベルをつけようと作用しないことや，平滑化項は隣接ボクセルまでの情報しか得られないため，さらにその隣以降に繋がるボクセルの CT 値が似たような値を取っていてもラベルが連続しているという情報を用いられていないことが挙げられる．高階ポテンシャルを用いると，隣接ボクセルだけではなく，さらにその隣以降に繋がるボクセルの CT 値を考慮してラベルを決めることが出来る．また， 3ボクセル以上の複数ボクセルのラベル値の関係を学習してポテンシャルに用いることも出来るので，より複雑なラベルの関係性を考慮出来る．従って，

高階ポテンシャルの導入はセグメンテーションにより良い影響をもたらすと考えられる．

(23)

図3.1: 肝臓が上手くセグメンテーションされていない例. 左はCT 画像，中央はGround Truth，右はデータ項，アトラス項，平滑化項の三つだけを用いたセグメンテーション結果を示す． CT 画像内の赤線およびGround Truth とセグメンテーション結果の白線がエッジを示す．肌色が肝臓，赤色が心臓を示す．セグメンテーション結果では肝臓左葉部分が上手くセグメンテーション出来ていない．

高階ポテンシャルを用いた場合，前章の (2.12)式は

E(L;X) =w_atlasE_atlas(L)+w_dataE_data(L;X)+w_smoothE_smooth(L;X) +w_higherE_higher(L;X).

(3.1) となる．

3.2 相互情報量によるクリーク配置の検討

提案する高階ポテンシャルは，データ項，アトラス項，平滑化項では捉えることが出来なかったボクセル間の関係を捉えられるものが良い．セグメンテーションでは，高階ポテンシャルのクリークに用いられる複数ボクセルのラベルによって，臓器内・臓器間の特徴といった何らかの関係を捉えられているクリーク配置が相応しい．

CT 画像内の総ボクセル数を N 個，エネルギーの階数を k 階とすると，高階エネルギーや非隣接のエネルギーを設計する際，クリークの配置はNP_k+1 通り存在し，セグメンテーションにエネルギーを与えるのに最適な配置をその中から選択するのは困難となる．

本研究では，最適な配置を選択する尺度として，相互情報量を導入する．2 変数の場合，相互情報量 I(X;Y)は事象 X, Y が起こる確率をそれぞれP(X), P(Y) とすると，

I(X;Y) =H(X) +H(Y)−H(X, Y) (3.2)

(24)

図 3.2: 相互情報量．相互情報量は 2 ボクセルのラベルの依存の度合いを示している．

H(X) =−P(X) lnP(X) (3.3)

H(Y) =−P(Y) lnP(Y) (3.4)

H(X, Y) = −P(X, Y) lnP(X, Y) (3.5) という計算式から求めることができる．本研究の場合，3.2のように，2ボクセルのうち，

一方のボクセルのラベルを X ，他方のボクセルのラベルを Y とした時に対応する．

相互情報量 I(X;Y) は X と Y の依存の度合いを示しており，X と Y が完全に独立している場合，相互情報量 I(X;Y)は 0 となる．相互情報量が高いほど，セグメンテーションにエネルギーを与えるクリークとして適していると考えられる．また，相互情報量は3 変数以上にも拡張ができるので，高階の場合の検討にも用いることができる．このようにして選択したクリークの配置を用いたエネルギーを設計し，セグメンテーションに利用する．

3.3 提案高階ポテンシャル

3.3.1 直角三角錐状に並んだクリークを用いた高階ポテンシャル

直角三角錐状に並んだクリークを用いた高階ポテンシャルを提案する．図3.3のような位置関係で並んだ 4 ボクセルを用いたクリークを定義する．この位置関係で並ぶボクセルのラベルの組の出現確率を学習してポテンシャルに用いる．

提案する高階ポテンシャルを E_tri4(L)，直角三角錐状に並んだボクセルの集合を C4 と

(25)

図 3.3: 直角三角錐状に並んだクリーク．この位置関係で並ぶボクセルのラベルの組の出現確率を学習する．

すると，提案ポテンシャルは

E_tri4(L) = ∑

(v0,v1,v2,v3)∈C4

f_v₀_v₁_v₂_v₃(L_v₀, L_v₁, L_v₂, L_v₃)

=− ∑

(v0,v1,v2,v3)∈C4

lnP(L_v₀, L_v₁, L_v₂, L_v₃).

(3.6)

のように表される．ここで，P(Lv0, Lv1, Lv2, Lv3)は，直角三角錐状に並んだボクセルのラベルの組の出現確率を示す．

3.3.2 離れたボクセルを用いたクリークのポテンシャル

離れたボクセル間のラベル値と距離，方向の関係を用いた高階ポテンシャルを提案する．図3.4のような位置関係で並んだボクセルを用いたクリークを定義する．この場合は高階ではないが， 1 階でもセグメンテーション精度改善の効果が期待される．1 階で上手くいけば，高階への拡張も可能だと思われる．

提案する高階ポテンシャルをE_apart(L)，基準となるボクセルv₀ とそのボクセルから離れたボクセルから成る組の集合をCA，ボクセル v₀ から |−→

d|離れたボクセルv_d=v₀+−→ d とすると，提案ポテンシャルは

E_apart(L) = ∑

(v0,vd)∈CA

f_v₀_v_d(L_v₀, L_v_d |−→ d)

=− ∑

(v0,vd)∈CA

lnP(L_v₀, L_v_d |−→ d).

(3.7)

のように表される．ここで，P(L_v₀, L_v_d |−→

d) は，離れた2ボクセルのラベルと方向と距離の組の出現確率を示す．

(26)

݀

⋱

଴ ଴

଴

଴ ଴

図 3.4: 離れたボクセルを用いたクリーク．この位置関係で並ぶボクセルのラベルと方向と距離の組の出現確率を学習する．

また，離れた 2ボクセルのうち，一方のボクセル v0 の座標を固定したポテンシャル E_apartF(L) = ∑

(v0,v_d)∈CA

f_v₀_v_d(L_v₀, L_v_d |−→ d , v₀)

=− ∑

(v0,v_d)∈CA

lnP(L_v₀, L_v_d |−→ d , v₀).

(3.8)

を提案する．一方のボクセルの座標を固定することで，エネルギーに位置の情報を加えられることが期待される．

3.3.3 高階平滑化項

従来の平滑化項では，隣接ボクセルのみの考慮，つまり 1 階であったが，その平滑化項を高階に拡張したものを提案する．図3.5のように， 2 階の場合は隣接 3ボクセルのラベルの関係ごとにエネルギーを変える．従来の平滑化項と同様に，隣り合うボクセルのラベルは激しく変動しないという事前知識をもとに設計する．隣接 3 ボクセルのラベルが全て同じ時のエネルギーを最も低くし，隣接 3 ボクセルのラベルが互いに異なる時のエネルギーを最も高くする．隣接 3ボクセルの中で隣り合う 2 ボクセルだけ同じラベルとなる場合のエネルギーを 2番目に低くし，隣接 3 ボクセルの中で真ん中のボクセルだけ異なるラベルとなる場合のエネルギーを 2 番目に高くする．

提案する高階ポテンシャルを E_HSL(L) ，隣接したボクセルの集合を C3 とすると，提

(27)

ݒ _଴ ݒ _ଵ ݒ _ଶ ݒ _଴ ݒ _ଵ ݒ _ଶ ݒ _଴ ݒ _ଵ ݒ _ଶ ݒ _଴ ݒ _ଵ ݒ _ଶ

図 3.5: 高階平滑化項のクリークとラベル値の組み合わせの例．ここでボクセルの色の違いはラベル値の違いを示す．並んだボクセルのラベル値によりエネルギーを変える．エネルギーの与え方は，ここで示す 4通りがある．

案ポテンシャルは

E_HSL(L) = ∑

(v0,v1,v2)∈C3

f_HSL(L_v₀, L_v₁, L_v₂) (3.9)

fHSL(Lv0, Lv1, Lv2) =











0 (L_v₀ =L_v₁ =L_v₂)

0.5 (L_v₀ =L_v₁ ̸=L_v₂andL_v₀ ̸=L_v₁ =L_v₂) 1 (L_v₀ =L_v₂ ̸=L_v₁)

2 (L_v₀ ̸=L_v₁ ̸=L_v₂)

(3.10)

のように表される．

さらに，上述した高階ポテンシャルに，各ボクセルの CT 値の差を反映させたポテンシャル

E_HSI(L;X) = ∑

(v0,v1,v2)∈C3

f_HSI(L_v₀, L_v₁, L_v₂, X_v₀, X_v₁, X_v₂) (3.11)

f_HSI(L_v₀, L_v₁, L_v₂, X_v₀, X_v₁, X_v₂)

=











0 (Lv0 =Lv1 =Lv2)

0.5

|Xv0−Xv2|+ 1 (L_v₀ =L_v₁ ≠ L_v₂andL_v₀ ≠ L_v₁ =L_v₂) 2

(|X_v₀ −X_v₁|+ 1) + (|X_v₁ −X_v₂|+ 1) (L_v₀ =L_v₂ ̸=L_v₁)

2 (L_v₀ ̸=L_v₁ ̸=L_v₂)

(3.12)

を提案する．ラベルが異なる場合に，異なるラベルをとるボクセル間の CT 値の差に応じてエネルギーを小さくすることで，異なるラベルを割り振られやすくしている．

また，隣接 3ボクセルが全部同じラベルとなる場合，隣接 3ボクセルのうち隣接2 ボクセルのみ同じラベルとなる場合，隣接 3 ボクセルのうち両端の2 ボクセルのみ同じラベルとなる場合，隣接 3ボクセルが全部異なるラベルとなる場合の出現頻度をそれぞれ

(28)

⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯

図 3.6: 離れた高階平滑化項のクリークとラベル値の組み合わせの例．ここでボクセルの色の違いはラベル値の違いを示す．並んだボクセルのラベル値によりエネルギーを変える．エネルギーの与え方は，ここで示す4 通りがある．

計算してエネルギーに用いる高階ポテンシャルを提案する．それぞれの場合の出現頻度を，C₁, C₂, C₃, C₄ とすると，提案ポテンシャルは

E_HSC(L) = ∑

(v0,v1,v2)∈C3

f_HSC(L_v₀, L_v₁, L_v₂) (3.13)

fHSC(Lv0, Lv1, Lv2) =











−ln C₁

C₁+C₂+C₃+C₄ (L_v₀ =L_v₁ =L_v₂)

−ln C₂

C₁+C₂+C₃+C₄ (L_v₀ =L_v₁ ̸=L_v₂andL_v₀ ≠ L_v₁ =L_v₂)

−ln C3

C₁+C₂+C₃+C₄ (L_v₀ =L_v₂ ̸=L_v₁)

−ln C₄

C₁+C₂+C₃+C₄ (L_v₀ ̸=L_v₁ ̸=L_v₂)

(3.14) と表すことが出来る．

さらに，各パターンの出現頻度を考慮した高階ポテンシャルのクリークを離して定義したポテンシャルを提案する．この場合，クリークは図3.6のように離して定義する．提案ポテンシャルを表す式は，式(3.13) ，式 (3.14)と同じになる．

3.3.4 CT 画像の輝度勾配を考慮した高階ポテンシャル

CT 画像の輝度勾配を用いた高階ポテンシャルを提案する．提案する高階ポテンシャルでは，CT 画像から抽出されたエッジを用いる．エッジとは CT 画像の輝度勾配の大きさが一定値以上の場所のことである．図3.1の Ground Truthで見られるように，エッジが存在するとき，その場所を境にラベル値が変化している．つまり，エッジを境に異なる臓器が存在することが分かる．エッジが有る場所を境にラベル値が異ならなけばならず，エッジが無い場所ではラベル値が連続しなければならない．図3.1において，セグメンテーション結果では， Ground Truth で肝臓と心臓の間に存在するエッジを境にしても，同じ心臓のラベルを付けてしまっている．これは，データ項，アトラス項，平滑化項だけでは，エッジを境に異なるラベルを付けるように作用していないためである．前章

(29)

の式 (2.11) で示したように，平滑化項は隣接する 2ボクセルに出来る限り同じラベルを割り振るように作用するが，2 ボクセルの CT 値の差が大きければ，与えるエネルギーを小さくすることによって 2 ボクセルに異なるラベルが割り振られやすくなるように作用するポテンシャルである．しかし，平滑化項はクリーク内でのエッジの存在を考慮していないため，エッジが有る場所でも同じラベルを割り振ることがある．そこで，提案する高階ポテンシャルでは，クリーク内でのエッジの有無によって与えるエネルギーを変化させるように設計することにより，平滑化項では用いられていないエッジの情報を取り入れる．また，高階を用いてクリークを設計することにより，エッジの周辺がエッジを境に異なるラベルが割り振られやすくする．

本項で提案する高階ポテンシャルを E_H(L)，高階ポテンシャルのクリークに用いる一直線上に並んだ n+ 1 個のボクセルの集合を Cn とする．図 3.7 は高階ポテンシャルに用いるクリーク配置の一例である．X¯ を画像 X のエッジ画像として，X(v) = 1¯ であれば，エッジ上にボクセル v は存在し， X(v) = 0¯ であればエッジがない場所にボクセル v が存在するとする．クリークC 内のエッジの有無を判定する関数は以下のように定義される．

g(C; ¯X) = 1−∏

v∈C

X(v).¯ (3.15)

また，提案するポテンシャルは

E_H(L) = ∑

C∈Cn

f_C(L_C), (3.16)

f_C(L_C) =







0 if g(C) = 0∧(∃l ∈ L^′,∀v ∈C, L_v =l), 1 otherwise,

(3.17)

と表され，ここでL^′ ⊂ L は，全ラベルの中から選択されたラベルの集合であり，このポテンシャルが影響するのは選択されたラベルの臓器だけである．

最終的には，式 (2.7)のエネルギーに，提案した高階ポテンシャルを加えたエネルギー E(L;X) =wAEA(L) +wDED(L;X) +wSES(L;X) +wHEH(L). (3.18) を用いてセグメンテーションを行う．ここで， w_H は高階項の重みであり， w_H = 0 の場合は，式 (2.7) のエネルギー関数と等しくなる．

L^′ ={Lliver} の場合， x 方向，つまり，体の左右方向に対して，肝臓は繋がる傾向に

(30)

図 3.7: x方向に並ぶクリーク． n+ 1 個の連続したボクセルから成る．

あるという事前知識を利用して，図3.7のような配置で設計した提案ポテンシャルを導入することにより，肝臓のラベルが右方向に伸びるように割り振られることが期待される．

L^′ = {L_liver, L_stomach, L_spleen, L_pancreas} の場合は，提案ポテンシャルを導入することにより，肝臓，胃，脾臓，膵臓の 4 つの臓器の領域が広がるように貼り付けられることが期待される．

3.3.5 CT 画像の輝度勾配を考慮して全臓器に作用する高階ポテンシ

ャル

前項で提案したポテンシャルを全臓器に作用するように拡張したポテンシャルを提案する．ここで提案するポテンシャルは，前項の提案ポテンシャルでL^′ =L の場合に該当するが，前項の提案ポテンシャルでは，各臓器がクリーク登録領域内に出現する確率を無視しているので，セグメンテーション結果に悪影響を及ぼしていた．そこで，本項で提案するポテンシャルでは，各臓器のクリーク登録領域内に出現する確率を考慮したポテンシャルを定義する．

本項で定義した提案ポテンシャルは三つある．一つ目は，

E_HA(L) = ∑

C∈Cn

f_C(L_C), (3.19)

f_C(L_C) =







−lnP(L_v)

max_L_v(−lnP(L_v)) if g(C) = 0∧(∀v ∈C, L_v =l, l̸= 0),

1 otherwise,

(3.20)

と定義する．ここで， P(L_v) はクリーク内の全ボクセルが同じラベルであった場合の，

各ラベルの出現確率を示す．EHA は，提案ポテンシャルによって返されるエネルギーを

2014 年度修士論文