逐次近似法による画像再構成を色々試してみた　1

CTの画像再構成手法には、大きく解析的手法と代数的手法が存在します。
解析的手法としてはCBP（Convolution back projection : コンボリューション補正逆投影法）や、FBP（Filtered back projection : フィルタ補正逆投影法）等が挙げられ、代数的手法としては逐次近似法が挙げられます。現在最も主流なのはFBP法（数学的にCBP法と等価）であり、短時間で高精度な画像が得られる手法として知られます。一方で、逐次近似法は繰り返し計算を必要とするため、解析的手法と比べ大量の計算時間が必要となりますが、少ない投影枚数でも高コントラストな画像を得ることができ、かつ外乱やアーチファクトの影響を受けにくいという性質を持ちます。近年では、コンピュータの高性能化が相まって、逐次近似法を利用したCT診断が増加傾向にあるように感じております。

また、我々はCTの性能向上（アーチファクト低減等）に関する研究を行っていますが、それらは全てFBP法に基づいたものでした。そこで、ここでは現在存在している様々な逐次近似法についてメモ書きし、それらを実装して性能を確認してみます。

 

逐次近似法は順問題と逆問題で構成されていて、フィードバックを与えながら順問題を逆問題を繰り返し計算していくことで、解に近づけていくものです。

CTの画像再構成における逐次近似法では、順問題として順投影、逆問題として逆投影が用いられます。

ここで、順投影・逆投影を定式化しておきます。入射X線強度を$I_{\rm in}$、透過X線強度を$I_{\rm out}$、2次元のX線吸収係数の分布を$\mu(x,y)$とすると、投影は次のように表されます（Beer-Lambertの法則）。

$$
I_{\rm in}(X,\theta) = I_{\rm out} \exp\left( – \displaystyle\int_{-\infty}^{\infty} \mu(x,y) dY \right)
$$

また、$I_{\rm in}$が既知であることから、投影データ$p$は上式を変形して次のように表されます。

$$
p(X,\theta) = \displaystyle\int_{-\infty}^{\infty} \mu(x,y) dY
$$

上式は$\mu(x,y)$を$Y$軸に沿って線積分した形になっていて、この変換はラドン変換とも呼ばれます。

一方、逆投影は被写体の方向に向かって直前状に投影データを加えていく操作であり、次式にように表されます。

$$
\ \ \ \ g(x,y) = \displaystyle\frac{1}{2\pi} \int_0^{2\pi} p(X,\theta) d\theta
$$

これで逐次近似を行う準備が整いました。それでは私の頭部断層画像（512px × 512px）を使って、実験していきます。

この断層画像に順投影処理を施すと、サイノグラムと呼ばれる次のようなデータが得られます。

今回は、512方向から撮影を行いました。このサイノグラムに対してFBP法を適用すると、次の再構成画像が得られます。

およそ1秒で再構成されました。若干放射状のようなモアレ状のようなボケがありますが、十分綺麗です。短時間で診断に十分な画質、これがFBPが最も使われる理由です。

また、今回はMAE（平均二乗誤差）、および再構成画像の画素値（最大と最小）を再構成精度の定量的指標に用いました。MAEの値は小さいほど良く、最大値は255、最小値は0に近いほど優れた定量性を持っているということになります。数値を見ると、定量性はあまり良くないようです。というより、FBP法では再構成の後キャリブレーションが必要だったかもしれませんが、気にせず進みます。

 

SIRT（Simultaneous Reconstruction Technique）法は、最も単純な逐次近似法による画像再構成手法の1つで、次の手順で表されます。

1. 適当に予想した画像を作成する
2. 予想画像に対して順投影処理を行う（順問題）
3. 手順2の結果と計測で得た投影データの差分を逆投影し、予想画像に加える
4. 手順2～3を、差分が十分少なくなるまで繰り返す

$k$回目の画像を$\mu(x,y)^{(k)} $、$k+1$回目の画像を$\mu(x,y)^{(k+1)} $としたとき、この繰り返しの式は次のように表されます。

$$
\mu(x,y)^{(k+1)} = \mu(x,y)^{(k)} + \displaystyle\frac{\alpha}{2\pi} \int_0^{2\pi} \left\{ p(X,\theta) -\int_{-\infty}^{\infty} \mu(x,y)^{(k)} dY \right\} d\theta
$$

5回、10回、50回繰り返した結果は次の通り。

最大最小を見ると、50回目ではそこそこの定量性を示せていますが、全体的にぼや～っとした画像になってしまっています。

 

SIRT法（加算型）との違いは、フィードバックの際に差分を利用するのではなく、比率を求めて乗算を利用するという点だけです。以下の式で表されます。

$$
\mu(x,y)^{(k+1)} = \mu(x,y)^{(k)} \times \displaystyle\frac{1}{2\pi} \int_0^{2\pi} \left\{ \frac{p(X,\theta)} {\int_{-\infty}^{\infty} \mu(x,y)^{(k)} dY} \right\} d\theta
$$

また、統計学的な理論に基づいているML-EM（Maximum Likelihood-Expectation Maximization Method）法というものがあるのですが、私はSIRT法（乗算型）との違いがイマイチわかりません。というより、等価であると思っております。そのため、グラフ等ではSIRT法（乗算型）をML-EM法と書いたりしていますが、同じものを指していることにご留意下さい。もし、これら2つは明確に分類されるべきという方がいらっしゃいましたら、是非ご教授頂きたい思いです。

ML-EM法は、係数行列を用いて次のように表されることが多いです。

$$
\lambda_j^{k+1} = \displaystyle\frac{\lambda_j^k}{ \displaystyle\sum_{i=1}^I C_{ij}} \sum_{i=1}^I \frac{y_i C_{ij}}{\displaystyle \sum_{m=1}^J C_{im}\lambda_m^k}
$$

5回、10回、50回繰り返した結果は次の通り。

MAEの値も、最大最小値も、SIRT法（加算型）と比べて優れた値を示しています。

 

SIRT法は全ての角度からの投影を行った後、差分を求めて逆投影をしていましたが、ART(Algebraic Reconstruction Technique）法では投影と逆投影を投影角度ごとに行い、画像を更新していく手法です。繰り返しの式は次のように表されます。

$$
\mu(x,y)^{(k+1)} = \mu(x,y)^{(k)} + \displaystyle \frac{\alpha}{2\pi} \int_{\theta_k} \left\{ p(X,\theta_k) -\int_{-\infty}^{\infty} \mu(x,y)^{(k)} dY \right\} d\theta
$$

ART法はSIRT法と比べ画像の更新回数が多く、画像再構成の収束を早めますが、反復計算において収束せず発散してしまう可能性が多少高くなります。一方、SIRT法では収束は遅くなりますが、値が発散する可能性は低くなります。

5回、10回、50回繰り返した結果は次の通り。

5回の段階でかなり綺麗な画像が得られていることから、収束がかなり早いことがわかります。定量性も良好です。

 

ART法で、フィードバックの際に比率を乗算を用いる手法です。次の式によって表されます。

$$
\mu(x,y)^{(k+1)} = \mu(x,y)^{(k)} \times \displaystyle \frac{1}{2\pi} \int_{\theta_k} \left\{ \frac{p(X,\theta_k)} {\int_{-\infty}^{\infty} \mu(x,y)^{(k)} dY} \right\} d\theta
$$

5回、10回、50回繰り返した結果は次の通り。

若干画像が他と比べて暗い印象ですが、収束は早く、定量性も良好です。

 

最適化に用いる評価関数の勾配を利用して、最適解を探索してく手法として、勾配法（Gradient method：GM）があります。式は後述の共役勾配法に含まれるので、省略します。

5回、10回、50回繰り返した結果は次の通り。

いまいちです。全体的にぼけた画像になっています。ただ、定量性はまずまずです。

 

最急降下法（Steepest descent method：SDM）は、勾配法の収束を早めるため、探索の方向をその時点で最も勾配が急な方に向ける手法です。式は後述の共役勾配法に（ほぼ）含まれるので、省略します。

5回、10回、50回繰り返した結果は次の通り。

勾配法と比べ、収束が早く、定量性も優れていることが確認できます。ただ、計算時間が約1.5倍に伸びています。

 

共役勾配法（Conjugate gradient method：CGM）とは、最急降下法の効率を更に向上させるために、共役勾配方向へ探索を進める手法です。評価関数の勾配分布が楕円形にひろがっており、楕円の中心が最適解と仮定した時の、楕円の中心方向が共役勾配方向に相当します。共役勾配法は、近年盛んに行われているCTやMRIにおいて少ないデータから画像再構成を行う研究（圧縮センシングによる画像再構成）において、使われることが多いです。

線形方程式
$$
C\vec{\mu} = \vec{y},\ \ \vec{g} = C^T (\vec{y} – C\vec{f})
$$
を共役勾配法によって解くときの更新式は、修正ベクトル$\vec{d}_k$，勾配ベクトル$\vec{g}_k$，2つのスカラー$\alpha_k, \beta_k$を用いて以下のように与えられます。
$$
\vec{\mu}_{k+1} = \vec{\mu}_k + \alpha_k \vec{g}_k
$$
$$
\alpha_k = \displaystyle\frac{\vec{d}_{k}^{\ T} \vec{g}_k}{\vec{d}_k^T C^T C\vec{d}_k} = \frac{\vec{g}_k^T \vec{g}_k}{\vec{d}_k^T C^T C\vec{d}_k}
$$
$$
\beta_k = -\displaystyle\frac{\vec{g}_{k+1}^T C^T C\vec{d}_k}{\vec{d}_k^T C^T C\vec{d}_k} = \frac{\vec{g}_k^T \vec{g}_k}{\vec{g}_{k+1}^T \vec{g}_{k+1}}
$$
$$
\vec{d}_0 = \vec{g}_0
$$
$$
\vec{d}_{k+1} = \vec{g}_{k+1} + \beta_k \vec{d}_k
$$
$$
\vec{g}_0 = C^TC(\vec{y} – C\vec{f}_0) = C^T C\vec{y}
$$
$$
\vec{g}_{k+1} = \vec{g}_k – \alpha_k C^T C \vec{d}_k
$$

5回、10回、50回繰り返した結果は次の通り。

収束速度、定量性ともに、最急降下法よりも優れていることがわかります。計算時間は最急降下法と同じくらいです。

 

それぞれの手法で256回繰り返しを行った際の、MAEの推移をグラフにプロットしました。ASIRTは加算型SIRT（Addition-type SIRT）、AARTは加算型ART（Addition-type ART）、MARTは乗算型ART（Multiplication-type ART）を意味します。

収束速度、MAE共に最も優れているのは、ART（乗算型）でした。ART（加算型）も、かなり優れています。共役勾配法は、ある瞬間において優れたMAEの値を示していますが、繰り返し回数を増やすことでMAEが上昇する結果になっています。

この結果だけ見ると、逐次近似法ではART法を使うのが良いように思えますが、上述のように、ART法では発散する危険が上昇するようです。また、共役勾配法は少ない投影枚数の際に使われることが多い等、それぞれの手法に一長一短があるようです。今回はそれらの差をしっかり確認することができませんでした。次回は、投影枚数を減らしてみたり、ART法では発散するケースを探してみたりして、各手法の性能を評価していきたいです。