前回の記事では，CTの投影がRadon変換によって与えられることを示し，
それをmatlab のコードにより実装しました。具体的には，投影領域に対して，ある角度で投影した際，以下のように与えられることを説明します。

ここでδはデルタ関数と呼ばれるもので，０になる部分で1，他の部分で0になるようなものです。

本記事で考えている題材は，いろいろな角度から投影が得られたとき，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

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 と変更した結果です。投影の回数が減ると，誤差が目に見えるようになってくることがわかります。

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