SPECT（ Single Photon Emission CT ）PET（ Positron Emission CT ） の原理
断層画像を得る方法
フィルタ重畳逆投影法
FBP （ Filtered Back Projectionj）
逐次近似再構成法 IterativeItti R
Reconstruction
tti
MLEM （ Maximun
M iLikelihood
Lik lih d Expectation
Et ti Maximization
M i i ti ）OSEM （ Ordered Subsets Expectation Maximization ）Convolution 重畳積分

(*)

FBP
FilteredB kBackProjection

Si l B kP j ti によるSPECT画像の作成
SimpleBackProjectionによるSPECT画像の作成

サイノグラムの各スライスの１次元配列は、
グ
各
３２方向の各々の角度から収集されたデ タ。
３２方向の各々の角度から収集されたデータ。
このデ タから、収集された各々の角度に傾いたこのデータから収集された各々の角度に傾いた２次元透視画像（ Pθ ）を作成する。

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

123I-IMP Brain SPECT 単純重ね合わせ像は回転中心近傍のカウントが持ち上がる。画像が不鮮明。

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

L ＝ Ｇ ／ｆｒ／なので

Ｇ＝Ｌ・ｆｒＧ ＝Ｌ・ｆｒの 意味は、 ２次元周波数空間上で、単純重ね合わせ画像をフーリエ変換した２次元データ Ｌ に、フィルタ関数 ｆｒ（ ｆｒは周波数空間上の原点からの距離）をかけると、正しい再構成画像をフーリエ変換したデータ Ｇになる。
Ｇ＝Ｌ・ｆｒに、畳み込みの定理 を用いると、以下のような実空間での計算に変換できる。
変換できる
この式を 逆フーリエ変換すると、
ｇ＝ｌ＊ｈ（ ｈ は フィルタｆｒ の 逆フーリエ変換 ）（ * は畳み込み演算 ）この式に、ｌ＝∫ Pθ dθ を代入すると、
ｇ＝ ∫ Pθ dθ ＊ｈ（ * は畳み込み演算 ）ｇ＝ ∫( Pθ ＊ｈ）dθ （ｈはθと独立した値なので交換可）
ｇ＝ ∫ Pθ dθ

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

Pθ = Pθ * ｈ（ * は畳み込み演算 ）フィルタ逆重畳画像再構成法
Filtered Back Projection (FBP)
サイノグラムの２次元透視画像 Pθ に 、実空間フィルタ h ( = frの逆フfrの逆フーリエ変換
リエ変換 ) を畳み込めば、を畳み込めば重ね合せると正確な断層像 ｇ になる２次元透視画像 Pθ を算出 きる算出できる。

ｈを畳み込んだ２次元透視画像 Pθ を単純に重ね合わせた画像 g は、
g = ∫ Pθ ｄθ

（ Filtered back projection ）Filtered back projection は、回転中心部ほど重ね合わせ回数が多くなる誤差がフィルタｈによって補正され、正しい再構成画像となる。
タ によ
補 され
し 再構成画像となる
ｈを畳み込んだ２次元透視画像 Pθ を単純に重ね合わせた画像 g =∫ Pθ ｄθ
（ Filtered back projection ）は、回転中心部ほど重ね合わせ回数が多くなる誤差がフィルタｈによって補正され、正しい再構成画像となる。

畳込みの定理 Convolution

畳み込みの定理
デ タ g をフーリエ変換して、データをフ リエ変換してその周波数空間成分 Ｇ に周波数空間 Rampフィルタ
Rフィルタ H をかけて逆フ リ 変換すると逆フーリエ変換すると、実空間で 実空間 Rampフィルタ
実空間で、実空間ルタ ｈ をｇ に畳み込みしたデータと同じになる。に畳み込みしたデ タと同じになる( G x H と ｇ * h は等価演算等価演算 )convolution.c の実行結果
半円形のサンプルデ タ g に半円形のサンプルデータ
実空間で Rampフィルタ
R
フ ルタ を畳み込む。
畳み込む畳み込み Convolution

l=g*h

サンプルデータ g[ ] に 実空間Rampフィルタ
h[ ] を畳み込んで、配列を畳 込ん
列 l [ ] に書き込む。
書き込む//--------- Convolution h[ ] into g[ ] --------------for( i=1; i<=256 ; 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( ii-jj >= 1 ) gh += g[ ii-jj ] * h[ j ]; }}l[i] = gh ;

実空間 Rampフィルタ h の積分値は １ なので、なので畳み込む前後のｇ と l（ = g*h） の積分値は同じ。

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

convolution_FFT.c
convolution
FFT cの実行結果

データ g の周波数成分 Gr、Gi に、周波数空間 Rampフィルタ H をかけて、逆フーリエ変換すると、逆変換す実空間で 実空間 Rampフィルタ
実空間で、実空間Rフィルタ ｈを ｇ に畳み込みしたデータと同じになる。 ( G x H と ｇ * h は等価演算 )Ｒａｍｐフィルタは理論的には正確な再構成フィルタだが、実際の臨床データに適用すると実際の臨床デタに適用すると ナイキスト周波数以上の高周波成分を不連続に遮断するために生じる高周波ノイズや 放射状ア チファクトが目立つ場合が多いため、アーチファクトが目立つ場合が多いため、臨床では高周波成分を抑制する工夫を施した再構成フィルタが用いられている。
用いられている
ＳＰＥＣＴで用いられる一般的な再構成フィルタとして、Shepp＆Loganフィルタがある。
Shepp＆Loganフィルタがある
周波数空間上で、ナイキスト周波数近傍の高周波成分を連続的に減衰させるように設計されており再構成画像に生じる高周波ノイズや放射状アーチファクトを抑制する効果をもつ。
抑制する効果をもつ

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

周波数空間上の６４ ６４画像用 Shepp
６４ｘ６４画像用
Sh& LoganLフィルタ

逆フーリエ変換された実空間 Shepp & Loganフィルタ
実空間フィルタの積分値は １

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

123I-IMP Brain SPECTFBP with Ramp filter

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

吸収補正後の 18F-FDG 脳PET サイノグラム

脳PETのサイノグラムデ タを並べ替えると脳PETのサイノグラムデータを並べ替えるとプロジェクション像

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

ＦＢＰ画像

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

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

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

再構成画像μの、画素 [ i ] [ j ] に対するサイノグラムλ[ yi ] [ yj ] への寄与率（検出確率）は、４次元配列
次元配列 C [ i ][ j ][ yii ][ yjj ] となる。となる
λ＝Σ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]C[i][j][ i][ j]OSEM は、 yj （サイノグラムの角度成分）の計算ループ
を間引を間引いてC （λ／（Σ（ ／（ C μ）））） ／ ΣCの値を求めて 次の推定画像 μの画素値を算出。の値を求めて、次の推定画像
μの画素値を算出
例えば、 yj が 0, 1, 2, 3, 4, 5, 6, 7, 8 の ９方向で、subsets を ３ に設定すれば、に設定すれば
まず yj = 0,
まず、yj0 3
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 計算結果
S b t 2
Subsetsk=0

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

k=4k = 10k = 20

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

再構成画像μの、各画素に対するサイノグラム λ への検出確率 C の分布を広くする。の分布を広くする 点広がり関数を加味。
点広がり関数を加味

検出確率 C の分布に 点広がり関数を加味すると、サイノグラム上でサイグラで広広がった分布がた分布
再構成画像上で 点 に収束するので分解能が向上し、ノイズが抑制される。

PVC （ Partial Volume Correction ）検出確率 C に ガウス分布 をかけて 再構成

MRP （Median Root Prior）画像再構成
ベイズ（Bayes）画像再構成法のひとつ
ベイズ（Bayes）画像再構成法のひとつ。
再構成式の中に条件式（先験確率 Prior ）を加える。

μｋ＋１／μμｋ =ΣC（λ／（ΣCμ））／
μ（ΣC + Prior ）Prior は、着目する画素値μは 着目する画素値μｋ と 周辺画素の中央値 M（median）の差が 小さくなるようにμμｋ＋１ を修正する。
Median Root Prior ＝ β （ μｋ ー M ） ／ M
M は、周辺画素（3x3画素など）の中央値 （median）。
β 、効果
βは、効果を調整する
す パラメータ （ 0 ＜ ββ＜ 1 ）。
（ β が ０ のときは OSEM と同じ ）18F-FDG 脳PETMRP
M matrix
i sizei
3x3x3

MRP 再構成法 は、画像輪郭を保ちながら 統計ノイズを抑制する。
画像の定量解析、統計解析に有効と考える。

MRPは、βを適切に設定すれば、画像の変動係数COV を抑制し、分解能をあまり劣化させない。

Alenius S, Ruotsalainen U. Eur J Nucl Med (1997) 24

MRP法は、繰り返し回数が多すぎても画素値の低下が少ない （画質が劣化しない）。

Alenius S, Ruotsalainen U. Eur J Nucl Med (1997) 24

逐次近似再構成で 画像を改善させる工夫
新たな先験確率 Prior の開発
検出確率 C の修飾
被検者のトランスミッションデ タを利用し被検者のトランスミッションデータを利用し被検者内で生じる散乱成分を補正する
検出確率分布 C を作成して画像を再構成する、など。
像を 構成す 、。