医用画像機器工学実習 断層画像CT（ Computed Tomography） を得る方法 フィルタ重畳逆投影法

Presentation on theme: "医用画像機器工学実習 断層画像CT（ Computed Tomography） を得る方法 フィルタ重畳逆投影法"— Presentation transcript:

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

2

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

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

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

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

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

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

9

10

11

12

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

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

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

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

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

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

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

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

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

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

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

24

25

26

27 畳込みの定理　Convolution

28

29

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

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

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

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

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

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

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

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

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

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

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

41 再構成画像μの、画素 [ 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

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

43

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

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

46 撮像した全方向について λ／（Σ 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

47 の値をかけて、次の推定画像 μｋ＋１ の画素値を算出。 μｋ＋１ ／μｋ ＝ Σ 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

48 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 を ３回繰り返した場合と同等の画像を得られる。

49 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ループ が何重も連続する。膨大な計算量。

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

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

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

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