PET検査の仕組み（第7回）「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$の結果です）。

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()