核医学機器工学概論
令和６年度の本講義は、Web講義です。
ホームページのスライドを見て自習して下さい。
Windowsプログラムの起動が可能な人は試して理解を深めてください。
出席票を毎回提出して下さい。
（講義後１週間以内に。）
試験は、11月18日 10:30 に、C304 で実施。

核医学機器工学概論
断層画像 CT（ Computed Tomography）
を得る方法
核医学断層画像のＰＥＴやＳＰＥＣＴも同じ方法
1. フィルタ重畳逆投影法
FBP （ Filtered Back Projection ）2. 逐次近似再構成法 Iterative Reconstruction
MLEM （Maximum Likelihood Expectation Maximization）
OSEM （ Ordered Subsets Expectation Maximization ）3. ディープ ラーニング法 Deep learning

PET/CTCTの世代分類
第１世代 Translate / Rotate （T / R） 方式
検出器は１個。X線は細く（ペンシルビーム）、X線管球が並進（translate）し、角度を変えて（回転 rotate）撮影。初期のCT。
第２世代 Translate / Rotate （T / R） 方式
検出器は10～20個。X線は１０度程の広がりのナロウファンビームが並進（translate）し、角度を変え（回転 rotate）撮影。回転角がファン角ごとに減り、第１世代より高速化。

第３世代 Rotate / Rotate （R / R） 方式
検出器は500個以上。X線は患者全体にあたるワイドファンビーム。 管球の並進は不要で、管球と検出器の回転で撮像できる。
現在のCTの主流であり、体軸方向へのビームを広げて、ワイドコーンビームによるマルチスライス、およびヘリカルスキャンで、高速、広範囲なCT撮影が可能となった。

フィルタ重畳逆投影法 FBP（Filtered Back Projection）
最も基本的な断層画像再構成法

このCTシミュレーションを、もう少し厳密に、256x256画素の画像で行うプログラムを作成してみた。

工学の領域では i は電流密度を表す場合が多いので、虚数単位は j で表す。

フーリエ変換とは、全ての関数はsin波とcos波の級数で表現できることを利用し、実空間の曲線や画像などのデータを周波数空間での周波数分布関数に変換する算術。フーリエ変換では、実空間関数での畳込み演算が、周波数分布関数で掛け算に変わり、計算が楽になる利点がある。

畳込み

実空間の関数の畳込みは周波数空間では掛算になるナイキスト周波数（Nyquist frequency）
医用画像はデジタル化して保存される。デジタル画像に含まれる最高周波数をナイキスト周波数という。ナイキスト周波数の波長は、その画像データの２画素長である。つまり１画素（pixel）の長さは、最高周波数の波長の半分（0.5 cycle/pixel）である。

このフィルタの積分値は１である点が再構成画像の定量性を保証する。

フィルタ重畳逆投影法（FBP）で計算された断層画像の画素値は、定量性が保たれている。

ゴッドフリー・ニューボルド・ハンスフィールド
Godfrey Newbold Hounsfield （1919 - 2004）

イギリスの電気技術者。1967年、コンピュータを用いたX線断層撮影技術（CT）の開発により、1979年にノーベル生理学・医学賞をアラン・コーマック とともに受賞した。
彼の名は現在のCT値の単位である
Hounsfield Unit (HU) として使用され、この値は -1000 HU を空気、0 HU を水とした吸収率の相対値で定義されている。

Windowsで動くプログラムをホームページにアップロードしました。
CT.zip
興味のある人は動かして、FBPの画像再構成を確認して下さい。
RAMP関数の働きを確認する PSF.exe
PAMP関数の実空間上の逆フーリエ変換データ RealRAMP256.txt
CTの畳込み再構成法を確認する CTFBP.exe
CTのプロジェクションデータ CTprojection
Zip解凍ソフト Ｌhaca で解凍し、プログラムの動作を確認してください。

（解凍ソフトLhacaをインストールして、デスクトップのLhacaのアイコン上へ、解凍したいzipファイルをドロップダウンする。）

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

Disp FBP Process? と出たらEnterキーを押す。

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

PETもFBPで断層画像を作成
脳PETのサイノグラムに Rampフィルタを重畳して重ね合わせ 18F-FDG PET脳画像

脳 ＦＤＧ PET の プロジェクションデータ。
PETの収集データは各角度から撮像された脳が並んでいる ３次元データ。（ CTと同じ ）。

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

MIBI 心筋SPECT再構成を 胆嚢のスライスで行うと胆嚢内の 99mTc-MIBI 停滞部位の再構成像が作成される（99mTc-MIBI は胆汁排泄あり）。
局所的に強いRI分布を示すスライスではフィルタ逆投影再構成像は、放射状アーチファクトが強く出ることが確認できる。 （99mTc-MIBI検査前は絶食の前処置（胆汁が少ない状態）が必要。）

123I-IMP Brain SPECT FBP with Shepp&Logan filter

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

畳込みの計算Convolution

畳込みの定理 Convolution

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

逐次近似画像再構成法で計算された断層画像の画素値は、定量性が保証されていない。 繰返し計算回数で画素値が変動する。

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

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

再構成画像μの、画素 [ i ] [ j ] に対するサイノグラムλ[ yi ] [ yj ] への寄与率（検出確率）は、４次元配列 C [ i ][ j ][ yi ][ yj ] となる。
λ＝ΣC μサイノグラム ＝ Σ（検出確率 x 再構成画像）

正確に記述するとi j

λ[ yi ] [ yj ] ＝ΣΣ C[ i ] [ j ][ yi ][ yj ] μk [ i ][ j ]μk [ i ][ j ] は、ｋ 番目の繰り返し計算後の画像測定したサイノグラム λ と 再構成画像 μ （初期値は全画素値１） について λ／（Σ C μ） を求める。
λ／（Σ C μ）＝ 真のサイノグラム ／ 画像μから推定されるサイノグラム

推定画像μの画素値が、真の値より大きすぎるとλ／（Σ C μ） は 1 未満 になる。
推定画像μの画素値が、真の値より小さすぎるとλ／（Σ C μ） は 1 以上 になる。

Σ C （λ／（Σ C μ）） ／ ΣC撮像した全方向について λ／（Σ C μ） の平均（検出確率 C をかけた加重平均）を求める。
正確に記述すると
yi y ji jΣ Σ C[i][j][yi][yj] （λ[yi][yj]／（ Σ Σ C[i][j][yi][yj] μk [i] [j] ））yi y j／ Σ Σ C[i][j][yi][yj]

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

μｋ＋１ [ i ][ j ]／μｋ [ i ][ j ] =yi y ji jΣ Σ C[i][j][yi][yj] （λ[yi][yj]／（ Σ Σ C[i][j][yi][yj] μk [i] [j] ））yi y j／ Σ Σ C[i][j][yi][yj]OSEM は、 yj （サイノグラムの角度成分）の計算ループを間引いて 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 プログラム 単純な加減乗除ばかりだが、//OSEM//------------------------------------------------------------------------------for(k=0;k<20;k++){forループ が何重も連続する。膨大な計算量。

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} // subfor(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 計算結果
Subsets 2k=0

繰り返し計算回数 ｋk=2k=4k = 10k = 20

サイノグラム （ 横から測定した全方向からのデータ ）から、確率の高い断面像を 逐次推定していく。

OSEMの初期画像は全て画素値１とする。

「寄与率」ボタンを押すと、逐次近似法における確率分布行列C[i][j][yi][yj]の準備を行う計算をする（膨大な行列を用意する作業で、やや時間がかかる）
画素の数（X方向128 x Y方向128）だけの寄与率曲線Cが観察できる「OSEM」ボタンを繰り返し押すと、逐次近似法が繰り返し実行される。（少し時間をあけて連打して下さい。）

平成２7年 国家試験

解答 ２

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

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

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

第３の断層画像再構成法
AI による Deep Learning を用いた方法
数十万通りのサイノグラムと断層画像との組み合わせをコンピュータに学習させ、力業でサイノグラムから断層画像を推定する。

CAE （Convolutional Auto Encoder ）Deep Learningによる画像生成のための学習モデルのひとつ。
高い特徴抽出力を持つConvolutionレイヤーを学習モデルに組込み、データを入力すると最終的に画像が出力されるよう設定する。

CT、SPECT、PETは人体の断層図を得るために、その断面の多方向からの透視像（サイノグラム）を収集。

AI による Deep Learning では、数十万通りのサイノグラムから、それぞれの断層像を学習させる。

Deep Learning による断層画像算出の利点は、１．断層画像を推定する時間が非常に短い。
２．FBP法で出るアーチファクトなどが出ない。
Deep Learning による断層画像の欠点は、定量性が保証できない。

金属義歯を含む上顎部のサイノグラム

FBP法ではDeep Learning
金属義歯周囲に では
金属アーチファクト 金属アーチファクトが出る。
が出ない。