トモグラフィー

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

トモグラフィー: tomography[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [注釈 1] [注釈 2] [注釈 3] [注釈 4] は、物理探査医療診断等で用いられる逆解析技術の一つ。日本語訳は、断層映像法または断層影像法である。

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

Figure1: CT撮影(人体)の様子[注釈 1]

概要[編集]

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

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

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

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

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

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

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

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

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


本節では、被写体の、座標(x,y)における吸収係数を\mu(x,y)とモデル化する。 そのうえで、

  1. 測定結果、即ち透過光によって得られた像のシリーズが、\mu(x,y)にラドン変換を施すことで得られた関数p(s,\theta)として表現される(モデル化される)こと
  2. 測定結果にラドン逆変換を施すことで、\mu(x,y)が復元されること。

を説明する。


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

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

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

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

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に沿った吸光度p_{l}と表すと、

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

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


\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座標系を角度θだけ回転した座標系である。このとき、回転行列の性質から、


\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を変数とする直線と見做せる。


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

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

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

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

が成り立つ。

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


平行ビーム光学系によりトモグラフ像の、復元[編集]

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

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

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

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

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

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

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

即ち、関数をフーリエ変換した後、フーリエ逆変換すれば、元の関数に戻る。 ここで、"\cdot"は、関数と関数の積(たんなる掛け算)を意味する。

ラドン逆変換[編集]

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

  • ラドン逆変換のステップ(1)
\tilde{p}(r,\theta)
={\int}_{s = -\infty}^{s = \infty} p(s,\theta)\exp(isr)ds

を計算する。

  • ラドン逆変換のステップ(2)
 {\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) >)drd\theta

を計算する。

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

証明第一段階[編集]

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

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

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

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

今、p(s,\theta)の定義式、即ち、

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

を、上式に代入すると、

\tilde{p}(r,\theta)
={\int}_{-\infty}^{\infty} p(s,\theta)\exp(isr)ds
={\int}_{-\infty}^{\infty} \{{\int}_{-\infty}^{\infty}\mu(s\cos \theta  -t\sin \theta,s\sin \theta + t\cos \theta )\,dt\} \exp(isr)ds
={\int}_{-\infty}^{\infty} {\int}_{-\infty}^{\infty}\mu(s\cos \theta  -t\sin \theta,s\sin \theta + t\cos \theta ) \exp(isr) \,dt ds

である。

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

{\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}
{\psi }_{\theta}(s,t)=
\begin{bmatrix} s(x,y)\\ t(x,y) \end{bmatrix} = \begin{bmatrix}
\cos \theta & \sin \theta \\
-\sin \theta & \cos \theta \\
\end{bmatrix}
\begin{bmatrix} s \\ t \end{bmatrix}

と定めると、明らかに、

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

である。さらに、


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

である。従って、

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

である。

さらに、上式を、


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

\exp(isr)
=\exp(i r(x\cos \theta + -y\sin \theta) ) 
=\exp(i <r(\cos \theta , -\sin \theta) | (x,y) >)

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

\tilde{p}(r,\theta)
={\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}| \,dx dy
={\int}_{-\infty}^{\infty} {\int}_{-\infty}^{\infty}\mu(x,y)) \exp(i <r(\cos \theta , -\sin \theta) | (x,y) >) \,dx dy


一方で、

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

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

\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) >) \,dx dy

従って、

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

が判る。

証明第二段階[編集]

(2) {\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) >)drd\theta
の証明

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

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

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


{\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) >)\frac{1}{r}drd\theta

=
{\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) >)\frac{1}{r}drd\theta

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

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


{\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) >)\frac{1}{r}drd\theta

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


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

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

であるため、

 {\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) >)drd\theta

を得る。

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

特に医療機器の場合には、図3のような扇形ビームを用いる場合が多い。
Figure3:扇型ビーム照射光学系によるトモグラフ像撮影原理。 ;

脚注・参考文献[編集]

脚注[編集]

  1. ^ a b 図1のような医療用CT撮影機器は、医療機器に該当するため薬事法の規制をうける。従って、それぞれの製品毎に必ず添付文書が必ず存在する(薬事法上の一般名称は、全身用X線CT診断装置)。添付文書は 医薬品医療機器総合機構のデータベース[1]から検索できる。本記事では、東芝メディカルシステムズ Asteion TSX-021Bの添付文書[2]と、日立メディコ製の全身用X線CT診断装置 SCENARIA[3] の添付文書[4]を参考にした。
  2. ^ いくつかの病院のホームページ上で、CTスキャン装置の写真を見ることができる。 一例→[5](麦島内科の例)
  3. ^ 図面の豊富なその他の特許として、例えば、次のようなものがある。 US6879657 [6]、 US6574296 [7] US6775346 [8] US7215734 [9]
  4. ^ A. M. Cormackの2論文は、AIPの重要論文とされている。[10]
  5. ^ プローブとするビームは、(主にX線)、磁場に加え、 電子線(平均自由行程が短いため、真空中に限る)変わり種としては、ミュオン(山の断層写真の撮影例がある)等がある。
  6. ^ 少なくとも医療機器の場合には、扇形ビームをベースとした方法がほとんどである。さらにメーカー独自の改良をが加わっている。最近の技術動向は、特許庁特許事務所等が作成したパテントマップ等からある程度解読可能である。 特許庁作成版としては平成15年版[11] あるいは平成23年版[12]の パテントマップが公開されている。本文では15年版を引用している(23年版は 過去の動向や歴史、基盤技術に関する分析が簡潔するため)本シリーズの前シリーズには、 「技術分野別特許マップ」[13] があり、基本特許の特定などの有益な情報が充実していたが、「技術分野別特許マップ」シリーズに比べ、15年版ですら そういった基礎的な情報は大幅にプアになっていて、産業統計に近い状況になっている。 尚、「技術分野別特許マップ」シリーズでは、CT関連を直接扱ったものはない。

参考文献[編集]

  1. ^ Avinash Kak & Malcolm Slaney (1988), Principles of Computerized Tomographic Imaging, IEEE Press, ISBN 0-87942-198-3. [14]
  2. ^ Excelによる画像再構成入門 (画像再構成シリーズ) [単行本] 篠原 広行 (著) , 坂口 和也 (著) , 橋本 雄幸 (著) 「Excelによる画像再構成入門 (画像再構成シリーズ)」医療科学社 (2007/9/6) [15]
  3. ^ 斎藤 恒雄 「画像処理アルゴリズム (アルゴリズム・シリーズ)」近代科学社(1993/2)
  4. ^ フーリエ解析(13): フーリエ変換の医療分野への応用例 [16]
  5. ^ 梅垣 寿春; 「情報数理の基礎―関数解析的展開 (Information & Computing) 」サイエンス社 (1993/07) ISBN-13: 978-4781907079
  6. ^ 河田 聡 (著), 南 茂夫 (著) ; 「科学計測のための画像データ処理―パソコン EWS活用による画像計測&処理技術 」CQ出版 (1994/04)
  7. ^ Johann Radon, Uber die Bestimmung von Funktionen durch ihre Integralwerte l?ngs gewisser Mannigfaltigkeiten, Computed tomography (Cincinnati, Ohio, 1982) Proc. Sympos. Appl. Math., vol. 27, Amer. Math. Soc., Providence, R.I., 1982, pp. 71?86 (German). MR 692055 (84f:01040)
  8. ^ A. M. Cormack;"Representation of a function by its line integrals, with some radiological applications" J. Appl. Phys. 34, 2722-2727 (1963) [17]
  9. ^ Hounsfield GN; Computerised transverse axial scanning (tomography) I. Description of system. Br J Radiol 46: 1016-1022, 1973.[18]
  10. ^ Hounsfield GNのノーベル賞講演 [19]
  11. ^ Godfrey Newbold Hounsfield US4,115,698 [20](Godfrey Newbold Hounsfieldによる特許)。 同特許のパテントファミリー等は、esp@cenet等より参照可[21]
  12. ^ 藤田保健衛生大学辻岡勝美研究室のページに、有益なリソースが多数ある。[22]
  13. ^ Steven W. Smith;"The Scientist and Engineer's Guide to Digital Signal Processing" [23]
  14. ^ 日立メディコ CTシステム関連技術 [24]
  15. ^ 特許庁総務部技術調査課  「医用画像診断装置に関する特許出願動向報告」 平成15年5月8日 [25]
  16. ^ 国際電気通信基礎技術研究所による特開平05-012438[26]
  17. ^ (新エネルギー・産業技術総合開発機構)即効型地域新生コンソーシアム研究開発 柔軟変形物ハンドリング用ビジョンチップの研究開発報告書 [27]および、立命館大学講義ノート [28]
  18. ^ MATLAB解説記事より[29]
  19. ^ [30]

関連項目[編集]