CTの原理②〜投影切断面定理とCT再構成の実装〜

CTの原理①～ラドン変換～ - IMACEL Academy -人工知能・画像解析の技術応用に向けて-|LPixel(エルピクセル)

CTの原理に関して詳しく説明します．今回はRadon変換についてです．

本記事で考えている題材は，いろいろな角度から投影が得られたとき，Objectを完全に復元できるのか？という問題です．

投影から復元する問題．

結論をいうと，出来ます．これには，投影切断面定理と呼ばれる以下の定理を用います．

投影切断面定理

投影切断面定理とはCTの投影とその領域の2次元フーリエ変換とを結びつける定理です．具体的には， CTの投影のフーリエ変換が，領域の2次元フーリエ変換の断面に一致するという定理です．領域の2次元フーリエ変換が分かれば，逆フーリエ変換をすることで元の領域を得ることができますので，投影から完全に復元できたことになります．

ここでは具体的な式変形は割愛させていただきます．興味のある方は調べてみてください．

Matlabによる実装

以下，Matlabのコードで実装してみました．
なお，Matlab の関数にはFilter Back Projection と呼ばれる処理を行う関数が入っており，上の定理に加え，周波数領域で空間フィルタをたたみ込むことで，高周波成分のノイズの拡大を防いでいます．

https://uk.mathworks.com/help/images/ref/iradon.html

1. まず画像をロードします．

P = phantom(128);
imshow(P)

Matlab コード1

Phantom 画像

2. 投影を計算します．

theta_0 = pi;  % 全方向からの投影
M = 100;  % 投影角度数．この場合，0からpi までで100方向から投影した．
theta_k = theta_0/M;
[Radon_mat,xp] = radon(P,(0:1:M-1)*theta_k * 180/pi);
plot(Radon_mat);

Matlab コード2

全ての投影をプロットしたもの．

3. 逆変換を計算します．

I1 = iradon(Radon_mat,(0:1:M-1)*theta_k * 180/pi);

%% plot
pcolor(I1);
hold on;
shading flat; shading interp; axis equal tight;
colormap(gray);

Matlab コード3

復元結果

元に戻りましたね！

十分な投影がなければ？

今回は，100方向から投影を行いました．では，投影の回数を減らすとどうなのでしょうか．現実は，CTの撮像をする際，放射線の被曝の問題があります．放射線の被曝を防ぐため，投影の回数をできるだけ減らすことが重要な課題になっているわけです．

これに関して数値実験を行ってみました．上のtheta_0 やMを変更するだけですぐチェックできます．

下の画像は上の設定したMを60,40,20,10 と変更した結果です．投影の回数が減ると，誤差が目に見えるようになってくることがわかります．

このような誤差のことを，アーチファクトと呼びます．

どうやって解決するの？

先ほど述べた通り，投影回数が少ない場合は多くの問題であり得ます．ここで用いられるのがスパースモデリングです．次の記事ではスパースモデリングを実装し，その威力を確かめてみます．

目次

ProjectionとObjectの関係