JP2006524059A - スパイラルｃｔ用の３ｐｉアルゴリズム - Google Patents

Abstract

２次元検出器を使用し３ＰＩアルゴリズムを用いたスパイラル走査される移動する対象物の画像を復元する方法及びシステムである。移動する対象物は、従来のシステムよりも最大約３倍ゆっくりした割合で走査することができる。好ましい実施形態では、本発明は、従来のスパイラル走査システムのものと同じ大きさの検出器アレイを使用したままで、スパイラル走査装置の中で寝台に載せた患者を最大約３倍ゆっくりと移動させることができる。

Description

本発明は、米国特許仮出願第６０／４３０，８０２号（出願日：２００２年１２月４日）の優先権の利益を主張するものである。また、本発明は、米国特許出願第１０／３８９，５３４号（出願日：２００３年３月１４日）の一部継続出願である。なお、米国特許出願第１０／３８９，５３４号は、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）の一部継続出願である、米国特許出願第１０／３８９，０９０号（出願日：２００３年３月１４日）の一部継続出願である。

本発明は、コンピュータ断層撮影（computer tomography：CT）に関し、特に、３ＰＩアルゴリズムを使用したスパイラル走査によって取得したデータから３次元画像を復元するプロセス、方法及びシステムに関する。

この３０年間、コンピュータ断層撮影（ＣＴ）は、スライス・バイ・スライス走査に基づく画像復元から、スパイラル走査による画像復元まで行われてきた。１９７０年代から１９８０年代までは、スライス・バイ・スライス走査が用いられた。この方法では、寝台に載せた患者をガントリーの中で少しずつ移動させるのと、ガントリーの回転とが交互に行われる。ガントリーの回転中、患者は静止しているので、患者の周りを回るＸ線源の軌道は円形になる。事前に選択された患者を横切るスライスは、そのような円形走査によって取得されたデータを使用して復元される。１９８０年代の中頃から今日まで、ＣＴでのデータ収集には、スパイラル型の走査が好ましい方法となった。スパイラル走査では、患者を載せた寝台は、寝台の周りを連続的に回転するガントリーの中を連続的に移動する。初めは、スパイラル走査は、データを１次元で受け取る１次元検出器を使用していた（検出器は一列に配置される）。その後、検出器が複数列（２列以上）で互いに隣接して配置される２次元検出器が導入された。ＣＴでは、特に２次元検出器を使用する場合に、画像の復元について重大な問題があった。以下では、２次元検出器で作成されるデータを、コーンビーム（cone-beam：CB）データ又はＣＢ投影と呼ぶ。

２次元検出器を使用したスパイラス走査により得られたデータからの３次元（別名：ボリュームメトリック）画像復元用としては、厳密アルゴリズムと近似アルゴリズムの２種類のアルゴリズムが知られている。厳密アルゴリズムと近似アルゴリズムは、どちらも既知の問題を有している。理想的な環境下では、厳密アルゴリズムは、精密な画像を復元することができる。したがって、厳密アルゴリズムには、非理想的な（すなわち、現実的な）環境下でも良質な画像を作成することが期待される。また一方、厳密アルゴリズムは、画像を復元するのに長時間かかることが知られている。また、厳密アルゴリズムは、その使用時にコンピュータの処理能力を多大に必要とする。これらのアルゴリズムは、大量のコーンビーム投影をメモリに保存する必要がある。さらに、ある厳密アルゴリズムでは、実施するために大きな検出アレイを必要としたり、走査する患者の大きさに制限があったりする。

近似アルゴリズムは、フィルタ補正逆投影（filtered back projection：FBP）の構造を有しているので、画像をとても効率良く作成することができる。また、厳密アルゴリズムよりも、使用するコンピュータの処理能力が少ない。また一方、近似アルゴリズムは、理想的な環境下でも、近似画像（精密な画像とは似ているものの、精密な画像とは異なっている）を作成することができる。特に、近似アルゴリズムは、画像における偽の特徴であるアーチファクトが作成され得る。ある環境下では、このようなアーチファクトはかなり強まる。

理論的には厳密であり推移不変（shift invariant）のＦＢＰ構造を有する前者のアルゴリズムは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）に開示されている。これらのアルゴリズムの短所は、１ＰＩモードで動作することである。このことは、検出器アレイが軸方向に大きい場合は、検出器すべてを使用するために、ガントリー内で患者を急いで平行移動させる必要があることを意味する。しかしながら、患者を急いで移動させることは、明白な理由により実用的ではない。とは言え、患者をゆっくりと移動させると、検出器の一部しか使用できない。この場合、Ｘ線照射の一部が無駄になる、打切りアーチファクト（discretization artifact）が強まる、ノイズ安定性が減少するなどの望ましくない結果をもたらす。

本発明の主な目的は、２次元検出器を使用したスパイラル方式により走査された対象物の画像を復元するための３ＰＩアルゴリズムを提供することにある。任意のボクセルで画像を復元するためには、これらのアルゴリズムは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許出願第６，５７４，２９９号）の１ＰＩのアルゴリズムよりも、スパイラルの断面が長いことが必要とされる。したがって、新規のアルゴリズムは、同じ大きさの検出器アレイを使用したままで、患者を約３倍ゆっくりと移動させることを可能にする。

本発明の第１の好ましい実施形態では、スパイラル走査された対象物の画像を復元するのに、５つの全体的なステップから成るプロセスを使用する。第１のステップでは、現在のＣＢ投影が測定される。次に、新規のアルゴリズムに従って、検出器上の線群が特定される。次に、隣接する投影間で導関数のコンピュータ処理が行われ、その後、選択された線群からの線に沿って、フィルタ処理によって導関数の畳み込みが行われる。次に、フィルタ処理されたデータを使用して逆投影を行うことにより、画像は更新される。最終的に、この前処理ステップは、全ての対象物が走査されるまで、各ＣＢ投影について繰り返される。この実施形態は、いくつかの（約２〜４）ＣＢ投影を同時にメモリに保存し、１つの線群を使用する。

第２の実施形態では、いくつかのＣＢ投影をメモリに保存するのと共に、異なる線群を利用することができる。

これらの実施形態を変更して、一度に、１つのＣＢ投影だけをコンピュータ・メモリに保存することが可能である。これは、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許出願第６，５７４，２９９号）（この参照により本発明に含まれるものとする）で行われたのと同じようにして実現することができる。

したがって、新規のアルゴリズムは、同じ大きさの検出器アレイを使用したままで、患者を約３倍ゆっくりと移動させることを可能にする。

本発明のさらなる目的及び利点は、以下の添付図面に概略的に図示した、現時点での好ましい実施形態についての詳細な説明によって明らかになるであろう。

本発明の開示された実施形態を詳細に説明する前に、本発明はここで説明される詳細によって限定されるものではなく、他の実施形態をとることもできることを理解されたい。また、ここで使用される専門用語は、説明目的のためのものであって、本発明を限定するためのものではない。

本発明は、米国仮出願第６０／３１２，８２７号（出願日：２００１年８月１６日）の優先権の利益を主張する米国特許出願第１０／１４３，１６０号（出願日：２００２年５月１０日、現在は米国特許第６，５７４，２９９号）（この参照により本発明に含まれるものとする）の一部継続出願である。

以下、本発明について詳細に説明する。

《理論背景》
まず、必要な表記について説明する。

と仮定する。ただし、ｈ＞０はスパイラル（らせん）であり、Ｕはスパイラルの完全に内側の開集合である。

Ｓ^２は、Ｒ^３における単位球面である。そして、

Ｄ_ｆ（ｙ，β）は、ｆのＣＢ変換である。

とすると、

ｊ＝１，２，…，は、平面

とＣとの交点を表すものとする。また、

である。
β∈Ｓ^２に関しては、β^⊥は大円｛α∈Ｓ^２：α・β＝０｝を表す。コンピュータ処理が必要なｆについて、任意のｘ∈Ｒ^３を決定する。ｆ（ｘ）をコンピュータ処理するために、後に決定される、有限の範囲のスパイラルの断面を使用する。とりあえず、今はＣ（ｘ）で示す。対応するパラメータ区間（parametric interval）はＩ（ｘ）で示される。Ｃ（ｘ）についての主な仮定は、下記の通りである。

〈性質Ｃ１〉（完全状態）
ｘを通る任意の平面は、少なくとも１点でＣと交差する。

反転公式の構成における重要な要素は加重関数ｎ（ｓ，ｘ，α），α∈β^⊥（ｓ，ｘ）である。関数ｎは、次のように理解される。ｘとαは平面П（ｘ，α）を決定し、加重関数ｎは、ｘの位置によって決まるｙ（ｓ）∈П（ｘ，α）ПＣに割り当てられる。従って、この解釈から、ｎ（ｓ，ｘ，α）＝ｎ（ｓ，ｘ，−α）であると推定する。ｎについての主な仮定は、下記の通りである。

〈性質Ｗ１〉（通常状態）

〈性質Ｗ２〉
ｎ（ｓ，ｘ，α）が任意の（ｓ，α）の近傍で局地的に一定である、有限の多くのＣ^１関数α_ｋ（ｓ，ｘ）∈β^⊥（ｓ，ｘ），ｓ∈Ｉ（ｘ）が存在する。ただし、ｓ∈Ｉ（ｘ）であり、α∈β^⊥（ｓ，ｘ），

である。

を示す。

Ｃｌ、Ｗ１、及びＷ２の仮定の下では、次の一般的な反転公式は、A. Katsevichの論文「A general scheme for constructing inversion algorithms for cone beam CT," International Journal of Mathematics and Mathematical Sciences, Vol. 21, pp. 1305-1321 (2003)」から導かれる。

ただし、θ_ｍ∈（０，π）は、φ（ｓ，ｘ，θ）が不連続になる点であり、ｃ_ｍ（ｓ，ｘ）はジャンプ値である。

ｎが適切に選択される限り、反転公式は必ずしもＦＢＰ型を必要としない。

《３ＰＩ線及びその性質》
あるｓ_０∈Ｒについて、Ｘ線源がｙ（ｓ_０）に固定されていると仮定する。検出器アレイはＸ線源と共に回転するので、検出器平面（detector plane）はｓ_０よって決まり、ＤＰ（ｓ_０）で示される。ＤＰ（ｓ_０）はスパイラルの軸と平行であり、Ｘ線源と反対側の点において円柱：ｙ_１ ^２＋ｙ_２ ^２＝Ｒ^２（（１）参照）に接すると仮定する。したがって、ｙ（ｓ_０）と検出器平面との間の距離は、２Ｒである（図３参照）。検出器平面内の座標は、次の通りである。ｄ_１軸がスパイラル軸と垂直であり、ｄ_２軸がスパイラル軸と平行であり、起点はｙ（ｓ_０）の投影と一致すると仮定する。スパイラルの上側及び下側の回転を検出器表面上に立体投影したのを図３に示す。これにより、次のパラメータ曲線が与えられる。

ただし、ρは、患者の支持体（support）の半径によって決まる。
ρ＝２ｃｏｓ^−１（ｒ／Ｒ）（（２）参照）。これらの曲線は、Γ_ｊ，ｊ＝±１，±２,…,で示される（図４参照）。

はｘの投影を示す。
（なお、

は、明細書中では、「ｘ´」とも表記する）
任意のＸ線源位置を、スパイラル上の全ての点ｙ（ｓ）に結び付ける。ただし、

である。

これにより、Ｓ_Ｕ ^３ＰＩ（ｓ_０）とＳ_L ^３ＰＩ（ｓ_０）で表される２つの表面が得られる。この２つの表面と、円柱：ｘ_１ ^２＋ｘ_２ ^２＝Ｒ^２とによって囲まれた領域は、Ｖ^３ＰＩ（ｓ_０）で示される。ｘは、固定されていると仮定する。ｓ_０が十分に小さければ、Ｓ_Ｕ ^３ＰＩ（ｓ_０）はｘよりも下側になる。ｓ_０が増大すると、２つの場合が有り得る。１つは、連続的照度（continuous illumination）として知られており、ｘはＳ_Ｕ ^３ＰＩ（ｓ_０）を通ってＶ^３ＰＩ（ｓ_０）に入り、Ｓ_L ^３ＰＩ（ｓ_０）を通ってＶ^３ＰＩ（ｓ_０）から出る。明らかに、上記の方法は、一意的な３ＰＩ区間［ｂ_０（ｘ），ｔ_０（ｘ）］と、対応する一意的な線Ｌ_０ ^３ＰＩ（ｘ）を与える。２番目の場合は、遮断照度（interrupted illumination）として知られており、ｘは、Ｖ^３ＰＩ（ｓ_０）に何回も出入りする。より正確には、ｘは、Ｓ_Ｕ ^３ＰＩ（ｓ_０）とＳ_Ｌ ^３ＰＩ（ｓ_０）の各表面と、正確に３回づつ交差する。したがって、上記の方法は、３つの３ＰＩ線を与える。Ｌ_ｉ ^３ＰＩ（ｘ），ｉ＝１，２，３。対応するパラメータの値は、ｂ_ｉ（ｘ），ｔ_ｉ（ｘ），ｉ＝１，２，３と表される。スパイラルは対称であるため、ｓ_０=ｂ_ｉ（ｘ）の場合、ｘはＶ^３ＰＩ（ｓ_０）から出る（に入る）。そして、ｓ_０=ｔ_ｉ（ｘ），ｉ＝１，２，３の場合、ｘはＶ^３ＰＩ（ｓ_０）に入る（から出る）。平面Ｘ_３＝０と見なす。この２つの場合の間の境界線を図５に示す。図５から、曲線は自分自身とは交差しておらず、
円：ｘ_１ ^２＋ｘ_２ ^２＜１を、Ｘ_１とＸ_３の領域に分割していることが分かる。Ｘ_１を示す中央の領域では、各ｘに対して１つの３ＰＩ線がある。Ｘ_１を示す外側の領域では、各ｘに対して３つの３ＰＩ線がある。

次の特性を確立できる。点ｘに対して３ＰＩが１つしかない場合は、

である。

点ｘに対して３つの３ＰＩがある場合は、

である。

検出器平面に関しては、方程式（１３ａ）及び（１３ｂ）は、ｓ∈Ｉ^３ＰＩ（ｘ）の時かつその時に限り、ｘ´はΓ_２とΓ_−２との間にあることを示す。したがって、１ＰＩの場合との類推から、この検出器上の領域を３ＰＩウインドウと呼ぶことができる。１ＰＩの場合では、ｘの１ＰＩ線の終点によって囲まれるパラメータ区間は、ｘの１ＰＩパラメータ区間と呼ばれる。同様に、Ｉ^３ＰＩ（ｘ）は、ｘの３ＰＩパラメータ区間と呼ばれる。

《補助構成》
まず、連続的照度が行われると仮定する。単位球面Ｓ^２の表面上にある次の２つの曲線を検討する。一方の曲線は、Ｌ_０Ｉ^３ＰＩ（ｘ）に対して垂直である全ての単位ベクトルから成り、Ａで示される。他方の曲線は、ベクトル

から成り、Ｔで示される。これらの曲線は、

により定義される球座標（θ_１，θ_２）を使用して、平面上で示されるので、これらの曲線をＳ^２上に直接的に示すのは、都合が良くない。αと−αは同一の平面を定義するので、θ_１を任意の長さπの区間に限定することができる。識別することにより、反対側の境界を「接着（glue）」することができる。

ｘ／Ｒ＝（０，０．２５，０）の典型的な状態を図６に示す。単位ベクトルとしてだけではなく、平面についても、Ｓ^２上の点について考えることは都合がよい。各α∈Ｓ^２は、法線ベクトルαを有するｘを通る一意的な平面と対応している。この対応は、反対方向のベクトルが特定できる場合は、１対１である。ＡとＴとで、Ｓ^２を複数の領域Ｄ_１，Ｄ_２,…,に分割する。固定されたｉの構成によって、Ｃ^３ＰＩ（ｘ）∩П（ｘ，α）における点の数は、全てのα∈Ｄ_ｉについて同じである。Ｄ_ｉがｋ交点（ＩＰｓ）を含む場合は、ｋＩＰ領域を呼ばれる。

次に、遮断照度が行われる場合を考える。同様の方法で、Ｓ^２の表面上に、いくつかの曲線を定義する。初めの３つの曲線は、それぞれ３つの３ＰＩ軸に垂直な単位ベクトルを考慮しながら取得される。それらは、Ａｋ、ｋ＝１，２，３と示される。曲線の第２の集合は、式（１５）中のｓを、区間［ｂ_１，ｂ_２］，［ｂ_３，ｔ_１］，［ｔ_２，ｔ_３］に限定することにより得られる。これらの曲線は、それぞれＴ_１２、Ｔ_３１、及びＴ_２３と示される。ｘ／Ｒ=（−０．５，０，０）についての典型的な状況を図７に示す。

ｘ_１及びｘ_２が互いに近接していれば、それらは同じダイアグラムを共有することは明らかである。したがって、数列を滑らかなに変化させることにより、あるダイアグラムは他のダイアグラムに変換される。そして対応する領域では、ＩＰの数と、Ｉ^３ＰＩでの部分区間におけるそれらの分配は、同じままである。本質的な変化は、３つの境界がＳ^２上の一点で互いに交差するという、「重大な事象」が起こるｘの近傍でのみ可能性がある。これらの点は、数的に発見することができる。その結果を図８に示す。これらの点から回転の中心までの最少距離は、

である。したがって、医療用に関心のある状況下では（ｒ_FOV５Ｒ）、図６に示した連続的照度と、図７に示した遮断照度との、２つの場合が起こり得る。

《加重関数ｎの構成》
任意のｓ∈Ｉ^３ＰＩ（ｘ）について、次の規則に従って、加重関数ｎ（ｓ，ｘ，α）、α∈β^⊥（ｓ，ｘ）を決定する。

・１ＩＰ領域では、唯一のＩＰは重量ｎ＝１を得る。

・３ＩＰ領域では、各ＩＰは重量ｎ＝１／３を得る。

・５ＩＰ領域では、２つのＩＰは重量ｎ＝０を得、残り全てのＩｐはｎ＝１／３を得る（図９及び図１０参照）。

《第１の復元アルゴリズム》
重量が求められると、アルゴリズムを構成する全ての要素が得られる。まず第１に、ｘの画像を復元するために、ｘ：Ｃ（ｘ）＝Ｃ^３ＰＩ（ｘ）とＩ（ｘ）＝Ｉ^３ＰＩ（ｘ）の３ＩＰスパイラル分割を行う（Section 1参照）。（６）や（７）から、各ｓ∈Ｉ^３ＰＩ（ｘ）について、フィルタ処理する方向は、

の切れ目を発見することにより決定されることが分かる。これらの切れ目の研究は、次の３つの線群をもたらす。

第１の群は、スパイラル接線に平行な線から成る。この群は、Ｌ_０と表される。この線と関連する係数を、図１１に示す。明らかに、任意のｘ´について、ｘ´を含む１つのＬ∈Ｌ_０だけがある。

図１２は、フィルタ処理する線の第２の群を示す。この群は、Γ_±１に接する線から成り、Ｌ_１と表される。１ＰＩウインドウ内の各ｘ´に対して、Ｌ_１から２つの線を引く。これらの線は、図１３に示した規則に従って決定する。また、この図は、関連した係数ｃ_ｍを示す。

図１４（左側の図）は、フィルタ処理したΓ_±１に接する線の群を示す。この群は、Ｌ_２と表される。任意のｘ´∈Ｆ_１∪Ｆ_２∪Ｆ_４∪Ｆ_５では、１つ以上のｘ´を含む線Ｌ∈Ｌ_２がある。一意的な線は、次の規則に基づいて決定される（図１４の右側の図を参照されたい）。ｘ´がｘ´∈Ｆ_１∪Ｆ_４である場合は、接点はｘ´の右側となる。また、ｘ´がｘ´∈Ｆ_２∪Ｆ_４である場合は、接点はｘ´の右側となる。全ての場合において、ｃ_ｍ＝１／３である。

図１５は、図１１〜１４に含まれる情報を要約しており、ｘ´となり得る全ての可能性がある場合と、全ての関連したフィルタ処理方向及び定数ｃ_ｍを示している。全ての場合において、フィルタ処理の方向は左から右である仮定する。この構成の結果として（ｃｆ．図１５）、一般的な復元式（７）、（８）中のθｍ（ｓ，ｘ）及びｃ_ｍ（ｓ，ｘ）は、β（ｓ，ｘ）だけによって、ｘによってのみ決まる。したがって、ｃ_ｍとα^⊥の理由で、ｘをβ（ｓ，ｘ）によって置き換えることができ、式（７）を次の形式で書き直すことができる。

また、ここで留意すべきは、フィルタ処理する線Ｌが定められると、その投影がＬに含まれるかつ上記した規則を満たす全てのｘは、同じフィルタ処理線Ｌを共有する（図１５参照）。したがって、（１９）は畳み込まれ、（１８）は逆投影となり、アルゴリズム（１８）・（１９）は畳み込みに基づくＦＢＰ型となる。このことは、米国特許第１０／１４３，１６０号（出願日：２００２年５月１０日）にも同様に見られる。逆投影法に関しては、他の方法及び技術を用いることもできる。さらに、例えば、Ｋａｔｓｅｖｉｃｈによる米国特許第６，５７４，２９９号を参照されたい（この参照により本発明に含まれるものとする）。

《改善された復元アルゴリズム》
Ψ（ｔ）をＲ上に定義される、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという性質を有する滑らかな関数と仮定する。Γ_±１∪Γ_±２：ｓ１，ｓ２，ｓ３と３つの交差点を有し、それらの点が

を満たす任意の既知の線Ｌ∈Ｌ_２´を求めることにより、新しい線群Ｌ_２´を定義する。ｓがｘ線現の現在位置であることを思い出されたい。この線Ｌ∈Ｌ´_２は、例えば、
ｓ_３，２π＜｜ｓ_３｜＜４πによりパラメータ化される。ｓ_１及びｓ_２の位置は、ｓ_３の位置によって決まる（図１６に示す）。

Ψの性質を利用して、それぞれのｘ´∈Ｆ_１∪Ｆ_２∪Ｆ_４∪Ｆ_５について、ｘ´を含む一意的なＬ∈Ｌ´_２を容易に決定することができる。また、
ｘ´→Ｌ_２ ^ｃｒであれば、ｓ_２，ｓ_３→２Δであり、ｓ_１→ｓである。同様に、
ｘ´→Ｌ_−２ ^ｃｒであれば、ｓ_２，ｓ_３→−２Δであり、ｓ_１→ｓである。

さらに重要なことには、ｘ´∈Ｆ_１∪Ｆ_２∪Ｆ_４∪Ｆ_５の場合は、フィルタ処理方向を２つから１つに減らすことができる。この付加的な利益は、検出器の使用法を改善する。したがって、図１５に２つの上側の図面は、図１７に示したダイアグラムと置き換えられる。ｘ´がＬ_０（即ち、ｘ´∈Ｆ_４∪Ｆ_５）の下側に現れる場合は、検出器平面の基点に対する対称により、図１７から得られる。

反転公式（１８）、（１９）の形式は、変更されない。第１の違いは、ｘ´∈Ｆ_１∪Ｆ_２∪Ｆ_４∪Ｆ_５の場合は、Ｍ（ｓ，β）＝１となることである（２ではなく）。第２の違いは、フィルタ処理の方向であり、ここでは式（２１）、（２２）で表される。

《復元アルゴリズムの一般的な概要》
図２は、本発明の基本処理ステップ１０、２０、３０、４０、５０の概要を示す。以下、各ステップについて説明する。

ステップ１０：現在のＣＢ（コーンビーム）投影を、コンピュータ・メモリにロードする。現在メモリに記憶されているＣＢ投影の中央点を、ｙ（ｓ_０）と仮定する。

ステップ２０：フィルタ処理する線群を求める。図１５を使用して、又は改善されたアルゴリズムの場合はさらに図１７を使用して、必要な線群を特定する。その後、各群から線の離散集合を選択する。

ステップ３０：フィルタ処理の準備。ステップ２０で選択された前記線上の点を、極角によりパラメータ化する。補間法を用いて、前記極角の離散値に対応する前記線上の点βにおけるＣＢデータ

の導関数を計算する。

ステップ４０：フィルタ処理。ステップ２０で特定された各線について、フィルタ１／sinγによってステップ３０で計算されたデータを畳み込む。

ステップ５０：逆投影。各復元点ｘについて、ステップ４０で求められたフィルタ処理されたデータを、方程式（１８）に従って逆投影する。その後、処理すべき新しいＣＢが無くなるまで、又は全ての必要な点ｘにおいて画像復元が完成するまで、ステップ１０に進む。

次に、図２に示した５つのステップ１０〜５０のアルゴリズムについて、詳細に説明する。

〈ステップ１０〉
現在のＣＢ（コーンビーム）投影を、コンピュータ・メモリにロードする。現在メモリに記憶されているＣＢ投影の中央点を、ｙ（ｓ_０）と仮定する。

〈ステップ２０〉
フィルタ処理する線群を求める。線群は、次のサブステップ２１、２２、及び２３により選択することができる。

ステップ２１：線群Ｌ_０から、スパイラル接線と平行であり、Γ_２とΓ_−２との間に位置する検出器の表面上における関心領域の投影を覆う線の等距離集合を選択する（図１１を参照）。

ステップ２２：線群Ｌ_１から、Γ_１とΓ_−１に接する線の離散集合を選択する（図１２参照）。Γ_１の接線における左端の点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。同様に、Γ_−１の接線における右端の点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。

ステップ２３：線群Ｌ_２から、Γ_２とΓ_−２に接する線の離散集合を選択する（図１４の左側の図を参照されたい）。どちらの場合も、接点は、検出器の表面上における関心領域の投影を超えて広がる必要はない。改善されたアルゴリズムの場合は、Γ_２とΓ_−２に接する線の代わりに、Γ_２とΓ_−２曲線上のｓ_３に対する離散（等距離）集合の値を選択し、方程式（２１）、（２２）を解くことにより、線Ｌ∈Ｌ_２´を決定する。両方の曲線上で、点ｓ_３は、検出器の表面上における関心領域の投影を超えて広がる必要はない。

〈ステップ３０〉
フィルタ処理の準備。

ステップ３１：ステップ２０で得られた線の集合から、線Ｌを決定する。

ステップ３２：ｙ（ｓ_０）とＬを通る平面内で、前記線上の点を極角γによりパラメータ化する。

ステップ３３：後にステップ４０で離散フィルタ処理に使用される、等距離値γ_ｊの離散集合を選択する。

ステップ３４：各γ_ｊについて、ｙ（ｓ_０）から、γ_ｋに対応するＬ上の点に向う単位ベクトルβ_ｊを求める。

ステップ３５：ｓ_０に近い少しの値ｑについてのＣＢ投影データＤ_ｆ（ｙ（ｑ），Θ）を使用して、全てのΘ=β_ｊについての導関数

を数的に求める。

ステップ３６：計算された導関数の値をコンピュータ・メモリに記憶させる。

ステップ３７：ステップ２０で特定されたすべての線Ｌについて、ステップ３１〜３６を繰り返す。このようにして、ｙ（ｓ_０）に位置するＸ線源に対応する処理されたＣＢ投影データを作成する。

〈ステップ４０〉
フィルタ処理。

ステップ４１：ステップ２０で特定された前記線群Ｌ_ｍの内の１つから、線Ｌを決定する。

ステップ４２：前記線に沿って、ステップ３０で計算された前記処理されたＣＢデータの値のＦＦＴを計算する。

ステップ４３：フィルタ１／sinγのＦＦＴを計算する。

ステップ４４：フィルタ１／sinγのＦＦＴ（ステップ４３の結果）と、前記処理されたＣＢデータの値のＦＦＴ（ステップ４２の結果）とを乗算する。

ステップ４５：ステップ４４の結果結果の逆ＦＦＴを行う。

ステップ４６：ステップ４５の結果をコンピュータ・メモリに記憶させる。

ステップ４７：前記線群の全ての線について、ステップ４１〜４６を繰り返す。このことにより、フィルタ処理されたＣＢデータΨ_ｍ（ｓ_０，β_ｊ）が得られる。ただし、ｍは、Ｌが選択される線群の番号を表す。ｍ＝０，１，２。

フィルタ処理は、当該技術分野では周知であり、例えば、タム（Tam）による米国特許第５，８８１，１２３号（この参照により本発明に含まれるものとする）に示され説明されているようにして実施される。

〈ステップ５０〉
逆投影。

ステップ５１：画像復元が求められている患者の内側の点である復元点ｘを求める。

ステップ５２：ｓ_０がＩ^３ＰＩ（ｘ）に属する場合は、前記フィルタ処理されたＣＢデータは、ｘにおける画像に影響を与える。そして、ステップ５３〜５８を実行する。ｓ_０がＩ^３ＰＩ（ｘ）区間の内側にない場合は、前記フィルタ処理されたＣＢデータはｘにおける画像復元に使用されない。この場合は、ステップ５１に戻り、他の復元点を選択する。

ステップ５３：検出器の平面ＤＰ（ｓ_０）上にｘの投影ｘ´を求める。そして、ｙ（ｓ_０）からｘに向う、単位ベクトルβ（ｓ_０，ｘ）を求める。

ステップ５４：図１１、図１３、及び図１４の右側の図を使用して、前記線群から線を特定する。また、前記投影ｘ´に近い前記線上の点を特定する。ｘ´がＬ_２ ^ｃｒよりも上側、又はＬ_２ ^−ｃｒよりも下側であれば、改善されたアルゴリズムを使用する場合は、前記投影ｘ´に近いフィルタ処理された線を見つけるために、方程式（２１）、（２２）、及び図１６を使用する。このことにより、β（ｓ_０，ｘ）に近いβ_ｊについての少しの値Ψ_ｍ（ｓ_０，β_ｊ）が得られる。

ステップ５５：補間法を使用して、前記したβ（ｓ_０，ｘ）に近いβ_ｊについての値Ψ_ｍ（ｓ_０，β_ｊ）の値から、Ψ_ｍ（ｓ_０，β（ｓ_０，ｘ））の値を推定する。

ステップ５６：Ψ_ｍ（ｓ_０，β（ｓ_０，ｘ））に−ｃ_ｍ（ｓ，β（ｓ_０，ｘ）／［４π^２｜ｘ−ｙ（ｓ_０）｜］を掛けることにより、前記フィルタ処理されたＣＢデータから、点ｘで復元される画像に対する寄与率（contribution）を計算する。適切な逆投影率ｃ_ｍは、図１１、図１３、及び図１４の右側の図を使用して選択される（概要については図１５を参照されたい）。ｘ´がＬ_２ ^ｃｒよりも上側、又はＬ_２ ^−ｃｒよりも下側であれば、改善されたアルゴリズムの場合は、適切な逆投影率ｃ_ｍを見つけるために図１７を使用する。

ステップ５７：方程式（１８）での積分の近似値を求めるために、前記寄与率を、点ｘにおいて復元される前記画像に加える。なお、前記画像は、予め選択された方法（例えば、台形法）に従って復元される
ステップ５８：ステップ５１に進み、異なる復元点ｘを選択する。全ての復元点ｘが処理されれば、ステップ５９に進む。

ステップ５９：ステップ１０へ進み、次のＣＢ投影を、コンピュータ・メモリに読み込む。

画像データは、画像復元処理が完了した全ての復元点ｘで表示することができる（すなわち、それらの点では、画像復元するためのそれに続くＣＢ投影は全て必要としない）。コンピュータ・メモリからは、画像復元処理が完了していない点での画像復元に不要な全てのＣＢ投影を捨て去る。アルゴリズムは、走査が終了したとき又は必要な点での画像復元が全て完了したときに、結論を出す。

例えば、検出器が３ＰＩアルゴリズムに必要なデータを取得しなかった場合は、欠けているデータを推定するために、様々な方法を採ることができる。この場合は、検出器は、各コーンビーム投影の測定データ及び推測データを取得する仮想の検出器となる。欠けているデータを正確に推測することができれば、３ＰＩアルゴリズムは正確に復元することができる。ここでは、依然として、欠けているデータが近似的に求められる現実の環境下における正確な復元について述べている。

以上、本発明の実施の形態について説明したが、本発明は前記した実施の形態のみに限定されるものではなく、本発明の技術的思想に基づく限りにおいて、種々の変形が可能である。

Ｘ線源及び検出器アレイを有する回転するガントリー内を移動する寝台に載せた患者の一般的な配置を示す。コーンビーム投影データ集合は、Ｘ線検出器によって受け取られる。そして、コンピュータにより画像復元処理が行われる。復元された画像はモニタに表示される。 本発明の基本処理ステップの概要を示す。 検出器平面上におけるスパイラルの立体投影を示す。 検出器平面上の様々な投影及び重要な線を示す。 境界曲線を示す。 連続的照度の場合を示す。ｘ／Ｒ＝（０，０．２５，０） 遮断照度の場合を示す。ｘ／Ｒ＝（−０．５，０，０） 様々な重要な事柄が生じる点である。 連続的照度の場合における５ＰＩ領域の内側の重量分布を示す。 遮断照度の場合における５ＰＩ領域の内側の重量分布を示す。 Ｌ_０に平行なフィルタ処理された線の群を示す。 Γ_±１に接するフィルタ処理された線の群を示す。 フィルタ処理された線を、ｘ´の位置によって選択する方法を示す。 Γ_±２に接するフィルタ処理された線の群を示す。 様々な場合における、フィルタ処理された線及び関連する定数ｃ_ｍを示す。左上の図：ｘ´∈Ｆ_１、右上の図：ｘ´∈Ｆ_２、中央の図：ｘ´∈Ｆ_３、及びΓ_１以上、左下の図：ｘ´∈Ｇ_１、右下の図：ｘ´∈Ｇ_２。 点ｓ_１，ｓ_２，ｓ_３が取り得る位置を示す。 様々な場合における、フィルタ処理された線及び関連する定数ｃ_ｍ。左の図：ｘ´∈Ｆ_１、右の図：ｘ´∈Ｆ_２、 図２に示したステップ２０に対応する、フィルタ処理する線群を求める３つのサブステップを表わすフローチャートである。 図２に示したステップ３０に対応する、フィルタ処理準備の７つのサブステップを表わすフローチャートである。 図２に示したステップ４０に対応する、フィルタ処理の７つのサブステップを表わすフローチャートである。 図２に示したステップ５０に対応する、逆投影の５つのサブステップを表わすフローチャートである。 図２に示したステップ５０に対応する、逆投影の後半部分の３つのサブステップを表わすフローチャートである。

Claims (26)

1. 少なくとも１つの検出器から得られたデータから画像を復元する方法であって、
   軸方向に、現在の放射源位置近傍におけるスパイラル４回転分の投影よりも幅広のコーンビーム投影を少なくとも１つ検出する検出器を少なくとも１つ使用して、対象物をスパイラル方式で走査するステップと、
   畳み込みに基づくＦＢＰ（フィルタ補正逆投影）アルゴリズムにより、走査した対象物の精密な画像を効率良く復元するステップとを含むことを特徴とする方法。
2. 請求項１に記載の方法であって、
   前記走査ステップは、前記検出器を使用して前記対象物の２次元コーンビーム（ＣＢ）投影データを取得することを特徴とする方法。
3. 請求項２に記載の方法であって、
   １Ｐ１アルゴリズムに必要とされる検出器と実質的に同様の検出器を使用するステップをさらに含むことを特徴とする方法。
4. 請求項２に記載の方法であって、
   前記走査ステップは、１Ｐ１アルゴリズムで使用されるコーンビーム投影よりも軸方向に幅広のコーンビーム投影を検出するステップをさらに含むことを特徴とする方法。
5. 請求項４に記載の方法であって、
   前記走査ステップは、１Ｐ１アルゴリズムで使用されるコーンビーム投影よりも軸方向に少なくとも３倍幅広のコーンビーム投影を検出するステップをさらに含むことを特徴とする方法。
6. 請求項１に記載の方法であって、前記対象物は人を含むことを特徴とする方法。
7. ２次元検出器を使用したスパイラル・コンピュータ断層撮影走査によって得られた精密な画像をコンピュータ処理する方法であって、
   （ａ）１Ｐ１アルゴリズムに必要とされる検出器よりも幅広い検出器から、現在の放射源位置近傍におけるスパイラル４回転分の投影の範囲をカバーするコーンビーム（ＣＢ）投影データを収集するステップと、
   （ｂ）コーンビーム投影と交差する平面П上の線群を特定するステップと、
   （ｃ）ＣＢ投影データを前処理するステップと、
   （ｄ）前記線に沿って、前記前処理したＣＢ投影データを畳み込みフィルタ処理するステップと、
   （ｅ）前記画像の前段階と称すべきものを作成すべく、前記フィルタ処理したデータを逆投影するステップと、
   （ｆ）前記対象物の精密な画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ），（ｄ）及び（ｅ）を繰り返すステップとを含むことを特徴とする方法。
8. 請求項７に記載の方法であって、前記走査は前記対象物のＸ線暴露を含むことを特徴とする方法。
9. 請求項７に記載の方法であって、前記ステップ（ａ）〜（ｆ）は３ＰＩアルゴリズムを含むことを特徴とする方法。
10. 検出器を使用したコンピュータ断層撮影走査によって得られた画像をコンピュータ処理する方法であって、
    （ａ）対象物の走査中に検出器からコーンビーム（ＣＢ）データを収集するステップと、
    （ｂ）前記コーンビームと交差する平面ＤＰ（ｓ）上における、次の３つの線群を特定するステップと、
    （ただし、ｓは、走査経路を表し、現在の放射源位置に相当するパラメータ値である）
    （ｂｉ）第１の線群は、現在の放射源位置における走査接線の方向である
    
    

    

    に対して平行である。
    （ｂｉｉ）第２の線群は、ｓ＜ｑ＜ｓ＋２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_１、及びｓ−２π＜ｑ＜ｓによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_−１に対して接する。
    （ただし、ｑは、走査経路に沿ったパラメータであり、投影される点を表す）
    （ｂｉｉｉ）第３の線群は、ｓ＋２π＜ｑ＜ｓ＋４πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_２、及びｓ−４π＜ｑ＜ｓ−２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_−２に対して接する。
    （ｃ）前記３つの群の前記線に沿って、前記データを前処理及び推移不変フィルタ処理するステップと、
    （ｄ）前記画像の前段階と称すべきもの作成すべく、前記フィルタ処理したデータを逆投影するステップと、
    （ｅ）前記対象物の画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ）及び（ｄ）を繰り返すステップとを含むことを特徴とする方法。
11. 請求項１０に記載の方法であって、
    前記前処理は、現在の放射源位置についてのＣＢデータの微分計算を含むことを特徴とする方法。
12. 請求項１０に記載の方法であって、
    前記推移不変フィルタ処理は、前記前処理されたデータを１／ｓｉｎγフィルタにより畳み込むことを特徴とする方法。
13. 請求項１０に記載の方法であって、
    前記第１の線群からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、
    
    

    

    に平行であり、Γ_２に接するＬ_２ ^ｃｒと、
    
    

    

    に平行であり、Γ_−２に接するＬ_２ ^−ｃｒとの間に位置する場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝２／３を乗じることを特徴とする方法。
14. 請求項１０に記載の方法であって、
    前記第１の線群の線からの前記フィルタ処理されたデータの逆投影は、
    ｘのＤＰ（ｓ）への投影が、Ｌ_２ ^ｃｒの上側に位置する、又はＬ_２ ^−ｃｒの下側に位置する場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝１／３を乗じることを特徴とする方法。
15. 請求項１０に記載の方法であって、
    前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
    ｘのＤＰ（ｓ）への投影が、Γ_１とΓ_−１との間に位置し、線がΓ_１∪Γ_−１と接する点がｘの１ＰＩパラメータ区間の内側にある場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝２／３を乗じることを特徴とする方法。
16. 請求項１０に記載の方法であって、
    前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
    ＤＰ（ｓ）へのｘの投影が、Γ_１とΓ_−１との間に位置し、線がΓ_１∪Γ_−１と接する点がｘの１ＰＩパラメータ区間の外側にある場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝−２／３を乗じることを特徴とする方法。
17. 請求項１０に記載の方法であって、
    前記第３の線群の線からの前記フィルタ処理されたデータの逆投影は、
    前記フィルタ処理されたデータに係数ｃ_ｍ＝１／３を乗じることを含むことを特徴とする方法。
18. 検出器を使用したコンピュータ断層撮影走査によって得られた画像をコンピュータ処理する方法であって、
    （ａ）対象物の走査中に検出器からコーンビーム・データを収集するステップと、
    （ｂ）前記コーンビームと交差する平面ＤＰ（ｓ）上における、次の３つの線群を特定するステップと、
    （ただし、ｓは、走査経路を表し、現在の放射源位置に相当するパラメータ値である）
    （ｂｉ）第１の線群は、現在の放射源位置における走査接線の方向である
    
    

    

    に対して平行である。
    （ｂｉｉ）第２の線群は、ｓ＜ｑ＜ｓ＋２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_１、及びｓ−２π＜ｑ＜ｓによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_−１に対して接する。
    （ただし、ｑは、走査経路に沿ったパラメータであり、投影される点を表す）
    （ｂｉｉｉ）第３の線群は、ｓ＋２π＜ｑ＜ｓ＋４πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_２、及びｓ−４π＜ｑ＜ｓ−２πによって定義される走査回転の平面ＤＰ（ｓ）上への投影Γ_−２に対して、少なくとも３つの交点ｓ_１，ｓ_２，ｓ_３を有する。
    （ｃ）前記３つの群の前記線に沿って、前記データを前処理及び推移不変フィルタ処理するステップと、
    （ｄ）前記画像の前段階と称すべきもの作成すべく、前記フィルタ処理したデータを逆投影するステップと、
    （ｅ）前記対象物の画像が完成するまで、前記ステップ（ａ），（ｂ），（ｃ）及び（ｄ）を繰り返すステップとを含むことを特徴とする方法。
19. 請求項１８に記載の方法であって、
    前記交点ｓ_１，ｓ_２，ｓ_３は、次の規則に従って決定されることを特徴とする方法。
    ｓ＋２π＜ｓ_３＜ｓ＋４πならば、ｓ_１−ｓ＝Ψ（ｓ_３−ｓ_２）
    ｓ−４π＜ｓ_３＜ｓ−２πならば、ｓ_３−ｓ_２＝Ψ（ｓ_１−ｓ）
    （ただし、Ψ（ｔ）は、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという性質を有する関数である）
20. 請求項１８に記載の方法であって、
    前記前処理は、現在の放射源位置に関するＣＢデータの微分計算を含むことを特徴とする方法。
21. 請求項１８に記載の方法であって、
    前記推移不変フィルタ処理は、前記前処理されたデータを１／ｓｉｎγフィルタにより畳み込むことを含むことを特徴とする方法。
22. 請求項１８に記載の方法であって、
    前記第１の線群からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、Ｌ_２ ^ｃｒとＬ_２ ^−ｃｒとの間に位置する場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝２／３を乗じることを特徴とする方法。
23. 請求項１８に記載の方法であって、
    前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
    ＤＰ（ｓ）へのｘの投影が、Γ_１とΓ_−１との間に位置し、線がΓ_１∪Γ_−１と接する点がｘの１ＰＩパラメータ区間の内側にある場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝２／３を乗じることを特徴とする方法。
24. 請求項１８に記載の方法であって、
    前記第２の線群の線からの前記フィルタ処理されたデータの逆投影は、
    ｘのＤＰ（ｓ）への投影が、Γ_１とΓ_−１との間に位置し、線がΓ_１∪Γ_−１と接する点がｘの１ＰＩパラメータ区間の外側にある場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝２／３を乗じることを特徴とする方法。
25. 請求項１８に記載の方法であって、
    前記第３の線群の線からの前記フィルタ処理されたデータの逆投影は、ｘのＤＰ（ｓ）への投影が、Ｌ_２ ^ｃｒ及びＬ_−２ ^ｃｒの上側に位置する場合に、前記フィルタ処理されたデータに係数ｃ_ｍ＝−２／３を乗じることを特徴とする方法。
26. コンピュータ断層撮影システムにおいて、対象物の走査に基づいて画像復元するのに使用される線群を特定する方法であって、
    （ｉ）現在の放射源位置ｓを決定するステップと、
    （ただし、ｓは、走査経路を表すパラメータである）
    （ｉｉ）コーンビーム投影と交差する平面ＤＰ（ｓ）を選ぶステップと、
    （ｉｉｉ）パラメータ値を用いて、次のｓ_１，ｓ_２，ｓ_３として表される走査経路上の３点を選択するステップと、
    （ｉｉｉａ）｜ｓ−ｓ_１｜＜２π
    （ｉｉｉｂ）２π＜ｓ_２−ｓ＜４π及び２π＜ｓ_３−ｓ＜４π、又は、−４π＜ｓ_２−ｓ＜−２π及び−４π＜ｓ_３−ｓ＜−２π
    （ｉｉｉｃ）２π＜ｓ_３−ｓ＜４πならばｓ_１−ｓ＝Ψ（ｓ_３−ｓ_２）、又は−４π＜ｓ_３−ｓ＜−２πならばｓ_３−ｓ_２＝Ψ（ｓ_１−ｓ）、
    （ただし、Ψ（ｔ）は、Ψ（０）＝０，Ψ´（ｔ）＞０，ｔ∈Ｒという関数である）
    （ｉｖ）前記３点をＤＰ（ｓ）上に投影するステップと、
    （ｖ）前記投影を通る境界線を引くステップとを含むことを特徴とする方法。
    

Images