招待論文

招待論文

^{医用画像論文特集}

医用画像に対する曲面モデルの非剛体位置合せのための統計形状モ デル

本谷 秀堅

^†a)

澤田 好秀

^†b)

Statistical Shape Model for Nonrigid Registration between Surface Model and Med- ical Images

Hidekata HONTANI

^†a)

and Yoshihide SAWADA

^†b)

あらまし 三次元画像中の臓器セグメンテーションを目的とした，Point Distribution Model (PDM)のSta- tistical Shape Model(SSM)の構築法について考察する．従来の線形なSSMには，部分空間モデルとMarkov

Random Fieldモデルがある．本論文では，これら従来のモデルが点間の条件付き独立性を正しく推定できない

可能性が高いことを指摘する．そして，正則化推定法を用いることにより，この独立性の推定精度を改善できる ことを述べる．その上で，graphical lassoを用いて構築したSSMを他のモデルと比較し，その推定性能が実際 に高くなっていることを示す．また，正則化の影響により対応する部分空間表現の次数が高くなることも実験的 に確認する．

キーワード 医用画像処理，非剛体位置合せ，画像セグメンテーション，統計形状モデル

1.

^{ま え が き}

三次元医用画像中の各臓器をセグメンテーションす ることは，診断支援システムや治療計画立案の支援シ ステムを実現するための重要な要素技術である．画像 中の対象領域をセグメンテーションする手法として，

数多くのアプローチが提案されている

[1]

〜

[3]

．これ らアプローチのうち，医用画像中の臓器セグメンテー ションに広く採用されているものは，臓器表面を表す 曲面モデルを用意し，この曲面モデルを画像へと位置 合せするアプローチである

[4]

．医用画像中の臓器は患 者ごとに形状も多様であり，また隣接臓器との接触や 病変などにより，臓器表面に画素値のコントラストが 存在しない場合も多い．これら要因のため，画素値の 違いやテクスチャの違いのみに基づいた手法では，精 度良く臓器領域をセグメンテーションすることは難し い．臓器表面の位置や形状の多様さを表す曲面モデル を用いることにより初めて，正確な自動セグメンテー

†名古屋工業大学，名古屋市

Nagoya Institute of Technology, Nagoya-shi, 466–8555 Japan a) E-mail: [email protected]

b) E-mail: [email protected]

ションが可能となる．

曲面モデルを画像へと位置合せする手法は，医用画 像処理

[5], [6]

だけではなく，レンジファインダーで取 得した三次元顔データの認識

[7], [8]

なども含め，様々 な応用を想定して多数提案されている．データの種類 や，位置合せに利用できる情報に依存して，用いるべ き位置合せ法も変化する．医用画像中の臓器セグメン テーションを応用とするとき，次に示す医用画像の性 質は，位置合せ法の良否を検討する際に重要である．

すなわち，

(1)

臓器形状は症例ごとに異なるが，

(2)

体 内における対象臓器のおおよその位置や姿勢，形状は 症例間で共通している．

(3)

レンジデータと異なり，臓 器表面の候補点を画像処理により抽出する必要があり，

(4)

隣接臓器との接触などにより，臓器表面候補点を 常に精度良く抽出できるとは限らない．

上記性質

(1)

により，剛体ではなく非剛体の位置合せ 法が臓器セグメンテーションには必須となる．非剛体位 置合せ法には，

Active Contour Model (ACM) [9]

や

Active Shape Model (ASM) [5], [10]

，

Nonrigid Iter-

ative Closest Point (Nonrigid ICP) [8]

など，様々な 手法が提案されている．これら手法は，曲面の形状の 統計を積極的に利用する手法と，曲面の滑らかさや剛

性などのみを利用する手法とに大別できる．性質に基 づけば，画像中の体の各部位をランドマークとして同 定することができれば，そのランドマークを指標とし て対象臓器を含む空間を正規化し，各臓器のおおよそ の位置や姿勢・形状を限定することが可能となる．例 えば肝臓を対象とする場合，肝臓を囲む肋骨の先端や 肺の下端をランドマークとして検出することができれ ば，それらランドマークに基づき肝臓を含む空間を正 規化し，肝臓のおおよその位置や姿勢，大きさを限定 できる．このことは，モデルの初期位置決定などに有 利な点であり，顔画像処理など，他の応用で常に期待 できる性質ではない．

性質

(3)

と

(4)

は，臓器セグメンテーションを難し くする要因になっている．仮に臓器形状がどんなに複 雑で，症例ごとに形状が大きく異なるものであったと しても，体内に対象臓器以外存在せず，臓器の内外で 画素値やテクスチャが異なるのであれば，臓器表面の 形状モデルを用いることなく，グラフカットなど画素 値に基づく従来法により臓器領域を正しくセグメン テーションできる．しかし実際には対象臓器以外から 多数の偽の臓器表面候補点が抽出される．また，隣接 臓器との接触などにより，画像中で臓器表面の一部が 頻繁に欠損する．対象臓器をセグメンテーションする には，臓器表面候補点のうち偽の点を除外し，欠損部 を補間しつつ，臓器表面を決定しなければならない．

先に述べたとおり，体型の正規化をすれば各臓器のお およその位置や大きさを決定できる．また，臓器形状 は症例ごとに異なるとはいえ，各臓器はそれぞれに固 有の形状特徴を有している．このため，偽の表面候補 点を除外したり，欠損部を補間したりする際には，対 象臓器の形状モデルを利用する方が，しない場合と 比べて位置合せの精度が高くなることを期待できる．

これが，対象臓器表面の統計形状モデル

(Statistical Shape Model (SSM)) [4], [10]

を利用する非剛体位置 合せ法が，臓器セグメンテーションでは広く採用され る理由である．本研究の目的は，この統計モデルの構 築法について考えることである．

曲面の表現には，

Point Distribution Model (PDM) [4], [10]

や

level set model [11], [12]

のほか，

spline

曲 面を用いる手法

[9]

や

current [13]

を用いる手法な どがある．本論文では，

PDM

における

SSM

の構 築法を提案する．

PDM

は曲面上に配置した

N

^点

{P

i

|i = 1 , 2 , · · · , N }

により，その曲面を表現する．

点

P

iの座標を

x

i

= [ x

i

, y

i

, z

i

]

^T で表すとき，曲面は

3 N

^{次元のベクトル}

x = [ x

^T1

, x

^T2

, · · · , x

^TN

]

^T により 表現される．

PDM

の非剛体位置合せは，与えられた データに基づいて，これら

3 N

個のパラメータを推 定することにほかならない．また，

PDM

の

SSM

は，

これら

3 N

個のパラメータの確率分布を表現するもの である．

SSM

は対象曲面の学習データ

M

^個の集合

{S

^l

|l = 1 , 2 , · · · , M }

^{を用いて構築される．}

PDM

のための代表的な

SSM

として，部分空間表 現と

Markov Random Field (MRF)

を挙げることが できる．前者は

ASM

で採用されている表現であり，

学習データに主成分分析を適用することにより構築さ れる．一方後者は点の座標間の関係を無向グラフによ り表現した，グラフィカルモデルである．

MRF

によ り

SSM

を表現する多くの手法

[14], [15]

においては，

隣接点間の統計は正規分布で表現される．このとき，

後に述べるとおり，

MRF

は対象曲面の形状の線形モ デルになっている．一方，前者の部分空間表現も曲面 形状の線形モデルである．ここで，同一の学習データ 集合を用いて，対象曲面の線形な

SSM

を，部分空間 表現と

MRF

表現のそれぞれに基づいて構築すること を考える．このとき得られる二つの

SSM

は，いずれ も同一の学習データから構築された線形モデルである にもかかわらず，曲面形状の表現としては異なるもの になり，位置合せの性能も同一にはならない．

本論文では，これら線形

SSM

の違いが，対応する グラフィカルモデルの構造に基づき見通しよく議論で きることを指摘する．その上で，グラフィカルモデル の構造を，

graphical lasso [16]

を用いて過学習を回避 しつつ，学習データに基づき決定する手法を述べる．

graphical lasso

は精度行列の疎性に基づく正則化によ り過学習を回避しつつ，精度行列を推定する手法であ る．正則化の強度を変化させると，グラフィカルモデ ルの構造も変化する．このグラフィカルモデルの構造 の変化により線形モデルの次数が変化することを示す．

また，肝臓を対象に，グラフィカルモデルの構造の変 化に対する位置合せの精度の変化を調査した結果を報 告する．グラフィカルモデルの構造変化に対する，部 分空間の次数の変化と肝臓に対する位置合せ精度の変 化は，本論文が初めて報告するものである．

2.

^線形な

SSM

のグラフィカルモデル表現 本 章 で は ，線 形 な

SSM

の グ ラ フィカ ル モ デ ル 表 現 に つ い て 確 認 す る ．

D

次 元 確 率 変 数

y =

[ y

1

, y

2

, · · · , y

D

]

^T の確率密度分布が，次式で示す正 規分布で表されるとする．

p ( y ) = 1 Z exp

− 1

2 ( y − µ )

^T

Λ( y − µ )

, (1)

ただし，

Z = (2 π )

^D/2

| Σ |

^{1/2 である．}

Λ

は精度行列 であり，共分散行列を

Σ

で表すとき，

Λ = Σ

⁻¹であ る．今

Λ

の

( i, j )

成分を

Λ

ijで表すことにし，非ゼロ の非対角成分に対応する要素のインデックスの集合を

Q = { ( i, j ) | Λ

_ij

= 0 }

^{で表す．右辺指数部の}

2

次の項 に注目すると，

Λ

_ij

= 0

のときのみ

2

変数の積

y

i

y

jが 現れることが分かる．このことから，上式の正規分布 は次のように

1

変数若しくは

2

変数の正規分布の積へ と分解できる．

p ( y ) =

D i=1

N ( y

i

)

(i,j)∈Q

N ( y

i

, y

j

) , (2)

ただし，

N ( · )

は

1

変数の正規分布を表し，

N ( ·, · )

は

2

変数正規分布のうち，

1

変数の正規分布二つの積 へと分解できないものを表す．今

D

個の確率変数

{y

i

|i = 1 , 2 , · · · , D}

から

y

iと

y

jを取り除いた

D − 2

個の確率変数の集合を

Y

\ijで表すことにする．式

(2)

より次に示す関係が成立することが分かる．

p ( y

i

, y

j

|Y

_\ij

) = p ( y

i

|Y

_\ij

) p ( y

j

|Y

_\ij

)

↔ Λ

_ij

= 0 . (3)

注目している

2

変数が，自分たち以外の確率変数の値 を固定したときに独立となるとき，それら

2

変数は 条件付き独立であるという．多変数正規分布において

Λ

ij

= 0

であることと，

y

iと

y

jが条件付き独立であ ることとは同値である．

多変数確率分布のグラフィカルモデル表現とは，各 変数に対応するノードを描き，注目する

2

変数が条件 付き独立ではないときに対応する二つのノードを辺で 結んだときに得られるグラフ表現のことである．正規 分布に従う多変数確率分布のグラフィカルモデルの構 造は，精度行列の非対角成分の非ゼロ成分の配置を表 現しており，各

2

変数間の条件付き独立性を示す．以 下，線形

SSM

の部分空間表現と通常の

MRF

表現の それぞれに対応するグラフィカルモデルの構造を示す．

ASM

で採用される部分空間表現は，縮退した共分 散行列をもつ正規分布に対応する．学習データ集合か ら計算された

D × D

の共分散行列

Σ

の固有値を降順

に並べたものを

λ

1

≥ λ

2

≥ · · · ≥ λ

Dで表し，それぞ れに対応する固有ベクトルを

u

iで表す．また，固有 ベクトルを並べた

D × D

の行列を次のように定める：

U = [ u

1

u

2

· · · u

D

]

．部分空間の次数を

d < D

^に定め るとき，この部分空間に対応する正規分布の精度行列

Λ

ASMは次式で表される．

Λ

ASM

= U

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣ 1 /λ

1

. . . 1 /λ

d

0 . . .

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦ U

^T

(4)

多くの場合固有ベクトルの全成分は非ゼロであり，上 式で定義される精度行列は，全成分が非ゼロになる．

このとき任意の

2

変数は条件付きに独立ではなく，グ ラフィカルモデルは各ノードが他の全てのノードと辺 で連結された，完全グラフになる．

一方，

PDM

の線形形状モデルを

MRF

で表現する 場合，多くの場合，測地距離の近い

2

点に対応する ノード同士のみが辺で連結される

[14], [17]

．このよう にして構築したモデルを，一般の

MRF

モデルと区別 するために，本論文では

GMRF (geodesic distance- based MRF model)

と呼ぶことにする．測地距離の 遠い

2

点は条件付きに独立であることが仮定されてお り，精度行列

Λ

GMRFは，ほとんどの成分がゼロの疎 な行列になる．線形な

GMRF

モデルは，共分散行列 が

Σ

GMRF

= Λ

⁻_GMRF¹ により表現される正規分布に対 応しており，

Σ

GMRFの固有値・固有ベクトルを求め れば部分空間表現を得ることもできる．

ASM

では，学習データより共分散行列を計算し，そ の固有ベクトルのうち，値の小さな固有値に対応する 固有ベクトルを消去することにより部分空間を構築す る．一方の

GMRF

では，グラフィカルモデルの構造 を学習データを解析する前に決定する．グラフィカル モデルの構造を決定することは，精度行列のうち，非 ゼロの要素の配置を決定することでもある．この配置 を決定した上で，学習データに基づき精度行列の非ゼ ロの要素の値を決定する．

ASM

と線形

GMRF

のそ れぞれで用いられるモデルの違いは，グラフィカルモ デルの構造の違い，若しくは線形部分空間表現の次元 や基底ベクトルの違いにより比較することができる．

3.

グラフィカルモデル上の推論による非剛 体位置合せ

本章では，

PDM (Point Distribution Model)

の

SSM (Statistical Shape Model)

の表現にグラフィカ ルモデルを採用したときの，

PDM

の位置合せアルゴ リズムを記す．

PDM

の

SSM

は，

PDM

を構成する各 点の座標を確率変数とする，多変数の確率モデルにな る．新たに入力された画像データに対する

SSM

の位 置合せは，その確率モデルを用いて，各点の座標値を 推定することにより行われる．この推定には，最ゆう 推定や

MAP

推定，ベイズ推定など，通常の統計的推 定の枠組みが利用できる．本研究ではベイズ推定の枠 組みを採用し，グラフィカルモデル上で，各点の座標 の周辺事後分布を推論する手法を採用する．最ゆう推 定や

MAP

推定など点推定により位置合せを行う場合 と異なり，周辺事後分布を推定することにより位置合 せ結果の確度も自己評価できるようになる．

SSM

が グラフィカルモデルにより陽に表現できていれば，そ れが共分散行列の固有ベクトルに基づき構築した部 分空間表現であっても，近接点同士を結んで構築した

MRF

表現であっても本手法により周辺事後分布を計 算することができる．以下，

PDM

の

SSM

のグラフィ カルモデル表現について述べた後，推論法を説明する．

3. 1 Statistical Shape Model

以下，

SSM

に従い変形する曲面をテンプレートと呼 び，画像より抽出された表面候補点の集合をターゲッ トと呼ぶ．テンプレートを構成する各点

P

iの座標を

x

i

( i = 1 , 2 , · · · , N )

で表し，テンプレート全体を

3 N

次元のベクトル

x = [ x

^T1

, x

2^T

, · · · , x

^TN

]

^Tで表す．入 力画像を

I = I ( x, y, z )

で表し，点

P

iを中心とする一 辺

L

ボクセルの立方体内部の画像を

I

iで表す．

I

iは 立方体内部のボクセル値を並べたもので，

L

^3次元の ベクトルである．ただし，各画像は体内の何らかのラ ンドマークに基づいて体型の正規化がなされているも のとする．正規化に用いたランドマークを利用するこ とにより，患者間で共通の座標系により

P

iの位置を 記述することができる．

本論文で扱う

SSM

は，テンプレート

x

^と画像

{I

i

}

の同時分布を次式により表現する．ここで

E

^は，グラ フィカルモデル上で辺により連結されている頂点のイ ンデックスの組の集合である．

p ({x

i

}, {I

i

}) =

N i=1

p (x

i

) p ( I

i

|x

i

)

(i,j)∈E

p (x

i

, x

j

) .

(5)

上式右辺の

p ( x

i

)

は各点の事前分布であり，各点

P

i

が正規化後の体内のどの位置に出現しやすいかを表し，

p ( x

i

, x

j

)

は点

x

iと

x

jの相対位置関係を表す．これ らを合わせた

p ( x ) =

i

p ( x

i

)

(i,j∈E)

p ( x

i

, x

j

)

に より，対象臓器表面の位置と形状の統計が表現されて いる．本論文では，

p ( x

i

)

と

p ( x

i

, x

j

)

のそれぞれを正 規分布で表現する．このとき，

p (x)

も正規分布に従 う．

p ( I

i

|x

i

)

は位置

x

iが与えられたときの，

P

i周辺 のアピアランスの事後分布を表している．本論文では，

この事後分布

p ( I

i

|x

i

)

も正規分布でモデル化する．右 辺に現れるこれら三つの正規分布のパラメータは，い ずれも学習データに基づき推定する．推定法について は後に示す．

上記モデルと

ASM

や

GMRF

モデルの違いを二点 指摘する．一つ目の相違点は，アピアランスに関する モデルの有無である．

ASM

や

GMRF

は，

(1)

画像 より対象候補点を抽出してターゲットとし，

(2)

この ターゲットに対してテンプレートを位置合せする．画 像からターゲットを抽出する際に用いるアピアランス モデルを，

ASM

や

GMRF

は含んでいない．

ASM

や

GMRF

は，エッジ抽出などにより対象候補点が抽出 済みであることを前提とし，位置合せの手続きのみを 提供する．一方，式

(5)

はアピアランスに関するモデ ル

p ( I

i

|x

i

)

を含んでおり，各点の位置

x

iのゆう度分 布の計算に利用できる．この計算により，各

P

iに対 応する対象候補点の分布が取得できるため，入力画像 アピアランスも考慮して

x

iの周辺事後分布を計算で きる．画像中の対象候補点の分布は，位置合せの精度 や確度に大きな影響を与える．例えば臓器表面の近傍 から多数の偽の候補点が抽出される部位があれば，そ の部位における位置合せの確度は低下する．式

(5)

は 偽候補点の分布も表現できるが，

ASM

や

GMRF

で は表現できない．

ASM

や

GMRF

の枠組みでは位置 合せの確度評価が困難なことが分かる．

二つ目の相違点は事前分布

p ( x

i

)

の有無である．統 計モデル構築の際，本論文では患者の体型の違いを 正規化するが，

ASM

や

GMRF

においては対象の位 置・姿勢・大きさの違いを

Procrustes Analysis [18]

などにより正規化する．体型正規化により，臓器の位 置・姿勢・大きさがどの程度正規化されるかは，体型

図1 ゆう度分布p(Ii|xi)の例 Fig. 1 Examples ofp(Ii|xi).

正規化の際に参照するランドマークの配置に依存する．

式

(5)

は形状の多様さと位置・姿勢・大きさの多様さ を区別しないモデルであり，一方

ASM

や

GMRF

は，

形状の多様さのみを表現し，対象の位置・姿勢・大き さの違いは正規化される．

ASM

や

GMRF

は位置に 関するモデルを含まないため，モデルの初期位置決定 に必要な情報を別途提供しなければならない．

3. 2

グラフィカルモデル上での周辺事後分布の推論 与えられた画像に対して，モデル

(5)

に基づいて，

PDM

の各点の座標

x

iの周辺事後分布

p ( x

i

|I )

を求め るアルゴリズムを記す．

p ( x

i

|I )

の期待値若しくは最 大値がテンプレートの位置合せ結果を示し，標準偏差 の大きさが推定の確度の低さを表す．アルゴリズムは 次のとおりである．

（

1

） 入力画像の体型を正規化する．

（

2

）

p ( I

i

|x

i

)

に基づき，

x

iのゆう度分布を求める．

（

3

）

p ( x

i

) p ( I

i

|x

i

)

を求め，パーティクルの集合に     より表現する．

（

4

）

MCMC

を利用し，

p (x

i

, x

j

) (( i, j ) ∈ E)

に     従ってパーティクルをサンプリングし直し，

    周辺事後分布を得る．

位置

x

iを固定するときの事後分布

p ( I

i

|x

i

)

を正規分 布で表現しても，画像

I

iを固定するときのゆう度分 布

p ( I

i

|x

i

)

は正規分布に従うとは限らない．実際図

1

に示すとおり，ゆう度分布はテンプレートの点ごとに 異なり，臓器表面に沿った非線形なものになる．この ため，

p ( x

i

) p ( I

i

|x

i

)

や

p ( x

i

|I )

は多数のパーティクル を用いて，ノンパラメトリックに表現することになる．

手続き

(3)

までにおいて，事前分布

p (x

i

)

とゆう度分 布

p ( I

i

|x

i

)

が計算される．いずれもテンプレートの 点

P

iごとに計算できるものであり，点間の相対位置 関係

p ( x

i

, x

j

)

は考慮されていない．すなわち，手続

き

(3)

までの処理結果には，対象臓器の形状モデルは 反映されていない．式

(5)

において形状は，

p (x

i

, x

j

) (( i, j ) ∈ E )

により表現されている．手続き

(4)

におい て形状モデル

p ( x

i

, x

j

)

が用いられ，最終的な周辺事 後分布が計算される．

上記手続き

(4)

の周辺事後分布の計算には

MCMC

が利用できる．ここでは，

MCMC

によるサンプリン グ法の代表例である，

Gibbs Sampling (GS)

による サンプリング法の概略を示す．

以下，各点

P

iの分布を表す

K

^{個のパーティクルの} 集合を

{z

^ki

|k = 1 , 2 , · · · , K}

により表現する．

z

^ki は パーティクルの位置を表す三次元ベクトルである．ま た，

Z

_\i^k∗により，各点

P

iに対応するパーティクルを 一つずつ含む集合

{z

^ki

|i = 1 , 2 , · · · , N }

^から

z

^ki^∗だけ を取り去ったものを表すことにする．

Z

_\i^k∗ は

N − 1

個のパーティクルの集合である．

N

個の点

{P

i

}

それ ぞれの分布を表すための

K

個ずつのパーティクルで の集合は次の手続きにより得ることができる．

-

正数

W

0と

W

1を定める．

W

1

= W

0

+ K

とする．

- {z

¹i

|i = 1 , 2 , · · · , N }

を

p ( x

i

) p ( I

i

|x

i

)

に従って サンプルする．

- for k = 2 , 3 , · · · , W

¹

- for i

^∗

= 1 , 2 , · · · , N

-

条件付き確率

p ( x

i^∗

|I, Z

_\i^k−∗¹

)

に従い

1

点 サンプルし，

z

^ki^∗として採用する．

-

各点

P

iに対して得られた

W

1

= W

0

+ K

^個 のサンプルの集合

{z

^ki

|k = 1 , 2 , · · · , W

1

}

^のう ち，

k ≤ W

0 のサンプルを削除し

K

個だけ残 す．右肩のインデックスが

1

から始まるように

z

^ki

← z

^k−W_i ^{0と修正し，終了する．}

条件付き確率

p ( x

i^∗

|I, Z

_\i^k−∗¹

)

は式

(5)

の

x

i

( i = i

^∗

)

に

z

^k−_i ¹を代入することにより得られる．初期段階で 取得されるサンプリングは，モデル式

(5)

に従わない 可能性があるため，最初の

W

0回までに取得されたサ ンプルは棄却する．適切な

W

0の値は初期サンプル

{z

¹i

|i = 1 , · · · , N}

の精度などに依存し，実験的に決 定される．各分布の表現に用いるサンプルの数

K

^は，

大きいほど分布の表現能が高くなるが，その分位置合 せに要する時間も長くなる．ただし，

K

^{を小さくしす} ぎると，分布を正しく表現できなくなり，位置合せの 性能が劣化する．適切な

K

の値も，位置合せの精度 と計算時間の双方を調べながら，実験的に決定する．

多峰性の複雑な分布からの効率良くサンプルを取得す

るためには，初期サンプルを変更しつつ並列でサンプ リングを行うことが有効である．

4.

^{統計モデルの構築}

式

(5)

に示した

SSM (Statistical Shape Model)

を 学習データに基づき構築する．本節では，その構築法 を述べる．具体的には，式

(5)

右辺の各分布と，グラ フィカルモデルにおいてノード間を連結する辺の集合

E

を決定する．先に述べたとおり，式

(5)

のモデルの 中で，臓器表面の形状は

p ( x

i

, x

j

) (( i, j ) ∈ E )

で表現 される．辺の集合

E

を決定することは，グラフィカル モデルの構造を決定することであり，臓器の形状の統 計を表現する上での重要なプロセスである．本節では，

式

(5)

右辺の各分布の推定法を述べ，グラフィカルモ デルの構造を

graphical lasso

を用いて決定する手続 きを記した後，グラフィカルモデルの構造の違いとモ デルにより表現される形状の多様性について考察する．

4. 1

確率分布の推定

M

症例の学習用画像それぞれの体型を正規化し，

画像の集合

{I

^l

|l = 1 , 2 , · · · , M }

^{を用意する．それ} ぞれの画像中の臓器表面を手作業で抽出し，学習用 の曲面の集合

{S

^l

|l = 1 , 2 , · · · , M }

を得る．

PDM

の

SSM

を構築するために，各曲面

S

^lに，

N

^{点の対応点}

{P

i^l

|i = 1 , 2 , · · · , N }

^{を生成する．各点}

P

i^lの座標を三 次元ベクトル

x

^liで表し，

S

^{l上の点の配置を}

3 N

^次元ベ クトル

x

^l

= [ x

^lT1

, x

2^lT

, · · · , x

^lTN

]

^Tで表す．対応点生成 法には

Entropy-based particle system [19]

や

Gen- eralized Multi-Dimensional Scaling (GMDS) [20]

な どが利用できる．各点

P

i^lを中心とした一辺

L

^ボクセ ルの立方体内部の画素値を並べ，

L

^{3 次元のベクトル}

I

i^lも生成する．

本研究では式

(5)

の右辺の各確率分布を正規分布 で表現する．このとき，学習データから推定すべき パ ラ メ ー タ は 分 布 の 平 均 と 共 分 散 行 列 で あ る ．点

P

i の事前分布

p ( x

i

)

の平均と共分散は，対応点の 集合

{x

^li

|l = 1 , 2 , · · · , M }

から計算する．また，局 所的なアピアランス

I

iの条件付き確率分布

p ( I

i

|x

i

)

の 平 均 と 共 分 散 は ，局 所 的 な ア ピ ア ラ ン ス の 集 合

{I

i^l

|l = 1 , 2 , · · · , M }

^{より計算する．}

相対位置関係の統計モデルを構築するには，分布

p ( x

i

, x

j

)

の推定だけではなく，グラフィカルモデル における辺の集合

E

も推定する必要がある．後者 の推定については次節で述べることとする．前者に ついては，

( i, j ) ∈ E

^を満たす

2

点

P

iと

P

jに対し

て，

x

^lij

= [ x

^lTi

, x

^lTj

]

^T

( l = 1 , 2 , · · · , M )

を用意する．

そして，

p ( x

i

, x

j

)

の平均と分散は学習データの集合

{x

^lij

|l = 1 , 2 , · · · , M }

より計算する．

4. 2

グラフィカルモデルの構造推定

PDM

の

SSM

をグラフィカルモデルで表現する場 合，グラフの各ノードは

PDM

を構成する各点の座標 を表現する．グラフィカルモデルの構造を決定するこ とは，どのノード間を辺で連結するかを決定し，

E

を 定めることである．

E

を決定することは，

2.

^でも述べ たとおり，

PDM

を構成する点間の条件付き独立性を 推定することに他ならない．

P

iと

P

jが条件付き独立 でないとき，そのときに限り，

P

iと

P

jは辺で連結さ れる．辺で連結された

2

点についてのみ，相対位置関 係を表す分布

p (x

i

, x

j

)

を推定し，式

(5)

右辺に組み 込むことになる．

変数間の条件付き独立性の有無は，精度行列の非対 角成分の有無により判定できる．先に線形モデルの例 として，

ASM

と

GMRF

モデルに言及した．

ASM

で は

SSM

のグラフィカルモデル表現は利用しないが，主 成分分析に基づいて構築した部分空間表現をグラフィ カルモデルへと変換すると，式

(4)

に示す精度行列

Λ

ASMに基づき，各点が他の全点と連結された完全グ ラフの構造をもつグラフィカルモデルが得られること になる．このとき，例えば肝臓の左葉端点と右葉端点 のように，測地距離の意味で互いに遠く離れた

2

点で あっても，条件付き従属になる．これは直感に反する．

隣接臓器との接触により局所的に頻繁に変形すること も考慮すると，臓器上の各点が，他の全点と条件付き で従属関係にあるとは考えにくい．一方の

GMRF

モ デルは測地距離の観点から近接する点間のみを辺で連 結するモデルであり，局所的な要因のみで形状の多様 性を表現するモデルである．近接点間が条件付きで従 属関係にあるとみなすモデルは，臓器表面の滑らかさ などに基づき支持できる．しかし，互いに離れた任意 の

2

点が常に条件付き独立か否かは自明ではない．

学習データに基づき精度行列を推定する場合，過学習 に注意する必要のあることが知られている

[16]

．

ASM

が採用する

SSM

は，通常どおり次式により算出され る共分散行列

Σ

の固有値・固有ベクトルに基づいて構 築される．次式は学習データ

{x

^l

|l = 1 , 2 , · · · , M}

^に 対する，共分散行列の最ゆう推定である．

Σ = 1 M

M l=1

( x

^l

− µ )( x

^l

− µ )

^T

, (6)

ただし，

μ = (

l

x

^l

) /M

は平均ベクトルである．最 ゆう推定の過学習を回避するアプローチには，ベイズ 推定や正則化法がある．前者は推定対象の事前分布を 導入する手法であり，後者は推定対象に関する拘束を 導入する手法である．いずれも推定対象に関する先見 知識を導入する点は共通している．本研究では，精度 行列を疎行列とする先見知識を導入する．すなわち，

臓器上の各点は，その近傍を含む少数の点とのみ条件 付きに従属であり，それ以外の点とは独立であると仮 定する．精度行列の疎性の仮定に基づき過学習を回 避しつつ精度行列を推定する正則化法に，

graphical lasso

がある

[16]

．推定すべき精度行列を

Λ

GLで表す．

graphical lasso

は次式に示すコスト関数を最小化する ことにより，

Λ

GLを推定する．

F (Λ) = − log det Λ + trace (ΣΛ) + ρ Λ

L1

. (7)

上式右辺の最初の

2

項が精度行列

Λ

の負の対数ゆう 度であり，最後の項，

Λ

L1，が精度行列の疎性に由 来する正則化項であり，

L0

のフロベニウスノルムを

L1

ノルムに緩和したものである．上式のコスト関数

F (Λ)

は凸関数であり，最小値を与える唯一の解をも つ．

Λ

GL

= arg min

Λ

F (Λ)

とおく．推定された精度 行列の非対角成分に基づいてグラフィカルモデルの辺 の集合を次のように定める．

E = { ( i, j ) | (Λ

GL

)

ij

= 0 }. (8)

式

(7)

右辺の

ρ

は正則化の強度を制御する正の係 数であり，

Λ

GLは

ρ

とともに変化する．

ρ = 0

のと きは疎性に関する正則化項が考慮されないため，得ら れる精度行列は式

(6)

を用いて普通に計算した共分散 行列の逆行列と一致する．このとき，先に述べたとお り，グラフィカルモデルは完全グラフになる．また，

ρ

を十分大きくすると非対角成分が全てゼロとなり，

式

(5)

のモデルが対象臓器の形状情報を表現しなくな る．

E = φ

となるため，グラフィカルモデルには辺が 一本も含まれないようになる．正則化の係数

ρ

^の値は，

実験的に決定する．ここでは，位置合せの精度が最良 となる係数を実験的に探索し，決定する．精度行列の 疎性の仮定を導入することにより過学習が抑制される ため，完全グラフの場合や

GMRF

モデルの場合より も位置合せの性能の向上を期待できる．このことを次 節における肝臓の位置合せを題材とする実験により確 認する．

4. 3

線形モデル構築の際の考え方

前節に述べた手法では，精度行列の疎性を仮定し，

グラフィカルモデルの構造を推定する．一方，

ASM

が採用する主成分分析に基づく手法では，共分散行列 の固有ベクトルのうち，小さな固有値に対応するもの を削除することにより部分空間表現を構築する．双方 共に線形なモデルであり，双方のモデルを部分空間で 表現することも，グラフィカルモデルで表現すること も可能である．同一の学習データから出発しても，異 なるアプローチをとることにより，異なるモデルが構 築される．

学習データに基づき統計モデルを構築するとき，構 築されたモデルの善し悪しを評価する尺度として，

generalization, specification, compactness

の三項目 を使うことができる

[21]

．

generalization

は学習デー タに含まれない未知の対象をどの程度表現できるか を示す指標であり，

specification

は対象以外をどの 程度表現から除外できているかを示す指標である．

compactness

はモデルを表現するパラメータ数であ る．

ASM

においてモデルを構築する際，部分空間の 次数に依存して位置合せの精度は変化する．これは，

部分空間の次数がモデルの

compactness

だけではな く，

generalization

や

specification

にも影響を与える からである．小さな固有値に対応する固有ベクトルを 表現から除外することは，学習データが含む雑音成 分などを削除することにより，モデルの自由度を抑制 し，

specification

を向上させることを期待した手続き である．一方，精度行列の疎性を正則化項として導入 することは，点間の条件付き独立性を強制することで あり，独立性が向上すればモデルの自由度が高くなる．

graphical lasso

に基づきグラフィカルモデルの構造を 決定することは，対応する部分空間の次数を高くし，

generalization

を向上させる手続きである．通常の部 分空間法と

graphical lasso

を利用するグラフィカルモ デルの構造決定法では，学習データからモデルを構築 するときによって立つ考え方が逆であることが分かる．

過学習は一般に学習データ数が少なすぎるときに起 こる．過学習を回避するために導入される先見知識 は，全ての対象に共通する性質を表現しなければなら ない．言い方を変えれば，先見知識を導入したせいで，

表現不可能となる対象があってはならない．このため，

例えば正則化により過学習を抑制するアプローチを採 用した場合，正則化項を重視すればするほど，モデル はより一般化され，

generalization

の性能が向上する．

ただし，先見知識により導入される性質は，表現すべ き対象以外のものにも共有されるはずである．対象以 外には一切共有されない性質の表現に成功すれば，先 見知識だけで理想のモデルが構築できる．このことか ら，先見知識を重視するほど，モデルの

specification

は劣化する．次章における実験では，式

(7)

中のパラ メータ

ρ

の変化に対する部分空間の次元の変化を確認 する．

5.

^{実 験}

X

線

CT

画像中の肝臓を対象とし，式

(7)

中のパ ラメータ

ρ

の変化に対する線形モデルの変化を調べ た実験結果を報告する．実験に用いた画像は医用画像

(CAD)

データベース

[22]

の三次元

X

線

CT

画像（門 脈相）であり，全て肝臓内に癌をもつ症例である．空 間分解能は

0 . 625 × 0 . 625 × 1 mm

³である．学習に用 いた症例数は

M = 18

であり，それとは別にテスト用 に画像を

18

症例用意した．スライス画像の肋骨と背 骨を手掛りに自動的に体型正規化を行った．また画素 値の分布もヒストグラムに基づいて正規化した．学習 用の画像とテスト用の画像の双方の肝臓表面を手作業 で抽出し，前者を統計モデルの構築に利用し，後者は 位置合せ結果の精度を評価するために利用した．

PDM

の表現に利用した点の数は

N = 1498

であり，学習の 際の対応点生成には

GMDS

を利用した．

5. 1

グラフィカルモデルの構造と部分空間の次数

graphical lasso

により推定された精度行列を

Λ

GL

= Λ

GL

( ρ )

で表す．引数

ρ

^は，式

(7)

において正則化項の 重みを制御するパラメータである．推定された精度行 列に対応する部分空間の次数を求めるため，対応する 共分散行列を

Σ

GL

( ρ ) = Λ

⁻_GL¹

( ρ )

として求める．ただ し，右辺の逆行列は一般化逆行列である．

Σ

GL

( ρ )

の 固有値・固有ベクトルを求め，固有値の分布の変化調 べた．図

2

に結果を示す．

ρ

^{の値が上から順に}

0.001

，

0.01

，

1.0

と大きくなるにつれ，固有値の分布は次第に 均一な分布へと変化している．各

ρ

の値において得ら れたグラフィカルモデルの辺の総数は，順に

2 . 3 ×10

⁵，

2 . 2 × 10

⁵，

1 . 4 × 10

⁵であった．正則化項の係数の増 加とともに辺の本数が減少し，グラフィカルモデルが 疎になるとともに，対応する部分空間の次元が増加し たことが分かる．図には掲載してないが，

ρ = 10

で は辺の本数はゼロとなり，全ての固有値が同一の値と なった．

図

3

に，係数

ρ

の変化に対する構造の変化の一部

(a)ρ= 0.001

(b)ρ= 0.01

(c)ρ= 1.0

図2 graphical lassoにおける正則化項の係数ρと，共 分散行列の固有値の分布の関係

Fig. 2 Relationship between the coefficient of the regularization term,ρ, for the graphical lasso, and the distribution of the eigen values.

を図示する．全ての辺は図示できないため，適当に 選択した

1

点と接続している辺のみを示す．図中に 緑色の四角形により選択した点を示し，赤線により

(a)ρ= 0.001

(b)ρ= 0.01

(c)ρ= 1.0

図3 graphical lassoにおける正則化項の係数ρと，グ ラフの構造の関係

Fig. 3 Relationship between the coefficient of the regularization term,ρ, for the graphical lasso, and the graph structure.

選択した点と条件付きに従属な点を示す．上から順に

ρ = 0 . 001

，

0 . 01

，

1 . 0

であり，図

2

のそれぞれと対 応している．係数の値の増加に伴い，接続されている 辺のうち遠方の点との接続数が減少し，近傍の点との 接続のみへと変化していることが分かる．グラフにお

表1 各係数ρで得られたグラフィカルモデルの辺の本数 とノードの次数

Table 1 The number of edges and the degree of nodes in the graphical model obained with each value ofρ.

ρ 0.001 0.01 0.1 1.0 10.0

# of edges (×10⁵) 2.3 2.2 2.0 1.4 0 Maximum degree 314 341 187 218 0 Mean degree 156 151 115 94 0 Minimum degree 108 101 65 12 0

いて各ノードに接続されている辺の本数は次数と呼ば れる．各

ρ

において得られたグラフの辺の本数と，全 ノード中の次数の最大値・平均値・最小値をそれぞれ 表

1

に示す．

ρ

の増加とともに，次数の平均値，最小 値も単調に減少していることが分かる．一方で，次数 の最大値は単調減少していないことも確認できる．こ のことは，

ρ

の値を変化させて得られたグラフ構造は 密な構造と疎な構造の間に包含関係がないことも示唆 している．すなわち，

ρ

の値を次第に大きくしていく とき，グラフィカルモデルから次第に辺が消えていく のではなく，辺の消失とともに新規の辺の追加が起こ る．式

(6)

により求めた共分散行列の逆行列を

Λ

で表 すとき，

Λ

の値がしきい値より小さい成分をゼロとす ることを考える．この手続きでも条件付き従属性の弱 い

2

変数を独立に扱うことが可能となる．このしきい 値処理による簡便な手法を採用するとき，しきい値を 次第に大きくして疎性を高くするとグラフィカルモデ ルから辺は順次消えていき，新規な辺の追加は起こら ず，密な構造が疎なグラフを包含する．共分散行列の ゆう度も考慮した式

(7)

の最小化により得られる精度 行列は，

Λ

の単純なしきい値処理とは異なることが分 かる．

5. 2

グラフィカルモデルの構造と臓器位置合せの 精度

テスト画像に対してモデルとの位置合せを行い，そ の精度を評価した．精度評価は，位置合せの結果と，

テスト画像

I

^jより手作業で抽出した肝臓表面

S

^jとの 距離に基づき行った．

画像

I

^{j上での点}

P

iの位置合せの誤差を

Δ

^j_iとする．

このとき，

Δ

^j_i は求めた周辺事後分布

p ( x

^j_i

|I

^j

)

の期待 値

x ˆ

^j_iと

S

^j上の最近傍の点

x ¯

^j_i との距離を表す．

Δ

^j_i

= x ˆ

^j_i

− x ¯

^j_i

.

また，

M

^{症例の画像に}

N

^点の

PDM

を位置合せした 結果得られる平均誤差

E

は次式で評価した．

図4 係数ρの変化に対する位置合せ誤差の変化 Fig. 4 Change of the registration error with respect

to the value ofρ.

図5 同一の学習データより構築した3種類の線形モデル の位置合せ精度の比較．赤：graphical lassoにより 構築．緑：隣接点間を連結して構築したGMRFモ デル．青：完全グラフ構造をもつモデル．横軸は症 例番号．縦軸が位置合せの誤差．

Fig. 5 Comparison of registration accuracy among three linear SSM constructed from an iden- tical set of training data. Red: Model con- structed by the graphical lasso, Green: GMRF model constructed by linking two neighboring points, Blue: A model represented by a com- plete graph. The horizontal axis indicates the image index, and the vertical one shows the registration error.

E = 1 M

1 N

M j=1

N i=1

(Δ

^j_i

)

²

.

ρ

の値を変えつつ位置合せの平均誤差

E

^{を評価した結} 果を図

4

に示す．図に示したとおり，

ρ = 0 . 001

のと きに位置合せの精度は最良であった．

次に，位置合せ精度が最良であった

ρ = 0 . 001

の ときに得られるグラフィカルモデルと，完全グラフ

( ρ = 0)

と，

GMRF

モデルの

3

種類のモデルの，

3

種 類のモデルの位置合せ精度を比較した．

GMRF

モデ

(a) Graphical Lasso (ρ= 0.001)

(b) GMRF model

(c) Complete Graph (ρ= 0)

図6 位置合せの平均誤差分布

Fig. 6 Distribution of the averaged error of registra- tion.

ルは測地距離がしきい値より小さな

2

点間を結ぶこと により構築した．それぞれのモデルの位置合せ誤差

E

を症例ごとに示したグラフを図

5

に示す．図中の赤が

graphical lasso

により得られた結果であり，緑と青が それぞれ

GMRF

モデルと完全グラフにより得られた 結果を示す．いずれの症例においても

graphical lasso

により構築したモデルの精度が最良であることが分か る．

p < 0 . 01

で統計的に有意に

graphical lasso

によ り構築したモデルが最良であることを確認した．

各点

P

iにおける位置合せの平均誤差

Δ

iは，各画像 で生じた誤差を画像間で平均することにより，次のよ うに評価できる：

Δ

_i

=

j

Δ

^j_i

/M

^{．この平均誤差の} 大きさに応じて，平均形状をもつ肝臓の各位置を着色

した結果を図

6

に示す．図は上から順に，

graphical lasso

により得られたモデル，

GMRF

モデル，完全グ ラフの平均誤差分布である．いずれの手法でも右葉端 点で推定精度が悪いが，

graphical lasso

のモデルによ る位置合せの精度が総じて高いことを視認できる．以 上の実験結果より著者らが得た考察は次のとおりであ る．すなわち，臓器

PDM

の

SSM

に線形モデルを採 用し，各点の周辺事後分布を推定することにより位置 合せをするのであれば，

(1)

まず学習用臓器データに 対応点を配置し，

(2)

適切な正則化法を導入するなど して精度行列を推定し，グラフィカルモデルの構造を 決定する．そして，

(3)

そのグラフィカルモデル上で 各点の座標を推定するのが良い．ただし，正則化項の 係数は交差検定などにより調整する必要がある．本論 文では正則化項に精度行列の疎性項を導入したが，真 の精度行列が本当に疎であるかどうかは議論の余地が あるかもしれない．今後の検討課題としたい．また，

ρ

の各値における

generalization

，

specification

などの 定量評価も今後の課題とする．

6.

む す び

Point Distribution Model (PDM)

の三次元画像 に対する非剛体位置合せに用いる

Statistical Shape Model (SSM)

の構築法について考察した．従来広く 用いられる

SSM

には，主成分分析により構築される 部分空間モデルと，近接点間の関係のみを表現する

GMRF

モデルがある．これらモデルの違いについて は，グラフィカルモデルの構造に着目することにより 統一的に議論できる．グラフィカルモデルは，条件付 きに独立ではない

2

辺数のみを辺で結ぶことにより得 られるモデルである．本論文では，従来の部分空間モ デルは共分散行列の過学習の影響により正確ではない 可能性があること，また

GMRF

モデルでは点間の条 件付き独立性を学習データに基づき推定してないこと を指摘した．

正規分布のグラフィカルモデル表現を構築する際に は，精度行列を推定すればよく，過学習を回避しつつ 精度行列を推定する手法として

graphical lasso

を採 用した．

graphical lasso

は精度行列の疎性を正則化項 として導入して推定する手法である．本論文では，過 学習を抑制しつつモデルを構築することがモデルの

generalization

の性能を向上させることにつながるこ とを指摘した．また，従来の主成分分析に基づく部分 空間モデルの構築法は，

specification

を向上させるこ

とを指向したものであり，

graphical lasso

を用いるモ デル構築法と考え方が逆であることを指摘した．その 上で，

graphical lasso

における正則化係数

ρ

^の変化に 対する部分空間の次数の変化，並びに肝臓の位置合せ 精度の変化を実験により調査した．その結果，正則化 の影響を強くするほどグラフの辺の本数は減少し，部 分空間の次数は高くなることを確認した．また，臓器 の位置合せ精度は，完全グラフ並びに

GMRF

モデル と比較して，

graphical lasso

を用いて構築したモデル の位置合せ精度が最良であることも確認した．

Graphical lasso

は精度行列が疎行列であることを 正則化項として導入し，精度行列を推定する手法で あった．臓器表面の各点が，その他全点と条件付き従 属であるとは考えにくい．例えば，

GMRF

が暗に仮定 しているように，臓器上の各点が，その近傍の点との み条件付き従属関係にあるのなら，精度行列は疎行列 となりそうである．ただし，統計モデルを構築する際，

体型などに基づき学習用の画像データが正規化される ことには注意を要する．

PDM

の統計モデルは，学習 用の臓器データ上に配置された対応点の座標に基づき 構築される．これら対応点の座標は，体型正規化に伴 い座標変換がなされている．このため，例えば，臓器 重心をそろえた上で拡大縮小するような操作が体型正 規化に含まれていれば，重心を挟んで対称に位置する 点には負の強い相関が生じ，条件付き従属となる可能 性がある．このことが，図

3

の結果が

ρ

を小さな値の ときに遠い点との間に条件付き従属関係がある理由で あろうと筆者らは考えている．グラフィカルモデルの 構造が，例えば臓器形状の決定の背後にある物理的な 要因，例えば生体膜の張力など，のみで決定されるわ けではなく，統計をとる際の正規化にも影響を受ける ということをここで指摘しておく．

グラフィカルモデルが完全グラフである場合，本来 は条件付きに独立な

2

点が，位置合せに際して互いに 影響を及ぼし合うことになる．また，測地距離のみに 基づいてモデルを構築する場合，測地距離の遠い

2

点 間の関係が一切表現されず，大局的な形状の表現能力 が低下すると予想できる．統計モデルを用いて推定を 行う場合，その推定精度はモデルが正確であるほど高 くなる．特に医用画像処理においては学習データを多 数そろえることが容易ではないため，先見知識を導入 するなどして，より正確なモデルを構築することは，

限られた学習データを有効に利用する観点からも有意 義である．