核医学機器工学概論 ２
断層画像CT（ Computed Tomography）を得る方法
1. フィルタ重畳逆投影法
FBP （ Filtered Back Projection ）2. 逐次近似再構成法 Iterative Reconstruction
MLEM （Maximum Likelihood Expectation Maximization）

OSEM （ Ordered Subsets Expectation Maximization ）プログラム PSF.exe の実行。フォルダ PSF 内 のPSF.exe をダブルクリック。
このテキストウィンドウ内をクリックしてから １を入力してSimple back projectioを実行。
次にプログラム PSF.exe を終了、再度実行。
２を入力してFiltered backprojectionを実行。

選択する再構成フィルタは、（real space filter は）フォルダ PSF 内 にあるRealRAMP256.ｔｘｔ を選択。

画像中心の１画素だけに画素値を与え、他を０に設定した２次元画像を作成。
その像を180度方向から横から透視したと想定した像 Pθを作成（１度ごと 合計１８０枚の線状画像）。
（スライドでは５度ごとの画像を表示）

１８０枚の線状画像 Ｐθ
を重ね合わせると、広がりをもつ分布が得られる。
点広がり関数 PSF(Point Spread Function)
PSF = 1/r

１８０枚の線状画像 Pθに、実空間 Ramp (RL) filter h を畳み込む （ Pθ＊ｈ ）。

（青く表示された画素はマイナスの値）

１８０枚の線状画像に
ＲＬ filter ｈ を畳み込みＰθ＊h を作成し、これらを重ね合わせると、広がりが消失し、１点の分布に戻る。
（青は マイナスの画素値）

PSF＊ｈ（ｒ） は 点に戻る

回転中心からの距離ｒに反比例した濃度に補正するフィルタ
１／ｒ を正確な断層像ｇに畳み込んだ像が lである。

式で表現すると ｌ＝ｇ＊(１／ｒ) となる。
ｌ、ｇ、１／ｒ のフーリエ変換を Ｌ、Ｇ、F(１／ｒ) と表現すると、畳み込みの定理より
Ｌ＝Ｇ・F(１／ｒ) となる。
２次元フーリエ変換の公式の極座標表現を用いると、（ｆｒは周波数空間上の原点からの距離 ）F(１／ｒ) ＝∫∫（１／ｒ）exp（－ｊ(2πｒｆr ））ｒｄｒｄθ ＝１／ｆｒ

これより

L ＝ Ｇ ／ｆｒ

なのでＧ＝Ｌ・ｆｒＧ ＝Ｌ・ｆｒの 意味は、 ２次元周波数空間上で、単純重ね合わせ画像をフーリエ変換した２次元データ Ｌ に、フィルタ関数 ｆｒ（ ｆｒは周波数空間上の原点からの距離）をかけると、正しい再構成画像をフーリエ変換したデータ Ｇになる。

Ｇ＝Ｌ・ｆｒに、畳み込みの定理 を用いると、以下のような実空間での計算に変換できる。
この式を 逆フーリエ変換すると、
ｇ＝ｌ＊ｈ（ ｈ は フィルタｆｒ の 逆フーリエ変換 ）この式に、ｌ＝∫ Pθ dθ を代入すると、
ｇ＝ ∫ Pθ dθ ＊ｈｇ＝ ∫( Pθ ＊ｈ）dθ （ｈはθと独立した値なので交換可）
ｇ＝ ∫ Pθ dθ

(Pθ ＝ Pθ ＊ｈ ）FBPの式Pθ に 実空間フィルタ h ( = frの逆フーリエ変換 ) を畳み込めば、重ね合せると正確な断層像 ｇ になる２次元透視画像 Pθ を算出できる。 これを Filtered Back Projection （FBP） という。
周波数空間での実際の計算においては、フィルタ Ｈ（ =ｆｒ）は常に正の値であり（絶対値）、さらにサンプリング定理より、ナイキスト周波数以上の成分を削除する必要があるので、周波数空間での再構成フィルタ Ｈは、
Ｈ ＝｜ｆｒ｜（ｆｒがナイキスト周波数未満の場合）

Ｈ ＝０（ｆｒがナイキスト周波数以上の場合）となる。 これをＲａｍｐフィルタという。

Ｒａｍｐフィルタを逆フーリエ変換して実空間Ｒａｍｐフィルタｈにしてから、実空間でPθにｈを畳み込む。

Pθ = Pθ * ｈ（ * は畳み込み演算 ）周波数空間 RL フィルタ ( １／ｆｒ ) を１次元逆フーリエ変換 して、実空間 RLフィルタ ｈ（ｒ） を作成。
フォルダ RAMP256 内の RAMP256.exeを実行。

Ramp フィルタの作成
Ramp.c 実行結果

プログラム CTFBP.exe の実行。フォルダ CT 内 のCTFBP.exe をダブルクリック。
このテキストウィンドウ内をクリックする。
選択するプロジェクション
データは、フォルダ CT 内 のCTprojection を選択。
１を入力してSimple back projectioを実行。
Disp FBP Process? と出たらEnterキーを押す。
Simple Back projection が実行される様子を観察。

求めたい正確な断層像 ｇ を算出するために、多方向の角度θ から ２次元透視像 Pθ を測定。

２次元透視像 Pθ は、サイノグラムの各行 θ の１次元データを ２次元に引き伸ばした画像。

例として胸部の１断面のサイノグラムから180度方向から横から透視したと想定した像 Pθを作成（ １度ごと、合計１８０枚の２次元透視画像 ）。
（ スライドでは５度ごとの画像を表示）

180度方向から透視した像 Pθを重ね合わせるとぼけて画像中心の値が持ち上がった像 l を得る。
正確な断層像 ｇ と l の関係式は、
ｌ=g*h

正確な断層像 g に 点広がり関数 h が 畳み込まれぼやけた断層像 l を得る と考える。

単純重ね合わせ再構成法 Simple Back Projection
収集された各々の角度に傾いた２次元透視画像（ Pθ ）を全部単純に重ねると再構成画像ができる。
(回転中心近傍の値が盛り上がった不正確な画像。)
スライスｊにおけるサイノグラムを求める。
サイノグラムの各スライスの１次元配列は、各々の角度から収集されたデータ。
サイノグラムの各スライスの１次元配列から、収集された各々の角度に傾いた２次元透視画像Pθを作成する。
Pθを単純に重ね合わせた画像をｌとするとｌ＝∫ Pθ dθ(Simple back projection)ｌ は、回転中心部ほど重ね合せ回数が多くなり、中心から距離が遠いほどカウントの低い像になる。

プログラム CTFBP.exe の再実行。
フォルダ CT 内 のCTFBP.exe をダブルクリック。
このテキストウィンドウ内をクリック。
選択するプロジェクション
データは、フォルダ CT 内 のCTprojection を選択。
２を入力してFiltered back projectioを実行。
選択するReal space filter はフォルダ CT 内 のRealRAMP256.txtを選択。
Disp FBP Process? と出たらEnterキーを押す。

１８０枚の２次元透視画像 Pθに実空間 Ramp (RL) filter h を畳み込む（ Pθ＊ｈ ）。
(青く表示された画素はマイナス値）

Pθ ＊ h を重ね合せると 正確な断層像 ｇ になる。
ｇ＝ ∫ Pθ dθ

(Pθ ＝ Pθ ＊ｈ ）FBPの式

畳み込みの定理
データ g をフーリエ変換して、その周波数空間成分 Ｇ に周波数空間 Rampフィルタ H をかけて逆フーリエ変換すると、実空間で、実空間 Rampフィルタ ｈ をｇ に畳み込みしたデータと同じになる。
( G x H と ｇ * h は等価演算 )データ数６４の１次元データ ｇ[0]～ｇ[63] と、右図に示すようなデータ数３２の実空間フィルタ
h[0]～h[31] がある （ h[0] が原点）。
データ ｇ に フィルタ h を畳み込んで配列 l [0]～ l [63] に書き込むプログラムを記述して下さい。
//--------- Convolution ｈ[ ] into ｇ[ ] ----------int
doublei, j;g[64] , h[32] , l[64] , gh ;for( i = 0 ; i <=63 ; i++ ){// i <64 と書いても同じ

gh = 0.0 ;for( j = 0; j <=31; j++ ){ if( i+j <=63 ) gh += g[ i+j ] * h[ j ] ; }for( j = 1; j <=31; j++ ){ if( i-j >= 0 ) gh += g[ i-j ] * h[ j ] ; }}l[i] = gh ;

ｇ に h を畳みこんだ結果を l とする。座標 i における ｌ の値 l[i] は、
ｌ[ i ] = g[ i ]*h[0] + g[i+1]*h[-1] + g[i+2]*h[-2] + g[i+3]*h[-3] + ・ ・ ・+ g[i-1]*h[ 1] + g[i-2]*h[ 2] + g[i-3]*h[ 3] + ・ ・ ・

ここで フィルタ h は偶関数（左右対称）なので、h[-j] = h[j] を代入し、ｌ[ i ] = g[ i ]*h[0] + g[i+1]*h[1] + g[i+2]*h[2] + g[i+3]*h[3] + ・ ・ ・+ g[i-1]*h[1] + g[i-2]*h[2] + g[i-3]*h[3] + ・ ・ ・
これを Ｃ言語で表す。 ｈ の要素数 j は h[0] から h[31] まで。
g の 要素数 は g[0] から g[63] までなので、g の 添字を表す i+j は63 以下、 i-j は 0 以上の場合だけ加算する条件式を入れる。
gh = 0.0 ;for( j = 0; j <=31; j++ ){ if( i+j <=63 ) gh += g[ i+j ] * h[ j ] ; }for( j = 1; j <=31; j++ ){ if( i-j >= 0 ) gh += g[ i-j ] * h[ j ] ; }l [ i ] = gh ;これを i が 0 から 63 まで繰り返すと、 l [0] から l [63] が計算される。

プログラム PETFBP.exe の実行。フォルダ PETFBP
内 の PETFBP.exe をダブルクリック。
このテキストウィンドウ内をクリックする。
選択するプロジェクション
データは、フォルダ PETFBP
内の PETsinogram を選択。
Select slice おすすめスライス
は、３６、３７、３８あたり。
２を入力してFiltered back projectioを実行。
選択フィルタは２種類。
RealRAMP256.txt
RealSheppLogan256.txt

プログラム PETFBP.exe
PET画像を Simple BackProjection、または FBP( Filtered BackProjection ) で再構成する。

PETsinpgram ファイル
PETの収集データは各スライスのサイノグラムが並んでいる ３次元データ。（ CTと同じ ）。

各スライスのサイノグラムが並んでいる３次元データを並べ替えれば、SPECTのプロジェクション像と同じ並びになる。 Sinogram[ i ][θ][ j ] = Projection [ i ][ j ][θ]PETFBP.c を実行して、Simple BackProjection とFBP の 再構成画像の違いを観察し、Ramp等のフィルタ関数の必要性を理解してください。

Ｒａｍｐフィルタは理論的には正確な再構成フィルタだが、実際の臨床データに適用すると ナイキスト周波数以上の高周波成分を不連続に遮断するために生じる高周波ノイズや 放射状アーチファクトが目立つ場合が多いため、臨床では高周波成分を抑制する工夫を施した再構成フィルタが用いられている。
ＳＰＥＣＴで用いられる一般的な再構成フィルタとして、Shepp＆Loganフィルタがある。
周波数空間上で、ナイキスト周波数近傍の高周波成分を連続的に減衰させるように設計されており再構成画像に生じる高周波ノイズや放射状アーチファクトを抑制する効果をもつ。

実空間フィルタを畳込んだ２次元透視画像（Pθ）を重ね合わせると フィルタ逆投影再構成像ができる。
フィルタの形状で、再構成画像の高周波成分が変る。

convolution.c の実行結果
半円形のサンプルデータ g に実空間で Rampフィルタ を畳み込む。

畳み込み Convolution

l=g*h

サンプルデータ g[ ] に 実空間Rampフィルタ
h[ ] を畳み込んで、配列 l [ ] に書き込む。
//--------- Convolution h[ ] into g[ ] --------------printf("¥n¥n Convolution Real space Ramp filter "); scanf("%c",&yn);

for( i=0; i<=255 ; i++ ) {gh = 0.0 ;for(j=0; j<=127; j++){ if( i+j < 256) gh += g[ i+j ] * h[ j ]; }for(j=1; j<=127; j++){ if( i-j >= 0 ) gh += g[ i-j ] * h[ j ]; }}l[i] = gh ;

周波数空間で Rampフィルタ を かける

convolution_FFT.cの実行結果

実空間 サンプルデータ ｇ をフーリエ変換して、周波数空間の 実数成分 Gr、虚数成分 Gi を表示。
実空間 サンプルデータ ｇ の積分値 が、周波数空間の実数成分 Gr の最初の成分（原点、直流成分）と同じ値になることを確認する。

周波数空間の 実数成分 Gr、虚数成分 Gi の二乗和の平方根（パワースペクトル G ）を表示。
実空間 サンプルデータ ｇ の積分値 が、周波数空間のパワースペクトル G の最初の成分（原点、直流成分）と同じ値になることを確認する。

//------ Convolution H[ ] into G[ ] ------for (i= 1; i<= 256; i++ ){ GHr[i] = Gr[i] * H[i] ; }for (i= 1; i<= 256; i++ ){ GHi[i] = Gi[i] * H[i] ; }//- Transform convolved GH into real space FFT1D( GHr, GHi, 256, -1 ) ;データ g の周波数成分 Gr、Gi に、周波数空間 Rampフィルタ H をかけて、逆フーリエ変換する。

データ g の周波数成分 Gr、Gi に、周波数空間 Rampフィルタ H をかけて、逆フーリエ変換すると、実空間で、実空間 Rampフィルタ ｈを ｇ に畳み込みしたデータと同じになる。 ( G x H と ｇ * h は等価演算 )

積分値が 元データ ｇ と同じことを確認する。
（カウント総和は保たれる。）