医師に、寄り添うAI。

EIRL - NEXT MEDICAL VISION EIRL - NEXT MEDICAL VISION

EN 資料ダウンロード 無料トライアル

製品

EIRL’s product lineup

資料ダウンロード 無料トライアル

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

EIRL Editor
Aug 27 2020

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

p(r,θ)=∫Df(x,y)δ(r−xcosθ−ysinθ)dxdy  (1)
図1. ProjectionとObjectの関係

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

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

図2. 投影から復元する問題

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

投影切断面定理

投影切断面定理とは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
図3. 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
図4. 全ての投影をプロット

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
図5. 復元結果

十分な投影がなければ?

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

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

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

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

図6. アーチファクト

AI技術

関連する記事

AI医療機器プログラムのバージョンアップ問題

Yuki Shimahara
Jan 28 2021

医用画像位置合わせの基礎⑧〜非線形変換を用いた画像位置合わせ〜

EIRL Editor
Dec 17 2020

医用画像位置合わせの基礎⑦ 〜Fiji pluginを用いたアフィン変換による画像位置合わせ〜

EIRL Editor
Dec 3 2020