JP4519974B2 - Ｘ線断層画像を再構成する方法及び装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明の分野は、計算機式断層撮影法に関し、より具体的には、Ｘ線医用計算機式断層撮影法に関する。
【０００２】
【従来の技術】
コーン・ビーム投影技術を用いた計算機式断層撮影（ＣＴ）方法では、典型的には、Ｘ線源を用いて円錐（コーン）の形状の頂点を形成する。Ｘ線源の走査軌道は、走査中にコーン・ビームの頂点がそこに沿って移動するような曲線となる。コーン・ビーム源から固定された距離に、１組の検出器が配置されている。現在では、コーン・ビーム投影と、これらの投影の３次元（３Ｄ）ラドン変換の１次導関数との間の関係を確立した方法が存在している。理論的には、各々の線源位置Ｓが、球上に位置するラドン・データ、即ち、Ｓのラドン・シェルを生成することができる。円軌道上の対応する各球は、物体のラドン空間を掃引するものとなり、これをトーラス（円環体）と呼ぶ。物体を包含する最小の球内で且つトーラスの外部にある領域は、シャドウ・ゾーンと呼ばれる。シャドウ・ゾーンは、円軌道からラドン・データが欠落している領域を表わす。
【０００３】
【発明が解決しようとする課題】
円−線軌道を用いる場合には、典型的には、線軌道データを用いてシャドウ・ゾーンを充填する。シャドウ・ゾーンは、円軌道によって生成される有効ラドン・データ全体の体積に比べると極く小さなものであるが、現状でのコーン・ビーム投影法は計算集約的である（すなわち、計算量が膨大になる）。従って、円−線軌道に基づいた画像再構成を簡単化する必要性が存在する。
【０００４】
【課題を解決するための手段】
本発明では、円−線軌道(circle-and-line orbit) に基づいてＸ線断層画像を再構成するシステムが提供される。本システムは、周波数領域において、円−線軌道データから別々に再構成された画像を結合する工程を含んでいる。詳しく述べると、この方法は、（１）シャドウ・ゾーンを周波数領域でシャドウ・コーンへ変換する工程と、（２）３次元フーリエ空間のシャドウ・コーンを１組の２次元フーリエ平面へ写像する工程と、（３）円軌道データのフーリエ変換から、シャドウ・コーンの投影としてマークされた領域内に位置するデータを除去すると共に、線軌道データのフーリエ変換からの同じデータで置き換える工程と、（４）各々の２次元フーリエ平面をそれぞれの２次元画像平面へ変換する工程と、（５）ホーン（horn）形状ボリュームを変換して格子(grid)形ボリュームへ戻す工程とを含んでいる。又、正フーリエ法（Direct Fourier Method、ＤＦＭ）を用いて線軌道データを再構成する補間法も提供される。
【０００５】
【発明の実施の形態】
図１は、本発明を説明する一実施例に全体的に従ったＣＴイメージング・システム１０を示す。システム１０内には、中央処理ユニット（ＣＰＵ）１６と、サンプリング装置１２と、その他の画像支援及び処理サブシステム(１４、１８、２０及び２２)とが含まれている。サンプリング装置１２は、Ｘ線源３４と画像検出アレイ３２とを、起動装置と共に含んでいる。画像支援装置(１４、１８及び２２)は、高速フーリエ・プロセッサ１４と、補間器１８と、メモリ２０と、表示装置２２とを含んでいる。
【０００６】
２次元（２Ｄ）断層像再構成において、フィルタ補正逆投影（ＦＢＰ）の１つの代替技術として正フーリエ法（ＤＦＭ；Direct Fourier Method）逆投影がある。正フーリエ法逆投影は、投影スライス定理に基づいており、この定理は、特定の角度における投影データの１次元（１Ｄ）フーリエ変換（ＦＴ）が同じ角度における物体の２次元ＦＴの断面に対応するものと仮定している。ＤＦＭは、線軌道に関連する従来の逆投影再構成の代替として用いられることになる。
【０００７】
【外１】

【０００８】
【数１】

【０００９】
【外２】

【００１０】
（例えば、Feldkamp等、及びHuによって与えられたようなもの）を用いて評価される。例えば、Feldkamp等は、Journal of the Optical Society of America、第１巻、第６号、１９８４年６月、第６１２頁・第６１９頁に記載された「実用コーン・ビーム・アルゴリズム（PRACTICAL CONE-BEAM ALGORITHM）」と題した論文においてコーン・ビーム・アルゴリズムを発表している。尚、本論文はここに参照されるべきものである。第２の要素は、第３項がラドン空間で評価される場合にのみ必要とされる。ここでは、第３項をフーリエ空間内で評価するので、第２項は無視される。
【００１１】
【外３】

【００１２】
は、円軌道によって回転軸の周りで掃引されるトーラス２４と、物体４０を包含する最小の球面とによって画定されるシャドウ・ゾーン２８を充填するために用いられる（図２）。ラドン空間内で画定されるこのシャドウ・ゾーン２８は、トーラス２４に比較すると小さいが、現在の逆投影法は計算集約的であり、このことは、例えば１９９７年１月３１日、GE社のMedical System事業部のThe Applied Science Laboratory刊のHui Huによる「円軌道並びに円−線軌道のためのコーン・ビーム再構成アルゴリズム（CONE BEAM RECONSTRUCTION ALGORITHM FOR THE CIRCULAR ORBIT AND A CIRCLE-AND-LINE ORBIT）」と題した論文において定義
【００１３】
【外４】

【００１４】
【数２】

【００１５】
ここで、
【００１６】
【数３】

【００１７】
また、ｗ（ｚ₀，Θ，L ）は、２L z₀ cos Θ＋ｚ₀ ² cos² Θ−ｄ² sin² Θ＞０の場合には、ｗ（ｚ₀，Θ，L ）＝１であり、その他の場合には、ｗ（ｚ₀，Θ，L ）＝０である。
また更に、ここで、Ｈ（ｚ₀ ，Θ，L ）は、ｚ＝ｚ₀ において得られるスナップショットから算出されるラドン・データの加重付き２次導関数である。物体空間
【００１８】
【外５】

【００１９】
ている。原点（０，０）及び（Ｙ₀，Ｚ₀）を直径の２端として用いて検出器平面上で円を描くと、円が切り取るすべての有効ラドン・データＨ（ｚ₀ ，Θ，L ）は、制限 L＝Ｙ₀ sin Θ＋Ｚ₀ cos Θを満たす。各々の点（Ｙ₀，Ｚ₀）毎に、加算された後に物体空間に逆投影されるべき対応するラドン・データの集合Ｈ（ｚ₀ ，Θ，L ）が存在する。従って、この処理の計算の複雑さは、Ｏ（Ｎ⁵ ）に近付く。
【００２０】
本発明は、物体４０（図２に示す）の線軌道走査から得られたデータを用いて計算効率のよいＣＴ再構成を決定する方法を記載する。第１に、ラドン空間におけるシャドウ・ゾーン２８を拡大して、周波数空間におけるシャドウ・コーン３０へ変換する（図３）。第２に、ＣＰＵ１６によって実行されるパッチ処理(patching)を、元の３次元周波数空間から多数の２次元周波数スライスへ単純化する。
【００２１】
第１段階は、３次元ラドン空間での半径線のフーリエ変換は、３次元フーリエ空間内の同じ半径線と等価であるという事実に基づく。３次元ラドン空間での半径線がシャドウ・ゾーン２８を切り取るならば、３次元フーリエ空間内の対応する半径線は、ラドン半径線における幾分かのラドン・データが欠落するために「汚染」される。シャドウ・ゾーン２８を切り取るすべての半径線がシャドウ・コーン３０を形成するので、３次元フーリエ空間における対応する半径線の集合も又、周波数欠落コーンを形成する。フーリエ空間で定義される式（１）の第１項の周波数欠落コーン（又は汚染されたコーン）は、式（１）の第３項のコーンによって置き換えられる。
【００２２】
第２段階は、線軌道が、（円軌道に起因する）汚染されたコーンを置き換えるのに十分な周波数データを支持することを示すものであり、この後に、２次元フーリエ空間での汚染されたコーンを置き換える手順が続く。
【００２３】
図６は、Ｘ線源３４及び検出器パネル３２を静止したままにして、物体４０をｚ軸（この軸は、円軌道の回転軸でもある）に沿って移動させることにより、線走査軌道を実現し得ることを示している。物体４０は、「ｚステップ」３６の増分に示す離散的な増分でｚ方向に移動するものとして示されている。複数の検出器３２から成る各々の垂直な列は、線源３４と共に、物体４０を通る平面を画定する。各々の垂直な平面内で、各々の検出器位置が特定の角度を画定し、各検出器によって収集されたデータは、ｚステップ３６の間隔でサンプリングされた物体の平行投影に対応する。図７は、角度θ_m ＝θ_max における物体４０のいくつかの平行投影を示しており、図８は、角度θ_m ＝θ_min における物体４０の平行投影を示している。これらの平行投影を用いて、２次元ＦＴ平面における対応する欠落した角セグメント(angular segment) を充填する。θ_max 及びθ_min は、シャドウ・ゾーン２８のデータ・ボリュームを置き換えるのに用いられるシャドウ・コーン３０（図３）のデータ・ボリュームが、シャドウ・ゾーンのデータ・ボリュームを包含するように選択される角度である。この条件は、図2に示すように、θ₀ ＝θ_m のときに満たされる。ここで、
θ₀＝ｓｉｎ^-1（Ｒ／ｄ）
であり、ここで、Ｒは物体４０を包含する最小の球の半径であり、ｄは円軌道の半径である。
【００２４】
３次元フーリエ空間内のシャドウ・コーン３０（図３）は、２次元シャドウ・コーンとして２次元フーリエ平面の集合へ投影することができる。これらの２次元フーリエ平面の各々の逆変換は、Ｘ線源及び検出器３２の列Ｎによって画定される画像平面に対応する（図６）。
【００２５】
先ず、Ｆ（ｕ，ｖ）を２次元関数ｆ（ｘ，ｙ）のフーリエ変換であるとする。ｘ＝ｘ₀ 上の１次元フーリエ・データは、次のようにして導かれる。
【００２６】
【数４】

【００２７】
同様に、３次元関数ｆ（ｘ，ｙ，ｚ）のフーリエ変換について、ｘ＝ｘ₀ における物体スライスの２次元フーリエ変換は、次のように表わされる。
【００２８】
【数５】

【００２９】
所望のｎ次元フーリエ空間が、拡大された（ｎ＋１）次元フーリエ空間の（ｎ＋１）次元と直交していない場合について述べる。例えば、所望の１次元フーリエ変換が、式ｘ cos θ＋ｙ sin θ−ｐ＝０によって定義される線上に位置しているならば、元の２次元座標（ｘ，ｙ）を新たな座標（ｘ′，ｙ′）へ角度θだけ回転させると適当である。Ｆ（ｕ′，ｖ′）をｆ（ｘ′，ｙ′）のフーリエ変換であるとする。すると、Ｆ（ｕ′，ｖ′）も又、Ｆ（ｕ，ｖ）を角度θだけ回転させることにより得られる。線ｘ cos θ＋ｙ sin θ−ｐ＝０に沿ったフーリエ変換は、次のように表わされる。
【００３０】
【数６】

【００３１】
同様に、ｘ cos ψ＋ｙ sin ψ−ｐ＝０によって画定される平面に沿った２次元フーリエ変換は、次の式によって与えられる。
【００３２】
【数７】

【００３３】
また、Ｆ（ｕ′，ｖ′，ｗ）は図４に示すように、ｗ座標を不変に保ったまま（ｕ，ｖ）座標を角度ψだけ回転させることにより得られる。回転後の座標系を図４及び図５に図形的に示しており、ここでは、新たな座標（ｕ′，ｖ′，ｗ）
【００３４】
【外６】

【００３５】
要約すると、３次元フーリエ空間においてラドン空間の面積を画定するシャドウ・コーン３０（図４）は、１次元逆フーリエ変換を適用することにより２次元フーリエ部分空間へ投影される。換言すると、３次元フーリエ空間のシャドウ・コーンを置き換えるのではなく、走査されている物体４０を一連の２次元フーリエ平面内に表わして、すべての個々の平面における投影後の２次元シャドウ・ディスク２６に置き換える。
【００３６】
以上の原理を適用するに当たって、次のアレイ及びスキャナのパラメータを仮定する。検出器アレイ３２（図６）は、Ｄ（Ｍ×Ｎ）という表現によって表わし、ここで、検出器アレイ３２はＭ行及びＮ列を有するものとする。円軌道の投影データをＰ_c （Ｍ×Ｎ×Ｖ_c ）という表現によって表わし、ここで、Ｖ_c は、ビューの総数である。線軌道の投影データは、Ｐ_L （Ｍ×Ｎ×Ｖ_L ）という表現によって表わし、ここでも又、Ｖ_L は、ビューの総数である。円軌道データから再構成される画像は、Ｉ_c （Ｋ×Ｋ×Ｋ）又はＩ_cs（Ｋ×Ｎ×Ｋ）という形態を有する。前者は立方ボリュームであり、後者はＸ線源及び検出器パネルによって画定される「ホーン形状」のボリュームである。線軌道データから再構成される画像は、Ｉ_L （Ｋ×Ｋ×Ｋ）又はＩ_Ls（Ｋ×Ｎ×Ｋ）という形態を有し、同様に、前者は立方ボリュームであり、後者は「ホーン形状」のボリュームである。
【００３７】
画像は、線軌道データから再構成される。第１段階として、線軌道投影データを修正して、図６、図７及び図８に図形的に示すように、複数の射線が線源３４から検出器アレイ３２へ走行する際の発散の効果を収容するようにする。検出器パネル上で検出される発散した射線は、検出器の値に次のような加重関数を乗算することにより修正される。
【００３８】
Ｗ（ｍ，ｎ）＝Ｄ／［Ｄ²＋（ｍ−ｍ₀）²ａ²＋（ｎ−ｎ₀）²ｂ²］^1/2
【００３９】
【外７】

【００４０】
検出器パネル３２の（ｍ₀ ，ｎ₀ ）へ投影されている。値ａ及びｂは、それぞれ行次元及び列次元での検出器ピッチである。検出器ピッチは、２つの隣接した検出器の中心と中心との間の距離を定義する標準的な定義である。
【００４１】
次に、ｎ番目の傾斜(angular) 垂直平面を再構成する。ｎ番目の傾斜垂直平面を再構成するためには、Ｐ_L 配列のｎ番目の列から投影データＰ_Ln（Ｍ×Ｖ_L ）を抽出する。ここで、ｍ番目の行は、固定されたビーム角度θ_m に対応している。
【００４２】
Ｐ_m の行要素は、この垂直スライスの物体中心が、ｖ＝Ｖ₁ ／２に再整列するようにシフトされる。Ｐ_Ln′を再整列後の投影配列であるとする。
【００４３】
物体空間に復帰するためには、２つの方法のうちいずれか一方を採用することができる。第１の選択肢としては、１次元ＦＦＴをＰ_Ln′の各々の行に適用する。そして、ＦＦＴデータの各々の集合を、ｍ＝１，．．．，Ｍとした場合の角度θ_m における半径線上に配置する。そして、逆ＦＦＴを用いて、フーリエ・データ集合を物体空間へ変換する。
【００４４】
第２の選択肢としては、すべてのＮ列についてのフィルタ補正後のＰ_Ln′を逆投影して、Ｘ線源３４及びＮ列の検出器列のＮによって画定されるＮ個の垂直平面を含むホーン形状ボリュームを得る。ｆ_n （ｘ，ｙ）をｎ番目の列に基づく再構成後の画像であるとし、Ｑ_n をフィルタ補正後の投影配列であるとすると、
【００４５】
【数８】

【００４６】
となる。
【００４７】
上の式をシフト後の行要素に適用した結果は、ホーン形状を有する立体配列Ｉ_Lsとして表わされる。
【００４８】
最終画像は、２つの可能な方法のうち一方によって再構成することができる。先ず、シャドウ・コーン３０を３次元フーリエ空間内で置き換えることができる。代替的には、シャドウ・コーン３０を２次元フーリエ空間内で置き換えた後に、ホーン形状ボリュームから立方ボリュームへ戻すように補間してもよい。
【００４９】
最初に、３次元空間でのシャドウ・コーン３０の置き換えについて述べる。第１段階として、Feldkampのアルゴリズムを用いて、円軌道データから、トーラス２４に関連した画像Ｉ_c を再構成する。次いで、Ｉ_c から３次元フーリエ変換Ｆ（Ｉ_c ）を得る。
【００５０】
ホーン形状ボリュームＩ_Lsを補間して、対応する立方ボリュームＩ_L を得る。Ｉ_L から３次元フーリエ変換Ｆ（Ｉ_L ）が得られる。
【００５１】
Ｆ（Ｉ_c ）のシャドウ・コーン部分は、Ｆ（Ｉ_L ）の対応部分で置き換えられ、修復されたフーリエ・データ配列を与える。次いで、修復後のフーリエ・データ配列に対して逆フーリエ変換を適用して、立方空間に復帰する。
【００５２】
次に、シャドウ・コーンを２次元フーリエ空間内で置き換える方法について述べる。第１段階として、この場合にもやはり、Feldkampのアルゴリズムを用いて、円軌道データから画像Ｉ_cを再構成する。次に、Ｉ_c を補間してＩ_csとする。ここで、Ｉ_s は、複数の傾斜垂直平面を画定する。Ｉ_Lsは、上述のように、線軌道データから得られる。
【００５３】
フーリエ変換Ｆ（Ｉ_cs）及びＦ（Ｉ_Ls）は、Ｉ_cs及びＩ_Lsから得られる。シャドウ・コーンは、Ｆ（Ｉ_cs）の各々の２次元平面内でＦ（Ｉ_Ls）において定義される対応データによって置換される。修復後のフーリエ・データの２次元集合の各々に対して逆フーリエ変換を適用する。最後に、ホーン形状ボリュームを補間して、立方ボリュームに戻す。
【００５４】
以上に述べたＤＦＭ法は、極座標から直交座標への補間を要求する。
【００５５】
上述のコーン・ビーム再構成は、円−線の軌道に基づいており、フーリエ領域の円軌道データの欠落したシャドウ・ゾーン２８を線軌道データによって置き換えることを必要とする。図６に示す線軌道配向の場合には、物体４０は、ｚステップ３６の増分として示す離散的な増分でｚ方向に移動するものとして示されている。複数の検出器３２から成る各々の垂直列は、線源３４と共に、物体４０を通る平面を画定する。各々の垂直な平面内で、各々の検出器位置が特定の角度を画定し、各検出器によって収集されたデータは、ｚステップの間隔でサンプリングされた物体の平行投影に対応する（図６、図７及び図８）。図７は、角度θ_m ＝θ_max における物体４０のいくつかの平行投影を示しており、図８は、角度θ_m ＝θ_min における物体４０の平行投影を示している。これらの平行投影を用いて、２次元ＦＴ平面における対応する欠落した角セグメントを充填する。θ_max 及びθ_min は、シャドウ・ゾーン２８のデータボリュームを置き換えるのに用いられるシャドウ・コーン３０（図３）のデータボリュームが、シャドウ・ゾーンのデータボリュームを包含するように選択される角度である。この条件は、図２に示すように、θ₀ ＝θ_m のときに満たされる。
【００５６】
１つのアプローチは、限定された角度ビューの集合から、平行投影を用いて、各々の傾斜平面毎にフィルタ補正逆投影（ＦＢＰ）画像を再構成し、次いで、２次元ＦＴを算出し、これを用いて欠落した角セグメントを置き換えるものである。
【００５７】
もう１つの代替的な方法は、極格子データを平行投影の１次元ＦＴから２次元ＦＴ領域内の直交格子へ直接変換するものである。この方法は、正フーリエ法逆投影に要求される補間に類似している。
【００５８】
線軌道では、各々の平面毎に、所定の数の投影（各々の垂直列におけるＭ個の検出器に等しい）を、θ_min からθ_max の角度範囲について取得する。線データのための補間手法は、コーン・ビーム再構成方法の計算に有用な以下に述べる補間法を展開するための方法論を提供する。
【００５９】
【外８】

【００６０】
【数９】

【００６１】
となる。
【００６２】
【外９】

【００６３】
【数１０】

【００６４】
【外１０】

【００６５】
【数１１】

【００６６】
となり、ここで、
ρ＞０、０≦φ＜２πである。
【００６７】
Ｍ（ρ，φ）が既知であると仮定すると、極座標でのこのＦＴデータは、表現Ｍ（ρ cosφ，ρ sinφ）＝Ｍ（ρ，φ）によってデカルト座標へ変換することができる。物体は、次の表現を用いて２次元逆ＦＴによって再構成される。
【００６８】
【数１２】

【００６９】
しかしながら、Ｍ（ρ，φ）は、すべての位置において既知であるわけではなく、離散的な点（ρ_j ，φ_k ）の有限な集合においてのみ既知であるに過ぎない。そこで、極座標点の既知の値から直交デカルト格子の全体に要求される値へ補間する問題となる。
【００７０】
この問題を解くために既存の定理が利用される。先ず、ｘ（ｔ）を、時間空間ｔにおいて周期Ｔを有する周期関数であるとする。ｃが定数であり、Ｋがｘ（ｔ）の標本数に関連した整数である場合に、ｘ（ｔ）が値［ｃ（１／Ｔ）］に角度に関して帯域限定されるならば、ｘ（ｔ）を次のように書くことができる。
【００７１】
【数１３】

【００７２】
ここで、Ｔ_s ＝Ｔ／（２Ｋ＋１）である。
【００７３】
この証明は、シャノンのサンプリング定理から、次のようにして導かれる。
【００７４】
【数１４】

【００７５】
ここで、データ配列の行数は、Ｎ＝２Ｋ＋１である。
【００７６】
角度補間は、関数Ｍ（ρ，φ）を用い、ここで、変数φは、２πの周期を有する。Ｍ（ρ，φ）は、角度に関してＫに帯域限定されている（即ち、フーリエ級数の係数が、｜ｎ｜＞Ｋについてｃ_n ＝０を満たす）ので、
【００７７】
【数１５】

【００７８】
となり、ここで、
σ（φ）＝［sin｛((2K＋1)/2)φ｝］／［（2K＋1）sin((1/2)φ)］（１５）
である。
【００７９】
次に、半径補間について記載する。Ａが、物体４０を包含する最小の球面の半径であるとして前に示した半径Ｒと等価である場合に、物体が直径２Ａに空間限定されているならば、物体の投影も又、直径２Ａに限定され、従って、シャノンのサンプリング定理を用いると、投影Ｐφ（ｕ）のＦＴを次の式として再構成することができる。
【００８０】
【数１６】

【００８１】
角度補間と半径補間とを組み合わせると、次の式が得られる。
【００８２】
【数１７】

【００８３】
Ｍ_c （ｕ_l ，ｖ_m ）は、デカルト空間全体にわたって次の式として算出される。
【００８４】
Ｍ_c（ｕ_l，ｖ_m）
＝Ｍ（ｕ_m ²＋ｖ_m ²）^1/2，cos ^-1［ｕ_m／（ｕ_m ²＋ｖ_m ²）^1/2］
補間誤差を減少させるために、以下に述べるようないくつかの段階を踏むことができる。先ず、空間内での不十分なサンプリングからの誤差を考察する。
【００８５】
線軌道の最中に、予め定められたｚステップ３６の間隔で平行投影データがサンプリングされる。投影は空間限定されるので、帯域限定されない。所定のｚステップ・ウィンドウ（例えば、２Ａ）にわたって投影をサンプリングすると、帯域幅を効率的に縮小すると共に、２次元ＦＴ領域にエイリアシング（折り返し）アーティファクトを導入することになる。結果として、ｚステップ３６は、実効帯域幅Ｂ＝１／（２＊ｚステップ）［点／サイクル］が投影の周波数内容を実質的に包含するように十分小さいものに選択される。
【００８６】
誤差のもう１つの原因は、半径周波数における不十分なサンプリングである。物体が、直径について１／（２＊Ａ）に空間限定されているならば、物体の投影も又、２Ａに空間限定され、このため、エイリアシングを回避するためにＦＴは少なくとも１／（２Ａ）［点／サイクル］の間隔で一様にサンプリングされる。
【００８７】
線データの場合には、実効帯域幅はＢと定義される。半径標本の数Ｎは、（２＊Ｂ）／Ｎ＝１／（Ｎ＊ｚステップ）≦１／（２＊Ａ）となるように選択される。又、２次元ＦＴ領域における格子間隔も、最大で１／２Ａ［点／サイクル］である。
【００８８】
誤差のもう１つの原因は、ビューの不十分な数である。非等方的で且つ中心がずれた物体は、角度に関して帯域限定されない。このため、投影の数を非常に大きくすることができる。しかしながら、実用上は、投影の数は、２次元ＦＴ領域において２つの隣接した半径線の間のサンプリングが、最大格子間隔１／２Ａ［点／サイクル］に近似的に等しくなるように選択される。
【００８９】
次に、２次元ＦＴ領域内での補間について述べる。線軌道データの場合には、２次元ＦＴ領域内で直交格子間隔１／２Ａ［点／サイクル］を用いる。コンピュータでの実現時には、典型的には、高速フーリエ変換（ＦＦＴ）を用いて各々の投影の１次元ＦＴを算出する。半径サンプリングを改善するために、ＦＦＴを適用する前に、投影データのゼロ充填を行う。本明細書では、ゼロ充填は、データ集合に対してゼロを付加して、ＦＦＴが行われた後のアーティファクトを最小化する工程を含んでいる。一旦、半径方向サンプリングが得られたら、最近接補間又は線形補間を用いて、新たな標本点の間の値を算出する。
【００９０】
前述のように、線軌道の場合には、小さな角セグメントについて複数の投影が要求される（即ち、角度方向は、半径サンプリングに比べてオーバサンプリングされる。）。１つの問題は、このデータを効率的に利用する仕方にある。最近接法のような単純な補間法は、利用可能なデータのすべてを利用しない。１つの代替的なアプローチは、行をグループ化し、次いで、以下に詳細に記載するような逆双１次補間法を利用して信号対ノイズ比（ＳＮＲ）を増大させるように各々のグループからの行を加算することにより、ビューの数をデシメート（decimate）するものである。
【００９１】
本発明では、逆双１次補間法を極格子フォーマットのデータに対して利用する。方位角方向での最適な補間のためには、式（１５）で定義されたσ（φ）関数を用いることができる。但し、この関数は、計算面では実用的でない。コンピュータでの実現のためには、補間関数における値をトランケート（切捨て）せねばならず、これにより、アーティファクトが生ずる。
【００９２】
コンピュータ計算の観点から、実現が最も簡単な補間アルゴリズムは、最近接（ＮＮ）アルゴリズムであり、このアルゴリズムでは、各々のピクセルに対し、各ピクセルに最も近接した標本の値を与える。ＮＮは、サンプリングされたデータを矩形関数によって畳み込むことと等価である。空間領域での矩形関数による畳み込みは、周波数領域において信号にｓｉｎｃ関数を乗算することと等価であるが、ｓｉｎｃ関数は顕著なサイドローブを有するため、不十分な低域通過フィルタである。
【００９３】
線形補間は、サンプリングされたデータの三角形ウィンドウによる畳み込みに相当する。この関数は、周波数領域においてＮＮよりも良好な低域通過フィルタ特性を有するｓｉｎｃ² 関数に対応する。
【００９４】
ＮＮ法は、単一の点を基盤として補間する。線形補間アルゴリズムは、２つの最も近接した点を基盤として補間する。補間に３つの点を用いると、１辺では２つの点が得られ、他辺では唯一の点が得られる。Ｂスプライン法は、各々の方向において２つの最も近接した点にわたる補間を可能にする。３次Ｂスプライン法は、ＮＮ法の矩形関数の４つの畳み込みを含む。３次Ｂスプラインは対称であるので、区間（０，２）において定義するだけでよい。数学的には、これらのＢスプラインは、次のように書くことができる。
【００９５】
ｆ（ｘ）＝ｘ³−ｘ²＋（４／６） ｘ∈（０，１）
ｆ（ｘ）＝−（ｘ³／６）＋ｘ²−２ｘ＋（８／６） ｘ∈（１，２）
上で定義した３次Ｂスプライン方程式は、シャドウ・コーン３０のデータの利用可能な集合の全体にわたって滑らかな３次多項式にフィットするように適合させることができる。
【００９６】
もう１つの補間関数は、σ（φ）のトランケートに関連している。前述したように、σ（φ）関数の角度方向における理想的な補間関数を見出すことが望ましい。しかしながら、σ（φ）関数は、かなりのエネルギを有しており、このため、長い距離にわたって急速に減衰せず、従って、空間領域での畳み込みとして容易に実現することができない。短い距離にわたってこの関数をトランケートすると自然であるが、トランケートは、上述のエネルギの殆どを無視する。更に、空間領域でのトランケートは、物体が空間領域から信号領域へ変換されるときに、リンギング・アーティファクトを生ずる。
【００９７】
双１次補間は、２次元の補間方法である。この補間法においては、ｆ（ｘ，ｙ）は、ｎ₁ Ｔ≦ｘ≦（ｎ₁ ＋１）Ｔ及びｎ₂ Ｔ≦ｙ≦（ｎ₂ ＋１）Ｔについて４つの最も近接した点におけるｆ（ｎ₁ ，ｎ₂ ）の線形結合によって評価される。双１次法における補間後の値ｆ（ｘ，ｙ）は、

となる。ここで、
Δ_x＝（ｘ−ｎ₁Ｔ）／Ｔ、及びΔ_y＝（ｙ−ｎ₂Ｔ）／Ｔ
である。
【００９８】
逆双１次補間は、逆双１次加重を介して極格子データを直交格子データへ変換する。即ち、各々の極格子要素は、４つの部分面積を正規化された加重として用いて各要素の４つの最も近接した直交格子に再分配される。この方法は、原データが各々の矩形方形に少なくとも１つの要素を有している限りにおいて、原データが規則的に分配されなくてもよいという意味で、ロバストなものである。このため、検出器アレイ３２（図６）の欠陥又は継ぎ目によって生ずる望ましくないデータ要素を無視することができる。この処理のもう１つの利点は、すべての良好なデータが完全に利用されることであり、この結果、よりよい信号対ノイズ比（ＳＮＲ）を提供する。
【００９９】
線軌道データ値ｘ、及び逆双１次補間に用いられる関連する格子点を図９に図形的に示す。データ値ｘは、次のようにその４つの最も近接した格子点に分配される。
【０１００】
ｇ_i,j＝ｇ_i,j＋ａ₁₁ｘ ｇ_i+1,j＝ｇ_i+1,j＋ａ₀₁ｘ
ｇ_i,j+1＝ｇ_i,j+1＋ａ₁₀ｘ ｇ_i+1,j+1＝ｇ_i+1,j+1＋ａ₀₀ｘ
各々の格子点は、積算された値ｇ_i,j を有しており、４つの隣接した正方形からの積算された加重後のｘの値を表わす。各々の格子点において積算された値を正規化するために、もう１つの変数ｗ_ijを保持する必要があり、この変数は、ピクセル（ｉ，ｊ）の双１次分布に関連した積算された加重（ａ_ij）で表わされる。実際には、この工程は、各々の格子点が積算された加重ｗ_ijで除算されるので、補間手順の終盤での正規化工程となる。
【０１０１】
２次元ＦＴ領域及び画像領域について、所載の補間法を用いて６組のシミュレーションを行った。シミュレーションの結果は、最近接（ＮＮ）補間法、線形補間法、３次Ｂスプライン補間法及びフィルタ補正逆投影（ＦＢＰ）補間法と比較した場合の逆双１次補間を示す。結果を表Ｉ・表VIに掲げる。
【０１０２】
最初の３つのシミュレーションについては、線形補間を半径方向に用いた。最近接補間法、線形補間法及び３次Ｂスプライン補間法は、方位角方向に用いた。２次元ＦＴ領域での角セグメントの再構成は、線軌道データを用いてシミュレートされた。関心のあるスライスは、図６の物体を通る中央のスライスとした。物体は３つの水平な円板から成り、ｚ（垂直）方向に尖端を有するようにした。
【０１０３】
又、トランケートσ（φ）補間法を方位角方向（８つの近接点を含めて）に用いた。線形補間は、半径方向に用いた。最近接補間及び線形補間は、３次Ｂスプライン補間よりも良好な性能の結果を生じた。しかしながら、３次Ｂスプラインは、トランケートσ（φ）補間における８つの点に対して４つの点しか用いない。従って、３次Ｂスプラインは、８点のトランケートσ（φ）補間よりも計算面でより好ましく、このため、トランケートσ（φ）補間は比較に含めなかった。
【０１０４】
極座標点に基づく双１次補間についても評価を行った。関心のあるスライスの２次元ＦＴを算出し、上述の方法を比較する基本的な基準として用いた。ＦＢＰ再構成の２次元ＦＴも又、比較に含めた。ＦＢＰ再構成は、２次元ＦＭＰ画像の再構成及びこの後に２次元ＦＴの算出を必要とする。しかしながら、極座標双１次補間は、画像領域へ戻す変換を必要としない。これにより、ＦＢＰ法を凌ぐ計算面の利点が得られる。
【０１０５】
各々の方法の性能を測定するために、様々な数のビュー及び半径標本について上述の方法を用いて補間後の直交格子値を算出した。誤差は、次のような実数２次元ＦＴによって算出した。則ち、ノルム１（全絶対差；表１）、ノルム２（全平方距離の平方根；表２）、及びノルム∞（最大絶対差；表３）である。尚、ノルム１、ノルム２及びノルム∞という用語は、統計分析において用いられる標準的な用語であることに留意されたい。
【０１０６】
【表１】

【０１０７】
【表２】

【０１０８】
【表３】

【０１０９】
又、補間された直交格子から、逆２次元ＦＴを用いて、対応する限定されたビュー画像も算出した（表４乃至表６）。誤差は、再構成画像と真の画像とで比較した。ＦＢＰ画像も比較に含めた。尚、線軌道の実際の実現時には、２次元ＦＴデータにパッチを当ててFeldkampの２次元ＦＴデータとする前に画像を再構成する必要はないことに留意されたい。
【０１１０】
【表４】

【０１１１】
【表５】

【０１１２】
【表６】

【０１１３】
半径サンプリングを増大させても、画像領域の誤差は減少しないことが観察されよう。しかしながら、半径サンプリングを減少させると、直交格子間隔が半径サンプリングと最早同程度でなくなるので、誤差が導入される。但し、ｚステップ３６の値は、帯域幅制限の故に固定されている。従って、各々の投影毎に、Ｎにｚステップを乗じたものが２ＡとなるようなＮ個のサンプルが利用可能となる。
【０１１４】
各々の例で、逆双１次補間は、２次元ＦＴ領域においても、画像領域においても、所載の他の補間法と比較して誤差が少ない。又、例えば、表VIのデータによって観察され得るように、画像誤差を大幅に増大させずに、２次元ＦＴ領域及び画像領域におけるビューの数を２５６から減少させることもできる。
【０１１５】
本発明に従ってＣＴ画像を再構成する様々な方法及び装置の特定の実施例を、本発明を実行し利用する方式を説明する目的で記載した。本発明の他の変形及び改変の実現並びにその様々な側面は当業者には明らかであり、本発明はここに記載した特定の実施例によって限定されないことを理解されたい。従って、本発明はここに開示し請求した基礎にある基本的な原理の要旨の範囲内にある任意の及びすべての改変、変形又は均等構成を網羅するものとする。
【図面の簡単な説明】
【図１】本発明のＸ線断層撮影イメージング・システムを示すブロック図である。
【図２】図１のシステムによって形成される円軌道の断面を示す模式図である。
【図３】図２のシャドウ・ゾーン及び対応するシャドウ・コーンを示す模式図である。
【図４】平面へ投影された図２のシャドウ・コーンの断面を示す模式図である。
【図５】平面へ投影された図２のシャドウ・コーンの断面を示す模式図である。
【図６】図１のシステムによってｚ方向に実行される漸増式処理工程を示す模式図である。
【図７】角度θ_max における平行投影データの集合（図６に示す検出器アレイのＭ番目の行から得られる）を示す模式図である。
【図８】角度θ_min における平行投影データの集合（図６に示す検出器アレイの１番目の行から得られる）を示す模式図である。
【図９】図１のシステムの逆双１次補間を示す模式図である。
【符号の説明】
１０ ＣＴイメージング・システム
１２ サンプリング装置
１４ 高速フーリエ変換プロセッサ
１６ 中央処理ユニット（ＣＰＵ）
１８ 補間器
２０ メモリ
２２ 表示装置
２４ トーラス
２６ シャドウ・ディスク
２７ 水平な平面
２８ シャドウ・ゾーン
３０ シャドウ・コーン
３２ 画像検出器アレイ
３４ Ｘ線源
３６ ｚ軸
４０ 物体

Claims (19)

1. 円軌道からのデータと線軌道からのデータとを含むＸ線断層画像を再構成する方法であって、前記円軌道が、３次元ラドン空間において、当該円軌道から３次元ラドン・データが欠落している部分の領域を表わすシャドウ・ゾーン（２８）を規定する、Ｘ線断層画像再構成方法において、
   円形の第１ホーン形状ボリューム内に、前記円軌道から得られる画像を表現する円軌道画像表現工程と、
   前記円形第１ホーン形状ボリュームに対応する線形の第２のホーン形状ボリューム内に、前記線軌道から得られる画像を表現する線軌道画像表現工程と、
   前記円形第１のホーン形状ボリュームの各々と前記線形第２ホーン形状ボリュームの各々を夫々横切る垂直平面（２７）にある２次元画像を複数の２次元フーリエ平面へ変換する２次元フーリエ平面変換工程と、
   前記３次元ラドン空間の前記シャドウ・ゾーン（２８）が前記３次元フーリエ空間にマッピングされることにより生ずる３次元シャドウ・コーン（３０）の前記複数の２次元フーリエ平面における投影（２４）を、前記円軌道データから、除去し、該投影（２４）を、前記線軌道データから得られた対応する投影により置き換えて、修復された複数の２次元フーリエ平面を形成する２次元フーリエ平面修復工程と、
   前記修復された複数の２次元フーリエ平面に対して逆２次元フーリエ変換を施して、結果として得られる第３のホーン形状ボリュームを表わすそれぞれの複数の２次元画像平面を得る２次元画像平面形成工程と、
   結果として得られた前記第３ホーン形状ボリュームを前記Ｘ線断層画像へ変換するＸ線断層画像変換工程と、
   を具備したことを特徴とするＸ線断層画像再構成方法。
2. 前記円軌道画像表現工程は、
   θ_０を３次元ラドン空間における前記シャドウ・ゾーン（２８）の限界を表わす角度θ_０として定義すると、コーン頂角２θ_０を有するシャドウ・コーンを選択する工程を含むことを特徴とする請求項１に記載のＸ線断層画像再構成方法。
3. 前記線軌道画像表現工程は、前記線軌道データから複数の垂直画像スライスを生成して、前記線形第２ホーン形状ボリュームを形成する垂直画像スライス生成工程を含んでいることを特徴とする請求項１に記載のＸ線断層画像再構成方法。
4. 前記垂直画像スライス生成工程は、前記線軌道データの投影から形成される複数の平行ビームを再構成する工程を含むことを特徴とする請求項３に記載のＸ線断層画像再構成方法。
5. さらに、前記垂直画像スライスをフィルタ補正逆投影法を用いて再構成するフィルタ補正逆投影法再構成工程を有することを特徴とする請求項４に記載のＸ線断層画像再構成方法。
6. さらに、前記複数の垂直画像スライスを正フーリエ法を用いて再構成する正フーリエ法再構成工程を有することを特徴とする請求項４に記載のＸ線断層画像再構成方法。
7. 前記正フーリエ法再構成工程は、逆双１次補間を用いて、データを極格子から直交格子へ変換する第１逆双１次補間工程を含んでいることを特徴とする請求項６に記載のＸ線断層画像再構成方法。
8. 前記第１逆双１次補間工程は、４つの最も近接した格子位置ｇ_ｉ，ｊ、ｇ_{ｉ＋１，ｊ}、ｇ_{ｉ，ｊ＋１}及びｇ_{ｉ＋１，ｊ＋１}に対して各々の線軌道フーリエ・データ値ｘに対して逆双１次補間を行う第２逆双１次補間工程を含むことを特徴とする請求項７に記載のＸ線断層画像再構成方法。
9. 前記第２逆双１次補間工程は、各々の前記線軌道データ値ｘのそれぞれの位置と、前記４つの最も近接した格子位置ｇ_ｉ，ｊ、ｇ_{ｉ＋１，ｊ}、ｇ_{ｉ，ｊ＋１}及びｇ_{ｉ＋１，ｊ＋１}の各々との間の正規化された部分面積ａ_００、ａ_０１、ａ_１０、ａ_１１を決定する部分面積演算工程を含むことを特徴とする請求項８に記載のＸ線断層画像再構成方法。
10. それぞれの前記線軌道データ値ｘに隣接した各々の格子点ｇ_ｉ，ｊを収集すると共に加重ｗ_ｉｊを割り当てる工程を更に含んでいることを特徴とする請求項９に記載のＸ線断層画像再構成方法。
11. 前記部分面積演算工程は、次の等式 ｇ_ｉ，ｊ＝ｇ_ｉ，ｊ＋ａ_１１ｘ ｇ_{ｉ＋１，ｊ}＝ｇ_{ｉ＋１，ｊ}＋ａ_０１ｘ ｇ_{ｉ，ｊ＋１}＝ｇ_{ｉ，ｊ＋１}＋ａ_１０ｘ ｇ_{ｉ＋１，ｊ＋１}＝ｇ_{ｉ＋１，ｊ＋１}＋ａ_００ｘに従って各々の前記格子位置を更新する工程を含んでいることを特徴とする請求項１０に記載のＸ線断層画像再構成方法。
12. それぞれのｇ_ｉ，ｊをｗ_ｉｊにより除算することにより各々の格子データｇ_ｉ，ｊを正規化する工程を更に含んでいることを特徴とする請求項１１に記載のＸ線断層画像再構成方法。
13. 円軌道からのデータと線軌道からのデータとを含むＸ線断層画像を再構成するＸ線断層画像再構成装置において、前記円軌道が、３次元ラドン空間において、当該円軌道から３次元ラドン・データが欠落している部分の領域を表わすシャドウ・ゾーン（２８）を規定する、Ｘ線断層画像再構成装置であって、
    円形の第１ホーン形状ボリューム内に、前記円軌道から得られる画像を表現する円軌道画像表現手段と、
    前記円形第１ホーン形状ボリュームに対応する線形の第２のホーン形状ボリューム内に、前記線軌道から得られる画像を表現する線軌道画像表現手段と、
    前記円形第１のホーン形状ボリュームの各々と前記線形第２ホーン形状ボリュームの各々を夫々横切る垂直平面（２７）にある２次元画像を複数の２次元フーリエ平面へ変換する２次元フーリエ平面変換手段と、
    前記３次元ラドン空間の前記シャドウ・ゾーン（２８）が前記３次元フーリエ空間にマッピングされることにより生ずる３次元シャドウ・コーン（３０）の前記複数の２次元フーリエ平面における投影（２４）を、前記円軌道データから、除去し、該投影（２４）を、前記線軌道データから得られた対応する投影により置き換えて、修復された複数の２次元フーリエ平面を形成する２次元フーリエ平面修復手段と、
    前記修復された複数の２次元フーリエ平面に対して逆２次元フーリエ変換を施して、結果として得られる第３のホーン形状ボリュームを表わすそれぞれの複数の２次元画像平面を得る２次元画像平面形成手段と、
    結果として得られた前記第３ホーン形状ボリュームを前記Ｘ線断層画像へ変換するＸ線断層画像変換手段と、
    を備えたことを特徴とするＸ線断層画像再構成装置。
14. 前記円軌道画像表現工程は、
    θ_０を３次元ラドン空間における前記シャドウ・ゾーン（２８）の限界を表わす角度θ_０として定義すると、コーン頂角２θ_０を有するシャドウ・コーンを選択する工程を含むことを特徴とする請求項１３に記載のＸ線断層画像再構成装置。
15. 前記線軌道画像表現工程は、前記線軌道データから複数の垂直画像スライスを生成して、前記線形第２ホーン形状ボリュームを形成する垂直画像スライス生成工程を含んでいることを特徴とする請求項１３に記載のＸ線断層画像再構成装置。
16. 前記垂直画像スライス生成工程は、前記線軌道データの投影から形成される複数の平行ビームを再構成する工程を含むことを特徴とする請求項１５に記載のＸ線断層画像再構成装置。
17. さらに、前記垂直画像スライスをフィルタ補正逆投影法を用いて再構成するフィルタ補正逆投影法再構成工程を有することを特徴とする請求項１６に記載のＸ線断層画像再構成装置。
18. さらに、前記複数の垂直画像スライスを正フーリエ法を用いて再構成する正フーリエ法再構成工程を有することを特徴とする請求項１６に記載のＸ線断層画像再構成装置。
19. 前記正フーリエ法再構成工程は、逆双１次補間を用いて、データを極格子から直交格子へ変換する第１逆双１次補間工程を含んでいることを特徴とする請求項１８に記載のＸ線断層画像再構成装置。

Images