■ フーリエ解析（19）：　フーリエ変換の医療分野への応用例
　　　　　　　　　　　　 ～ CTスキャンの原理と体験 ～ （JavaScript版）

　フーリエ変換は、複雑な関数を周波数成分に分解してより簡単に記述することができるため、音や光、振動、コンピュータグラフィックス、医療など幅広い分野で用いられています。

　ここでは医療分野への応用の１例として、CTスキャンの原理についておおまかに説明し、それに基づいて作成したプログラムを通してCTスキャンによる断面画像再構成を体験してみましょう。

　CT（Computed Tomography：コンピュータ断層撮影）は、X線などの放射線を用いた測定とコンピュータ処理により、物体の内部画像（断面図）を再構成する技術です。
　最近の画像処理技術の向上により、単に「断面」に限定して用いる検査方法ではなく、3次元グラフィックスとして表示されることも多くなってきています。

　CTにより人体内部の断面図を得るCTスキャンは、短時間でほとんど苦痛なく検査ができ、しかも多くの情報を得ることができることから、超音波検査と並び最もよく用いられる画像検査法のひとつです。

　下図の「元画像」の欄の人体を模した円形（と内部）の上下に、X線源と検出器があります。　CTスキャンではX線源から人体に向けて照射されたX線は人体内部を通過し、一部吸収されて減衰した後、反対側の検出器で受けます。
　人体内部の各位置でのX線吸収率の違いにより、検出値の分布（下図中央の図）が変化します。　そして、X線源と検出器の組を少しずつ回転させて撮影を繰り返します（下図左側「元画像」左右の矢印付近をクリックすると、X線源と検出器の組が回転し、それに応じて中央の分布図が変化する様子を確認できます）。

---

形状の数	分割数	ｽｷｬﾝ範囲	ｽｷｬﾝ刻み
K	N	θmax（度）	Δθ（度）

　このようにして得られた情報をコンピュータ処理することによって、対象物中のそれぞれの位置でのX線の吸収量を、黒から白に至る輝度（明るさ）として表示したものがＣＴ画像です。

　CT画像では、骨の部分などはX線の吸収が大きいので白く見え、空気など吸収が少ない部分は黒く見えます。また、これらの中間の吸収を示す部分（水や筋肉など）は灰色に見えます。

-----------------------------------------
　CTスキャンの数学的な原理は以下のとおりです。
-----------------------------------------

　今、人体の画像（減衰率分布）を f(x, y) とし、ある角度θでの検出値の分布（X軸への投影データ）を R(θ, X)とすると、
　　R(θ, X) = ∫_-∞^∞ f(x, y) dY

となります。　ここで、X-Y はｘｙ座標系をθだけ回転させた座標系です。
　　X = xcosθ + ysinθ
　　Y = -xsinθ + ycosθ
R(θ, X)を f(x, y) のラドン（Radon）変換と呼びます。

　色々な角度θに対するR(θ, X)から f(x, y)を求めることができれば、人体の画像（断層写真）が得られることになります。　このために、以下のようにフーリエ変換を利用します。

　R(θ, X)をフーリエ変換したものをG(θ, r)すると、（dXdY = dxdy）
　　G(θ, r) = ∫_-∞^∞ R(θ, X)e ^-jrX dX

　　　　　　　= ∫_-∞^∞ ｛∫_-∞^∞ f(x, y) dY ｝e ^-jrX dX

　　　　　　　= ∫_-∞^∞∫_-∞^∞ f(x, y) e ^{-jr(xcosθ + ysinθ)} dxdy

　　　　　　　= F (rcosθ, rsinθ)

　ここで、F(u, v) は f(x, y) の２次元フーリエ変換：
　　F(u, v) = ∫_-∞^∞∫_-∞^∞ f (x, y)e^{-j (ux + vy)} dxdy

であり、任意の方向θへの投影データR(θ, X)のフーリエ変換G(θ, r)は、同じ方向に沿ったF(u, v)の値に等しくなることがわかります。

　従って、色々な方向のG(θ, r)からF(u, v)の全体が得られ、これを２次元フーリエ逆変換すれば f(x, y) が求められることになります。

　　f (x, y) = 1/(4π²)∫_-∞^∞∫_-∞^∞ F(u, v)e ^{j (ux + vy)} dudv

　ここで、極座標： u = rcosθ, v = rsinθを導入すると、（dxdy = rdrdθ）
　　f (x, y) = 1/(4π²)∫₀^2π∫₀^∞ F (rcosθ, rsinθ)e ^{jr(xcosθ + ysinθ)} rdrdθ

　　　　　　= 1/(4π²)∫₀^2π∫₀^∞ G(θ, r)e ^{jr(xcosθ + ysinθ)} rdrdθ

　　　　　　= 1/(4π²)∫₀^2π｛∫₀^∞ G(θ, r) r e ^{jr(xcosθ + ysinθ)} dr ｝dθ

　この式は、投影データ（ラドン変換）R(θ, X)のフーリエ変換 G(θ, r)に r を掛けたものをフーリエ逆変換し、それをあらゆる角度θについて求めて積分すれば（平均値を取れば）元画像 f(x, y) が得られることを示しています。

　上のアプリを使って、CTスキャンの処理の様子を体験してみましょう。

・形状の数：K を指定すると、様々な円や長方形が合わせてK個表示されます。
　　実際には内部の様子はこの段階では不明です。　これをCTで検出（再構成）する訳です。

・「計算開始」を押すと、各回転角度でのスキャニング（データのX軸への投影）が行われます。
・計算が終了すると、再構成画像 ～投影データから元画像を復元したもの～ が表示されます。

・「分割数：N」は画像の解像度を指定します。
・「データセット」は形状の数：Kを変化させずに、大きさ・位置を変えます。

・「CT連続スキャン」はX線源と検出器の組を連続的に回転させて、投影データの変化を表示します。
　　左側の図中で矢印部をクリックすると、X線源と検出器の組が左右に回転します。
・「計算／スキャン中止」で、計算あるいは連続スキャンを中止します。

・再構成画面上にマウスがあると、その点での元画像(f0)/再構成画像(f1)の画素値を表示します。

　このアプリの処理手順は概略次のとおりです。
（１）指定された形状の数：K に対して、乱数を利用して位置・大きさの異なる円/長方形を配置します。
　　　（画面左側）

（２）回転角度θを少しずつ変化させて、
　　（２－１）投影データ R(θ, X)を作成し、表示します（ラドン変換、画面中央）。
　　（２－２）これをフーリエ変換（FFT）します　－＞ G(θ, r)。
　　（２－３）G(θ, r)に r を掛けたものをフーリエ逆変換（IFFT）します。

（３）（２－３）の結果を積算し、正規化します。

（４）これを画像（再構成画像）として表示します（画面右側）。