医学  ·  2018/01/12

PET検査の仕組み（第4回）「画像再構成③：FBP法」

（著）山たー

　第3回では画像再構成を試みましたがうまくいきませんでした。今回はフィルタ補正逆投影(filtered back projection: FBP)法と呼ばれる方法を用いて画像の再構成をすることを考えます。とは言え、特別にFBP法と呼ばれてはいるのですが、実はFBP法は変数変換を行った2次元逆フーリエ変換と同じです。なので数学的には前回述べた画像再構成法と同じことです。

 

フィルタ補正逆投影法

$F(u,v)$の2次元逆フーリエ変換を極座標で表すと、$dudv=|J|d\omega d\theta, J= \left[ \begin{array}{cc} \cos\theta&-\omega\sin\theta\\ \sin\theta&\omega\cos\theta \end{array} \right]=\omega $ により、

\begin{align*} f(x,y)&=\frac{1}{4\pi^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)\exp[j(ux+vy)]dudv\\ &=\frac{1}{4\pi^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(\omega\cos\theta,\omega\sin\theta)\exp[j\omega(x\cos\theta+y\sin\theta)]|J|d\omega d\theta\\ &=\frac{1}{4\pi^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(\omega,\theta)\exp[j\omega(x\cos\theta+y\sin\theta)]|\omega|d\omega d\theta\ \ (\because 投影切断面定理) \end{align*}

となります。よってFBP法では次の手順を踏みます。

[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)$を2次元逆フーリエ変換する.

$$ f(x,y)=\frac{1}{4\pi^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(\omega,\theta)\exp[j\omega(x\cos\theta+y\sin\theta)]|\omega|d\omega d\theta $$

　Step2は投影データの空間で高周波強調を行うフィルタ（ランプフィルタ（下図））を作用させ、逆フーリエ変換を行うということを意味しています。もしランプフィルタを除いて計算すると、画像がぼやけてしまいます。ランプフィルタはJacobianに基づくので、人為的に新たな処理を加えているというわけではないのですが、結果として合理的な処理をしていることが分かります。これは大変興味深い点です。

Illustration

FBP法による画像の再構成

　前回と同じテストデータを用いてFBP法により画像の再構成をすると、次図のようになります。

Illustration

明瞭な像が再構成できていることが分かります。このように、変数変換をした2次元フーリエ逆変換を用いればデータの欠損は無くなるのです。

※実はライブラリを使っているので、あまり説得力のない説明である。

 

ソースコード

　thetaNumは投影の方向の数であるので、これを小さくすれば復元度は下がります。なお、このコードはscikit-imageの"Radon transform - skimage v0.14dev docs""

[http://scikit-image.org/docs/dev/auto_examples/transform/plot_radon_transform.html]

を参考にしました。 

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

import numpy as np
import matplotlib.pyplot as plt

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

#画像のインポート
image = imread(data_dir + "/phantom.png", as_grey=True)

imgNum=max(image.shape)
thetaNum=500
theta = np.linspace(0., 180., thetaNum, endpoint=False)

#ラドン変換
sinogram = radon(image, theta=theta, circle=True)

#元画像とサイノグラムの表示
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("Radon transform\n(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による画像の復元(filterのデフォルトはランプフィルタ. filter=Noneでフィルタなしに)
reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)

#復元画像と元の画像の誤差
error = 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\nFiltered back projection")
ax3.imshow(reconstruction_fbp, cmap=plt.cm.Greys)
ax4.set_title("Reconstruction error\nFiltered back projection")
ax4.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
plt.show()

参考文献

・日本医用画像工学会(2012)『医用画像工学ハンドブック』, 日本医用画像工学会
・"Radon transform - skimage v0.14dev docs""

[http://scikit-image.org/docs/dev/auto_examples/transform/plot_radon_transform.html]