PET検査の仕組み（第3回）「画像再構成②：逆Radon変換」

　第2回では2次元画像$f(x,y)$から、その投影データ$p(r,\theta)$をRadon変換によって求める方法について述べました。また、X線CT, PET, SPECTではそもそもとして投影データはRadon変換によって得られるのではなく、放射線と生体との作用によって得られるのだということも述べました。よってすべきことは投影データ$p(r,\theta)$を元の2次元画像$f(x,y)$に対応付ける変換, すなわち逆Radon変換を求めることです。今回はこの逆Radon変換について説明します。

　逆Radon変換を導出するためには投影切断面定理という定理が必要です。そのため、先ずは、この定理の導出をすることにします。 $f(x,y)$の2次元フーリエ変換は次式で定義されます。 $$ F(u,v)=\int_{-\infty}^{\infty}\!\!dx\int_{-\infty}^{\infty}\!\!dy\ f(x,y)\exp[-j(ux+vy)] $$ ただし、$j$は虚数単位、$u,v$は角周波数です（なんで$i$を虚数単位としないのかというと、電気電子系においては$i$を電流を表す文字に使っているためです。そのため、アルファベットで次の$j$を虚数単位にします)。次に、$u=\omega\cos\theta, v=\omega\sin\theta$として直交座標から極座標に座標変換すると、 $$ F(\omega\cos\theta,\omega\sin\theta)=\int_{-\infty}^{\infty}\!\!dx\int_{-\infty}^{\infty}\!\!dy\ f(x,y)\exp[-j\omega(x\cos\theta+y\sin\theta)] $$ となります。一方、投影データ$p(r, \theta)$の$r$に関する1次元フーリエ変換$P(\omega,\theta)$を計算すると、 \begin{align*} P(\omega,\theta)&=\int_{-\infty}^{\infty}\!\!dr\ p(r,\theta)\exp(-j\omega r)\\ &=\int_{-\infty}^{\infty}\!\!dx\int_{-\infty}^{\infty}\!\!dy\int_{-\infty}^{\infty}\!\!dr\ f(x,y)\delta(r-x\cos\theta-y\sin\theta)\exp[-j\omega r]\\ &=\int_{-\infty}^{\infty}\!\!dx\int_{-\infty}^{\infty}\!\!dy\ f(x,y)\exp[-j\omega(x\cos\theta+y\sin\theta)] \end{align*} となります。よって $$ F(\omega\cos\theta,\omega\sin\theta)=P(\omega,\theta) $$ が成り立ちます。これを投影切断面定理と呼びます。変換がいっぱい出てきてややこしいですが、図にすると以下のようになります。

　投影切断面定理により逆Radon変換が求められます。2次元画像$f(x,y)$を求めるために逆Radon変換として次の手順を踏めば良いわけです。

[Step1]
投影データ$p(r,\theta)$の$r$に関する1次元フーリエ変換$P(\omega,\theta)$を次式によって求める. $$P(\omega,\theta)=\int_{-\infty}^{\infty}p(r,\theta)\exp(-j\omega r)dr$$
[Step2]
$P(\omega,\theta)$を放射状に配置し, 極座標から直交座標へ$u=\omega\cos\theta, v=\omega\sin\theta$により座標変換を行う。こうして$F(\omega\cos\theta,\omega\sin\theta)=F(u,v)$を求める.
[Step3]
$F(u,v)$に対して2次元フーリエ逆変換 $$ f(x,y)=\frac{1}{4\pi^2}\int_{-\infty}^{\infty}\!\!du\int_{-\infty}^{\infty}\!\!dv\ F(u,v)\exp[j(ux+vy)] $$ を行い, $f(x,y)$を求める。

　見ての通りですが、画像が再構成できていません。もちろん私のプログラミングが悪いという話も大いにありますが、こうなる原因は上記の逆Radon変換の手順のStep2にあります。投影データは連続的な$\theta$についてではなく離散的な$\theta$について取っており、極座標から直交座標に変換する際に、一部$F(u,v)$のデータが欠けてしまうのです。ゆえにStep2とStep3を分けず、まとめて処理することを次の第4回で考えます。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

#scikit-imageをインポート
from skimage.io import imread
from skimage import data_dir
from skimage.transform import radon
import cv2

#phantom imageは頭部の断層画像っぽい画像のこと
image = imread(data_dir + "/phantom.png", as_grey=True)

#データのパラメータ
thetaNum=500
imageNum= max(image.shape)
imageNum2nd=imageNum/2
theta = np.linspace(0., 180.,thetaNum, endpoint=False)
thetaRad=np.radians(theta)
sinogram = radon(image, theta=theta, circle=True)

# 信号を生成
G=np.zeros((imageNum,thetaNum),dtype=np.complex)

# 高速フーリエ変換
for i in range(thetaNum):
    G[:,i] = np.fft.fft(sinogram[:,i])

#実部と虚部に分解し,shift
G_real =np.fft.fftshift(np.real(G))
G_imag= np.fft.fftshift(np.imag(G))
#原画像の表示
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
ax1.set_title("FFT real")
ax1.set_ylabel("Projection angle (deg)")
ax1.imshow(G_real, cmap=plt.cm.Greys)
ax2.set_title("FFT imagin")
ax2.set_ylabel("Projection angle (deg)")
ax2.imshow(G_imag, cmap=plt.cm.Greys)
fig.tight_layout()
plt.show()

#変数変換後のフーリエ変換
F_real=np.zeros((imageNum,imageNum))
F_imag=np.zeros((imageNum,imageNum))

#極座標変換
for i in range(imageNum):
    for j in range(thetaNum):
        x=np.clip(int(round((i-imageNum2nd)*np.cos(thetaRad[j])+imageNum2nd)),
                  0,imageNum-1)
        y=np.clip(int(round((i-imageNum2nd)*np.sin(thetaRad[j])+imageNum2nd)),
                  0,imageNum-1)
        F_real[x][y]=G_real[i][j]
        F_imag[x][y]=G_imag[i][j]

#途中経過表示
fig2, (ax3, ax4) = plt.subplots(1, 2, figsize=(8, 4.5))
ax3.set_title("FFT2 real")
ax3.imshow(F_real, cmap=plt.cm.Greys_r)
ax4.set_title("FFT2 imagin")
ax4.imshow(F_imag, cmap=plt.cm.Greys_r)
fig2.tight_layout()
plt.show()

#Gaussian平滑化をする場合
"""
F_real =cv2.blur(F_real,(3,3))
F_imag= cv2.blur(F_imag,(3,3))

fig2, (ax3, ax4) = plt.subplots(1, 2, figsize=(8, 4.5))
ax3.set_title("FFT2 real Blur")
ax3.set_xlabel("Projection angle (deg)")
ax3.imshow(F_real, cmap=plt.cm.Greys_r)
ax4.set_title("FFT2 imagin Blur")
ax4.set_xlabel("Projection angle (deg)")
ax4.imshow(F_imag, cmap=plt.cm.Greys_r)
fig2.tight_layout()
plt.show()
"""

F=np.zeros((imageNum,imageNum), dtype=np.complex)
F=F_real+1j*F_imag

#2次元フーリエ逆変換
f_ishift = np.fft.ifftshift(F)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

#結果表示
fig, (ax3, ax4) = plt.subplots(1, 2, figsize=(8, 4.5))
ax3.set_title("Original")
ax3.imshow(image, cmap=plt.cm.Greys)
ax4.set_title("IFFT")
ax4.imshow(img_back, cmap=plt.cm.Greys)
fig.tight_layout()
plt.show()