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

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

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

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

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

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

プログラム PSF.exe の実行。フォルダ PSF 内 の PSF.exe をダブルクリック。 このテキストウィンドウ内を クリックしてから　１を入力して Simple back projectioを実行。 次にプログラム　PSF.exe を 終了、再度実行。 ２を入力して Filtered backprojectionを実行。 選択する再構成フィルタは、 （real space filter は） フォルダ PSF 内 にある RealRAMP256.ｔｘｔ を選択。

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

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

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

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

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

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

回転中心からの距離ｒに反比例した濃度に補正するフィルタ １／ｒ を正確な断層像ｇに畳み込んだ像が l　である。 式で表現すると ｌ＝ｇ＊(１／ｒ) となる。 ｌ、ｇ、１／ｒ のフーリエ変換を Ｌ、Ｇ、F(１／ｒ) と 表現すると、畳み込みの定理より Ｌ＝Ｇ・F(１／ｒ) となる。 ２次元フーリエ変換の公式の極座標表現を用いると、 （ｆｒは周波数空間上の原点からの距離 ）　 F(１／ｒ) ＝ ∫∫（１／ｒ）exp（－ｊ(2πｒｆr ））ｒｄｒｄθ ＝１／ｆｒ　　 これより L ＝ Ｇ ／ｆｒ なので　Ｇ＝Ｌ・ｆｒ　

Pθ に 実空間フィルタ h ( = frの逆フーリエ変換 )を畳込めば、 この式に、ｌ＝∫ Pθ dθ を代入すると、 ｇ＝ ∫ Pθ dθ ＊ｈ　 ｇ＝ ∫( Pθ ＊ｈ）dθ （ｈはθと独立した値なので交換可） ｇ＝ ∫ Pθ dθ　 (　Pθ ＝ Pθ ＊ｈ ）　FBPの式　　 Pθ に 実空間フィルタ h ( = frの逆フーリエ変換 )を畳込めば、 重ね合せると正確な断層像 ｇ になる ２次元透視画像 Pθ を算出できる。 これを Filtered Back Projection （FBP） という。

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

Pθ ＊ h を重ね合せると 正確な断層像 ｇ になる。 ｇ＝ ∫ Pθ dθ　 (　Pθ ＝ Pθ ＊ｈ ）　FBPの式

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

畳込みの定理　Convolution

メタル アーチファクト Metal Artifact 体内に金属がある場合をシミュレーションして 特定の部位のサイノグラム値を高くすると ＣＴ値が局所的に高い部位の周辺に 線状のアーチファクトが生じる。

フォルダ CT 内 のCTFBP.exe をダブルクリック。 このウィンドウ内をクリック。 選択するプロジェクション データは、フォルダ CT 内 の Stomach_metal_projection.txt ２を入力して Filtered back projectioを実行。 選択するReal space filter は フォルダ CT 内 の RealRAMP256.txtを選択。 HU Centerを５０、 HU Widthを５０程度で Metal artifact の状態を観察。

脳PETのサイノグラムに Rampフィルタを 重畳して重ね合わせ　　　　　ＦＢＰ画像

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

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

各スライスのサイノグラムが並んでいる３次元データを並べ替えれば、SPECTのプロジェクション像と同じ並び になる。　Sinogram[ i ][θ][ j ] = Projection [ i ][ j ][θ]

PETFBP.c を実行して、Simple BackProjection と FBP の 再構成画像の違いを観察し、 Ramp等のフィルタ関数の必要性を理解してください。

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

逐次近似法 サイノグラム λ[ yi ] [ yj ] 再構成画像 μ[ i ] [ j ] 4次元の変数に よる繰り返し計算

再構成画像μの、画素 [128] [10] に対する サイノグラム λ[ yi ] [ yj ] への寄与率（検出確率）

再構成画像μの、画素 [128] [128] に対する サイノグラム λ[ yi ] [ yj ] への寄与率（検出確率）

再構成画像μの、画素 [ i ] [ j ] に対する サイノグラムλ[ yi ] [ yj ] への寄与率（検出確率）は、 ４次元配列 C [ i ][ j ][ yi ][ yj ] となる。 λ＝ΣC μ サイノグラム ＝ Σ（検出確率 x 再構成画像） 正確に記述すると λ[ yi ] [ yj ] ＝ΣΣ C[ i ] [ j ][ yi ][ yj ] μk [ i ][ j ] 　　μk [ i ][ j ]　は、ｋ 番目の繰り返し計算後の画像 i j

プログラム PETOSEM.exe の実行。 フォルダ PET_OSEM 内 の OSEM.exe をダブルクリック。 寄与率をクリック。曲線が出るまで待つ。 yi, yj の変化スライド。　　寄与率曲線の変化を確認 サイノグラムはPETsinogramを選択。 空白枠内に再構成するスライスを入力。 46 あたりを入力してスライス選択を押す。 OSEMを押すたびに繰り返し再構成を実行。 繰り返し再構成回数が表示されるのを確認。

OSEM 計算結果 繰り返し回数を多くするほど 画像が鮮明化することを確認する。 Subsets 2 繰り返し計算回数　ｋ 　　k = 0 k = 2 k = 4 k = 10 k = 20 サイノグラム （ 横から測定した全方向からのデータ ） から、確率の高い断面像を 逐次推定していく。

測定したサイノグラム λ と 再構成画像 μ （初期値は 全画素値１） について λ／（Σ C μ） を求める。 λ／（Σ C μ） ＝ 真のサイノグラム ／ 画像μから推定されるサイノグラム 推定画像μの画素値が、真の値より大きすぎると λ／（Σ C μ） は 1 未満 になる。 推定画像μの画素値が、真の値より小さすぎると λ／（Σ C μ） は 1 以上 になる。

撮像した全方向について λ／（Σ C μ） の平均 （検出確率 C をかけた加重平均）を求めている。 Σ C （λ／（Σ C μ）） ／ ΣC 撮像した全方向について λ／（Σ C μ） の平均 （検出確率 C をかけた加重平均）を求めている。　 正確に記述すると Σ Σ C[i][j][yi][yj] （ λ[yi][yj]／（ΣΣC[i][j][yi][yj] μk [i] [j] ） ） ／ Σ Σ C[i][j][yi][yj] 　　 この式の値は配列（ 要素数は i x j ）　　　　　　　　　　　　　　　 yi yj i j yi yj

の値をかけて、次の推定画像 μｋ＋１ の画素値を算出。 μｋ＋１ ／μｋ ＝ Σ C （λ／（Σ C μ）） ／ ΣC k 番目の再構成画像μｋ の 各画素ごとに Σ C （λ／（Σ C μ）） ／ ΣC　 の値をかけて、次の推定画像 μｋ＋１　の画素値を算出。 μｋ＋１　／μｋ　＝　Σ C （λ／（Σ C μ）） ／ ΣC 　　　　　　 逐次近似再構成法 MLEM、OSEM の式 正確に記述すると μｋ＋１　[ i ][ j ]／μｋ [ i ][ j ] = ΣΣ C[i][j][yi][yj] （λ[yi][yj]／（ΣΣC[i][j][yi][yj] μk [i] [j] ）） ／ ΣΣC[i][j][yi][yj] 　 yi y j i j yi y j

OSEM は、 ｊｙ （サイノグラムの角度成分）の計算ループ を間引いて　 C （λ／（Σ C μ）） ／ ΣC　 の値を求めて、次の推定画像 μの画素値を算出。 例えば、 yj が 0, 1, 2, 3, 4, 5, 6, 7, 8 の ９方向で、 subsets を ３ に設定すれば、 まず、yj = 0, 3, 6 の値で μｋ　を計算する。 次に、yj = 1, 4, 7 の値で μｋ　を基に μｋ+1 を計算する。 更に、yj = 2, 5, 8 の値で μｋ+1を基に μｋ+2 を計算する。 計算量は MLEM の １回繰り返しと同量だが、 MLEM を ３回繰り返した場合と同等の画像を得られる。

OSEM プログラム 単純な加減乗除ばかりだが、 forループ が何重も連続する。膨大な計算量。 //------------------------------------------------------------------------------- // 　OSEM for(k=0;k<20;k++){ 　for(sub=0; sub<8; sub++){ s1 = sub - 2*(int)((double)sub/2.0) ; s2 = 1-s1; 　　　　for(j=0;j<192;j++){ for(i=0;i<192;i++){ S_YC_CM[i][j] = SC[i][j] = 0.0; }} 　　　　　　for(j=0;j<192;j++){ printf("\n j= %d ", j); for(i=0;i<192;i++){ 　for(yj=sub; yj<32; yj+=8){ for(yi=CZL[j][i][yj][0]; yi<=CZL[j][i][yj][1];yi++){ CM=0.0; for(jj=0;jj<192;jj++){ for(ii=CZM[yj][yi][jj][0];ii<=CZM[yj][yi][jj][1];ii++){ CM += C[ii][jj][yi][yj] * M[ii][jj][k][s1]; }} S_YC_CM[i][j] += Yi[yi][yj] * C[i][j][yi][yj] / CM ; SC[i][j] += C[i][j][yi][yj]; }} // yi, yj }} // i, j for(j=0;j<192;j++){ for(i=0;i<192;i++){ if(SC[i][j]>0.) M[i][j][k][s2] = M[i][j][k][s1] * S_YC_CM[i][j] / SC[i][j] ; }}　// j, i 　} // sub for(j=0;j<192;j++){ for(i=0;i<192;i++){ M[i][j][k+1][s2] = M[i][j][k][s2]; }}　// j, i Disp_M(k,s2); printf("\n\nNext iteration ? ");scanf("%c",&yn); if(yn=='n')break; } // k OSEM プログラム 単純な加減乗除ばかりだが、 forループ が何重も連続する。膨大な計算量。

平成２7年　国家試験　　　　解答　２　

断層画像の投影データ Pθ の値の合計は、 どの角度 θ でも、だいたい同じ程度の値になるはず。 逐次近似法を計算するには、初めに適当な初期値 が、断層画像の画素に入っていなければならない。 そこで、まず右方向への透視データの合計を算出。 それを断層画像（この問題では ２ｘ２ の画素数）の すべてに同じ画素値が入っているという初期条件を考える。

断層画像の画素すべてに画素値５が入っていると いう初期条件で、右方向への透視データを算出する。 正しい投影データPθ と、初期条件での投影データの差分を算出。 その差を、透視した画素の数で割って、 それぞれ透視した画素の画素値に加える。

次に、透視の向きθを下方向にして、同様の計算。 正しい下向き投影データPθ と、 計算途中の断層 画像の下向き投影データの差分を算出。 その差を、透視した画素の数で割って、 それぞれ透視した画素の画素値に加える。 このように、逐次、投影データと整合する断層画像を 計算する方法が、逐次近似画像再構成法。