“Hfi¡.book

Transcription

1 MEDICAL IMAGING TECHNOLOGY Vol.23 No.1 January 特集論文 / 3D / 4D イメージングの最前線 逐次近似法を用いた CT 画像再構成法の考え方と驚異 Iterative Methods for Tomographic Image Reconstruction: Foundations and Surprizing Examples 工藤博幸 * Hiroyuki KUDO 要旨 CT の画像再構成法には解析的手法と代数的手法があり, 近年 SPECT や PET では代数的手法の実用化が進んでいる. 代数的手法は反復法を用いて再構成画像を求めるのが一般的であることから, 逐次近似法と呼ばれることが多い. 本論文では, 逐次近似法を用いた CT 画像再構成の考え方を平易に解説し, これにより優れた画像再構成が可能となる実例を紹介する. キーワード : 画像再構成, 逐次近似法,CT,SPECT,PET, アンジオグラフィー Two different strategies (analytical methods and algebraic methods) exist to perform tomographic image reconstruction. The algebraic methods are often called iterative methods because they compute reconstructed images by using iterative methods. This paper explains mathematical principles of iterative reconstruction methods, and show surprizing examples which demonstrate that the iterative methods can provide significantly better images compared to the analytical methods. Key words: Image reconstruction, Iterative method, CT, SPECT, PET, Angiography Med Imag Tech 23(1): 23-29, はじめに CT の画像再構成法としては, フィルタ補正逆投影 (FBP) 法に代表される解析的手法と代数的手法が知られている 1. 代数的手法は一般に反復法を用いて再構成画像を求めることから, 逐次近似法と呼ぶのが慣例となっている. 解析的手法と逐次近似法の優劣については古くから多様な意見があるが, 近年逐次近似法を上手に使うと解析的手法では困難と思われる優れた画像再構成が可能になることがわかってきた 1 ~ 15. とくに,SPECT や PET の核医学イメージングでは OS-EM(Ordered Subset Expectation Maximization) 法と呼ばれる逐次近似法の実用化 * 筑波大学大学院システム情報工学研究科コンピュータサイエンス専攻 つくば市天王台 :Department of Computer Science, Graduate School of Systems and Information Engineering, University of Tsukuba. kudo@cs.tsukuba.ac.jp 論文受付 :2004 年 8 月 24 日最終稿受付 :2005 年 1 月 19 日 が進んでいる 3 ~ 8. また, 実用的かどうかは別問題として, 逐次近似法による画像再構成法を導く過程には数理計画 数値解析 統計学などの観点から数学的な奥深いおもしろさがある. 本論文では, 逐次近似法を用いた CT 画像再構成法の考え方を平易に解説し, 著者らの画像再構成に関する研究の中から逐次近似法により優れた再構成画像が得られる実例を紹介する. 2. 不足した投影データからの画像再構成 1) 不足した投影データからの画像再構成不足した投影データから十分な画質のCT 画像を再構成することができれば,CT の応用範囲の拡大や新しい測定方式のイメージング装置の実現などが期待できる. たとえば, 不足した投影データからの画像再構成問題には,(1) 少数方向投影データからの画像再構成,(2) 投影角度範囲に制限がある画像再構成,(3) 投影データにトランケーションがある場合の画像再構成, などがある. 本章では, 不足した投影データか

2 24 Med Imag Tech Vol.23 No.1 January 2005 らの画像再構成を逐次近似法により解く一般的方法について解説する. 画像の画素値を一列に並べた画像ベクトルを x = ( x 1,, x N ) T, 測定によって得られた投影データを一列に並べたベクトルをy = ( y 1,, y M ) T とする. ただし,N は画像の画素数,M は測定された投影データの標本点数である. また, x とy を関係付ける投影演算を M N 係数行列 A={a ij } で表す. このとき, 画像 再構成は線形方程式 Ax = yを解いてy からx を求める問題として定式化できる. 不足した投影データからの画像再構成とは, 投影データ数 M が不足して画素数 N より少 ない (M < N) 投影データy に雑音が入っていないという設定の問題である. この場合, 式の数 M が未知数の数 N より少ないためAx = yを満足する画像ベクトルx は無数に存在し, 測定データy の情報のみで十分な画質の画像再構成を行うことは困難である. そこで, 対象とする被写体に関する事前情報から再構成画像のよさを評価する凸関数 f( x) を構成しておく. ただし, 評価関数 f( x) は再構成画像がよいほど値が小さくなるように設計する. そして, 画像再構成を次の線形等式制約条件付きの最適化問題として定式化する. (1) 式 (1) の定性的意味は, Ax = yを満たす測定デー タに矛盾がないx の中から評価関数 f( x) を最小とするものを選び出す と考えることができる. 式 (1) の形に定式化された画像再構成法として, 歴史的に有名なものに加法的 ART(Additive Algebraic Reconstruction Technique) と乗法的 ART (Multiplicative ART) の 2 つの逐次近似法がある 1, 2. 加法的 ART と乗法的 ART の具体的な反復式は, 次のように表される. (2) ただし,k は反復回数,i は内部反復回数である. 式 (2) の反復式の特徴は 1 回の解の更新に投影データの 1 つの標本値 y i のみを使う点であり, この構造を持つ反復法を Row Action 型と呼ぶ 1, 2.Row Action 型の反復法は, 実装が容易で投影行列 A が疎行列の場合メモリが少なくてすみ画像再構成で好んで使われる. 加法的 ART と乗法 的 ART の反復を初期ベクトルを ( 0) T x N ) は, x ( 0) x ( 0 ) = ( 1,, > 0 として実行したときの評価関数 f( x) (3) であることが知られている 2. すなわち, 加法的 ART の評価関数はx の L 2 ノルム, 乗法的 ART の評価関数はx の Shannon エントロピーである. 不足した投影データから, どの程度の画質の再構成画像が得られるかは, 評価関数 f( x) の選び方によって決まる. よく使われる一般的な評価関数として再構成画像 x の L 2 ノルムや滑らかさがある. しかし, これらのあまりに一般的すぎる評価関数では良好な再構成画像を得ることは困難である. 十分な画質の再構成画像を得るためには, イメージングの各状況における被写体の直接的な性質や再構成画像の利用のされ方など, 強力に解を拘束する事前情報を用いて f( x) を構成することが必要である. 有効な評価関数の構成法に関する一般論は存在せず, 個別事例にならざるを得ない. 以降では, f( x) を上手に選べば測定条件が驚くほど悪い状況でも十分な画質の再構成画像が得られることを, 著者らの最近の研究例から紹介する. 2) アンジオグラフィーにおける少数方向投影からの画像再構成 10, 11 一番目の例として, アンジオグラフィー装置で撮影した少数方向 (4 ~ 20) のコーンビーム投影データから 3 次元血管像を再構成する問題を取り上げる. アンジオグラフィーにおける血管の撮影では, 造影剤の注入前後でのサブトラクションにより血管以外の背景を消す操作が行

3 Med Imag Tech Vol.23 No.1 January われる. また, サブトラクションを行わない場合でも, 造影剤の吸収率は背景と比較して大きいので血管のみが濃く写っている. このような血管画像のもっとも顕著な特徴とは, ( 血管部に相当する ) 大きな値を持つ画素はごく少数でほとんどの画素は値が零 ( か零に近い ) ことである. この性質を造影血管画像の疎 (Sparse) な性質と呼ぶ. 著者らの研究 10, 11 では, 疎な画像を抽出するf( x) として次の関数を採用した. (4) から, 肺の 3 次元血管像を再構成した結果を示す. 3) トポロジー情報を利用した画像再構成 12 二番目の例として,PET の透過型スキャンによる吸収マップ再構成への適用例を示す.PET 装置では, 体内における γ 線の吸収を補正して定量性のある再構成画像を得るため, 弱い γ 線源を用いて透過型スキャンを行う 4. 透過型スキャンは PET 検査の多くの時間を占め γ 線源を回転させる複雑なメカが必要であり,( 可能であ 式 (4) はx の L 1 ノルムに他ならない. そして, 画像再構成を次の最適化問題として定式化した. (5) なぜ L 1 ノルムで疎な性質を評価したことになる かを,N = 2, M = 1( 未知数 2 個で式数 1 個 ) の簡単な場合について Fig. 1 を用いて説明する. この場合, 測定データを表す方程式 a 11 x 1 +a 12 x 2 =y 1 は 2 次元平面の直線となり, 疎な解はこの直線と座標軸 (x 1 軸か x 2 軸のどちらか ) が交差する点となる. 加法的 ART の評価関数 L 2 ノルムの等値線は Fig. 1(a) のように同心円になるため, 加法的 ART では原点から直線に降ろした垂線の足 x opt を解として選び出してしまいこれは疎でない解である. これに対して,L 1 ノルムの等値線は Fig. 1(b) のように同心の正方形になるため, 座 Fig. 1 L1 norm cost function can pick up a sparse solution but L2 norm cost function cannot. 標軸との交点 x opt を解として厳密に選び出すことができる. 式 (5) の解を求める Row Action 型の反復法は文献 10, 11 で示されている.Fig. 2 にシミュレーション実験の結果を示す. 角度 0, 45, 90, 135 の 4 方向からコーンビーム投影データを測定する状況を想定し, 加法的 ART と著者らの方法の比較を行った. 加法的 ART の L 2 ノルムでは疎な解を上手に抽出できず再構成誤差がかなり大きいのに対し, 著者らの方法による再構成画像はほぼ完全である.Fig. 3 に, 米国 Marquette 大学で開発された小動物用マイクロCT 装置で測定したラットの肺の 9 方向投影データ Fig. 2 Reconstructed image for the blood-vessel phantom from 4 simulated cone-beam projections taken from 0, 45, 90, 135 view angles.

4 26 Med Imag Tech Vol.23 No.1 January 2005 複数の領域から構成され各領域の濃度値は一定で既知 トポロジー( 領域数と領域間の接続関係 ) は既知 領域の境界線( 形 ) のみが未知であると仮定し, 次の評価関数 f( x) を採用した. Fig. 3 Reconstructed blood vessel image for a rat lung from 9 cone-beam projections measured with the micro CT system developed at Marquette university. (6) ただし,C, V( ) はクリーク ( 隣接画素対 ) の集合とクリークポテンシャルで, この部分で領域間の境界線の滑らかさを評価している.Fig. 4 に, PET 透過型スキャン実データへの適用例を示す. 解の計算は極小解を避けるため確率的緩和法 (Stochastic Relaxation) を拡張した方法により行った 12.4 方向の投影データから画像再構成を行い, 透過型スキャンの画像再構成で一般的に使われている OS-CONVEX 法, 著者らの方法, 著者らの方法からトポロジー条件を外したものの 3 つの手法を比較した. 3. 統計雑音の抑制を目的とした画像再構成 Fig. 4 Reconstructed attenuation map from 4 projections of PET transmission data. れば ) 透過型スキャンなしで計算のみで吸収補正を行いたい要望がある. 著者らの研究 12 では, 少数方向の投影データや透過型スキャンを行わない放射型スキャンの測定データのみの情報で, 吸収マップを再構成することを試みた. ここで採用した評価関数 f( x) は吸収マップのトポロジーを巧みに利用したものである. トポロジーとは画像の領域数や領域間の接続関係のことである. 吸収マップは 1) 統計雑音の抑制を目的とした画像再構成 SPECT や PET の核医学イメージングでは統計雑音の影響が大きく, その抑制を目的とした画像再構成法の開発が重要である. このような画像再構成法は 1982 年に Shepp と Vardi により初めて研究され 6,1994 年に Hudson と Larkin により OS-EM 法と呼ばれる収束の高速化を図った手法が提案され 7, 商用の SPECT や PET 装置で実用になった. このタイプの手法は, 雑音の統計的性質を利用することから統計的手法と呼ばれることも多い. 近年では, 核医学のみではなく X 線 CT への適用も検討され始めている. 本章では, 統計雑音の抑制を目的とした画像再構成法の考え方について解説する. 画像の画素値を一列に並べた画像ベクトルを x= ( x 1,, x N ) T, 測定によって得られた検出器データを一列に並べたベクトルを b= ( b 1,, b M ) T とする. ただし, N は画像の画素数,M は測定された検出器データの標本点数である. また, x と投影データを関

5 Med Imag Tech Vol.23 No.1 January 係付ける投影演算を M N 係数行列 A ={a ij } で表す. このとき, 以下のような設定の問題を考える. 測定データ数 Mが画素数 Nより多く測定データが冗長 (M > N) 測定データb に多くの統計雑音が入っている 雑音の統計モデルが既知このタイプの問題では, 雑音の統計的性質を利用していかに統計雑音を抑制するかが主要な論点であり, 最尤推定により以下のように定式化できる. 放射線を用いたイメージングにおける測定データb i は, 放射型と透過型おのおのの場合について次の統計モデルで表すことができる. (7) ただし,d i は透過型 CT において放射線源から放出される γ 線光子数の平均,Poisson(t) は平均 t のポアソン分布を表す. 式 (7) にポアソン分布の具体的な表現を代入し画像 x に条件付けた測定データb の ( 負の ) 対数尤度を計算すると, 次のようになる. (8) したがって, 最尤推定に基づく画像再構成は次の最適化問題として定式化される. (9) ただし, x 0 は再構成画像の画素値が非負の物理量であることを表す制約条件で, 非負条件と呼ばれる. また, 評価関数 f( x) は凸関数となるので, 画像再構成は非負制約条件付きの凸関数の最小化問題に帰着される. この問題の解を求める反復法の構成については膨大な数の研究がある. 好んで使われる高速に収束する反復法としては, ブロック反復法 (OS-EM 法, Row Action Maximum Likelihood Algorithm (RAMLA) 13 ), 座標降下 (Coordinate Descent) 法, 前処理付き共 役勾配法などがある 5. 以下では, 放射型 CT で実用化されている OS-EM 法について述べる. まず, 測定データ b 1,..., b M を L 個の部分集合 ( ブロック )S 1,..., S L に分割しておく. そして, 次式により反復計算を行う. (10) ただし,k は反復回数,l は内部反復回数である. この反復式の特徴は 1 回の更新に 1 つのブロックの測定データのみを使う点であり, この構造を持つ反復法をブロック反復法と呼ぶ. 収束が速く好んで用いられるブロック分割の仕方は, たとえば投影データ方向数が8でブロック数が L = 4 の場合 (11) のように, 逐次する 2 回の更新で使われる測定データが持つ情報の独立性が大きくなるものである. 実際の SPECT や PET 装置における測定データの標本化は冗長性が大きい 3 次元データ収集モードの PET でも高々 M 100N 程度であり, 測定データの冗長性だけでは十分な統計雑音の軽減が達成できない. したがって実用的には, 上述の枠組みに統計雑音と解像度の関係のトレードオフを制御する平滑化を組み合わせて使うことが必要である. たとえば,FBP 法ではランプフィルタの周波数特性がその役割を担っている. 逐次近似法に平滑化の枠組みを導入するには, いくつかの異なる方法がある. たとえば, 評価関数に画像の滑らかさを評価するコストを加える事後確率最大 (Maximum A Posteriori, MAP) 推定 8, 反復の最後に再構成画像に平滑化を施すポストスムーシング 14 が代表的なものである. 2) 画像再構成の実例国内外の学会に参加すると, 前節の定式化に基づいて構成した OS-EM 法などを用いると解析

6 28 Med Imag Tech Vol.23 No.1 January 2005 的手法と比較して良好な画像が得られることを示す強烈な実例が出てくることも, 珍しくなくなってきた. ここでは, 著者らの PET 画像再構成に関する実験例を紹介する.Fig. 5 に 5 分スキャンの FDG 検査の測定データを統計的手法と FBP 法で再構成した結果を示す. 両者とも画像の平滑化は行っていないので画像の解像度はほぼ同じであるが, 統計的手法の方が統計雑音の影響がより少ないことがわかる.Fig. 6 に 3 分スキャンの透過型データの再構成画像を示す. この場合, 統計的手法の優位性はより顕著である. 透過型スキャンにおいて FBP 法を用いると統計雑音の影響が大きいのは, 直感的には以下の理由によると考えられる.FBP 法を適用する際には対数演算 y i = ln(d i /b i ) により検出器データ b i を線積分値 y i に変換する必要があるが, 弱い γ 線を用いると分母の b i が零に近い小さな数になり対数演算の際の統計雑音の増大はきわめて大きいものとなる. 逐次近似法では実測データとの距離は対数尤度の空間で評価しており, 上述の問題は緩和されている. 3) 超高速に真の解へ収束する逐次近似法 14, 15 実用化されている手法は EM 法 6 の収束速度を改善した OS-EM 法 7 であるが, この方法では Limit Cycle と呼ばれる評価関数を最小とする真の解とは異なる偽の解に収束する問題点がある 5. この問題点を改善するため, 最近 Tanaka らは DRAMA(Dynamic RAMLA) という手法を考案し 14, 李らは非線形計画問題の双対性に基づく手法を考案した 15.Limit Cycle を回避するトリックは 2 つの方法でまったく異なるが, どちらの方法も OS-EM 法の問題点を克服し 1 ~ 2 回の反復で十分な画質の再構成画像を得ることができる. 長年逐次近似法の欠点とされてきた計算量の問題は, ほぼ克服されたと言っても過言ではない. 4. おわりに 逐次近似法を用いた画像再構成法で良好な画像が得られる理由は, 一言で述べると 非線形の画像再構成法が作りやすい ということである. 実際に, 本論文で述べた画像再構成法は加法的 ART を除いてすべて非線形である. もちろん, 解析的手法で絶対に非線形の画像再構成法が作れないわけでないが, 解析的手法はフーリエ変換 コンボリューション 微積分などの線形演算の組み合わせで逆演算を実現するという制約があるため, どうしても非線形の画像再構成法が作りにくい. Fig. 5 Reconstructed emission image from 5 min scan PET emission data. Fig. 6 Reconstructed attenuation map from 3 min scan PET transmission data.

7 Med Imag Tech Vol.23 No.1 January 文献 1 Herman GT : Image reconstruction from projections. Academic Press (Orlando), Censor Y: Finite series-expansion reconstruction methods. Proc IEEE 71: , Leahy RM, Qi J: Statistical approaches in quantitative positron emission tomography. Statistics and Computing 10: , 北村圭司 :PET におけるデータ補正と画像再構成. Med Imag Tech 19: , 工藤博幸 :SPECT / PET 画像再構成のための高速な反復法.Med Imag Tech 18: 40-45, Shepp LA, Vardi Y: Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging 1: , Hudson HM, Larkin RS: Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging 13: , Lange K, Fessler J: Globally convergent algorithms for maximum a posteriori transmission tomography. IEEE Trans Image Process 4: , Ohura N, Ogawa K, Kunieda E: Fast vascular reconstruction with MAP-EM method from few projections. Conference Record of 1999 IEEE Med Imaging Conf: , Li M, Yang H, Kudo H: An accurate iterative reconstruction algorithm for sparse objects: application to 3D blood vessel reconstruction from a limited number of projections. Phys Med Biol 47: , Li M, Kudo H, Hu J et al: Improved iterative algorithm for sparse object reconstruction and its performance evaluation with micro CT data. IEEE Trans Nucl Sci 51: , 工藤博幸, 中村宏貴 : トポロジー拘束条件付きラベリング法を用いた吸収マップ再構成. 電子情報通信学会論文誌 (D-II) J85-D-II: , Obi T, Matej S, Lewitt RM et al: 2.5-D simultaneous multislice reconstruction by series expansion methods from Fourier-rebinned PET data. IEEE Trans Med Imaging 19: , Tanaka E, Kudo H: Subset-dependent relaxation in block-iterative algorithms for image reconstruction in emission tomography. Phys Med Biol 48: , 李美花, 工藤博幸 : 非線形計画問題の双対性を用いた統計的 PET 画像再構成. 電子情報通信学会論文誌 (D-II) J87-D-II: 62-70, 2004 工藤博幸 ( くどうひろゆき ) 1985 年東北大 工 通信卒.1990 年同大大学院博士課程了. 同年同大助手. 現在, 筑波大 システム情報 CS 専攻 助教授. 工博.CT を中心とした医用イメージング, 画像処理の研究に従事.IEEE, 電子情報通信学会各会員. * * *