<4D F736F F D F385F322089E6919C8DC48D5C90AC82CC8AEE B92808E9F8BDF8E CC8CB4979D815B>

Transcription

1 画像再構成の基礎 - 逐次近似法の原理 - Basic of Image Reconstruction -Fundamentals of iterative method- 首都大学東京篠原広行 Shinohara Hiroyuki はじめに解析的画像再構成法のフィルタ補正逆投影法 (FBP 法 ) の特徴は, 回の計算で解 ( 再構成像 ) が求まることである. 一方, 反復計算によって解を求める画像再構成法は逐次近似法と呼ばれる. 逐次近似法ははじめに初期画像を仮定し, この画像から計算で作成した投影 ( 順投影 ) と実測投影との整合性を反復計算によって高めていく. 逐次近似法は計算時間を多く必要とするが, 統計雑音の性質, 装置の分解能, 被写体の滑らかさなどの事前情報などを式中に組み込める柔軟性を有し, コンピュータの高速化 高度化に伴いその重要性が増している. さらに, 近年, 情報科学では圧縮センシングの理論が発展し, このことも逐次近似法が注目される要因になっている. 本稿では, 逐次近似法の基礎, 応用例について述べ, 最後に, 逐次近似法を用いた圧縮センシングによる少数投影からの画像再構成についても紹介する.. 最小二乗法逐次近似法は被写体と投影の関係を線形システムとして捉え方程式を解く方法である -3). 被写体断面内の物理量 ( 原画像 )f を矩形画像とし, 辺の画素数を N, 総画素数を J (= N N), 投影角度数 M, 角度あたりの投影数 N, 総投影数を I (= M N), 原画像と投影の関係を表す係数行列を C (= J I), 投影を y とするときこれらの関係は以下で表される. C f C f C f y J J C f C f C f y J J () C f C f C f y I I IJ J I 本稿ではベクトルを f や y のようにイタリック体の太字で表している.() 式の行列表示は C f y () となる.C は行列,f と p は列 ( 縦 ) ベクトルである.Fig. は () 式の原画像, 係数行列, 投影の関係を示す. 係数行列には計測過程の様々な要因を組み込むことができシステム行列あるいは投影行列とも呼ばれる. 本稿では,Fig. のように i 番目の X 線が画素 j を横切る長さ lij を検出確率 Cij, 検出確率の行列を係数

2 行列 C としている. 係数行列 C は次元 ( 行列の行または列の数 ) が大きいこと, また, つの投影に沿って X 線が横切る画素は全画素の一部のため, 行列の要素にゼロを多く含む疎行列 ( スパース行列 ) となることが特徴である. 最小二乗法は次式の評価関数 Q を最小にする f を求める. Q( f ) C f y ( C f y) T ( C f y ) (3) ここで はベクトルの L ノルムの 乗を表し, ベクトルの成分の 乗和であ る ( 下付の が L ノルムの を上付きの が 乗和の を表す ).T は行列の行と列を入れ換える転置を表す. 評価関数 Q をベクトル f の成分 ( f, f,, f n ) について偏微分し 0 と置くと次式が得られる. Q( f ) T C ( C f y) 0 (4) f これから原画像の近似値は次式で与えられる. ( T T f CC) C y (5) この解き方の つに Fig. 3 の特異値分解を利用する方法がある. 特異値分解 (singular value decomposition: SVD) によって C を 3 つの行列に分解する 4). それらの行列から擬似逆行列が得られる.C の分解は C T UWV (6) で表され, 行列 U と V は正規直交行列, 行列 W は対角行列である. これから擬似逆行列 C + は次式で与えられる. C V W U T (7) Fig. 4 は 8 8 画素の画像について特異値分解による再構成像を示す 5). 比較に,FBP 法,ML-EM 法,OS-EM 法の再構成像も載せている. 雑音がない場合, 特異値分解による再構成像は分解能に優れ, FBP 法に見られる線状のアーチファクトが認められない. 再構成像の視覚的な印象に一致し平均絶対誤差 ( 画素値の差の絶対値 / 総画素数 ) も他の方法に比べ小さい. このように擬似逆行列を求めることで, 最小二乗法に基づく画像再構成を特異値分解による行列演算で実行することができる. 一方, 反復計算によって画像と投影の整合性を一致させるには, 評価関数の最大化ではなく最小化が目的なので, (4) 式の勾配の逆方向に ( 誤差が減る方向に ) 解を探索する. そこで Q( f ) T g C ( yc f) (8) f と置くと g は誤差の勾配ベクトルを表し反復計算は次式で表される.

3 T fk fk kgk, gk C ( ycf k) (9) 本節では転置行列に上付き文字を使用しているので, 後述の統計的画像再構成と異なり, 反復の添字を下付文字で表している. f k は反復前の画像, f k は反復後の画像, k は勾配ベクトルに掛ける係数であり, この係数の違いによって最急降下法や共役勾配法など一般に勾配法と呼ばれる最適化法がある.Fig. 5 は反復回数に関係なく勾配ベクトルの係数を k と定数にした場合について, 画像再構成の過程を示す.(c) の f 0 は初期画像ですべての画素の値を にしている. (d) は (c) の画像と係数行列から計算で作成した順投影を示す.. 統計的逐次近似法統計的逐次近似法の ML-EM (maximum likelihood-expectation maximization) 法は次式で表される 6). yc k I k j i ij j I J i k C ij Cimm i m ここで, 検出確率 Cij は j 番目の画素から発生した光子が i 番目の検出器に入射する確率を表す.k は繰り返し回数,j は画像の画素番号 ( 総画素数は J),i は投影の検出器番号 ( 総検出器数は I) を表す. k j と k j はそれぞれ k 回目と k+ 回目の画像の画素値,yi は実測の投影である. この式を計算の手順に沿って分解して考えると以下のようになる 7,8). k 回目の画像から k 回目の投影を計算で作成する ( 順投影 ). y J k k i Cimm m (0) () k 回目の順投影と実測の投影との比を求める. k y ' y / y () i i i 3 その比を逆投影する. ' y ' C j I i ij i Cij i I 4 k 回目の画像に逆投影した画像を掛けて k+ 回目の画像に更新する. k k j j j (3) ' (4) この更新の様子を Fig. 6 に示す. 画像の初期値 ( 入力 ) (0) j はすべて値が の画 像である. 初期値としては 正の値であること という制限のもとで, 一般には 3

4 画像全体に一様分布を仮定するか, 画像の内接円内のみに一様分布を仮定する. 以上の計算手順を繰り返すことによって, 画像 λ はその投影が実測の投影に近づいていく.Shepp-Logan ファントムの投影から ML-EM 法で再構成した画像を Fig. 7 の 行に示す. 繰り返しの回数が 回,0 回,50 回の画像を並べて表示している. 繰り返しの回数を重ねるごとに原画像に近づいていく様子が見られる. 統計的画像再構成法は核医学において広く臨床応用されている. 係数行列に計測過程の要因を組み込むことで, その影響を補正し画像再構成することができる. これが ML-EM 法の柔軟性を高くしている. SPECT,PET の画像再構成に 画素 j から出た光子が検出器 i に到達するまでの減弱の割合 ( 吸収補正 ) 9,0) コリメータによって画素 j から出た光子が距離に依存して広がりを持つ割合 ( 深さに依存する分解能補正 ) 9,) 3 放射性薬剤の体内残存分布による散乱成分の組み込み ( 散乱補正 ) ) を考慮した係数行列が使われている. OS-EM (ordered subset EM) 法は投影をいくつかの組 ( サブセット ) に分割しておき, このサブセットに属する投影だけで, 順投影, 逆投影, 比較, 更新を行い, それをサブセットごとに繰り返す方法である 3). すべてのサブセットで画像の更新を行った時点で, 回の繰り返しとしている.Fig. 8 に投影角度数 6( 両矢印の線は つの投影としている ) におけるサブセットの分割例を示す. サブセットを としたときが ML-EM 法に相当する. また, サブセットの数を投影角度数に等しくしたときは, 加減と乗除の違いはあるが, 後述の加算型代数的方法 (AART 法 ) の考え方と同じになる. 通常サブセットは 8 や 6 などが使われる. OS-EM 法では, サブセットに分けることによって 回の繰り返しで画像を更新する回数が多くなり, 結果として速く収束する. 画像の更新回数 =( サブセット数 ) ( 繰り返し回数 ) の関係が成り立ち, 一般にこの更新回数が同じであれば, ほぼ同様な再構成像が得られる. サブセットの数や使用する順序は, なるべく離れた角度の投影ごとにサブセットを構成するのがよいといわれている. サブセットを 8 としたときの 8 回の更新の様子を Fig. 7 の 行に示す. サブセット 8 の画像が OS-EM 法での 回の繰り返し画像となる. 3. ベイズ型逐次近似法画像の滑らかさなどの事前情報を統計的逐次近似方法に組み込んだものを MAP (maximum a posteriori) 法あるいはベイズ型画像再構成法という. 事前情報の確率 ( 事前確率 ) を Pp(λ) とし, 画像 λ が得られたときに投影 y が得られる条 4

5 件付き確率を P(y λ) とすると, 投影 y を得たときの λ の事後確率 P(λ y) はベイズの定理から次式で表される. Py ( ) Pp ( ) P( y) P( y ) Pp ( ) (5) Py ( ) MAP 法は近傍の画素間で画素値の変化が少ないなど, 画像の事前確率 Pp(λ) を (5) 式に与え, そのときの P(λ y) を最大にするような λ を求める. これを次式で表す. Py Pp Py Pp ˆ arg max ln ( ) ( ) arg max ln ( ) ln ( ) (6) 透過型 CT の最尤推定による画像再構成として, 次式の Ordered-Convex 法 (OSC 法 ) が報告されている 4). ( kq, ) ( kq, ) Rl ( i, ) j lij Bie Ti ( kq, ) ( kq, ) i josc, j ( kq, ) ( kq, ) Rl ( i, ) lrl ij ( i, ) Be i i ここで, j は画素 j の線減弱係数, l ij は i 番目の X 線が画素 j を通過する長さ, (7) R( l, ) l (8) i ij j jii は画像から検出器 i への順投影を表す.Bi は被写体が存在しないときの入射 X 線強度,Ti は被写体が存在するときの透過強度である. 外部線源として 68 Ge(48 MBq) を装着した PET 装置を用い透過型 CT 画像再構成の実験を行った例を紹介する. 透過型 CT 画像再構成はギブス分布に基づく事前確率にメディアンフィルタを用いた MRP-OSC 法の他,MSRP-OSC 法 5) で行った. MSRP-OSC 法は MRP-OSC 法におけるメディアンフィルタ P med med ( j m j) ( ) amed exp j m j (9) を事前確率としたものと, 透過型 CT 像をセグメント化し解剖学的情報としての線減弱係数分布を事前確率 P seg seg ( j S j ) ( ) aseg exp j S j (0) に用いる. これらの つの事前確率を互いに独立と仮定することで Pu ( ) P ( ) P ( ) () med MSRP-OSC 法は次式で表される. seg 5

6 ( kq, ) ( kq, ) j, MSRPOSC ( kq, ) ( kq, ) j, OSC j M j j Sj med seg M j Sj MSRP-OSC 法は つの事前確率を用いたベイズ型画像再構成法である. Table は PET による放射濃度の測定値を示す.Fig. 9(a) に示す胸部ファントムの肝臓部位に 8.0 kbq/ml の 8 F- FDG を封入し長時間撮像を行い, 統計雑音の少ない透過型 CT 像による線減弱係数分布を作成した. この分布を用い吸収補正した PET 再構成像から求めた放射濃度は 7.90 kbq/ml であり, この値を比較の基準値とした (MAC 法 ). 次に, 計数値が長時間撮像の /50 に相当する透過データを用い, 透過型 CT 像 ((b) (e)) を作成し PET の吸収補正を行った.MSE は (a) との平均二乗誤差を表す. SAC 法は透過型 CT 像をセグメント化し, 肺野, 軟部組織, 骨などの領域ごとに一定の線減弱係数を与え吸収補正する方法である. MSRP-OSC 法による放射能度は,7.4 ± 0.96 kbq/ml と基準値に近くかつ MRP- OSC 法に比べ標準偏差も小さい. 以上のことから PET によって高い正確さで放射濃度の測定を行えることが明らかになった. () 4. 3 次元 OS-EM 法による分解能補正 Fig. 0 のように SPECT はコリメータの幾何学的な構造によって, 検出器からの距離が遠くなる程分解能が低下する. 実測の投影はこのような分解能の影響を受けているので, 係数行列にこの影響を含めれば画像から計算する投影と実測の投影との整合性を高めることができる. 3 次元被写体に対する 次元検出器の応答関数として次式の 次元ガウス関数を仮定し, 分解能の影響を OS-EM 法の式中に組み込める. ( r / ( d)) hxyd (,, ) e (3) ( d) r x tx y ty ( ) ( ) (4) ここで,(tx, ty) はガンマ線が検出器に入射する位置,(x, y) は検出器上の座標, r は距離を表す. 検出器から被写体のボクセルまでの距離 d に比例してガウス関数の半値幅が直線的に変化するとし,FWHM = a d + c, ( d) 0.45FWHM と置き, 分解能の影響を係数行列に含めた 3 次元 OS-EM 法が報告されている ). 画素の長さ 0.3 cm の 画素の立方体の中央スライスに, 直径 5 画素の球を中心から 0,±5.0,±0.0,±5.0 cm の位置に置いた数値ファントムを FBP 法と分解能補正付き 3 次元 OS-EM(OS-EM with DRC) で画像再構成した.Fig. に示すように FBP 法の再構成像は中心では円となるが, 中心から離れるに従い楕円になり分解能の影響が顕著に現れている. 一方, 分解能補正付 6

7 き 3 次元 OS-EM 法の再構成像は位置に関係なくほぼ円に近い. プロファイルにおいても FBP 法の半値幅は位置に依存して変化するのに対し, 分解能補正付き 3 次元 OS-EM 法の半値幅は変化が僅かである.Fig. は 3 方向について中心からの距離と半値幅の関係を示す. 5. 圧縮センシングによる少数投影からの画像再構成原画像は 画素で 投影角度あたりのデータ数 ( 投影数 ) を 56 とすると, 画像再構成には標本化定理から 40 の投影角度数を必要とする. 圧縮センシングによって原画像が疎の性質を持つ画像 ( 非ゼロの画素数が少ない画像のことで英語では sparsified image と呼ばれる. 本稿ではスパース画像ということにする.) に変換される場合には, 少数投影から画像再構成を行える 6). これまで 画素の Shepp-Logan ファントムの画像再構成には投影角度数 56 を用いたが, 本節では僅か投影角度数 6 から画像再構成を行えることを示す. この魔法のような話しには L ノルムが重要な役割を果たす. 画像再構成に画像の勾配を正則化に用いた最小二乗法の評価関数は次式で表される. Q( f ) C f y f (5) 右辺の第 項はこれまでの式と同じであり, 第 項の f は画像の勾配を示し, f は勾配についての L ノルム ( 勾配の大きさの総和 ) を示す. はこの L ノルムに掛ける重み係数である. トータルバリエーション ( 全変動 ) ノルム (TV ノルム ) とは勾配の大きさの L ノルムとして次式で表される. f f TV f ( x, y) dxdy dxdy x y (6) したがって,TV ノルムの最小化は勾配画像 ( 勾配の大きさを表す画像 ) の L ノルムの最小化と同じことである. (6) 式の TV ノルムを (5) 式の第 項に用いると評価関数 Q は次式で表される. Q( f ) C f y TV( f ) (7) (7) 式の評価関数 Q を小さくするように f を求めるには, 画像と投影との整合性の指標である第 項の誤差を小さくするとともに, 第 項も小さくする必要がある.TV ノルムの正則化は再構成像における勾配の変動の総和を最小にする働きがあり, 強度が一定であるべき領域を滑らかにする一方でエッジを保存し, 雑音やギブスリンギングアーチファクトに対しては抑制する. その結果, 少ない投影角度数でもエッジを保存しアーチファクトの少ない再構成像を得ることができる. そのためには を適切に選ぶ必要があり, を必要以上に大きくす 7

8 ることは画像の微小な強度変化など本来の構造を失わせる結果となる. 代数的方法の SIRT 法,7) は, 順投影と実測投影との差をすべての投影角度について平均し画像の更新に用いる. 一方,AART 法は投影角度ごとに差を反復画像に反映させる. 実測投影との比較に差を用いる AART 法は次式で表される. J k yi C n N imf m k k m f j f j C J ij (8) in C im m ここで,N は 投影角度あたりのデータ数,J は画素の総数, m( m J) は画像の画素を表す添字,i は投影を表す添字,n は比較を行うサイノグラムの縦方向 ( 角度方向 ) の座標を示す.Fig. 3 に 80 度について投影角度数 6 における FBP 法 ((a)-(c)) と AART 法 ((d)-(f)) の再構成像を示す. 投影にはガウス雑音を仮定し, その大きさを 0 db,30 db,40 db と変えている.0 db は投影の平均値の約 %,30 db は約 0.%, 40 db は約 0.0% の雑音レベルである. Fig. 4 は (a) Shepp-Logan ファントム,(b) 勾配画像を示す. このファントムはそれぞれの楕円の値が一定なのでそれらの領域内の勾配はゼロとなる. すなわち勾配画像にはゼロの領域が多い. その結果,(7) 式の TV ノルムは小さな値になる.(c) -(e) は 0.の TV 正則化を用いた AART 法 (AART_TV 法と略 ) の再構成像を示す 8). なお,AART 法と AART_TV 法では, 画像の内接円の外はゼロ, および被写体は負の値をもたないとする正則化を付けている.0 db の雑音では FBP 法の再構成像は内部構造が不明である.AART 法の再構成像は構造をそれなりに再現しているが, アーチファクトが顕著である. 一方, 圧縮センシングによる AART_TV 法の再構成像はファントムに近い画像となっている. 圧縮センシングで重要なことは, 原画像が図 4(b) のようにスパース画像に変換可能なことである. スパース画像は原画像よりもゼロの画素を多く含み, 値として意味のある画素が少なくなる. このような性質を持つスパース画像を利用することで, 標本化定理よりも少ない投影角度数から画像再構成を行うことができる. 医用画像は値がほぼ一様な複数のセグメントからなるとみなせるので圧縮センシングの対象になる. 本稿ではスパース画像への変換に TV を使用したが, 他にウエーブレット変換 9) による方法がある. 圧縮センシングによる画像再構成は,X 線 CT では被ばくの低減に,MRI では検査時間の短縮に繋がる可能性があり今後の研究が注目される. 8

9 Figure Captions Fig. 線形連立方程式で表される原画像, 係数行列, 投影の関係 Fig. 画像と投影の関係を表す係数行列 Fig. 3 係数行列 C の特異値分解と擬似逆行列 C + Fig. 4 特異値分解による再構成像 Fig. 5 反復計算による最小二乗法の解法 Fig. 6 ML-EM 法の計算過程 Fig. 7 ML-EM 法と OS-EM 法の再構成像 Fig. 8 投影角度数 6 におけるサブセットの分割例 Fig. 9 MSRP-OSC 法による胸部ファントムの線減弱係数分布 (Table,Fig. 9 は Sakaguchi K, et al.: Ann Nucl Med : 69-79, 008 5) からの転載許可を日本核医学会から得た.) Fig. 0 次元ガウス関数で近似した検出器応答関数 Fig. 分解能補正付き 3D OS-EM 法の再構成像 Fig. 分解能補正付き 3D OS-EM 法の半値幅 (Fig., は Yokoi T, et al.: Ann Nucl Med 6: -89, 00 ) からの転載許可を日本核医学会から得た.) Fig. 3 少数投影からの FBP 法と代数的方法の再構成像 Fig. 4 圧縮センシングによる少数投影からの再構成像 9

10 参考文献 ) Kak AC, Slaney M: Principles of Computerized Tomographic Imaging. IEEE press, New York, 988: ) Gordon R, Bender R, Herman GT: Algebraic reconstruction technique (ART) for three dimensional electron microscopy and X-ray tomography. J Theor Biol 9: , 97 3) 篠原広行, 中世古和真, 坂口和也, 他 : 逐次近似画像再構成の基礎. 医療科学社, 東京, 03: ) 河田聡, 南茂雄編 : 科学計測のための画像データ処理.CQ 出版, 東京, 994: -6. 5) 寺松翔太, 篠原広行, 橋本雄幸 : 特異値分解による CT 画像再構成の検討. 日本放射線技術学会雑誌 67: 073, 0 6) Lange K, Carson R: EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomgr 8: , 984 7) 橋本雄幸, 篠原広行 :C 言語による画像再構成の基礎. 医療科学社, 東京, 006: ) 篠原広行, 坂口和也, 橋本雄幸 :Excel による画像再構成入門. 医療科学社, 東京, 007: ) 橋本雄幸, 横井孝司, 篠原広行 :SPECT 画像再構成の基礎. 医療科学社, 東京, 006: -6, ) Daisaki H, Shinohara H, Terauchi T, et al.: Multi-bed-position acquisition technique for deep inspiration breath-hold PET/CT: a preliminary result for pulmonary lesions. Ann Nucl Med 4: 79-88, 00 ) Yokoi T, Shinohara H, Onishi H: Performance evaluation of OSEM reconstruction algorithm incorporating three-dimensional distance-dependent resolution compensation for brain SPET: A Simulation study. Ann Nucl Med 6: -8, 00 ) Yokoi T, Shinohara H, Takaki A: Improvement of signal-to-noise ratio using iterative reconstruction in a 99mTc-ECD split-dose injection protocol. Eur J Nucl Med and Molecular Imaging 30: 5-33, 003 3) Hudson HM, Larkin RS: Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging 3: , 994 4) Lange K, Fessler JA: Globally convergent algorithms for maximum a posteriori transmission tomography. IEEE Trans Med Imaging 4: 430-4,995 5) Sakaguchi K, Shinohara H, Hashimoto T, et al.: An iterative reconstruction using median root prior and anatomical prior from segmented μ map for count-limited 0

11 transmission data in PET imaging. Ann Nucl Med : 69-79, 008 6) 工藤博幸, イサムラシド : 圧縮センシングを用いた少数投影データからの CT 画像再構成. 映像情報 Medical 43: , 0 7) 橘篤志, 橋本雄幸, 坂口和也, 他 :k 空間データの極座標変換を用いた MRI 画像再構成法. 日本放射線技術学会雑誌 68: 43-4, 0 8) 篠原広行, 小畠隆行, 橋本雄幸 : 断層映像法の基礎第 4 回圧縮センシングによる少数投影からの画像再構成. 断層映像研究会誌, 印刷中. 9) 篠原広行, 伊藤猛, 橋本雄幸 : 医用画像位置合わせの基礎. 医療科学社, 東京,0: 59-8.