トモグラフィー

トモグラフィー（英: tomography）^{[1][2][3][4][5][6][7][8][9][10][11][12][注釈 1][注釈 2]}は、物理探査、医療診断等で用いられる逆解析技術の一つ。日本語訳は、断層映像法または断層影像法である。コンピュータを用いて処理することで、画像を構成する技術はコンピュータ断層撮影と呼ばれる。

その多くは、対象領域を取り囲む形で、走査線（線源と検出器）を配置し、内部の物性（音速、比抵抗、音響インピーダンス、密度など）の分布を調べる技術である。評価したい対象物によって、X線CT、地震波トモグラフィー、海洋音響トモグラフィーなどと呼ばれている。

Figure1: CT撮影(人体)の様子^{[注釈 3]}。

概要

本記事では、トモグラフ像の撮影と、復元について、原理と装置構成を説明する。トモグラフ像の撮影方法には、主に、平行ビーム光学系を用いる方法（図2参照）と、扇形ビーム（ファンビーム）光学系（図3参照）と円錐ビーム(コーンビーム)を用いる方法がある^{[注釈 4][注釈 5]}。

画像再構成アルゴリズム

CT画像再構成法は解析的再構成法、代数的再構成法、統計的再構成法に大別され、逆投影法は解析的再構成法に分類され、逐次近似画像再構成法は代数的再構成法と統計的再構成法に分類される^[1][13]。 これまでCT画像再構成法の主流はフィルタ補正逆投影法（filtered back projection：FBP法）であったが、近年では画像ノイズ低減効果やアーチファクト低減効果が期待される逐次近似画像再構成法（iterative reconstruction：IR法）が増えつつある^[1][13]。

逆投影法（Back projection）

逆投影法では1回の計算で解(再構成像)が求まる^[13]。

逐次近似画像再構成法（Iterative reconstruction）

逐次近似法は最初に初期画像を仮定してこの画像から計算で作成した投影(順投影)と実測投影との整合性を反復計算によって高めていく手法で、反復計算により計算時間を多く必要とするものの、コンピュータの高速化に伴って統計雑音の性質、装置の分解能、被写体の滑らかさなどの事前情報などを式中に組み込める柔軟性や近年発展の著しい圧縮センシング理論を取り入れることにより徐々に増えつつある^[13][14]。

平行ビーム光学系を用いたトモグラフ像の撮影と復元

トモグラフィーの数学的な基礎はラドン変換と、ラドン逆変換である^[1]。ラドン変換は、トモグラフィーの基本原理であるばかりでなく、例えば、ハフ変換等にも応用される^[15][16][17]。応用範囲の広い数学的手法であるが、ここでは、トモグラフィーのモデル化という観点に重きをおいて説明する。

Figure2: 平行ビーム照射光学系によるトモグラフ像撮影原理；
被写体と、透過光との角度を、${\displaystyle \theta }$とするような、平行ビーム照射光学系を考える。この光学系による投影像は、ある種の線積分の結果と見做すことが出来る。前記の線積分は、被写体をビームが貫通した際に生じた減衰量(attenuation)を表している。図中の符号はそれぞれ、以下の通り (1)被写体, (2)平行ビーム光源, (3)スクリーン, (4)透過光, (5)平行ビーム光源、スクリーンの軌道, (6)平行ビーム光源、スクリーンの軌道の中心, (7)透影像(一次元画像;${\displaystyle p(s,\theta )}$をそれぞれあらわす。

関数${\displaystyle \mu (x,y)}$のラドン変換は、以下の式で与えられる。

${\displaystyle p(s,\theta )=-{\int }_{-\infty }^{\infty }\mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )\,dt}$

即ち、 「${\displaystyle \mu (x,y)}$のラドン変換の ${\displaystyle (\theta ,s)}$での値 ${\displaystyle p(s,\theta )}$」は、 「関数${\displaystyle \mu (x,y)}$の直線${\displaystyle {l}_{[\theta ,s]}(t)}$に沿う線積分の値」である。但し、${\displaystyle {l}_{[\theta ,s]}(t)}$は、

${\displaystyle {l}_{[\theta ,s]}(t)={\begin{bmatrix}s\cos \theta -t\sin \theta \\s\sin \theta +t\cos \theta \\\end{bmatrix}}}$

で定まる直線（tについての直線）である。ここで、上式をtについての直線とみなす際には、 ${\displaystyle \theta }$,sは、固定されているものと考えるが、その際、${\displaystyle \theta }$は前記の直線の傾き角を表し、sは、前記の曲線と原点との間の距離を表すことに注意されたい。

本節では、座標(x,y)における被写体の吸収係数を${\displaystyle \mu (x,y)}$とおく。そのうえで、

1. 吸収係数の位置依存性 ${\displaystyle \mu (x,y)}$ にラドン変換を施すことで、測定結果、即ち透過光によって得られた像 ${\displaystyle p(s,\theta )}$ が得られる(モデル化される)こと
2. 測定結果にラドン逆変換を施すことで、${\displaystyle \mu (x,y)}$が復元されること。

を説明する。

平行ビーム光学系によるトモグラフ像撮影のモデル化

被写体を光線が透過した際に、透過光がどれだけ減衰するかを考えることで、上記のラドン変換が導出される。以下、その導出を行う。

ラドン変換を考える際、光線は幾何光学的な光を考える。即ち、光線は、極めて直進性がよく、吸収はされるが、回折や散乱をしないと考え、さらに反射もしないと考えてよいとする。例えばX線を、人体に透過させる場合には、このように考えて差しさわりない。幾何光学において、光線は直線で表される。光線の軌跡が、x-y断面上の直線${\displaystyle l}$で表される場合について考える。

吸光が、ランベルト・ベールの法則^[1]に従うとすると、 前記光線の入射強度を${\displaystyle {I}_{0}}$、透過後の強度を${\displaystyle I}$表記したとき、

${\displaystyle I=I_{0}\exp \left({-\int \mu (x,y)\,dl}\right)=I_{0}\exp \left({-{\int }_{-\infty }^{\infty }\mu (l(t))\,|{\dot {l}}(t)|dt}\right)}$

が成り立つ。従って、光線lに沿った吸光度を${\displaystyle p_{l}}$と表すと、

${\displaystyle p_{l}=\ln(I/I_{0})=-\int \mu (x,y)\,dl=-{\int }_{-\infty }^{\infty }\mu (l(t))\,|{\dot {l}}(t)|dt}$

次に、x-y平面に対し、角度θをなす光束を考える。この光束の像について考察しよう。新たに${\displaystyle s-t}$座標系を、

${\displaystyle {\begin{bmatrix}s\\t\end{bmatrix}}={\begin{bmatrix}\cos \theta &\sin \theta \\-\sin \theta &\cos \theta \\\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}}$

により定義する。即ち、s-t座標系は、x-y座標系を角度θだけ回転した座標系である。このとき、回転行列の性質から、

${\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}={\begin{bmatrix}\cos \theta &-\sin \theta \\\sin \theta &\cos \theta \\\end{bmatrix}}{\begin{bmatrix}s\\t\end{bmatrix}}={\begin{bmatrix}s\cos \theta -t\sin \theta \\s\sin \theta +t\cos \theta \\\end{bmatrix}}}$

となる。今、上式において、sとθを固定すると、上式は、tを変数とする直線と見做せる。

${\displaystyle {\begin{bmatrix}x\\y\end{bmatrix}}=t{\begin{bmatrix}-\sin \theta \\\cos \theta \\\end{bmatrix}}+{\begin{bmatrix}s\cos \theta \\s\sin \theta \\\end{bmatrix}}}$

のように書くと、より直線らしく見えるであろう。

即ち、x-y平面に対し、角度θをなす光束は、以下の${\displaystyle {l}_{[\theta ,s]}(t)}$で定まる直線を、すべてのsにわたって集めてきたものと考えられる。

${\displaystyle {l}_{[\theta ,s]}(t)=t{\begin{bmatrix}-\sin \theta \\\cos \theta \\\end{bmatrix}}+{\begin{bmatrix}s\cos \theta \\s\sin \theta \\\end{bmatrix}}}$

さらに、光線${\displaystyle {l}_{[\theta ,s]}(t)}$による吸光量を、${\displaystyle p(s,\theta )}$と書くと、

${\displaystyle p(s,\theta )=-{\int }_{-\infty }^{\infty }\mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )\,dt}$

が成り立つ。

以上から、x軸に対し、角度θをなす平行光束による透過像（一次元透過像）のプロファイルが、ラドン変換によって与えられることが判った。この、${\displaystyle p(s,\theta )}$が、CT撮影により測定される測定データである。

平行ビーム光学系によるトモグラフ像の復元

二次元フーリエ変換についての補足

ラドン逆変換とは、実測された${\displaystyle p(s,\theta )}$から、${\displaystyle \mu (x,y)}$を復元する作業のことを指すが、これを説明するためには2変数関数のフーリエ変換について知っておく必要があるので、簡単に復習する。

まず、${\displaystyle \mu (x,y)}$のフーリエ変換とは、

${\displaystyle {\hat {\mu }}(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (x,y)\cdot exp(i(ux+vy))dxdy}$

である。${\displaystyle {\hat {\mu }}(u,v)}$のことを、 ${\displaystyle {\mathcal {F\mu }}(u,v)}$と書く場合もある。 ここで、"${\displaystyle \cdot }$"は、関数と関数の積（単なる掛け算）を意味する。

先の${\displaystyle {\hat {\mu }}(u,v)}$と、${\displaystyle \mu (x,y)}$に対し、以下の等式が成立する。これを、フーリエ逆変換と呼ぶ。

${\displaystyle {\mu }(x,y)={\frac {1}{4{\pi }^{2}}}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\hat {\mu }}(u,v)\cdot exp(i(xu+yv))dudv}$

即ち、関数をフーリエ変換した後、フーリエ逆変換すれば、元の関数に戻る。 ここで、"${\displaystyle \cdot }$"は、関数と関数の積（単なる掛け算）を意味する。

ラドン逆変換

実測された${\displaystyle p(s,\theta )}$から、今度は${\displaystyle \mu (x,y)}$を復元することを考える。この復元操作は、数学的に、 以下の2ステップで行われる。

• ラドン逆変換のステップ（1）

${\displaystyle {\tilde {p}}(r,\theta )={\int }_{s=-\infty }^{s=\infty }p(s,\theta )\exp(isr)ds}$

を計算する。

• ラドン逆変換のステップ（2）

${\displaystyle {\mu }(x,y)={\frac {1}{4{\pi }^{2}}}{\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\tilde {p}}(r,\theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$

を計算する。

上記2ステップの計算を実施することを、ラドン逆変換という。以下、上記ステップにて、μ(x,y)が復元されることを示す。

証明第一段階

（1）${\displaystyle {\tilde {p}}(r,\theta )={\hat {\mu }}(r\cos \theta ,-r\sin \theta )}$ の証明

${\displaystyle p(s,\theta )}$を、変数sについてフーリエ変換（一変数関数としてフーリエ変換）したものを、${\displaystyle {\tilde {p}}(r,\theta )}$とする。即ち、

${\displaystyle {\tilde {p}}(r,\theta )={\int }_{s=-\infty }^{s=\infty }p(s,\theta )\exp(isr)ds}$

とする。これは、上記ステップ（1）の変換に他ならない。 当たり前のことだが、この${\displaystyle {\tilde {p}}(r,\theta )}$は、${\displaystyle {\hat {p}}(u,v)}$とは別物である。そもそも定義が異なる。

今、${\displaystyle p(s,\theta )}$の定義式、即ち、

${\displaystyle p(s,\theta )={\int }_{t=-\infty }^{t=\infty }\mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )\,dt}$

を、上式に代入すると、

${\displaystyle {\tilde {p}}(r,\theta )={\int }_{-\infty }^{\infty }p(s,\theta )\exp(isr)ds}$

${\displaystyle ={\int }_{-\infty }^{\infty }\left\{{\int }_{-\infty }^{\infty }\mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )\,dt\right\}\exp(isr)ds}$

${\displaystyle ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )\exp(isr)\,dtds}$

である。

今、2変数ベクトル値関数${\displaystyle {\varphi }_{\theta }(x,y)}$と、${\displaystyle {\psi }_{\theta }(s,t)}$を、それぞれ、

${\displaystyle {\varphi }_{\theta }(x,y)={\begin{bmatrix}s(x,y)\\t(x,y)\end{bmatrix}}={\begin{bmatrix}\cos \theta &-\sin \theta \\\sin \theta &\cos \theta \\\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}}$

${\displaystyle {\psi }_{\theta }(s,t)={\begin{bmatrix}x(s,t)\\y(s,t)\end{bmatrix}}={\begin{bmatrix}\cos \theta &\sin \theta \\-\sin \theta &\cos \theta \\\end{bmatrix}}{\begin{bmatrix}s\\t\end{bmatrix}}}$

と定めると、明らかに、

${\displaystyle {\psi }_{\theta }({\varphi }_{\theta }(x,y))=(x,y)}$

である。さらに、

${\displaystyle \mu (s\cos \theta -t\sin \theta ,s\sin \theta +t\cos \theta )=\mu ({\varphi }_{\theta }(s,t))}$

である。従って、

${\displaystyle {\tilde {p}}(r,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (({\varphi }_{\theta }(s,t)))\exp(isr)\,dtds}$

である。

さらに、上式を、

${\displaystyle dtds=|J{\varphi }_{\theta }|dxdy=dxdy}$

${\displaystyle \exp(isr)=\exp(ir(x\cos \theta -y\sin \theta ))=\exp(i<r(\cos \theta ,-\sin \theta )|(x,y)>)}$

に注意して積分の変数変換を施すと、

${\displaystyle {\tilde {p}}(r,\theta )}$

${\displaystyle ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu ({\psi }_{\theta }({\varphi }_{\theta })(x,y)))\exp(i<r(\cos \theta ,-\sin \theta )|(x,y)>)|J{\varphi }_{\theta }|\,dxdy}$

${\displaystyle ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (x,y))\exp(i<r(\cos \theta ,-\sin \theta )|(x,y)>)\,dxdy}$

一方で、

${\displaystyle {\hat {\mu }}(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (x,y)\exp(i(ux+vy))dxdy}$

の、(u,v) に　${\displaystyle (r\cos \theta ,-r\sin \theta )}$を代入すると、

${\displaystyle {\hat {\mu }}(r\cos \theta ,-r\sin \theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mu (x,y)\exp(i<r(\cos \theta ,-\sin \theta )|(x,y)>)\,dxdy}$

従って、

${\displaystyle {\tilde {p}}(r,\theta )={\hat {\mu }}(r\cos \theta ,-r\sin \theta )}$

が判る。

証明第二段階

（2）${\displaystyle {\mu }(x,y)={\frac {1}{4{\pi }^{2}}}{\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\tilde {p}}(r,\theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$の証明

第一段階の結論、すなわち、

${\displaystyle {\tilde {p}}(r,\theta )={\hat {\mu }}(r\cos \theta ,-r\sin \theta )}$

より、以下の等式が、任意の(x,y)に対して成り立つ。

${\displaystyle {\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\tilde {p}}(r,\theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$

${\displaystyle ={\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\hat {\mu }}(r\cos \theta ,-r\sin \theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$

上式の右辺に、以下の変数変換:${\displaystyle \xi (r,\theta )={\begin{bmatrix}u(r,\theta )\\v(r,\theta )\end{bmatrix}}=r{\begin{bmatrix}\cos \theta \\-\sin \theta \end{bmatrix}}}$

を施すと、 積分の変数変換の公式から、

${\displaystyle {\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\hat {\mu }}(r\cos \theta ,-r\sin \theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$

${\displaystyle ={\int }_{u=-\infty }^{u=\infty }{\int }_{v=-\infty }^{v=\infty }{\hat {\mu }}(u,v)\exp(i<(u,v)|(x,y)>)dudv}$

一方、二次元のフーリエ逆変換を考えると、

${\displaystyle {\mu }(x,y)={\frac {1}{4{\pi }^{2}}}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\hat {\mu }}(u,v)\exp(i(xu+yv))dudv}$

であるため、

${\displaystyle {\mu }(x,y)={\frac {1}{4{\pi }^{2}}}{\int }_{r=0}^{r=\infty }{\int }_{\theta =0}^{\theta =2\pi }{\tilde {p}}(r,\theta )\exp(i<(r\cos \theta ,-r\sin \theta )|(x,y)>)rdrd\theta }$

を得る。

扇形ビーム光学系によるトモグラフ像の測定と復元

特に医療機器の場合には、図3のような扇形ビーム^[1]を用いる場合が多い。

円錐形ビーム光学系によるトモグラフ像の測定と復元

→詳細は「コーンビームCT」を参照

円錐状のビームをフラットパネルディテクタのような2次元に配置された検出素子で検出する^[1]。1998年にソニーがコーンビームCT の先駆けとなる大視野3次元X線CTの開発に成功した^[18]。その後、コーンビームCTやマルチスライスCTで使用される。

トモグラフィーの種類

• ラジオグラフィー - 放射線を用いる手法。
  • X線トモグラフィー（X線撮影）
  • 中性子トモグラフィー（中性子イメージング） - 中性子を利用して内部構造を可視化する^[19]。
  • ミューオントモグラフィー（ミュオグラフィ） - 1kmの岩盤でも透過する宇宙から降り注ぐミュー粒子を使用して内部構造を調べる手法。火山のマグマの上昇やピラミッドの内部構造の調査^[20]や原子炉の炉心の調査に使用される^[21][22][23]。
• 地震波トモグラフィー - 地震波の伝播時間を用いて、地球内部の3次元速度構造を求める手法。
• 海洋音響トモグラフィー - 海洋中の音波の伝播時間を測定して、音波の伝わった海洋の内部の状態を調べる手法。
• 弾性波トモグラフィ - 弾性波を使用して地質構造を可視化する^[24]。
• 比抵抗トモグラフィ - 地中の電気伝導度を元に地下構造を可視化する^[24]。
• 音響トモグラフィ - 高周波数の音響波を用いて、地盤や構造物の内部構造、岩盤やコンクリート内部の亀裂、樹木内部の空洞を非破壊で可視化する^[25]。
• 光コヒーレンストモグラフィ - 細胞レベルの高速・高分解能で深さ5mm未満の断層画像化に使用される^[26]。
• 電磁誘導トモグラフィ - 電磁波の伝播経路の伝播特性から内部構造を可視化する^[27]。
• 電気インピーダンストモグラフィ - 電気インピーダンスを元に内部構造を可視化する^[28][29]。
• 電気容量トモグラフィ - 電気容量を元に内部構造を可視化する^[29][30]。
• 磁気誘導トモグラフィ - 分析対象内に送信コイルによって誘導電流を生じさせてそれによって発生した磁場を検出することで対象の内部構造を可視化する^[31]。
• 熱音響トモグラフィ - 熱音響トモグラフィは加熱によって生じた音響を検出して内部構造を可視化する^[32]。
• テラヘルツ波トモグラフィ - テラヘルツ波(周波数100GHz～ 10THz、波長30um～3mm)を利用して内部構造を可視化する^[33]。
• 量子トモグラフィ - 未知の量子的対象を完全に同定する手法の総称であり、量子情報処理実験において重要な役割を担っている^[34]。
• 電離層トモグラフィ - 電離層内での電波の伝播を利用して構造を可視化する^[35]。
• 水蒸気トモグラフィ - 多数のGPS観測点でグローバル・ポジショニング・システム(GPS)からの電波を受信する事で大気圏を伝播する電波の水蒸気による衛星視線方向の大気遅延量を利用して大気中の水蒸気の空間分布を算出する^[36][37]。詳細はGPS気象学を参照。