基礎CT値

核医学装置にも X線CTが装備された機器が増えつつあり、CT値 （HU : HounsfieldUnit ）についての知識を確認しておきたい。CT値は、組織の密度に比例する値である。

CT断層像の画素値の基になる値は、体内の各組織の線減弱係数 μｔ だが、透過するX線のエネルギーによって変化する。その変化を抑制し、さらに臨床的な理解度を容易にするために μｔ に比例した値が下記の式で算出され、CTの画素値に使われる。
CT値 = 1000 x （ μt - μw ） ／ μw
μw : 水のX線吸収係数 （線減弱係数）
μt : 組織のX線吸収係数 （線減弱係数）
理解が困難そうな式だが、具体的な数値を代入してみると非常にわかりやすい。
空気のCT値 は -1000
1000 x （ μ air - μw ） ／ μw = -1000 （HU）
厳密には空気には質量があり線減弱係数 μ air は０ではないが、水や人体組織と比べると極めて小さい値なので、CT値を計算する場合は μ air = ０ とする。
水のCT値 は 0 （ 比重１が 0 HU ）1000 x （ μw - μw ） ／ μw = 0 （HU）
水の２倍の密度のCT値は 1000 （ 比重２が 1000 HU ）1000 x （ 2 μw - μw ） ／ μw = 1000 （HU）
人体の中では歯が最も高密度で、歯の密度は約２でCT値は約1000である。

CT値は体内の密度（比重）を表す数値である。
肺野のCT値が約-800で、密度は0.2、肺気腫の症例では、肺野が黒っぽく見えるが、CT値が約-900で、密度が0.1程度に低下していることを示す。

血液のCT値は約60だが、これは密度（比重）が1.06であることを示す。
造影剤で白く見える血液のCT値は、約１00から200を示すが、造影剤によって血液の密度が1.1から1.2程度に重くなっていることを示す。

基礎

フィルタ重畳逆投影法

フィルタ重畳逆投影法（FBP; Filtered Back Projection）
CTやSPECT、PETなどの人体断層画像を得るための基本的な計算法はフィルタ重畳逆投影法である。
例として、5x5画素の断層像gを得るために45度おきに4方向から透視（投影）データを収集したとする。
断層像gの正解は、中心に画素値１があるだけで周囲の画素値は全て0とする。

断層像g を4方向から投影して得るデータを求めると、左図①のようにプロジェクションデータ1から4が収集される。
収集されたプロジェクションデータから断層像gを得るために、まず単純に4方向から得た投影データを左図②のように5x5画素の行列内に逆投影する。断層像は、中心に最大の画素値を得るが、その周囲にも0ではない画素値が分布する。
これが単純重ね合わせ（Simple BackProjection）で、本当は中心に画素値１があるだけなのに、周囲にも画素値が発生して点像がぼけている。
単純重ね合わせでは点が点像に戻らない。それを補正するため投影デー
タにフィルタｈをかける。ここでは簡単に要素数3個の行列(-1,3,-1)をｈとする。投影データの値に行列ｈをかける。 このようにスカラー値に行列を畳込むように掛ける演算を畳込み、重畳（convolution）という。 行列ｈを畳込まれた各プロジェクションを左図③に示す。これらを4方向から5x5画素の行列内に逆投影する。これをフィルタ重畳逆投影（FBP；FilteredBack Projection）という。

基礎

フィルタ重畳逆投影法２行列ｈを畳込まれた4方向の投影データ（演算記号＊は畳込みを表す）を、左図④のように5x5画素の行列内に逆投影すると、単純重ね合わせよりも、最大値を示す中心画素の近傍の画素値が低く抑えられ明瞭な点像が得られ、本来の断層像gに近い画像になる。
つまりフィルタ重畳逆投影法は、点像をぼかさずに点像に戻す事ができるので、横方向からの投影データから正しい断層画像を算出できる。
画像中心の１画素だけ値が１で、他は全て0 の256x256画素の画像を作成し、これを画像ｇとする。その像を1度毎に180度方向から投影した投影像Pθを作成する。左図では5度毎のPθを表示。
180枚の投影像Pθから単純重ね合わせ像を作ると、画像中心の１画素が最大値を示し、その近傍に中心からの距離 r に反比例する値が示される。
つまり１点の画素値が 1/r の関数に広がっている。
この単純重ね合わせ像を画像ｌとする。

１/ｒとは、１点の画素を、画素からの距離ｒに反比例した値の分布に変換する関数と解釈すると、画像ｇに関数１/ｒを畳込んだ像が l である。
式で表現すると ｌ＝ｇ ＊ (１/ｒ) となる。
l、ｇ、1/ｒのフーリエ変換を Ｌ、Ｇ、F(１/ｒ) と表現すると、畳み込みの定理より
Ｌ＝Ｇ・F(１/ｒ) となる。
ここで １/ｒのフーリエ変換を計算すると、 F(１/ｒ) ＝ １/ fr となる。
（ fr はフーリエ空間での周波数を表す。）
したがって、Ｌ＝Ｇ / fr 、 Ｇ ＝ Ｌ・fr
フーリエ空間（周波数空間）で、画像ｌのフーリエ変換像Ｌにフィルタ frをかけると、画像ｇのフーリエ変換Ｇが得られる。
この周波数空間での計算は、以下の式を用いて実空間で済ませられる。

g = l ＊ h （ ＊ は畳込み演算）（h は fr の逆フーリエ変換）

基礎

フーリエ変換
フーリエ変換 Fourier transform
実空間の曲線や画像データなどを周波数空間のデータに変換する演算。実空間データに含まれる cos 成分と sin 成分を各周波数ごとに算出する。
オイラーの公式 exp( j θ) = cosθ + j sinθ を用いて、以下の式で示される。
曲線 g(t) の周波数分布を G(f) とすると、G(f) = ∫ g(t) exp(-j・2πft) dt （ j は虚数単位）

逆フーリエ変換も類似した式で表される。

g(t) = ∫ G(f) exp( j・2πft) df
高速フーリエ変換 FFT Fast Fourier Transform
フーリエ変換、逆フーリエを行う際に使用されるアルゴリズム、プログラム。
フーリエ変換の公式は、波形データが連続曲線と想定して exp(-j・2πft) という三角関数のオイラー式表記をしているが、実際に扱う数値はデジタルデータであり離散的なので、 exp(-j・2πft) の項を Wnk のように変形して高速に計算できる巧みなアルゴリズムである。
１/ｒのフーリエ変換が １/ fr になる理由
２次元フーリエ変換の公式の極座標表現を用いると（ rとθ は実空間の原点からの距離と偏角、frは周波数空間の原点からの距離）関数 g(r) のフーリエ変換は ∫ ∫ g(r) exp(-j・2π r fr ) r dr dθ （ヤコビアンに注意）
g(r) = 1/r とすると ∫ ∫(1/r) exp(-j・2π r fr ) r dr dθ = ∫ ∫ exp(-j・2π r fr ) dr dθ = 1/fr

周波数空間を理解することは難しい印象を持つが、周波数空間や周波数分布を表現しているものは一般にも結構ある。
たとえば楽譜は、音楽の周波数分布、曲のフーリエ変換とも解釈できる。五線譜の下側の音符は低音、低周波数成分を表し、上側の音符は高音、高周波成分を示している。
虹は白色光に含まれる色の周波数分布を示している。雨粒中で屈折しやすい高周波の紫色成分から、屈折しにくい低周波の赤色成分まで、周波数の順序で色が並んでいる。

基礎

畳込みの定理畳込みの定理実空間での畳込み演算は、周波数空間では掛け算になる、という定理である。
l (x) に h(x) を畳込んで g(x) を算出する式は、
g(x) = ∫ l(n) h(x-n) dn これを g = l＊h と表し、畳込み（convolution）という。
h(x-n) のフーリエ変換を H(f) とすると、
h(x-n) = ∫ H(f) exp( j・2πf(x-n) ) df （逆フーリエ変換の公式を利用）
これを、畳込みの式 g(x) = ∫ l(n) h(x-n) dn に代入し（つまり g = l＊h である）
L(f) = ∫ l(n) exp( - j・2πfn) dn （フーリエ変換の公式）
g(x) = ∫ G(f) exp( j・2πfn) df （逆フーリエ変換の公式）を利用すると、
g(x) = ∫ l(n) ∫ H(f) exp( j・2πf(x-n) ) df dn= ∫ l(n) ∫ H(f) exp( j・2πfx) exp(- j・2πfn) ) df dn= ∫ { ∫ l(n) exp( - j・2πfn) dn } H(f) exp( j・2πfx) df= ∫ L(f) H(f) exp( j・2πfx) df= ∫ G(f) exp( j・2πfn) df
これよりG(f) = L(f) H(f) である。
したがって g = l＊h であれば、G(f) = L(f) H(f) である。
投影像 Pθ を単純重ね合わせした画像ｌを２次元フーリエ変換し、周波数空間でフィルタ fr を乗算したデータを２次元逆フーリエ変換すると、正確な断層画像gが得られる。
畳込みの定理を使うと、これらの周波数空間での処理が実空間で簡略化される。
投影像 Pθ にフィルタh （ fr の１次元逆フーリエ変換）を畳込み、それを重ね合わせると正確な断層画像g が得られる。

基礎

Rampフィルタ
周波数空間でのフィルタ frとは、どのような形をしているのか。１次元周波数空間の縦軸の周波数成分が fr の値を持つフィルタであるが、横軸が周波数軸すなわち fr である。つまり周波数空間で、正比例のような形状になる。

周波数にはマイナスの値は無いが、高速フーリエ変換（FFT）のアルゴリズムの都合で見かけ上、マイナスの周波数データが存在する。縦軸の周波数成分にはマイナスは無いので、周波数がマイナス側の成分は絶対値となり、原点を中心に左右対称な絶対値 |fr| の関数として表現される。
さらにサンプリング定理によって、ナイキスト周波数（0.5 cycle / pixel）以上の高周波成分はノイズなので除去する。
よって、このフィルタは上図のようになり、これをRampフィルタという。
周波数空間の RampフィルタをH と表現する。
H を１ 次元逆フーリエ変換すると、実空間でのRampフィルタ h が左図のように算出される。

画像中心の１画素だけ値が１で、他は全て0 の256x256画素の画像ｇを1度毎に180度方向から透過投影した投影像Pθを作成する。次に全ての Pθ にRampフィルタhを畳込む（式で表すと Pθ ＊h）（上図で青色の部位はマイナスの画素値）（全てのPθ
＊hの断面が hと同じ形状をしている）。さらに、全ての Pθ ＊h を重ね合わせて画像ｌを算出する。
式で表すと ｌ＝∫ (Pθ ＊h) dθ
画像ｌは中心の１画素だけ値が１で、周囲は殆ど0 に近い値を示し、ほぼ画像ｇに戻っている。
Pθの単純重ね合わせ像では、画像中心の１画素が最大値を示し、その近傍に中心からの距離 r に反比例する1/r の関数に画素値が広がっていたが、Rampフィルタhを畳込むと、その広がりが補正されて元の画像が再現される。このRampフィルタhの機能によって断層画像を算出することができる。

基礎

ナイキスト周波数、統計雑音CT装置やガンマカメラの総合分解能が 半値幅4mmのとき、収集マトリクスサイズは（収集画像の１画素の１辺の長さ）何mmに設定するのが適切か。
正解は2mmであるが、これはデジタルデータのサンプリング定理に基づいて求められる。
総合分解能の半分以下の画素長で撮像しないと、せっかく細かい収集をしたのに画像をぼかして記録することになり、カメラの性能が無駄になる。
デジタル画像に含まれる最高周波数の波は、１周期分が２画素で表される。
つまり１画素には1/2個の最高周波数の波が入っており、その最高周波数をナイキスト（Nyquist）周波数という。
サンプリング定理だけを考慮すると、さらに小さい画素長で収集すればよいのではと考えるが、画素が細かすぎると１画素内のカウントが減り、統計雑音が増加して画像の信号/雑音比（S/N 比）が悪くなる。
したがってカメラの総合分解能のちょうど半分の長さが、最も適切な収集マトリクスサイズになる。
ナイキスト周波数の理論は、アナログデータをデジタル変換（AD変換）する際のサンプリング間隔を考えると理解しやすい。
１波長を8分割（信号の8倍の周波数でサンプリング）した場合は、かなりアナログ波形の形状を保ったサンプリング（標本化）が出来ている。
では、１波長を2分割したサンプリング（信号の2倍の周波数でサンプリング）ではアナログデータの波長と振幅の値は、ぎりぎりでデジタルデータに記録されている。
ところが、アナログデータの半波長を少しでも超えたサンプリング（信号の周波数の2倍より低い周波数でサンプリング）では、アナログデータの波長も振幅もわからない標本化が行われている。
つまり、データの最高周波数（ナイキスト周波数）をデジタルデータに記録（標本化）するためには、その2倍の周波数で標本化する必要がある。
統計雑音（ホワイトノイズ）の式
Nカウントの画素値に含まれる統計雑音は√Nである。
これは正規分布の公式から導かれる式である。
一般的に核医学画像の統計雑音は10%以下が良い。
そのためには最低でも画素値が100カウント必要である（√ 100 = 10）。カウントが多いほど S/N比は良くなる。

基礎

画像再構成フィルタ
CT装置が収集した各々の角度における投影画像（ Pθ ）から断層画像を再構成する方法を考える。
スライスｊでのサイノグラムが求められたとする。
サイノグラムの各スライスの１次元配列は、各々の角度から収集された投影データである。
サイノグラムの各スライスの１次元配列から、各々の角度に傾いた２次元投影画像Pθを作成する。
２次元投影画像 Pθ は、サイノグラムの各 θ の１次元データを ２次元に引き伸ばした像である。
Pθを単純に重ね合わせた画像をｌとすると
ｌ＝∫ Pθ dθ (Simple back projection)
画像ｌ は、回転中心部ほど重ね合せ回数が多くなり、中心から距離が遠いほど値の低い像になる。
フィルタ重畳逆投影法（FBP）
回転中心からの距離ｒに反比例した濃度に補正す
るフィルタ 1/r を正確な断層像ｇに畳み込んだ像が lである。式で表現すると ｌ＝ｇ＊(1/r) となる。
ｌ、ｇ、1/r のフーリエ変換を Ｌ、Ｇ、F(1/r) と表現すると、畳み込みの定理より Ｌ＝Ｇ・F(1/r)となる。ここで F(1/r) ＝ 1/fr なので L ＝ Ｇ/fr
Ｇ＝Ｌ・fr この式を逆フーリエ変換して g = l ＊ h

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

基礎

Shepp & Logan フィルタ、
プレフィルタRampフィルタは理論的には正確な再構成フィルタで、X線CTやPETなど投影データの画質が良好な場合は適切な断層画像を算出する。しかしSPECTなど投影データのカウントが低く統計雑音が多い場合は、高周波成分のS/N比の悪さが影響して、左図の心筋SPECTのように放射状アーチファクトが目立つ断層画像を再構成する。
そこで、SPECT画像再構成には、S/N比の良くない高周波成分を抑制する効果を加えたShepp & Loganフィルタという画像再構成フィルタが用いられる。

周波数空間でナイキスト周波数に近い高周波成分の増幅を抑えることで、高周波雑音を低減させ、放射状アーチファクトが少ない断層像を再構成する。

S/N比の良くないSPECT投影（プロジェクション）画像の高周波成分を抑制する方法として、再構成を行う前に、投影画像に高周波成分抑制フィルタをかける方法がある。これを前処理フィルタ、プレフィルタという。
良く使用される前処理フィルタは、バターワースフィルタ（Butterworth）である。ナイキスト周波数から中間周波数成分を左図のような連続関数で除去する。
遮断周波数を決めるだけで使い易いフィルタである。
凝った前処理フィルタとして、ウィナーフィルタ（Wiener）がある。通信技術にも使われているフィル
タで、理論的に正確な雑音除去フィルタである。
画像情報を多く含む中間周波数成分を増幅しながら、低周波成分の抑制と、ナイキスト周波数近くの高周波成分も抑制する効果をもつが、パラメータが多くて各症例に適切なウィナーフィルタを設定するのは臨床的にはかなり困難である。

逐次近似画像再構成法

基礎

逐次近似再構成法

Iterative Reconstruction

PETやSPECT装置が収集した各々の角度における投影画像（ Pθ ）から断層画像を再構成する別の方法として確率的手法を用いた逐次近似画像再構成法がある。
フィルタ重畳逆投影法と比較し、投影データのカウントが少なくても放射状アーチファクトが少なく、画質の良好な断層画像を算出できる再構成だが、短所として、計算時間が長い、繰り返し回数を決める理論がない（何回繰り返し計算をするか主観的に決める場合が多い）。したがって、再構成画像の定量性もフィルタ重畳逆投影法と比較し、不確かである。
MLEM （Maximum Likelihood Expectation Maximization）という方法が基本的な逐次近似再構成法であるが、計算量が多く時間がかかるので、これを改良したOSEM （Ordered Subsets Expectation Maximization ）という計算時間を短縮できる方法が使用されている。
角度yj における座標yi の投影データ（サイノグラム）をλ [yi][yj] と示す。
サイノグラムから再構成される断層画像を μ[i][j] とする。
逐次近似画像再構成法は、これらの4個の変数 yi、yj、i、j による4次元の繰り返し演算である。
再構成画像μ の、画素 [ i ] [ j ] に対するサイノグラム λ[ yi ] [ yj ] への寄与率C（検出確率）は 4次元配列 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 以上 になる。

MLEM、OSEM基礎

撮像した全方向について λ / Σ C μ） の平均（検出確率 C をかけた加重平均）を求める。式で表すと、

Σ C （λ / Σ C μ）/ ΣC この式の値は要素数が i x j の配列なので正確に記述すると、
yi yji jyi yjΣ Σ C[i][j][yi][yj] ( λ[yi][yj] / ( Σ Σ C[i][j][yi][yj] μk [i] [j] ) ) / Σ Σ C[i][j][yi][yj]k 番目の再構成画像μk の 各画素ごとに Σ C （λ /（Σ C μ））/ ΣCの値をかけて、次の推定画像 μk+1 の画素値を算出する。式で表すと、
μk+1 / μk ＝ Σ C （λ /（Σ C μ））/ ΣC
正確に記述すると

逐次近似再構成法 MLEM の式

yi yji jΣ Σ C[i][j][yi][yj] ( λ[yi][yj] / ( Σ Σ C[i][j][yi][yj] μk [i] [j] ) )μk+1 [ i ][ j ] / μk [ i ][ j ] ＝yi yjΣ Σ 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 の値で μk を計算する。
次に、yj = 1, 4, 7 の値で μk を基に μk+1を計算する。
更に、yj = 2, 5, 8 の値で μk+1を基にμk+2 を計算する。
計算量は MLEM の １回繰り返しと同量だが、MLEM を ３回繰り返した場合とほぼ同等の画像を得られるので、OSEMでは計算量、計算時間が 1/3 になる。
OSEMは、サイノグラムを全部使わないで断層画像を推定していくが、逐次異なる角度からのサイノグラムを使うのでプロジェクションデータは全部有効に利用している上で、計算の分量が減らせるアルゴリズムである。

基礎

PVC、MRP

逐次近似画像再構成法は、プログラムの中に再構成画像の画質を改善させるために、様々な分布関数やアルゴリズムを加えられる長所がある。
PVC （ Partial Volume Correction ）検出確率 C の分布に 点広がり関数（ガウス分布）を加味すると、サイノグラム上で 広がった分布が再構成画像上で 1点 に収束するので分解能が向上しノイズが抑制される。

MRP （Median Root Prior）
画像再構成MRP 再構成法 は、画像輪郭を保ちながら 統計ノイズを抑制する。
画像の定量解析、統計解析に適する画質を導出する。ベイズ（Bayes）画像再構成法のひとつ。 再構成式の中に条件式（先験確率 Prior ）を加える。

μ k+1 / μk = Σ C ( λ / ΣCμ ) ) / ( ΣC + Prior )Prior は、着目する画素値 μk と 周辺画素の中央値 M（median）の差が 小さくなるように μ k+1 を修正する。
Median Root Prior＝β （ μk ー M ） ／ M

M は、周辺画素（3x3画素など）の中央値 （median）。β は効果を調整する パラメータ （ 0 < β< 1 ）。（ β が 0 の場合は一般のOSEM と同じ 。）
メディアンフィルタの効果を加えながら、サイノグラム（横から投影した全方向からのプロジェクションデータ）から、断層画像を逐次推定していく。