第8章 逐次近似再構成

再構成は、連立一次方程式を解く問題として定式化できます。画像を $N^2$ 個の画素値を並べたベクトル $x$、サイノグラムを $M$ 個（投影数 × 検出器数）の測定値ベクトル $p$ とすると、順投影は系統行列と呼ばれる巨大な疎行列 $A$（要素 $a_{ij}$ はレイ $i$ が画素 $j$ を横切る重み）による線形写像になり、再構成は $Ax = p$ を $x$ について解くことに帰着します。1 パスで解析的に逆変換する FBP に対して、こちらは方程式系を反復で解く見方です。投影数が少なければ方程式は劣決定（未知数 > 独立な式）になり、測定にノイズが乗れば厳密に整合する解は存在しません。そこで登場するのが、初期画像から出発して測定との食い違いを少しずつ減らしていく反復解法です。これを逐次近似再構成と呼びます。

$$A\mathbf{x} = \mathbf{p}$$

最も古い逐次近似法が ART（Algebraic Reconstruction Technique）で、数学的には Kaczmarz 法として知られています。方程式を 1 本（= 1 レイ）ずつ取り出し、現在の画像 $x$ をその方程式が定める超平面 $a_i \cdot x = p_i$ の上へ射影する操作を、全レイについて巡回します。更新式の分子 $p_i - a_i \cdot x$ は、そのレイの測定値と現在画像の順投影との食い違い（残差）です。これをレイの通り道に塗り戻して食い違いを打ち消します。緩和係数 $\lambda$ は 1 回の更新で超平面へどこまで近づくかを決め、小さくするほどノイズが平均化されて挙動が安定します。ノイズのない整合した方程式系なら、この射影の繰り返しが解に収束することが証明されています。

$$\mathbf{x} \;\leftarrow\; \mathbf{x} + \lambda\,\frac{p_i - \mathbf{a}_i\!\cdot\!\mathbf{x}}{\lVert \mathbf{a}_i \rVert^2}\,\mathbf{a}_i$$

ART は 1 レイごとに画像を更新するため初期収束が速い一方、個々のレイに含まれるノイズの影響を受けやすいという性質があります。SIRT（Simultaneous Iterative Reconstruction Technique）は全レイの残差をまとめて計算し、レイごとの行和 $R$ と画素ごとの列和 $C$（感度）で正規化した補正 $A^{\mathsf T}$ を一括で適用します。1 反復あたりの変化は小さいものの、多数のレイを平均するためノイズの影響を抑えられます。SART は両者の中間に位置します。1 つの投影角度（ビュー）分のレイをまとめて SIRT と同じ式で逐次更新し、収束速度と安定性のバランスを取ります。

$$\mathbf{x} \;\leftarrow\; \mathbf{x} + \lambda\, C^{-1} A^{\mathsf T} R^{-1}\,(\mathbf{p} - A\mathbf{x}), \qquad R = \operatorname{diag}\Bigl(\sum_j a_{ij}\Bigr),\quad C = \operatorname{diag}\Bigl(\sum_i a_{ij}\Bigr)$$

ここまでは幾何学的な方程式解法でしたが、測定の統計モデルから出発する方法もあります。光子計数を直接扱う放射型トモグラフィでは、測定値をポアソン分布とみなし、尤度を最大にする画像を求める MLEM（Maximum Likelihood Expectation Maximization）が使われます。測定と順投影の比 $p/(Ax)$ を逆投影し、感度で正規化して現在画像へ掛ける乗法型の更新なので、正の初期値から始めれば非負性が保たれます。OSEM（Ordered Subsets EM）は投影をサブセットに分けて更新を早める方法で、PET/SPECT では標準的です。

$$x_j \;\leftarrow\; \frac{x_j}{\sum_i a_{ij}} \sum_i a_{ij}\, \frac{p_i}{(A\mathbf{x})_i}$$

FBP と逐次近似は、投影数、ノイズ、欠損の有無、計算時間に応じて使い分けます。投影数が十分でノイズが小さい条件では、1 パスで処理できる FBP が高速です。一方、投影数が少ない場合、線量が低い場合、データに欠損がある場合には、FBP の $|\omega|$ フィルタがアンダーサンプリングとノイズを増幅し、画質が低下します。逐次近似は測定データとの整合を反復的に改善するため、このような条件でもアーチファクトを抑えられます。ただし計算量は FBP より大きくなります。臨床 CT では、被ばく低減を目的として逐次近似再構成とその発展手法が利用されています。

シミュレーション：逐次近似の反復を動かす

手法を選び、「1反復」で 1 ステップずつ、「▶ 連続実行」で自動的に反復を進めます。ゼロ画像から再構成像が立ち上がっていく様子を観察してください。ART は数反復で像が現れ、SIRT はゆっくり滑らかに、OSEM は MLEM よりサブセット数倍ほど速く収束します。手法・$\lambda$・サブセット数・投影数・ノイズを変更すると反復は最初からやり直しになります（MLEM/OSEM は非負画像を前提とするため、手法切替時の初期化は必須です）。表示ウィンドウは真値ファントムの窓に固定してあり、RMSE は FOV 中心 90% 円内で真値と比較した値です。

シミュレーション：収束の観察

左のグラフは、上のシミュレータで実行中の手法の相対残差 $\|Ax - p\|/\|p\|$ を反復ごとにプロットします（縦軸は対数）。残差は「現在の画像を順投影した結果が測定とどれだけ食い違うか」の指標です。右の比較モードでは、現在と同一の条件（投影数・ノイズ・$\lambda$・サブセット数）で全 5 手法を 10 反復ずつ実行し、残差曲線を重ねて表示します。ART の初期収束の速さ、SIRT の穏やかさ、サブセットによる OSEM の加速を見比べてください。ノイズがある場合、残差はノイズ水準より下がらなくなります。それ以上下げると、ノイズへの過適合になるからです。

シミュレーション：FBP との使い分け

厳しい条件のプリセットでは、そのサイノグラムから FBP（Ram-Lak）と、選択中の手法を 10 反復した結果を並べます。少数投影（24 本）では FBP に強い放射状ストリークが現れるのに対し、逐次近似ではストリークが抑えられます。強ノイズ（$I_0 = 3\times10^3$）では Ram-Lak FBP がノイズを増幅し、SIRT や OSEM はより滑らかな像になります。手法や $\lambda$ を変えると結果はクリアされます。

実務での使い分け

再構成は連立一次方程式 $Ax = p$ の求解として定式化できます。ART（レイごとの射影）、SIRT/SART（正規化した同時更新）、MLEM/OSEM（ポアソン統計に基づく乗法更新）は、いずれも残差を減らしながら解へ近づく方法です。少数投影や低線量でも FBP より安定する場合がありますが、計算量は大きく、反復回数、緩和係数、サブセット数によって画質が変わります。

参考文献

• Andersen AH, Kak AC. Simultaneous algebraic reconstruction technique (SART). Ultrasonic Imaging 6, 81–94 (1984).
• Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Transactions on Medical Imaging 1, 113–122 (1982).
• Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Transactions on Medical Imaging 13, 601–609 (1994).
• Thibault JB et al. A three-dimensional statistical approach to improved image quality for multislice helical CT. Medical Physics 34, 4526–4544 (2007).