数理工学センターコラム
「物を壊すことなく内部を「見る」には？」

CTの仕組みを数式で表現するため、まず、次のような設定を考えます（図１）。

観察対象となる物体の断面が $xy$ 平面上にあるとし、その中の点 $(x,y)$ におけるX線吸収係数の分布を $f(x,y)$とします。 $f(x,y)$ はほぼ通常の物質の密度に比例するので、$f(x,y)$ が分かれば物体の内部構造が分かると考えられます。次に、この物体に、原点からの（符号付き）距離が $t$ で $x$ 軸からのなす角が $\theta$ の方向の法線を持つ直線 $\ell(t,\theta)$ に沿って、$X$ 線を照射するとします（図２）。光源でのX線の強度が $I_0$ であるとき、物体を通過して検出器で観測されるX線の強度 $I(t)$ は、

$ \begin{align} { I(t)=I_{0}\exp\{-P_{f}(\ell(t, \theta))\} } \tag{1}\ \end{align} $

と表されます。ここで

$ \begin{align} { P_{f}(\ell(t,\theta))=\int_{-\infty}^{\infty}f(t\cos\theta -s\sin\theta, t\sin\theta+s\cos\theta)ds } \tag{2}\ \end{align} $

です。よって、実際の観測（撮影）結果から $\log \{ \frac{I(t)}{I_0} \}$ の値を求めることで、各直線 $\ell(t,\theta)$ に対する $P_f(\ell(t,\theta))$ の値が求められます。

式（２）は、与えられた $x$ と $y$ の関数 $f(x,y)$ に対し、$t$ と $\theta$ の関数 $P_f(\ell(t,\theta))$ を定める積分変換を表しています。そこでこの $P_f(\ell(t,\theta))$ を $Rf(t,\theta)=P_f(\ell(t,\theta))$ のような記号で表すと、記号 $R$ は、関数 $f$ に関数 $Rf$ を対応させるという変換 $f \mapsto Rf$ を表します。この変換をラドン変換と呼びます。よって、CTを実現するためには、各 $t$ と $\theta$ に対して（つまり各直線 $\ell(t,\theta)$ に対して）関数値 $Rf(t,\theta)$ が与えられたときに関数 $f$ を求めよ、という問題を解くことになります。

まず、式（１）および（２）の導出を示しておきます。直線 $\ell(t,\theta)$ 上の点 $(x,y)$ に対し、図２のように点 $H$ からの（符号付き）距離を $s$ とおくと、

$x=t\cos\theta - s\sin\theta, y=t\sin\theta + s\cos\theta$

となります。点 $(x,y)$ でのX線強度を $I$ とし、これの直線 $\ell(t,\theta)$ に沿った変化を考えます。点 $(x,y)$ から $\ell(t,\theta)$ に沿って微小変位 $\Delta s$ の範囲を通過する際の、X線の強度変化を $\Delta I$ とおくと、
 
$\displaystyle\frac{\Delta I}{I}=-f(x,y)\Delta s=-f(t\cos\theta - s\sin\theta,\ t\sin\theta+s\cos\theta)\Delta s$

となります。これはランベルトの吸収則と呼ばれます。これを $s$ で積分することで式（１）および（２）が得られます。

さて、前述の問題を解くためには、ラドン変換の逆変換$g \mapsto R^{-1}g$ を求めればよいことになります$（R^{-1}$は、数の-１乗ではなく、逆変換を表す記号です）。この逆変換 $R^{-1}$ を表す式はラドンによって1917年に発見されており、やはりある積分を用いて書けることが知られています。

しかし既に述べたように、積分の式をそのまま用いて計算を行うことは殆どの場合に困難です。そこで、式（２）を計算が容易な形で近似することを考えます。ここでは以下に示す図３のように、物体を平面上の等間隔格子で近似し、各格子内ではX線吸収係数の分布 $f(x,y)$ は一定値と見なす近似を行います。式で書けば、各格子の正方形を $B_{ij}$ で表すとき、

$ \begin{align} { f(x,y)\approx\displaystyle\sum_{i,j}f_{ij}\chi_{ij}(x,y) } \tag{3}\ \end{align} $

という近似を考えます（$\approx$は近似を表します）。ここで $\chi_{ij}(x,y)$ は

$\chi_{ij}(x,y)=\left\{\begin{array}{l} 1 ((x,y)\in B_{ij})\\ 0 ((x,y)\not\in B_{ij}) \end{array}\right.$

という関数です。格子の間隔を $0$ に近づけていけば、近似の精度は上がっていきます。
式（３）の右辺を式（２）に代入することで

$ \begin{align} { \eqalign{ Rf(\displaystyle t,\theta) &= P_{f}(\ell(t,\theta))\approx\int_{-\infty}^{\infty}\displaystyle\sum_{i,j}f_{ij}\chi_{ij}(t\cos\theta-s\sin\theta, t\sin\theta+s\cos\theta)ds \\ &= \sum_{i,j}f_{ij}\int_{-\infty}^{\infty}\chi_{ij}(t\cos\theta-s\sin\theta,t\sin\theta+s\cos\theta)ds \\ &= \sum_{i,j}f_{ij}a_{ij}(t,\theta) } } \tag{4}\ \end{align} $

という近似が得られます。ここで $a_{ij}(t,\theta)$ は、正方形 $B_{ij}$ が直線 $\ell(t,\theta)$ から切り取る線分の長さです（図３参照）。
さらに、ここまでは任意の $t$ と $\theta$ の組 $(t,\theta)$ について $Rf(t,\theta)$ の値が観測されると考えましたが、実際には、無限個の $(t,\theta)$ についての観測結果を得ることは不可能です。そこで、有限個の $(t,\theta)$ について観測結果が得られるとして、それらを $(t_k,\theta_k )$　$(k=1,…,N)$ とおきます。これらに対して得られた観測値を $c_k=Rf(t_k,\theta_k)$　$(k=1,…,N)$ とおき、さらに式⑷を等式に置き換えた式を考えれば

$ \begin{align} { c_{\kappa}=\displaystyle \sum_{i,j}f_{ij}a_{ij}(t_{k},\theta_{k})　　(k=1,\ \ldots,N) } \tag{5}\ \end{align} $

となり、関数 $f$ の近似値 $f_{ij}$ に関する連立一次方程式が得られます。これを解けば関数 $f$ が近似的に求められ、CTが実現されます。しかし実際には、方程式（５）がきちんと解けるものになるにはどのような直線 $\ell(t,\theta)$ を選べば十分かという問題や、実際にコンピューターで方程式（５）の解を計算するにはどのような計算法がよいかという問題など、考えるべき問題は数多くあります。また、今回ここで紹介した方法以外にも、例えばフーリエ（Fourier）解析を用いた方法などもあります。詳細については、参考文献[２]の第１章およびその中の参考文献などを参照してもらえればと思います。