· 

PET検査の仕組み(第7回)「PETシミュレーション」

(著)山たー

 X線CTの場合は外部から光子(放射線)の方向をコントロールできるが、PETやSPECTの場合は光子(放射線, γ線)が放射される方向はランダムです。果たして、ランダムな方向に光子が放出される場合においても、画像の再構成が可能なのでしょうか。もちろん、可能だからPET装置は成り立っているのですが、少し疑問だったため、シミュレーションをすることにしました。

 

シミュレーションにおける仮定

 シミュレーションにおいては次のことを仮定します。


 ・1つの画素につき、光子は画素値に比例する数が放出される。すなわち、比例定数を$N_p$として$N_p\times$(画素値)[個]の光子対が放出されるとする。
 ・放出された光子は2次元平面上を移動する。
 ・2つの光子は正確に正反対($180^\circ$)に放出される。
 ・光子の放出角度はいずれかの検出器において検出される。そのため角度は検出器の数に依存する。
 ・偶発同時計数、散乱同時計数による誤差は存在しない。
 ・被写体による光子の吸収はない。
 ・陽電子飛程は存在しない(=0m)。核種が存在する場所で対消滅が生じる。

 この条件を第6回で用いた記号を用いて表すと、$C(r,\theta)=S(r,\theta)=0, A(r,\theta)=D(r,\theta)=1$となり, $P(r,\theta)=E(r,\theta)$となっています。つまり、エミッションデータが真の線和に等しいという仮定を行っているのです。

 実装する上では、1つの光子が放出され、どの検出器に入るかが検出する前から分かっているとします。また、光子は画素ごとに順番に放出します。今回は偶発同時計数が存在しないと仮定しているので、光子の放射される順番は像の形成に影響しません。

  テストに用いた画像は前回までと同様、Shepp-Loganの頭部ファントムです。計算の流れは次の図のようになります(なお、ここで用いている画像は検出器数400, $N_p=2000$の結果です)。

シミュレーションの結果

 結果は次の図のようになりました。Radon変換を用いた投影を用いた場合ほどではないのですが、ある程度の像を再構成することができています。当然ですが、検出器の数と検出光子数が増えるほど画像は鮮明になります。ただし、これは誤差のない理想的な環境を想定したものであり、実際には検出光子数が増えると偶発同時計数が増加してしまいます。

ソースコード

 今回用いたプログラムのコードを以下に記します。

import numpy as np
import matplotlib.pyplot as plt
import random

from skimage.io import imread
from skimage import data_dir
from skimage.transform import iradon

#phantomイメージの読み込み。2次元配列で保存。0~1の値を取る.
image = imread(data_dir + "/phantom.png", as_grey=True)

#画像サイズ
imgNum=max(image.shape)
imgNum2nd=int(imgNum/2) #imgNum/2

#thetaを等間隔に設定
thetaNum=180 #光子の検出器の数(全体の検出器の数はこの2倍になる)
theta = np.linspace(0., 180., thetaNum, endpoint=False)

sinogram = np.zeros((imgNum,thetaNum)) #サイノグラムの初期化
photonNum=50; #放射される光子の個数の比例定数

#光子をいずれかの方向に放射
for i in range(imgNum):
    #繰り返し回数を表示
    if((i+1)\%100==0):
        print(i+1)
    for j in range(imgNum):
        #光子数は画素値×photonNum
        for k in range(int(image[i][j]*photonNum)):
            angle_num=random.randint(0,thetaNum-1) #放射される方向は一様分布に従う
            angle=np.radians(theta[angle_num]) #角度をradianに変換
            #座標をずらす(中央を原点にする)
            x=i-imgNum2nd
            y=j-imgNum2nd
            #原点を通り角度angleの直線との符号付き距離(この数式は結果から調整している)
            d=int(x*np.sin(angle)+y*np.cos(angle))
            #範囲に制限をつける
            const_d=np.clip(d+imgNum2nd,0,imgNum-1)
            sinogram[imgNum-const_d-1][thetaNum-angle_num-1]+=1

#中央のみ平滑化(これがないと中央にノイズが入る)
for i in range(thetaNum):
    sinogram[imgNum2nd-1][i]=(sinogram[imgNum2nd][i]+sinogram[imgNum2nd-2][i])/2

#途中経過の表示
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
#原画像の表示
ax1.set_title("Original")
ax1.imshow(image, cmap=plt.cm.Greys)
#得られたサイノグラムの描画
ax2.set_title("Sinogram")
ax2.set_xlabel("Projection angle (deg)")
ax2.set_ylabel("Projection position (pixels)")
ax2.imshow(sinogram, cmap=plt.cm.hot,
           extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')
fig.tight_layout()
plt.show()

#FBPによる画像の復元
reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)

#復元画像を最大値1の配列にする(正規化する).最小値まで0にすると精度が落ちる
arrayMax = reconstruction_fbp.max()
normalized_reconstruction_fbp=reconstruction_fbp.astype(float)/arrayMax.astype(float)

#復元画像と元の画像の誤差
error = normalized_reconstruction_fbp - image

#平均二乗誤差
print('FBP rms reconstruction error: \%.3g' \% np.sqrt(np.mean(error**2)))

#結果画像表示
imkwargs = dict(vmin=-0.2, vmax=0.2)
fig2, (ax3, ax4) = plt.subplots(1, 2, figsize=(8, 4.5),
                               sharex=True, sharey=True,
                               subplot_kw={'adjustable': 'box-forced'})
ax3.set_title("Reconstruction Filtered back projection\n
              Number of photons for value of pixels: "+str(photonNum)+"\n
              Number of detectors:"+str(thetaNum))
ax3.imshow(normalized_reconstruction_fbp, cmap=plt.cm.Greys)
ax4.set_title("Reconstruction error\nFiltered back projection")
ax4.imshow(error, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()

参考文献

・日本医用画像工学会(2012)『医用画像工学ハンドブック』, 日本医用画像工学会