1
画像再構成の基礎2－逐次近似法の原理－
Basic of Image Reconstruction 2-Fundamentals of iterative method-
首都大学東京 篠原広行 Shinohara Hiroyuki はじめに 解析的画像再構成法のフィルタ補正逆投影法（FBP 法）の特徴は，1 回の計算 で解（再構成像）が求まることである．一方，反復計算によって解を求める画像 再構成法は逐次近似法と呼ばれる．逐次近似法ははじめに初期画像を仮定し，こ の画像から計算で作成した投影（順投影）と実測投影との整合性を反復計算によ って高めていく．逐次近似法は計算時間を多く必要とするが，統計雑音の性質， 装置の分解能，被写体の滑らかさなどの事前情報などを式中に組み込める柔軟 性を有し，コンピュータの高速化・高度化に伴いその重要性が増している．さら に，近年，情報科学では圧縮センシングの理論が発展し，このことも逐次近似法 が注目される要因になっている．本稿では，逐次近似法の基礎，応用例について 述べ，最後に，逐次近似法を用いた圧縮センシングによる少数投影からの画像再 構成についても紹介する． 1. 最小二乗法 逐次近似法は被写体と投影の関係を線形システムとして捉え方程式を解く方 法である1-3)．被写体断面内の物理量（原画像）f を矩形画像とし，1 辺の画素数 をN，総画素数を J (= N × N)，投影角度数 M，角度あたりの投影数 N，総投影 数をI (= M × N) ，原画像と投影の関係を表す係数行列を C (= J ×I)，投影を y とするときこれらの関係は以下で表される． 11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 J J J J I I IJ J I C f C f C f y C f C f C f y C f C f C f y (1) 本稿ではベクトルを f や y のようにイタリック体の太字で表している．(1) 式の 行列表示は C f y (2) となる．C は行列，f と p は列（縦）ベクトルである．Fig. 1 は (1) 式の原画像， 係数行列，投影の関係を示す．係数行列には計測過程の様々な要因を組み込むこ とができシステム行列あるいは投影行列とも呼ばれる．本稿では，Fig. 2 のよう にi 番目の X 線が画素 j を横切る長さ lijを検出確率Cij，検出確率の行列を係数
2 行列C としている．係数行列 C は次元（行列の行または列の数）が大きいこと， また，1 つの投影に沿って X 線が横切る画素は全画素の一部のため，行列の要 素にゼロを多く含む疎行列（スパース行列）となることが特徴である． 最小二乗法は次式の評価関数Q を最小にする f を求める． 2 2 ( ) ( ) (T ) Q f C f y C f y C f y (3) ここで 2 2はベクトルの L2ノルムの 2 乗を表し，ベクトルの成分の 2 乗和であ る（下付の2 が L2ノルムの2 を上付きの 2 が 2 乗和の 2 を表す）．T は行列の行 と列を入れ換える転置を表す．評価関数Q をベクトル f の成分 ( , , , )f f1 2 fn に ついて偏微分し0 と置くと次式が得られる． ( ) 2 T( ) 0 Q C C f f y f (4) これから原画像の近似値は次式で与えられる． 1 (C CT ) CT f y (5) この解き方の 1 つに Fig. 3 の特異値分解を利用する方法がある．特異値分解 （singular value decomposition: SVD）によって C を 3 つの行列に分解する4)．そ
れらの行列から擬似逆行列が得られる．C の分解は T C U W V (6) で表され，行列U と V は正規直交行列，行列 W は対角行列である．これから擬 似逆行列C+ は次式で与えられる． 1 T C V W U (7) Fig. 4 は 128×128 画素の画像について特異値分解による再構成像を示す 5)．比 較に，FBP 法，ML-EM 法，OS-EM 法の再構成像も載せている．雑音がない場 合，特異値分解による再構成像は分解能に優れ， FBP 法に見られる線状のアー チファクトが認められない．再構成像の視覚的な印象に一致し平均絶対誤差（画 素値の差の絶対値／総画素数）も他の方法に比べ小さい．このように擬似逆行列 を求めることで，最小二乗法に基づく画像再構成を特異値分解による行列演算 で実行することができる． 一方，反復計算によって画像と投影の整合性を一致させるには，評価関数の最 大化ではなく最小化が目的なので， (4) 式の勾配の逆方向に（誤差が減る方向 に）解を探索する．そこで 1 ( ) ( ) 2 T Q C C f g y f f (8) と置くと g は誤差の勾配ベクトルを表し反復計算は次式で表される．
3 1 , T( ) k k k k k C C k f f g g y f (9) 本節では転置行列に上付き文字を使用しているので，後述の統計的画像再構成 と異なり，反復の添字を下付文字で表している． f は反復前の画像，k fk1は反 復後の画像，kは勾配ベクトルに掛ける係数であり，この係数の違いによって 最急降下法や共役勾配法など一般に勾配法と呼ばれる最適化法がある．Fig. 5 は 反復回数に関係なく勾配ベクトルの係数をk と定数にした場合について，画 像再構成の過程を示す．(c) の f は初期画像ですべての画素の値を 1 にしている．0 (d) は (c) の画像と係数行列から計算で作成した順投影を示す． 2. 統計的逐次近似法
統計的逐次近似法のML-EM (maximum likelihood-expectation maximization) 法 は次式で表される6)． 1 kj I i ij k j I J i k ij im m i m y C C C
１
１ １ (10) ここで，検出確率Cijは j 番目の画素から発生した光子が i 番目の検出器に入射 する確率を表す．k は繰り返し回数，j は画像の画素番号（総画素数は J），i は投 影の検出器番号（総検出器数はI）を表す． k j とkj1はそれぞれk 回目と k+1 回 目の画像の画素値，yiは実測の投影である．この式を計算の手順に沿って分解し て考えると以下のようになる7,8)． ① k 回目の画像から k 回目の投影を計算で作成する（順投影）． 1 J k k i im m m y C
(11) ② k 回目の順投影と実測の投影との比を求める． ' / k i i i y y y (12) ③ その比を逆投影する． 1 1 1 ' I ' j I i ij i ij i y C C
(13) ④ k 回目の画像に逆投影した画像を掛けて k+1 回目の画像に更新する． 1 ' k k j j j (14) この更新の様子をFig. 6 に示す．画像の初期値（入力） (0) j はすべて値が1 の画 像である．初期値としては「正の値であること」という制限のもとで，一般には4 画像全体に一様分布を仮定するか，画像の内接円内のみに一様分布を仮定する． 以上の計算手順を繰り返すことによって，画像 λ はその投影が実測の投影に近 づいていく．Shepp-Logan ファントムの投影から ML-EM 法で再構成した画像を Fig. 7 の 1 行に示す．繰り返しの回数が 1 回，10 回，50 回の画像を並べて表示 している．繰り返しの回数を重ねるごとに原画像に近づいていく様子が見られ る． 統計的画像再構成法は核医学において広く臨床応用されている．係数行列に 計測過程の要因を組み込むことで，その影響を補正し画像再構成することがで きる．これがML-EM 法の柔軟性を高くしている． SPECT，PET の画像再構成 に ① 画素 j から出た光子が検出器 i に到達するまでの減弱の割合（吸収補正） 9,10) ② コリメータによって画素 j から出た光子が距離に依存して広がりを持つ 割合（深さに依存する分解能補正）9,11) ③ 放射性薬剤の体内残存分布による散乱成分の組み込み（散乱補正）12) を考慮した係数行列が使われている．
OS-EM (ordered subset EM) 法は投影をいくつかの組（サブセット）に分割し ておき，このサブセットに属する投影だけで，順投影，逆投影，比較，更新を行 い，それをサブセットごとに繰り返す方法である 13)．すべてのサブセットで画 像の更新を行った時点で，1 回の繰り返しとしている．Fig. 8 に投影角度数 16（両 矢印の線は 2 つの投影としている）におけるサブセットの分割例を示す．サブ セットを1 としたときが ML-EM 法に相当する．また，サブセットの数を投影角 度数に等しくしたときは，加減と乗除の違いはあるが，後述の加算型代数的方法 （AART 法）の考え方と同じになる．通常サブセットは 8 や 16 などが使われる． OS-EM 法では，サブセットに分けることによって 1 回の繰り返しで画像を更新 する回数が多くなり，結果として速く収束する．画像の更新回数=（サブセット 数）×（繰り返し回数）の関係が成り立ち，一般にこの更新回数が同じであれば， ほぼ同様な再構成像が得られる．サブセットの数や使用する順序は，なるべく離 れた角度の投影ごとにサブセットを構成するのがよいといわれている．サブセ ットを8 としたときの 8 回の更新の様子を Fig. 7 の 2 行に示す．サブセット 8 の 画像がOS-EM 法での 1 回の繰り返し画像となる． 3. ベイズ型逐次近似法 画像の滑らかさなどの事前情報を統計的逐次近似方法に組み込んだものを MAP (maximum a posteriori) 法あるいはベイズ型画像再構成法という．事前情報 の確率（事前確率）をPp(λ) とし，画像 λ が得られたときに投影 y が得られる条
5 件付き確率をP(y | λ)とすると，投影 y を得たときの λ の事後確率 P(λ | y) はベイ ズの定理から次式で表される． ( | ) ( ) ( | ) ( | ) ( ) ( ) p p P y P P y P y P P y (15) MAP 法は近傍の画素間で画素値の変化が少ないなど，画像の事前確率 Pp(λ) を (15) 式に与え，そのときの P(λ | y) を最大にするような λ を求める．これを次式 で表す．
ˆ arg max ln ( | ) ( ) P y Pp arg max ln ( | ) ln P y Pp( )
(16) 透過型CT の最尤推定による画像再構成として，次式の Ordered-Convex 法（OSC 法）が報告されている14)． ( , ) ( , ) ( , ) ( , ) ( , 1) ( , ) , ( , ) ( , ) ( , ) k q i k q i k q R l j ij i i k q k q i j j OSC k q R l ij i i i l B e T l R l B e
(17) ここで，jは画素j の線減弱係数，l は i 番目の X 線が画素 j を通過する長さ， ij ( , ) i i ij j j I R l l
(18) は画像から検出器i への順投影を表す．Biは被写体が存在しないときの入射X 線 強度，Tiは被写体が存在するときの透過強度である．外部線源として68Ge（148 MBq）を装着した PET 装置を用い透過型 CT 画像再構成の実験を行った例を紹 介する．透過型CT 画像再構成はギブス分布に基づく事前確率にメディアンフィ ルタを用いた MRP-OSC 法の他，MSRP-OSC 法 15)で行った． MSRP-OSC 法はMRP-OSC 法におけるメディアンフィルタ 2 ( ) ( ) exp 2 j j med med med j j m P a m (19) を事前確率としたものと，透過型CT 像をセグメント化し解剖学的情報としての 線減弱係数分布を事前確率 2 ( ) ( ) exp 2 seg j j seg seg j j S P a S (20) に用いる．これらの2 つの事前確率を互いに独立と仮定することで ( )P u Pmed( ) Pseg( ) (21) MSRP-OSC 法は次式で表される．
6 ( , 1) ( , 1) , ( , ) ( , ) , 1 1 k q k q j MSRP OSC k q k q j OSC j j j j med seg j j M S M S (22) MSRP-OSC 法は 2 つの事前確率を用いたベイズ型画像再構成法である． Table 1 はPET による放射濃度の測定値を示す．Fig. 9(a) に示す胸部ファントムの肝臓 部位に28.0 kBq/ml の18F- FDG を封入し長時間撮像を行い，統計雑音の少ない透
過型CT 像による線減弱係数分布を作成した．この分布を用い吸収補正した PET 再構成像から求めた放射濃度は 27.90 kBq/ml であり，この値を比較の基準値と した（MAC 法）．次に，計数値が長時間撮像の 1/250 に相当する透過データを用 い，透過型CT 像（(b) – (e)）を作成し PET の吸収補正を行った．MSE は (a)と の平均二乗誤差を表す． SAC 法は透過型 CT 像をセグメント化し，肺野，軟部 組織，骨などの領域ごとに一定の線減弱係数を与え吸収補正する方法である． MSRP-OSC 法による放射能度は，27.42 ± 0.96 kBq/ml と基準値に近くかつ MRP-OSC 法に比べ標準偏差も小さい．以上のことから PET によって高い正確さで放 射濃度の測定を行えることが明らかになった． 4. 3 次元 OS-EM 法 による分解能補正 Fig. 10 のように SPECT はコリメータの幾何学的な構造によって，検出器から の距離が遠くなる程分解能が低下する．実測の投影はこのような分解能の影響 を受けているので，係数行列にこの影響を含めれば画像から計算する投影と実 測の投影との整合性を高めることができる． 3 次元被写体に対する 2 次元検出 器の応答関数として次式の2 次元ガウス関数を仮定し，分解能の影響を OS-EM 法の式中に組み込める． 2 2 ( / 2 ( )) 2 1 ( , , ) 2 ( ) r d h x y d e d (23) 2 2 ( ) ( ) r x tx y ty (24) ここで，(tx, ty) はガンマ線が検出器に入射する位置，(x, y) は検出器上の座標， r は距離を表す．検出器から被写体のボクセルまでの距離 d に比例してガウス関 数の半値幅が直線的に変化するとし，FWHM = a × d + c，( ) 0.425 FWHMd と 置き，分解能の影響を係数行列に含めた3 次元 OS-EM 法が報告されている11)． 1 画素の長さ 0.31 cm の 128×128×128 画素の立方体の中央スライスに，直径 5 画素の球を中心から0，±5.0，±10.0，±15.0 cm の位置に置いた数値ファント ムをFBP 法と分解能補正付き 3 次元 OS-EM（OS-EM with DRC）で画像再構成 した．Fig. 11 に示すように FBP 法の再構成像は中心では円となるが，中心から 離れるに従い楕円になり分解能の影響が顕著に現れている．一方，分解能補正付
7 き 3 次元 OS-EM 法の再構成像は位置に関係なくほぼ円に近い．プロファイル においてもFBP 法の半値幅は位置に依存して変化するのに対し，分解能補正付 き3 次元 OS-EM 法の半値幅は変化が僅かである．Fig. 12 は 3 方向について中心 からの距離と半値幅の関係を示す． 5．圧縮センシングによる少数投影からの画像再構成 原画像は256×256 画素で 1 投影角度あたりのデータ数（投影数）を 256 とす ると，画像再構成には標本化定理から 402 の投影角度数を必要とする．圧縮セ ンシングによって原画像が疎の性質を持つ画像（非ゼロの画素数が少ない画像 のことで英語では sparsified image と呼ばれる．本稿ではスパース画像というこ とにする．）に変換される場合には，少数投影から画像再構成を行える16)．これ まで256×256 画素の Shepp-Logan ファントムの画像再構成には投影角度数 256 を用いたが，本節では僅か投影角度数 16 から画像再構成を行えることを示す． この魔法のような話しにはL1ノルムが重要な役割を果たす．画像再構成に画像 の勾配を正則化に用いた最小二乗法の評価関数は次式で表される． 2 2 1 ( ) Q f C f y f (25) 右辺の第1 項はこれまでの式と同じであり，第 2 項の f は画像の勾配を示し， 1 f は勾配についての L1ノルム（勾配の大きさの総和）を示す． はこの L1 ノルムに掛ける重み係数である．トータルバリエーション（全変動）ノルム（TV ノルム）とは勾配の大きさのL1ノルムとして次式で表される． 2 2 1 T V f x y dxdy( , ) f f dxdy x y
(26) したがって，TV ノルムの最小化は勾配画像（勾配の大きさを表す画像）の L1ノ ルムの最小化と同じことである． (26) 式の TV ノルムを(25) 式の第 2 項に用い ると評価関数Q は次式で表される． 2 2 1 ( ) ( ) Q f C f y T V f (27) (27) 式の評価関数 Q を小さくするように f を求めるには，画像と投影との整 合性の指標である第 1 項の誤差を小さくするとともに，第 2 項も小さくする必 要がある．TV ノルムの正則化は再構成像における勾配の変動の総和を最小にす る働きがあり，強度が一定であるべき領域を滑らかにする一方でエッジを保存 し，雑音やギブスリンギングアーチファクトに対しては抑制する．その結果，少 ない投影角度数でもエッジを保存しアーチファクトの少ない再構成像を得るこ とができる．そのためには を適切に選ぶ必要があり， を必要以上に大きくす8 ることは画像の微小な強度変化など本来の構造を失わせる結果となる．代数的 方法のSIRT 法 1,17)は，順投影と実測投影との差をすべての投影角度について平 均し画像の更新に用いる．一方，AART 法は投影角度ごとに差を反復画像に反映 させる．実測投影との比較に差を用いるAART 法は次式で表される． 1 1 1 J k i im m n N k k m j j J ij i n im m y C f f f C C
(28) ここで，N は 1 投影角度あたりのデータ数，J は画素の総数，m(1 m J)は画 像の画素を表す添字，i は投影を表す添字，n は比較を行うサイノグラムの縦方 向（角度方向）の座標を示す．Fig. 13 に 180 度について投影角度数 16 における FBP 法（(a)-(c)）と AART 法（(d)-(f)）の再構成像を示す．投影にはガウス雑音 を仮定し，その大きさを20 dB，30 dB，40 dB と変えている．20 dB は投影の平 均値の約1%，30 dB は約 0.1%， 40 dB は約 0.01%の雑音レベルである． Fig. 14 は(a) Shepp-Logan ファントム，(b) 勾配画像を示す．このファントムはそれぞれ の楕円の値が一定なのでそれらの領域内の勾配はゼロとなる．すなわち勾配画 像にはゼロの領域が多い．その結果，(27) 式の TV ノルムは小さな値になる．(c) -(e) は 0.1のTV 正則化を用いた AART 法（AART_TV 法と略）の再構成像を 示す18)．なお，AART 法と AART_TV 法では，画像の内接円の外はゼロ，および 被写体は負の値をもたないとする正則化を付けている．20 dB の雑音では FBP 法 の再構成像は内部構造が不明である．AART 法の再構成像は構造をそれなりに 再現しているが，アーチファクトが顕著である．一方，圧縮センシングによる AART_TV 法の再構成像はファントムに近い画像となっている． 圧縮センシングで重要なことは，原画像が図 14(b) のようにスパース画像に 変換可能なことである．スパース画像は原画像よりもゼロの画素を多く含み，値 として意味のある画素が少なくなる．このような性質を持つスパース画像を利 用することで，標本化定理よりも少ない投影角度数から画像再構成を行うこと ができる．医用画像は値がほぼ一様な複数のセグメントからなるとみなせるの で圧縮センシングの対象になる．本稿ではスパース画像への変換にTV を使用し たが，他にウエーブレット変換 19)による方法がある．圧縮センシングによる画 像再構成は，X 線 CT では被ばくの低減に，MRI では検査時間の短縮に繋がる可 能性があり今後の研究が注目される．9 Figure Captions Fig. 1 線形連立方程式で表される原画像，係数行列，投影の関係 Fig. 2 画像と投影の関係を表す係数行列 Fig. 3 係数行列 C の特異値分解と擬似逆行列 C + Fig. 4 特異値分解による再構成像 Fig. 5 反復計算による最小二乗法の解法 Fig. 6 ML-EM 法の計算過程
Fig. 7 ML-EM 法と OS-EM 法の再構成像
Fig. 8 投影角度数 16 におけるサブセットの分割例
Fig. 9 MSRP-OSC 法による胸部ファントムの線減弱係数分布（Table 1，Fig. 9 は Sakaguchi K, et al.: Ann Nucl Med 22: 269-279, 200815) からの転載許可を日本核医
学会から得た．）
Fig. 10 2 次元ガウス関数で近似した検出器応答関数 Fig. 11 分解能補正付き 3D OS-EM 法の再構成像
Fig. 12 分解能補正付き 3D OS-EM 法の半値幅（Fig. 11，12 は Yokoi T, et al.: Ann Nucl Med 16: 11-189, 200211) からの転載許可を日本核医学会から得た．）
Fig. 13 少数投影からの FBP 法と代数的方法の再構成像 Fig. 14 圧縮センシングによる少数投影からの再構成像
10 参考文献
1) Kak AC, Slaney M: Principles of Computerized Tomographic Imaging. IEEE press, New York, 1988: 275-293.
2) Gordon R, Bender R, Herman GT: Algebraic reconstruction technique (ART) for three dimensional electron microscopy and X-ray tomography. J Theor Biol 29: 470-481, 1971 3) 篠原広行，中世古和真，坂口和也，他：逐次近似画像再構成の基礎．医療科 学社，東京， 2013: 68-94. 4) 河田 聡，南 茂雄 編：科学計測のための画像データ処理． CQ 出版，東京， 1994: 121-126. 5) 寺松翔太，篠原広行，橋本雄幸：特異値分解による CT 画像再構成の検討. 日 本放射線技術学会雑誌67: 1073, 2011
6) Lange K, Carson R: EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomgr 8: 306-316, 1984
7) 橋本雄幸，篠原広行：C 言語による画像再構成の基礎．医療科学社，東京， 2006: 257-261. 8) 篠原広行，坂口和也，橋本雄幸：Excel による画像再構成入門．医療科学社， 東京， 2007: 95-123. 9) 橋本雄幸，横井孝司，篠原広行：SPECT 画像再構成の基礎．医療科学社，東 京， 2006: 122-126, 175-182.
10) 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 24: 179-188, 2010
11) 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 16: 11-18, 2002 12) 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: 1125-1133, 2003
13) Hudson HM, Larkin RS: Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging 13: 601-609, 1994
14) Lange K, Fessler JA: Globally convergent algorithms for maximum a posteriori transmission tomography. IEEE Trans Med Imaging 4: 1430-14,1995
15) Sakaguchi K, Shinohara H, Hashimoto T, et al.: An iterative reconstruction using median root prior and anatomical prior from segmented μ map for count-limited
11
transmission data in PET imaging. Ann Nucl Med 22: 269-279, 2008
16) 工藤博幸，イサム ラシド：圧縮センシングを用いた少数投影データからの CT 画像再構成．映像情報 Medical 43: 1093-1099, 2011 17) 橘 篤志，橋本雄幸，坂口和也，他：k 空間データの極座標変換を用いた MRI 画像再構成法．日本放射線技術学会雑誌 68: 413-422, 2012 18) 篠原広行，小畠隆行，橋本雄幸：断層映像法の基礎 第 41 回 圧縮センシ ングによる少数投影からの画像再構成．断層映像研究会誌，印刷中． 19) 篠原広行，伊藤 猛，橋本雄幸：医用画像位置合わせの基礎．医療科学社， 東京，2011: 259-281.
