JP2008541982A - コーンビームct用高速再構成アルゴリズム - Google Patents
コーンビームct用高速再構成アルゴリズム Download PDFInfo
- Publication number
- JP2008541982A JP2008541982A JP2008515352A JP2008515352A JP2008541982A JP 2008541982 A JP2008541982 A JP 2008541982A JP 2008515352 A JP2008515352 A JP 2008515352A JP 2008515352 A JP2008515352 A JP 2008515352A JP 2008541982 A JP2008541982 A JP 2008541982A
- Authority
- JP
- Japan
- Prior art keywords
- grid
- interest
- volume
- detector
- plane
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Withdrawn
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/421—Filtered back projection [FBP]
Abstract
コーンビームコンピュータトモグラフィの高速再構成アルゴリズムが提示される。アルゴリズムは、フィルタリングされる逆投影のタイプであり、検出器面上で規定され一様にサンプリングされる関数と、関連付けられるフーリエ面内の極グリッド上でサンプリングされるこれらの関数のフーリエ変換との間で切り替えを行うために非一様な高速フーリエ変換を使用する。
Description
本発明は、コンピュータトモグラフィ(CT)の分野に関する。特に、本発明は、対象のコーンビーム投影の組から、この対象の関心ボリューム(volume of interest、関心体積)の3次元画像を再構成する方法、画像処理装置、検査装置、コンピュータ可読媒体及びプログラム要素に関する。
コーンビームCTにおいて、X線源が、対象の周囲のX線源軌道に沿って移動する間、対象の関心ボリューム(VOI)のX線コーンビーム投影の組が、取得される。対象のVOIの3次元(3D)画像は、このコーンビーム投影の組から再構成される。投影は、コーンビームCTスキャナによって取得される。このようなコーンビームCTスキャナは、点のようなX線源及び広い領域のX線検出器を備える。一般に、検出器は、平坦であり、小さい検出器素子の2次元(2D)アレイに細分される。このような平坦な検出器は、当然ながら、3D空間における2D平面、つまり検出器面を規定し、かかる検出器面は、検出器にしっかりと取り付けられており、検出器が移動するとともに移動する。コーンビームは、線源から出て、検出器に入るX線によって、形成される。線源が固定の位置にあるときに得られる検出器の2Dアレイの読み取り値は、生のコーンビーム投影を構成する。検出器素子によって読み取られた値は、線源から出てこの検出器素子に入るX線ビームレットの強度を表す。前処理ステップにおいて、検出器の読み取り値は、対象内のX線減衰係数の線積分に変換される。積分の線は、線源の位置及び元の値を読み取った検出器素子によって決定される。以下、特に言及されない限り、検出器読み取り値のこの前処理が理解され、「コーンビーム投影」なる語は、知られている積分のラインに沿った線積分の2Dアレイをさす。取得されたコーンビーム投影の組は、線源軌道に沿って線源位置に関してサンプリングされ、各々のサンプリングされる線源位置ごとに、関連付けられるコーンビーム投影が、検出器面内の検出器素子の位置に関してサンプリングされることが留意されるべきである。概念的には、指定された線源位置に関連付けられるコーンビーム投影は、この線源位置に関連付けられる検出器面上に規定される2つの変数の関数のサンプリングされた評価として観察されることができる。画像は、画像処理装置によってコーンビーム投影の組から再構成される。画像処理装置は、再構成アルゴリズムを実現するコンピュータプログラムを実行するコンピュータである。再構成された画像は、対象のVOI内のX線減衰係数の3Dマップのサンプリングされた評価を提供する。このようなマップは、例えば医用診断又は材料試験のような分野において貴重な情報を提供することができる(コーンビームCTの医療用途では、対象は患者である)。
VOI内の再構成されたマップと真のマップとの間の類似性は、さまざまな要求の充足に依存する。(前処理された)検出器読み取り値が、必要とされる線積分の正確な評価であり、線源軌道に沿った隣り合う線源位置の間の及び検出器面内の隣り合う検出器素子の間のサンプリング距離が、十分に小さいことが要求されることは明らかである。更に、少なくとも対象の理想化された、サンプリングされない、エラーの無いコーンビーム投影は、良好な精度をもって真のマップを決定すべきである。これは、投影が切頭されず、線源軌道が、スキャン中、対象のVOIを含んだボリュームに関して完全である場合にあてはまることが知られている。
この文脈において、投影されている対象が、投射コーンビームによって完全にカバーされない場合、コーンビーム投影は、切頭されるといわれる。
更に、あるボリュームと交差する各々の面が線源軌道とも交差する場合、線源軌道は、このボリュームに関して完全であるといわれる。以下、「ボリュームVに関して完全である」という句は、通常、「完全である」と略され、特定されないボリュームが、文脈から推測されることができる。ある真のボリュームに関して完全であるべき線源軌道は、非平面でなければならないことが留意される。一般的な線源軌道は、円である。この線源軌道は、単一軸を中心とする回転運動により実現するのが簡単であるが、平面であるので完全ではない。
2つの言及される条件が、真のマップを決定するのに十分であるが、絶対に必要というわけではない。例えば、線源軌道が円である場合、円の中心の周りのボリュームにおいて許容できる再構成を得ることがなお可能でありうる。
コーンビームCTスキャナは、さまざまなやり方で実現されることができる。例えば、X線源及びX線検出器は、対象の周りで動かされることが可能なCアームの両端部に取り付けられてもよい。米国特許第6,582,120号明細書に記載されるように、2つの回転運動の適当な組み合わせによって、完全な線源軌道を実現することが可能である。線源及び検出器は、回転ガントリ上に取り付けられてもよい。このようにして、環状の線源軌道を実現することが可能であるが、ガントリが傾斜しうる場合、完全な線源軌道が実現されることもできる。完全な線源軌道を実現する能力をもつガントリの別のタイプは、米国特許第5,124,914号明細書に開示されている。
Cアームに基づくスキャナは、通常、平坦な検出器を具備する。ガントリに基づくスキャナは、平坦な又は非平坦な検出器を有することができる。平坦な検出器は、それに関連付けられる自然な検出器面を有し、通常、小さい検出器素子の2Dアレイに細分される。このような場合、コーンビーム投影の取得された組の各々のコーンビーム投影は、検出器面内の等距離のデカルトグリッド上でサンプリングされる。検出器が平坦でない場合であっても、仮想検出器面を検出器と関連付けることが可能であるとともに、各々のコーンビーム投影がこの仮想検出器面内でサンプリングされるように見せかけることが可能である。このような場合、サンプリング点のグリッドは、平面であるが、概して等距離のデカルトグリッドではない。
コーンビームCTスキャナは、再構成アルゴリズムを実行することによって画像を再構成する画像処理装置を有する。更に、コーンビームCTスキャナは、再構成された画像を表示するためのビューイングコンソールを有する。最後に、コーンビームCTスキャナは、オペレータがスキャナを操作することを可能にする手段を有する。
データが十分良好にマップを決定するという条件で、対象のVOI内のX線減衰係数の所望の3Dマップを再構成することができるさまざまな再構成アルゴリズムが提案されている。
これらの再構成アルゴリズムのうちの多くは、データ取得プロセスの数学的モデルを逆にすることによって見い出される「連続するバージョン」を有する。モデルは、入力関数として対象のX線減衰の3Dマップを受け入れ、コーンビーム投影の組を表す出力関数を決定する、公式又は公式の組からなる。再構成アルゴリズムの連続するバージョンは、このモデルを反転させる:コーンビーム投影の組を表す入力関数が与えられる場合、連続するバージョンは、X線減衰係数の3Dマップを表す別の関数を決定する公式又は公式の組からなる。コンピュータ上で再構成アルゴリズムを実現するために、連続するバージョンの離散化されたバージョンが必要とされる。このような離散化されたバージョンを案出するために、「連続する」関数は、離散的な又はサンプリングされた近似と置き換えられ、公式又は公式の組は、1又は複数の適当な数値アルゴリズムと置き換えられる。最後に、結果的に得られた再構成アルゴリズムの離散的なバージョンは、コーンビームCTスキャナの画像処理装置によって実行されることが可能なコンピュータプログラムとして公式化される。再構成アルゴリズムの離散的なバージョンは、それ自体の「離散的なエラー」を再構成された画像に取り入れる。離散化プロセスは、さまざまなサンプリンググリッドを伴い、かかるグリッドのパターン及び密度は、離散化エラーを小さく保つために適切に選択される必要がある。
各々のコーンビーム投影のサンプリンググリッドは、平面であり、等距離のデカルトグリッドであることが望ましい。検出器が平坦である場合、平面性は自然に得られ、検出器が平坦でない場合、仮想検出器面を採り入れることによって平面性が達成可能である。検出器が平坦である場合であっても、検出器の検知領域を含む面と異なる仮想検出器面を採り入れることができる。平面グリッドが等距離のデカルトグリッドでない場合、必要に応じて、コーンビーム投影を等距離のデカルトグリッド上にリサンプリングすることが可能である。以下、このような任意の種類のリサンプリングが、前処理ステップの一部として考えられる。リサンプリングプロセスは、取得されたコーンビーム投影の精度を低下させることがあるので、可能であれば、リサンプリングを回避することが望ましいことがある。以下、サンプリンググリッドは平面であるものとする。「検出器面」なる語は、仮想検出器面が使用される場合には、その仮想検出器面を意味する。検出器面のサンプリンググリッドは、等距離のデカルトグリッドであってもよく又はそうでなくてもよい。
再構成アルゴリズムが高速であることは、常に望ましく、コーンビームCTの多くのアプリケーションにとって必須である。
再構成アルゴリズムの1つの例は、Defrise及びClackによるコーンビームフィルタード逆投影アルゴリズムであり、かかるアルゴリズムは、M. Defrise及びR. Clackによる「A cone-beam reconstruction algorithm using shift-variant filtering and cone-beam backprojection」(IEEE Trans. Med. Imag., vol. 13, no. 1, pp. 186-195, 1994)に記載されており、参照によって本願明細書に盛り込まれるものする。Defrise及びClackのアルゴリズムは、好ましいアルゴリズム構造を有する:各々の投影は、「フィルタリングされ」、「逆投影され」、後になって捨てられてもよい。Defrise及びClackのアルゴリズムのフィルタリングステップは、やや複雑であり、さまざまな中間関数の順方向及び逆方向の2Dラドン変換を計算することを伴う。これらの中間関数のうちのいくつかの或る偏微分を得ることも必要である。逆投影ステップは、多くのCT再構成アルゴリズムの一般的な部分である。線源軌道が完全でない場合であっても、可能な範囲で、Defrise及びClackアルゴリズムは合理的な結果を与える。切頭された投影は、邪魔なアーチファクトを生じさせることがありうる。しかしながら、拡張された投影が、VOIより少し大きい対象の切頭されていない投影のようにほぼ見えるように、再構成の前に切頭された投影を拡張することによって、実質的にこれらのアーチファクトを低減することが可能である。このような方法は、米国特許第6,542,573号明細書に開示されている。
他方、その複雑なフィルタリングステップは、Defrise及びClackアルゴリズムを、計算的に非常に遅くする。原則として、フィルタリングステップを速めるために、いわゆるフーリエ方法を使用することが可能である。重要なことは、ラドン変換及びフーリエ変換の間の良く知られている関係である。この関係は、いわゆる投影定理によって言い換えられ、かかる投影定理は、2つの変数の関数の2Dラドン変換の行ごとの1Dフーリエ変換が、2Dフーリエ空間における半径方向線に沿ったこの関数の2Dフーリエ変換に等しいことを述べている。更に、前述の偏微分は、フーリエ変換に関するいわゆる微分定理によって、フーリエ空間において実現されることができる。更に、出力関数又は入力関数のいずれも等距離のデカルトグリッドにサンプリングされない、2Dフーリエ変換及びそれらの逆変換を計算するための高速計算方法が必要とされている。適当な態様で、これらの考えと、高速な2Dの順方向及び逆方向の「離散フーリエ変換」とを組み合わせることは、「フーリエフィルタリングされた」逆投影アルゴリズムをもたらす。
1つのこのようなフーリエ変換された逆投影アルゴリズムは、C. Jacobsonによる論文「Fourier methods in 3D reconstruction from cone-beam data」(Linkoping University, Sweden, 1996, pp. 112-113)に示されている。Jacobsonアルゴリズムは、必要な順方向及び逆方向の高速2D離散フーリエ変換のために、いわゆるリノグラム(linogram)技法を使用する。Linogram技法は、例えばP. R. Edholm、G. T. Herman及びD. A. Robertsによる論文「Image reconstruction from linograms: Implementation and evaluation」(IEEE Trans. Medical Imaging, vol. 7, pp. 239-246, 1988)及びM. Magnussonによる論文「Linogram and Other Direct Fourier Methods for Tomographic Reconstruction」(Dissertation No. 320, Linkoping University, S-58183 Linkoping, Sweden, 1993)に記載されている。リノグラム技法は、2Dフーリエ空間において特別なグリッドの使用を必要とし、そのグリッドのグリッド点は、同心正方形上に位置付けられる。これは、そのようなグリッドが、必要な2D離散フーリエ変換を計算するための相対的に高速な手段として、いわゆるチャープz変換を呼び出すことを可能にするからである。チャープz変換は、例えばMagnussonの論文に記述されている。残念なことに、同心正方形のグリッドは、或る角度異方性及びラジアル異方性を、再構成された画像に取り込む。結果として生じる離散化エラーを低減するために、同心正方形グリッドのグリッド点の相対的に高い密度を選択することが強いられ、このことは、グリッド点の数及びゆえに計算時間を増加させる。
Defrise及びClackのアルゴリズムの利点を維持しながら、Jacobsonのアルゴリズムの不利益を回避するコーンビームCTのための再構成アルゴリズムを有することが望ましい。
本発明の例示的な実施例によれば、線源軌道に沿った複数の線源位置から取得される対象のコーンビーム投影の組から、対象のVOIの3D画像を再構成する方法であって、フーリエ空間における同心円グリッドに基づいて、投影の組をフィルタリングして、フィルタリングされた投影データセットを生成するステップと、フィルタリングされた投影データセットから関心ボリュームを再構成して、関心ボリュームの再構成された画像を生成するステップと、を含む方法が提供されることができる。
本発明のこの例示的な実施例によれば、2Dフーリエ空間におけるグリッドが使用され、そのグリッド点は、正方形ではなく同心円上に位置する。
本発明の別の例示的な実施例によれば、必要とされる順方向及び逆方向の2D離散フーリエ変換は、いわゆる非一様な高速フーリエ変換を使用して実現される。
「非一様な高速フーリエ変換」(nonuniform FFT)なる語は、出力グリッド、入力グリッド又はそれらの両方が、等距離のデカルトグリッドでないとき、(三角和によって表される)離散フーリエ変換を計算するために使用されうる方法のクラスの総称としてここで使用される。関係するグリッドタイプに依存して、一様−非一様なFFT、非一様−一様なFFT、及び非一様−非一様なFFTの間で区別がなされることが可能である。すべての3つのケースにおいて、フーリエ変換又は逆フーリエ変換のいずれが計算されることができるかに依存して、「順方向のバージョン」及び「逆方向のバージョン」がある。非一様なFFTは、例えばA. Fessler及びB. P. Suttonによる論文「Nonuniform fast Fourier transforms using min-max interpolation」(IEEE Trans. Signal Processing, vol. 14, no. 3, pp 560-574, 2003)に記述されている。この論文は、密接に関連する方法、特に「グリッドを作る方法」のクラスについて他の文献にも言及している。関係する両方のグリッドが、等距離のデカルトグリッドである場合、一様−一様なFFTとしてここに示される古典的な高速フーリエ変換を使用することができる。
本発明の別の例示的な実施例によれば、対象のコーンビーム投影の投影の組からこの対象のVOIを再構成するための画像処理装置が、提供されることができる。画像処理装置は、計算ユニット及びデータを記憶するためのメモリを有する。計算ユニットは、本発明の例示的な実施例による方法を実施するように適応される。
本発明の別の例示的な実施例によれば、対象の投影の組から画像を再構成するための検査装置であって、上述の方法ステップを実施するように適応される計算ユニットを有する検査装置が、提供されることができる。
本発明の別の例示的な実施例によれば、検査装置は更に、電磁放射線源及び広領域の検出器を有し、放射線源は、放射線源軌道に沿った線源位置の組から関心対象に電磁放射線を発するように適応され、検出器は、それが対象のVOIのコーンビーム投影を測定するように配置され、適応される。
更に、本発明の別の例示的な実施例によれば、検査装置は、手荷物検査装置、医用診断機器、材料試験装置及び材料科学解析装置を含むグループのうちの1つとして構成される。本発明の応用分野は、手荷物検査でありうる。これは、本発明の規定された機能が、手荷物アイテムの中身の安全な、信頼できる、高速な解析を可能にしうるので、それにより疑わしい中身を検出するとともに、このような手荷物アイテム内の材料のタイプさえも決定することを可能にする。
本発明の別の例示的な実施例によれば、検査装置によって関心対象の投影の組から画像を再構成するためのコンピュータプログラムであって、プロセッサによって実行されるとき、上述の方法ステップを実行するように適応されるコンピュータプログラムが記憶されているコンピュータ可読媒体が提供されることができる。
本発明は更に、関心対象の投影データセットから画像を再構成するためのプログラム要素であって、プロセッサによって実行されるとき、上述の方法ステップを実行するように適応されるプログラム要素に関する。プログラム要素は、好適には、データプロセッサの作業メモリへロードされることができる。こうして、データプロセッサは、本発明の方法の例示的な実施例を実行する能力を備えることができる。コンピュータプログラムは、例えばC++のような任意の適切なプログラミング言語で書かれることができ、例えばCD−ROMのようなコンピュータ可読媒体に記憶されることができる。更に、コンピュータプログラムは、例えばワールドワイドウェブのようなネットワークから入手可能でありえ、かかるネットワークから、画像処理装置、プロセッサ又は任意の適切なコンピュータにダウンロードされることができる。
本発明のこれら及び他の見地は、以下に記述される実施例から明らかになり、それらの実施例を参照して解明される。
本発明の例示的な実施例は、添付の図面を参照して以下に記述される。
コーンビームCTスキャナ
例示的なコーンビームCTスキャナの概略図が、図1に示されている。X線源100及び広い検知領域を有する平坦な検出器(フラットディテクタ)101が、Cアーム102の両端部に取り付けられている。Cアーム102は、曲がったレールである「スリーブ」103によって保持される。Cアームは、スリーブ103内を摺動することができ、それにより、Cアームの軸を中心に「ロール運動」を実施する。スリーブ103は、回転ジョイント部を介してLアーム104に取り付けられ、このジョイント部の軸を中心に「プロペラ運動」を実施することができる。Lアーム104は、別の回転ジョイント部を介して天井に取り付けられ、このジョイント部の軸を中心に回転を実施することができる。さまざまな回転運動は、サーボモータによって達成される。3つの回転運動の軸及びコーンビーム軸は、単一の固定点、すなわちコーンビームCTスキャナの「アイソセンタ」105において常に交わる。線源軌道に沿ってすべてのコーンビームによって投影されるボリュームが、アイソセンタ周囲にある。この「投影ボリューム」(VOP)の形状及び大きさは、検出器の形状及び大きさ並びに線源軌道に依存する。図1において、球110は、VOP内に収まる最大のアイソセントリックな球を示している。イメージングされるべき対象(例えば患者又は手荷物アイテム)は、対象のVOIがVOPを満たすように、テーブル111に配置される。対象が十分に小さい場合、その対象は完全にVOP内に収まる;対象が十分に小さくない場合、その対象は完全にはVOPに収まらない。従って、VOPは、VOIの大きさを制限する。
例示的なコーンビームCTスキャナの概略図が、図1に示されている。X線源100及び広い検知領域を有する平坦な検出器(フラットディテクタ)101が、Cアーム102の両端部に取り付けられている。Cアーム102は、曲がったレールである「スリーブ」103によって保持される。Cアームは、スリーブ103内を摺動することができ、それにより、Cアームの軸を中心に「ロール運動」を実施する。スリーブ103は、回転ジョイント部を介してLアーム104に取り付けられ、このジョイント部の軸を中心に「プロペラ運動」を実施することができる。Lアーム104は、別の回転ジョイント部を介して天井に取り付けられ、このジョイント部の軸を中心に回転を実施することができる。さまざまな回転運動は、サーボモータによって達成される。3つの回転運動の軸及びコーンビーム軸は、単一の固定点、すなわちコーンビームCTスキャナの「アイソセンタ」105において常に交わる。線源軌道に沿ってすべてのコーンビームによって投影されるボリュームが、アイソセンタ周囲にある。この「投影ボリューム」(VOP)の形状及び大きさは、検出器の形状及び大きさ並びに線源軌道に依存する。図1において、球110は、VOP内に収まる最大のアイソセントリックな球を示している。イメージングされるべき対象(例えば患者又は手荷物アイテム)は、対象のVOIがVOPを満たすように、テーブル111に配置される。対象が十分に小さい場合、その対象は完全にVOP内に収まる;対象が十分に小さくない場合、その対象は完全にはVOPに収まらない。従って、VOPは、VOIの大きさを制限する。
さまざまな回転運動が、制御ユニット120によって制御される。Cアーム角度、スリーブ角度及びLアーム角度の各々の三つ組が、X線源の位置を規定する。時間と共にこれらの角度を変えることによって、線源は、定められた線源軌道に沿って動かされることができる。Cアームの他端の検出器は、対応する運動を行う。線源軌道は、アイソセントリックな球の表面に制限される。
図2は、検出器101の検知領域を示している。検知領域は、矩形であり、等しい大きさの検出器素子の2Dアレイに細分される。検出器102の出力は、データ要素の対応するアレイである。或る線源位置に関連付けられるデータのアレイは、生のコーンビーム投影を構成する。線源が定められた線源軌道に沿って移動する間、生のコーンビーム投影の組が取得される。これらの生のコーンビーム投影に関連付けられる線源位置は知られている。生コーンビーム投影は、前処理装置130へ転送され、前処理装置130において前処理される。前処理されたコーンビーム投影は、画像処理装置140に供給される。画像処理装置140は、再構成アルゴリズムを実現するコンピュータプログラムを実行するコンピュータである。投影が、切頭されており、再構成の前に拡張されるべきである場合、画像処理装置140は更に、投影の拡張をも実施する。再構成された画像は、コンソール150のモニタ151に表示される。コンソール151は、コーンビームCTスキャナを制御するためにも使用され、オペレータにユーザインタフェースを提示する。
コーンビームCTスキャナは、例えばX線管用高電圧発生器、X線管用冷却手段及び電気ケーブルのような他の補助的装置によって補われる。このような補助的装置は、図に示されていない。
数学的変換
本発明の例示的な実施例の詳細な説明は、良く知られているいくつかの数学的変換を利用する。これらの変換は、このサブセクションにおいて規定される。これらの変換のいくつかの良く知られている特性についても記述される。
本発明の例示的な実施例の詳細な説明は、良く知られているいくつかの数学的変換を利用する。これらの変換は、このサブセクションにおいて規定される。これらの変換のいくつかの良く知られている特性についても記述される。
関数g:Rd→Cのフーリエ変換は、
によって規定される。関数h:R→Cのヒルベルト変換は、
(F1H1h)(σ)=i sign(σ)(F1h)(σ)
によって規定される。微分演算子D1は、
(F1D1h)(σ)=i 2πσ(F1h)(σ)
によって規定される。慣例により、F1、F1 −1、H1又はD1が、2つ以上の変数の関数に適用される場合、それらは、最初の変数に働く。関数g:R2→Cのラドン変換は、
によって規定される。関数g:R×[0,π]→Rの逆投影は、
によって規定される。反転公式2πR2 −1=−B2D1H1は、恒等式2πR2 −1H1=B2D1を与える。R2の投影定理は、
(F1R2g)(σ,θ)=(F2g)(σcosθ,σsinθ)
を述べる。R2のシフト定理は、
(R2ga,b)(s,θ)=(R2g)(s−a cosθ−b sinθ,θ)
を述べる。上式にて、
ga,b(u,v)=g(u−a,v−b)
である。慣例により、F2、F2 −1、R2、R2 −1又はB2が、3つ以上の変数の関数に適用される場合、それらは、最初の2つの変数に働く。
によって規定される。関数h:R→Cのヒルベルト変換は、
(F1H1h)(σ)=i sign(σ)(F1h)(σ)
によって規定される。微分演算子D1は、
(F1D1h)(σ)=i 2πσ(F1h)(σ)
によって規定される。慣例により、F1、F1 −1、H1又はD1が、2つ以上の変数の関数に適用される場合、それらは、最初の変数に働く。関数g:R2→Cのラドン変換は、
によって規定される。関数g:R×[0,π]→Rの逆投影は、
によって規定される。反転公式2πR2 −1=−B2D1H1は、恒等式2πR2 −1H1=B2D1を与える。R2の投影定理は、
(F1R2g)(σ,θ)=(F2g)(σcosθ,σsinθ)
を述べる。R2のシフト定理は、
(R2ga,b)(s,θ)=(R2g)(s−a cosθ−b sinθ,θ)
を述べる。上式にて、
ga,b(u,v)=g(u−a,v−b)
である。慣例により、F2、F2 −1、R2、R2 −1又はB2が、3つ以上の変数の関数に適用される場合、それらは、最初の2つの変数に働く。
高速な一様及び非一様な離散フーリエ変換
このサブセクションは、離散フーリエ変換及び非一様な高速フーリエ変換について背景情報を提供する。
このサブセクションは、離散フーリエ変換及び非一様な高速フーリエ変換について背景情報を提供する。
{xk∈Rd:k∈K}及び{yj∈Rd:j∈J}を、Rdにおける2つの任意の有限グリッドとする。ある関数g:Rd→Cのサンプルg(xk),k∈Kが与えられており、関数、
のサンプル、
を計算しようとしているものとする。この問題を解決するための一般的な方法は、適切な直交係数ckに関する三角和、
によって、
を近似し、
を、
の近似として得ることでる。点広がり関数(PSF)、
に関して、
であることに留意されたい。サンプル、
は、サンプルg(xk),k∈Kの離散フーリエ変換(DFT)を構成する。サンプル、
からサンプルg(xk),k∈Kを計算する逆の問題は、
の三角和によってgを近似し、サンプルg*(xk),k∈Kを計算することによって、解決されることができ、これは、サンプル、
の逆DFTを構成する。PSF、
についてg*=p*gであることに留意されたい。(逆)DFTの入力グリッドが、等距離のデカルトグリッドであり、出力グリッドが、このようなグリッドでない場合、一様−非一様変換(U−NU)DFTという。同様に、残りのケースを、NU−U DFT、NU−NU DFT及びU−U DFTと称する。4つのすべてのケースにおいて、順方向(Fdgの計算)又は逆方向(
の計算)の計算が行われることが可能である。入力グリッドが、a∈R+ d及びK∈Z+ dのとき、式{ka:−K/2≦k≦K/2}の等距離のデカルトグリッドであり、出力グリッドが、b=1/(Ka)のとき、共役のグリッド{nb:−K/2≦n<K/2}である場合、U−U DFTは、本質的に古典的なDFTである。(ka及び1/(Ka)のような表現は、成分に関して理解される。)更に、Kに対する適度な制約下において、U−U DFTは、ここでU−U FFTと呼ばれる古典的な高速フーリエ変換(FFT)を使用して、非常に効率的に計算されることができる。NU−U FFT及びU−NU FFTとここで呼ばれる、NU−U DFT及びU−NU DFTを計算するための高速アルゴリズムが、より最近になって利用できるようになっている。これらのアルゴリズムは更に、NU−NU FFTを組み立てるために使用されることができる。例えば、順方向NU−NU FFTは、順方向NU−U FFT、逆方向U−U FFT及び順方向U−NU FFTをつなぎ合わせることによって得られる。非一様なFFTなる語は、総称的に、U−NU FFT、NU−U FFT及びNU−NU FFTについて使用される。フーリエ変換を離散フーリエ変換と区別するために、前者は「連続するフーリエ変換」と称されることがある。
のサンプル、
を計算しようとしているものとする。この問題を解決するための一般的な方法は、適切な直交係数ckに関する三角和、
によって、
を近似し、
を、
の近似として得ることでる。点広がり関数(PSF)、
に関して、
であることに留意されたい。サンプル、
は、サンプルg(xk),k∈Kの離散フーリエ変換(DFT)を構成する。サンプル、
からサンプルg(xk),k∈Kを計算する逆の問題は、
の三角和によってgを近似し、サンプルg*(xk),k∈Kを計算することによって、解決されることができ、これは、サンプル、
の逆DFTを構成する。PSF、
についてg*=p*gであることに留意されたい。(逆)DFTの入力グリッドが、等距離のデカルトグリッドであり、出力グリッドが、このようなグリッドでない場合、一様−非一様変換(U−NU)DFTという。同様に、残りのケースを、NU−U DFT、NU−NU DFT及びU−U DFTと称する。4つのすべてのケースにおいて、順方向(Fdgの計算)又は逆方向(
の計算)の計算が行われることが可能である。入力グリッドが、a∈R+ d及びK∈Z+ dのとき、式{ka:−K/2≦k≦K/2}の等距離のデカルトグリッドであり、出力グリッドが、b=1/(Ka)のとき、共役のグリッド{nb:−K/2≦n<K/2}である場合、U−U DFTは、本質的に古典的なDFTである。(ka及び1/(Ka)のような表現は、成分に関して理解される。)更に、Kに対する適度な制約下において、U−U DFTは、ここでU−U FFTと呼ばれる古典的な高速フーリエ変換(FFT)を使用して、非常に効率的に計算されることができる。NU−U FFT及びU−NU FFTとここで呼ばれる、NU−U DFT及びU−NU DFTを計算するための高速アルゴリズムが、より最近になって利用できるようになっている。これらのアルゴリズムは更に、NU−NU FFTを組み立てるために使用されることができる。例えば、順方向NU−NU FFTは、順方向NU−U FFT、逆方向U−U FFT及び順方向U−NU FFTをつなぎ合わせることによって得られる。非一様なFFTなる語は、総称的に、U−NU FFT、NU−U FFT及びNU−NU FFTについて使用される。フーリエ変換を離散フーリエ変換と区別するために、前者は「連続するフーリエ変換」と称されることがある。
すべての上述の一様及び非一様なFFTにより、入力データ又は出力データが実数である場合、計算時間をほぼ半分にすることが可能である。
再構成問題
線源軌道は、有限数の滑らかなセグメントからなるべきである。説明を簡単にするため、線源軌道が、単一セグメントのみからなるものとする。
線源軌道は、有限数の滑らかなセグメントからなるべきである。説明を簡単にするため、線源軌道が、単一セグメントのみからなるものとする。
右手デカルトx−y−z座標系を、対象テーブルに結び付けることによって、物理空間におけるポイントとR3におけるポイントとの間の一致が確立される。この座標系は、固定座標系と称される。線源軌道は、滑らかなマッピングa:Λ→R3によって表現される。ここで、Λ={λmin,λmax}⊂Rは、間隔である。λmin≦λ0<・・・<λK−1≦λmaxであって、線源が、位置a(λk),0≦k<Kにあるとき、コーンビーム投影が取得される。
検出器の検知領域は、幅W及び高さHの矩形である。検出器素子のアレイは、I列及びJ行を有する。説明を簡単にするため、I及びJは、等しいものとする(必要に応じて、ダミー行又は列を加えることによってI又はJを等しくすることも可能である)。
検出器の検知領域は、3D空間における2D面、つまり検出器面に属し、かかる検出器面は、検出器が動くと共に動く。(より一般的には、仮想検出器面が認められる。u軸及びv軸が、検出器素子によって形成される行及び列と平行であり、w軸が、検出器の前の半空間を指すように、右手デカルトu−v−w座標系が検出器に結びつけられる。この検出器座標系の起点は、検出器素子の中心点が、座標(ui,vi,0)=(iΔu,jΔv,0),−I/2≦i<I/2,−J/2≦j<J/2を得るように、選択される。ここで、(Δu,Δv)=(W/I,H/J)は、検出器素子の大きさである。(−I/2≦i<I/2,−J/2≦j<J/2のような同様のインデックスレンジでもよい。)起点のこの選択は、図2に示されている。検知領域のポイントのu座標及びv座標は、矩形D0⊂R2を構成する。より一般的には、この組は、λに依存し、D0に代わってD0(λ)と書く。
固定座標系において、検出器座標系の起点の軌道は、マッピングd:Λ→R3によって表される。u軸及びv軸の向きは、2つの他のマッピング、
によって表される。マッピング、
及び組D0は、取得ジオメトリを規定する。取得ジオメトリは、図3に示されており、知られているものとする。線源ポイントa(λ)に関連付けられる検出器面は、
によって与えられる。あらゆるp∈Δ(λ)は、d(λ)を中心とする検出器座標(u,v)に関して、
の形で書くことができる。各々のλ∈Λについて、c(λ)∈R3を、a(λ)のΔ(λ)への正射影であるものとする。あらゆるp∈Δ(λ)は、c(λ)を中心とする検出器座標、
に関して、
の形で書くことができる。c(λ)のd(λ)を中心とする検出器座標は、取得ジオメトリから計算可能であり、(uc(λ),vc(λ))によって示される。
によって表される。マッピング、
及び組D0は、取得ジオメトリを規定する。取得ジオメトリは、図3に示されており、知られているものとする。線源ポイントa(λ)に関連付けられる検出器面は、
によって与えられる。あらゆるp∈Δ(λ)は、d(λ)を中心とする検出器座標(u,v)に関して、
の形で書くことができる。各々のλ∈Λについて、c(λ)∈R3を、a(λ)のΔ(λ)への正射影であるものとする。あらゆるp∈Δ(λ)は、c(λ)を中心とする検出器座標、
に関して、
の形で書くことができる。c(λ)のd(λ)を中心とする検出器座標は、取得ジオメトリから計算可能であり、(uc(λ),vc(λ))によって示される。
対象のX線減衰係数は、まだ知られていない関数f:R3→Rによって表される。対象の(理想的な)コーンビーム投影は、関数、
によって表される。ここで、
は、a(λ)から
までを指す単位ベクトルである。gの最初の2つの変数は、d(λ)を中心とする検出器座標であることに留意されたい。「測定された」ドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}に対するgallの制約は、gmによって示される。更に、g:R2×Λ→Rは、ドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}においてgmに適合し、このドメインの外側のgallを評価する関数とする。例えば、gは、米国特許第6,542,573号明細書に記述される方法を使用して、測定されたデータgmを拡張することによって得ることが可能である。測定されたドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}においてg=gall=gmであることに留意されたい。λ∈Λ及び(u,v)∈R2\D0(λ)であって、gall(u,v,λ)≠0である場合、(理想的な)投影は切頭されているといわれる。データ取得プロセス及び取得された生データののちの前処理は、サンプルg(ui,vj,λk),−I/2≦i<I/2,−J/2≦j<J/2,0≦k<K(の評価)を提供する。gのこのサンプリングされたバージョン並びに取得ジオメトリ及びサンプリング点(ui,vj,λk)の知識を使用して、VOPのいくつかのサブセットVにおいてfを再構成することが望まれる。
によって表される。ここで、
は、a(λ)から
までを指す単位ベクトルである。gの最初の2つの変数は、d(λ)を中心とする検出器座標であることに留意されたい。「測定された」ドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}に対するgallの制約は、gmによって示される。更に、g:R2×Λ→Rは、ドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}においてgmに適合し、このドメインの外側のgallを評価する関数とする。例えば、gは、米国特許第6,542,573号明細書に記述される方法を使用して、測定されたデータgmを拡張することによって得ることが可能である。測定されたドメイン{(u,v,λ):λ∈Λ,(u,v)∈D0(λ)}においてg=gall=gmであることに留意されたい。λ∈Λ及び(u,v)∈R2\D0(λ)であって、gall(u,v,λ)≠0である場合、(理想的な)投影は切頭されているといわれる。データ取得プロセス及び取得された生データののちの前処理は、サンプルg(ui,vj,λk),−I/2≦i<I/2,−J/2≦j<J/2,0≦k<K(の評価)を提供する。gのこのサンプリングされたバージョン並びに取得ジオメトリ及びサンプリング点(ui,vj,λk)の知識を使用して、VOPのいくつかのサブセットVにおいてfを再構成することが望まれる。
再構成アルゴリズム
線源軌道がVに関して完全であり、投影が切頭されていない場合、fは、知られている正確なコーンビーム再構成アルゴリズムのうちの1つによって、データ内のエラーつまりデータのサンプリングされた性質によって引き起こされるエラー及び再構成アルゴリズムによって導入されるエラーに耐えて、Vにおいて再構成されることができる。1つの候補の再構成アルゴリズムは、ここでCBFBPと称されるDefrise及びClackのコーンビームフィルタード逆投影アルゴリズムである。残念なことに、CBFBPのフィルタリングステップは、計算的に非常に遅い。以下、この欠点を克服するCBFBPの変更されたバージョンが提示される。これは、変更されたバージョンが一様及び非一様なFFTを使用して実現されることができるように、フィルタリングステップを再公式化することによって達成される。このコーンビームフーリエフィルタード逆投影アルゴリズムは、CBFFBPと称される。
線源軌道がVに関して完全であり、投影が切頭されていない場合、fは、知られている正確なコーンビーム再構成アルゴリズムのうちの1つによって、データ内のエラーつまりデータのサンプリングされた性質によって引き起こされるエラー及び再構成アルゴリズムによって導入されるエラーに耐えて、Vにおいて再構成されることができる。1つの候補の再構成アルゴリズムは、ここでCBFBPと称されるDefrise及びClackのコーンビームフィルタード逆投影アルゴリズムである。残念なことに、CBFBPのフィルタリングステップは、計算的に非常に遅い。以下、この欠点を克服するCBFBPの変更されたバージョンが提示される。これは、変更されたバージョンが一様及び非一様なFFTを使用して実現されることができるように、フィルタリングステップを再公式化することによって達成される。このコーンビームフーリエフィルタード逆投影アルゴリズムは、CBFFBPと称される。
元のバージョンのように、変更されたバージョンは、冗長性補償関数(RCF)を伴う。RCFは、線源ポイントを含み、その線源ポイントに関連付けられる検出器面と交差する、面の組において規定される。Pがこのような面であり、a(λ)がその面内の線源ポイントであるとする。PのΔ(λ)との交差は、線であり、この線のまさに1つのポイントが、c(λ)に最も近い。この最も近いポイントは、c(λ)を中心とする極の検出器座標、
に関して、
の形で書くことができる。このようにして、RCFは、三つ組パラメータ、
の関数、
になる。RCFは、或る正規化条件を満たさなければならないが、線源軌道によってユニークに決定されるわけではない。適切なRCFが、Defrise及びClackの引用された論文に提案された。このRCFは、F. Noo、M. Defrise及びR. Clackによる論文「Direct reconstruction of cone-beam data acquired with a vertex path containing a circle」(IEEE Trans. Image Processing, vol. 40, no. 4, pp. 1092-1101, August 1998)に記述される方法を使用して事前に計算されることができる。c(λ)を中心とするデータ、
に関して作用する元のバージョンとは異なり、変更されたバージョンは、d(λ)を中心とするデータg(u,v,λ)に関して作用する。アルゴリズムは更に、線源軌道の微分a'(λ)、線源と検出器との間の距離w(λ)=‖a(λ)−c(λ)‖、及びベクトル、
を伴う。ここで、
である。CBFFBPの連続するバージョンは、以下の通りである:
に関して、
の形で書くことができる。このようにして、RCFは、三つ組パラメータ、
の関数、
になる。RCFは、或る正規化条件を満たさなければならないが、線源軌道によってユニークに決定されるわけではない。適切なRCFが、Defrise及びClackの引用された論文に提案された。このRCFは、F. Noo、M. Defrise及びR. Clackによる論文「Direct reconstruction of cone-beam data acquired with a vertex path containing a circle」(IEEE Trans. Image Processing, vol. 40, no. 4, pp. 1092-1101, August 1998)に記述される方法を使用して事前に計算されることができる。c(λ)を中心とするデータ、
に関して作用する元のバージョンとは異なり、変更されたバージョンは、d(λ)を中心とするデータg(u,v,λ)に関して作用する。アルゴリズムは更に、線源軌道の微分a'(λ)、線源と検出器との間の距離w(λ)=‖a(λ)−c(λ)‖、及びベクトル、
を伴う。ここで、
である。CBFFBPの連続するバージョンは、以下の通りである:
1.フィルタリング:各々のλ∈Λについて、以下の関数を計算する:
2.逆投影:以下を計算する:
上式にて、u(r,λ)及びv(r,λ)は、a(λ)からΔ(λ)までのrの透視投影のd(λ)を中心とする検出器座標である。
2.逆投影:以下を計算する:
上式にて、u(r,λ)及びv(r,λ)は、a(λ)からΔ(λ)までのrの透視投影のd(λ)を中心とする検出器座標である。
CBFFBPの連続するバージョン及びCBFBPが等価であることを見い出すのは容易である:R2に関する投影定理及びD1のフーリエ表現を使用して、上式(3a)−(5)は、
h2=R2g1,h3=D1h2 (11)
と等価であることが分かる。
h2=R2g1,h3=D1h2 (11)
と等価であることが分かる。
H1のフーリエ表現及びR2に関する投影定理を使用して、更に、上式(7)−(9b)が、
h5=2πH1h4,g6=R2 −1h5 (12)
と等価であることが分かる。
h5=2πH1h4,g6=R2 −1h5 (12)
と等価であることが分かる。
こうして、式(3a)−(5)の代わりに式(11)を使い、式(7)−(9b)の代わりに式(12)を使うことによって、ここでCBFBP2によって示される等価なアルゴリズムに至る。恒等式2πR2 −1H1=B2D1及びR2に関するシフト定理を使用して、最後に、CBFBP2が、元のバージョンCBFBPと等価であることは確認することは単純明快である。実際に、これらの3つのアルゴリズムの等価性は、そのu軸及びv軸が検出器面内にある限り、検出器座標系の起点及び向きから独立している。
CBFFBPのフィルタリングステップを数値的に実現するため、関係のあるさまざまな関数について、適当なサンプリンググリッドを決める必要がある。入力データ及びすべての中間関数は、第3の変数としてλに依存する。この変数について、グリッド{λk:0≦k<K}が採用される。最初の2つの変数は、考慮中の関数に依存する。関数g、g1及びg6については、再構成問題のサブセクションにて取り込まれたグリッドが採用される。d(λ)を中心とするデータ(g)は、直接、このグリッド上で得られる。(c(λ)を中心とするデータ、
は一般にこのグリッド上で得られない)。関数h3及びh4の場合、
の形の標準の「サイノグラム」グリッドが採用される。ここで、
であり、例えばNrad及び3Nrad/2の間のあるNangについて、
NangΔθ=π
である。関数、
について、共役のグリッド、
を使用する。ここで、σμ=μΔσ,Δσ=1/(NradΔs)である。関数、
について、関連付けられる極グリッド、
が使用される。このような極グリッドの一例は、図4に示されている。グリッド点は、半径方向の線及び同心円の交点に位置する。
は一般にこのグリッド上で得られない)。関数h3及びh4の場合、
の形の標準の「サイノグラム」グリッドが採用される。ここで、
であり、例えばNrad及び3Nrad/2の間のあるNangについて、
NangΔθ=π
である。関数、
について、共役のグリッド、
を使用する。ここで、σμ=μΔσ,Δσ=1/(NradΔs)である。関数、
について、関連付けられる極グリッド、
が使用される。このような極グリッドの一例は、図4に示されている。グリッド点は、半径方向の線及び同心円の交点に位置する。
フィルタリングステップは、各々のλk,0≦k<Kについて実施される。サブステップ(2)、(4)、(6)及び(8)の数値的な実現は、直接的である。残りのサブステップについては、高速な一様及び非一様なDFTを用いる。具体的には、式(3a)及び(3b)を介してg1から、
を計算するために、一定の直交係数cij=ΔuΔvに関して2次元の順方向U−NU DFT/FFTが使用される。式(9a)及び(9b)を介して、
からg6を計算するために、直交係数、
に関して2D逆方向NU−U DFT/FFTが使用される。再構成問題サブセクションにおいて選ばれた特別な検出器座標系は、特に効率的な実現を可能にする。サブステップ(5)及び(7)について、一定の直交係数Δσ及びΔsに関する複数の1次元逆方向及び順方向U−U DFT/FFTが使用される。さまざまなDFTに関連付けられるサンプリンググリッドが十分に微細な場合、結果として得られるPSFは、振る舞いが良く、わずかなエイリアシング及び切頭エラーのみを生じさせる。すべてのこれらのFFTについて、入力又は出力データは実数でありえ、これは、計算の負担をほぼ半分に低減することを可能にする。
を計算するために、一定の直交係数cij=ΔuΔvに関して2次元の順方向U−NU DFT/FFTが使用される。式(9a)及び(9b)を介して、
からg6を計算するために、直交係数、
に関して2D逆方向NU−U DFT/FFTが使用される。再構成問題サブセクションにおいて選ばれた特別な検出器座標系は、特に効率的な実現を可能にする。サブステップ(5)及び(7)について、一定の直交係数Δσ及びΔsに関する複数の1次元逆方向及び順方向U−U DFT/FFTが使用される。さまざまなDFTに関連付けられるサンプリンググリッドが十分に微細な場合、結果として得られるPSFは、振る舞いが良く、わずかなエイリアシング及び切頭エラーのみを生じさせる。すべてのこれらのFFTについて、入力又は出力データは実数でありえ、これは、計算の負担をほぼ半分に低減することを可能にする。
関数h3及びh4についてのサイノグラムグリッドに代わって、関数、
について「リノグラム」グリッド及び対応する共役グリッドを使用したくなるかもしれない。(「リノグラム」の概念についてEdholm他の引用される論文又はMagnussonの引用される論文を参照されたい。)関数、
は、同心円ではなく同心正方形上でサンプリングされる。このような同心正方形グリッドは、図5に示されている。順方向U−NU DFTの直交係数は変化しない。逆方向NU−U DFTの直交係数は変更されなければならないが、適切な係数が見つけられることができる。しかしながら、逆方向NU−U DFTのPSFが十分に良く振る舞うように、正方形の対角線に沿ったそのグリッド点の密度が、同心円の場合のグリッド点の半径の密度と同じくらい高くなるように、同心正方形グリッドのパラメータを調整する必要がある。これは、同心正方形グリッドのグリッド点の総数を増加させ、ゆえに関連付けられる非一様なFFTの計算の負担を増加させる。代替例として、チャープz変換を使用して通常のリノグラムのやり方で同じ非一様なDFTを計算することができる。しかしながら、演算数は、これらのアルゴリズムが、それらの非一様なFFTの対応するものより速くないことを明らかにする。
について「リノグラム」グリッド及び対応する共役グリッドを使用したくなるかもしれない。(「リノグラム」の概念についてEdholm他の引用される論文又はMagnussonの引用される論文を参照されたい。)関数、
は、同心円ではなく同心正方形上でサンプリングされる。このような同心正方形グリッドは、図5に示されている。順方向U−NU DFTの直交係数は変化しない。逆方向NU−U DFTの直交係数は変更されなければならないが、適切な係数が見つけられることができる。しかしながら、逆方向NU−U DFTのPSFが十分に良く振る舞うように、正方形の対角線に沿ったそのグリッド点の密度が、同心円の場合のグリッド点の半径の密度と同じくらい高くなるように、同心正方形グリッドのパラメータを調整する必要がある。これは、同心正方形グリッドのグリッド点の総数を増加させ、ゆえに関連付けられる非一様なFFTの計算の負担を増加させる。代替例として、チャープz変換を使用して通常のリノグラムのやり方で同じ非一様なDFTを計算することができる。しかしながら、演算数は、これらのアルゴリズムが、それらの非一様なFFTの対応するものより速くないことを明らかにする。
図6は、結果として得られるCBFFBPの離散的なバージョンのフローチャートを示している。
投影が切頭される場合、再構成の前に、好ましくは米国特許第6,542,573号明細書に開示される方法を使用することによって、投影を拡張することが望ましい。
線源軌道が円である場合、CBFBPの連続するバージョンが、FDKアルゴリズムの連続するバージョンになることは良く知られている。FDKアルゴリズムは、環状の線源軌道の標準コーンビームCT再構成アルゴリズムであり、L. A. Feldkamp、L.C. Davis及びW.J. Kressによる「Practical cone-beam algorithm」(J. Opt. Soc. Amer. A, vol. 1, no. 6, pp. 612-619, 1984)に記載されている。CBFFBPの連続するバージョン及びCBFBPは等価であるので、対応する記述が、CBFFBPについても当てはまる。
デモンストレーション
アルゴリズムCBFFBPは、パーソナルコンピュータ上で「C」プログラミング言語で実現され、シミュレートされたデータを用いてテストされた。FDKアルゴリズムの匹敵する実現との比較も行われた。2つの線源軌道が試みられた。これらのうちの1つは、
によって規定される「ねじれた円」atc:[0,2π]→R3であった。ここで、rsrc=810mmであり、α=15π/180の場合θ(λ)=αcos2λである。図7(a)、(b)及び(c)は、それぞれ、x軸、y軸及びz軸に沿って見られる、この線源軌道を示している。これらの図に示される小さい球は、25cmの直径を有し、VOP内の最大のアイソセントリックな球を示す。ねじれた円は、一定スピードのプロペラ運動を正弦波のロール運動と組み合わせることによって、図1に表される例示的なコーンビームCTスキャナによって実現されることができる。(プロペラジョイント部は、制限されない角度レンジを支援しなければならない。)他の線源軌道は、α=0にセットすることによってatcから得られる通常の円ac(λ)=rsrc(sinλ,0,cosλ)である。
アルゴリズムCBFFBPは、パーソナルコンピュータ上で「C」プログラミング言語で実現され、シミュレートされたデータを用いてテストされた。FDKアルゴリズムの匹敵する実現との比較も行われた。2つの線源軌道が試みられた。これらのうちの1つは、
によって規定される「ねじれた円」atc:[0,2π]→R3であった。ここで、rsrc=810mmであり、α=15π/180の場合θ(λ)=αcos2λである。図7(a)、(b)及び(c)は、それぞれ、x軸、y軸及びz軸に沿って見られる、この線源軌道を示している。これらの図に示される小さい球は、25cmの直径を有し、VOP内の最大のアイソセントリックな球を示す。ねじれた円は、一定スピードのプロペラ運動を正弦波のロール運動と組み合わせることによって、図1に表される例示的なコーンビームCTスキャナによって実現されることができる。(プロペラジョイント部は、制限されない角度レンジを支援しなければならない。)他の線源軌道は、α=0にセットすることによってatcから得られる通常の円ac(λ)=rsrc(sinλ,0,cosλ)である。
シミュレートされた検出器は、幅及び高さW=H=380mmの正方形の検知領域を有し、大きさ(Δu,Δv)=(W/I,H/J)のIxJ=512x512個の検出器素子に細分された。検出器座標系は、再構成問題サブセクションに記述されたように選ばれた。ベクトル、
は、全てのλ∈[0,2π]について面y=0と平行だった。c(λ)のd(λ)を中心とする検出器座標は、uc(λ)=−Δu/2及びvc=−Δv/2であった。線源と検出器との間の距離は、w(λ)=1195mmであった。
は、全てのλ∈[0,2π]について面y=0と平行だった。c(λ)のd(λ)を中心とする検出器座標は、uc(λ)=−Δu/2及びvc=−Δv/2であった。線源と検出器との間の距離は、w(λ)=1195mmであった。
シミュレートされる被検体は、面y=0がファントムの横断面になり、面x=0が矢状面になるように、固定座標系の起点に置かれる、変更されたFORBILD頭部ファントムであった。FORBILD頭部ファントムの詳述は、インターネットサイトhttp://www.imp.uni-erlangen.de/forbild/deutsch/results/head/head.htmlにおいて入手可能である。ファントムは、図7に示す小さい球にちょうど収まり、そのコーンビーム投影は、切頭されない。シミュレーションプログラムが、λに関して等距離に間隔を置かれたファントムのK=801コーンビーム投影を計算するために使用された。
画像は、各々の大きさが(0.5mm)3である5123個の3次元ボクセルの中心を持つ3次元ボリュームにおいて再構成された。VOPの外側の再構成された画像の部分は、ゼロにセットされる。アルゴリズムCBFFBPは、両方の線源軌道に関して使用され、FDKアルゴリズムは、円のみに関して使用された。
図8(a)−(f)は、変更されたFORBILD頭部ファントムの再構成されたスライス画像を示している。図8(a)、(c)及び(e)は、横断スライスを示しており、図8(b)、(d)及び(f)は、矢状スライスを示している。図8(a)及び(b)は、円の場合にFDKによって得られるスライスを示している。図8(c)及び(d)は、CBFFBPを用いて得られる対応するスライスを示している。図8(e)及び(f)は、ねじれた円の場合にCBFFBPを用いて得られる再構成されたスライスを示している。全てのケースにおいて、グレースケールウィンドウは、レンジ[0HU,100HU]に対応する。
円の場合、CBFFBP及びFDKによって再構成される画像は、ほぼ同等である。2つのアルゴリズムの連続するバージョンは等価であるので、これは、驚くべきことではない。線源軌道を含む面から外れた矢状スライスの領域には、わずかな強度低下及びいくらかな幾何学的ひずみも見られる。これらのアーチファクトは、良く知られており、環状の線源軌道の完全性が欠如する結果である。ねじれた円の場合にCBFFBPによって生成される横断スライスは、有利に、円の場合にCBFFBP又はFDKによって生成されるものと同等である。しかしながら、ねじれた円の完全性のため、CBFFBPによって生成される矢状スライスは、円の場合にCBFFBP又はFDKによって生成される矢状スライスによって示されるアーチファクトはない。上述の例において、CBFFBPのフィルタリングステップは、FDKアルゴリズムの同様の実現より6.9倍遅く、CBFBP2の同様の実現より約70倍速かった。CBFFBP及びFDKの両方の場合、逆投影ステップに必要な計算時間は、フィルタリングステップに必要な計算時間をはるかに上回った。
本発明の先に述べた例示的な実施例において、離散フーリエ変換(三角和)が、連続するフーリエ変換を近似するために使用され、一様又は非一様なFFTが、これらの離散フーリエ変換を計算するために使用された。原則として、連続するフーリエ変換を近似し、計算する他の方法も使用されることができる。本発明は、2Dフーリエ空間における同心円グリッドの使用によって特徴付けられる。この同心円グリッドは、角度及び半径において等距離であるように選択されることが勧められるが、これは必ずしも必要でない。
先に述べた例示的な実施例において、コーンビーム投影は、等距離のデカルトグリッド上でサンプリングされ、順方向U−NU FFTは、2Dフーリエ空間における同心円グリッド上で、
を計算するために使用された。投影が、等距離のデカルトグリッド上でサンプリングされない場合、このようなグリッド上にそれらをリサンプリングすることが可能であり、又は順方向U−NU FFTの代わりに順方向NU−NU FFTが使用されることができる。
を計算するために使用された。投影が、等距離のデカルトグリッド上でサンプリングされない場合、このようなグリッド上にそれらをリサンプリングすることが可能であり、又は順方向U−NU FFTの代わりに順方向NU−NU FFTが使用されることができる。
図9は、本発明による方法の例示的な実施例を実行するための本発明による画像処理装置の例示的な実施例を示している。図9に示される画像処理装置400は、データを記憶するためのメモリ402に接続される中央処理ユニット(CPU)又はデータプロセッサ401を有する。データプロセッサ401は、複数の入出力ネットワーク又はCT装置のような診断装置に接続されることができる。データプロセッサ401は、データプロセッサ401にて計算され又は処理された情報又は画像を表示するための、例えばコンピュータモニタのような表示装置403に接続される。オペレータ又はユーザは、図9に示されていないキーボード404及び/又は他の出力装置を介して、データプロセッサ401と対話することができる。更に、バスシステム405を介して、データプロセッサ401を、例えば関心対象の動きを監視する動きモニタに接続することが可能である。例えば、患者の肺がイメージングされるようとする場合、動きセンサは、呼気センサでありうる。代替として、心臓がイメージングされようとする場合、動きセンサは、ECGマシンでありうる。
「有する、含む」なる語は、他の構成要素又はステップを除外せず、単数形の構成要素は、その複数の存在を除外せず、単一のプロセッサ又はシステムが、請求項に詳述されるいくつかの手段又はユニットの機能を果たすことができることが留意されるべきである。
請求項におけるいかなる参照符号も、請求項の範囲を制限するものとして解釈されるべきではないことが留意されるべきである。
Claims (8)
- 対象周囲の線源軌道に沿って複数の線源位置において取得される該対象のコーンビーム投影の組から、該対象の関心ボリュームの3D画像を再構成する方法であって、
フーリエ空間における同心円グリッドに基づいて、前記投影の組をフィルタリングして、フィルタリングされた投影データセットを生成するステップと、
前記フィルタリングされた投影データセットから、前記関心ボリュームを再構成して、前記関心ボリュームの再構成された画像を生成するステップと、
を含む方法。 - 前記関心ボリュームの再構成が、式
の離散的な近似による、前記フィルタリングされた投影データセットの逆投影を含み、上式にて、
rは、デカルト固定座標系を基準にした位置ベクトルであり、
frec(r)は、位置rにおける再構成された画像であり、
Vは、関心ボリュームであり、
Λは、境界付けされた間隔であり、
a:Λ→R3は、線源軌道であり、
w(λ)は、線源ポイントa(λ)と該線源ポイントに関連付けられる検出器面Δ(λ)との間の距離であり、
(u(r,λ),v(r,λ))は、前記検出器面Δ(λ)のデカルト検出器座標系を基準とした、a(λ)から前記検出器面Δ(λ)までのrの透視投影の座標であり、
(uc(λ),vc(λ))は、a(λ)の前記検出器面Δ(λ)までの正射影の検出器座標であり、
g(u,v,λ)が、前記線源が位置a(λ)にあるときの検出器検知領域の検出器座標の組である、取得されたコーンビーム投影gm(u,v,λ)に等しい関数であって、(u,v)∈D0(λ)のとき、g6(u,v,λ)は、
(a)u−v面のグリッドにおいて、
を計算するステップと、
(b)ξ−η面の極グリッドにおいて、
を計算するステップであって、前記極グリッドのグリッド点が、半径方向の線及び同心円の交点に位置する、ステップと、
(c)
によって、σ−μ面において、
を規定するステップと、
(d)σ−μ面のグリッドにおいて、
を計算するステップと、
(e)s−μ面のグリッドにおいて、
を計算するステップと、
(f)
として、s−μ面のグリッドにおいて、
を計算するステップと、
(g)σ−μ面のグリッドにおいて、
を計算するステップと、
(h)σ−μ面のグリッドにおいて、
を計算するステップと、
(i)
によって、ξ−η面において、
を規定するステップと、
(j)u−v面のグリッドにおいて、
を計算するステップと、
の離散的なバージョンを実行することによって得られる、フィルタリングされたコーンビーム投影である、請求項1に記載の方法。 - 非一様な高速フーリエ変換が、前記ステップ(b)及び(j)について使用される、請求項2に記載の方法。
- 対象周囲の線源軌道に沿って複数の線源位置に沿って取得される、該対象のコーンビーム投影の組から、該対象の関心ボリュームの3D画像を再構成する画像処理装置であって、
多次元データセットを記憶するメモリと、
フーリエ空間における同心円グリッドに基づいて、前記投影の組をフィルタリングして、フィルタリングされた投影データセットを生成し、前記フィルタリングされた投影データセットから、前記関心ボリュームを再構成して、前記関心ボリュームの再構成された画像を生成するように適応される計算ユニットと、
を有する画像処理装置。 - 対象のコーンビーム投影の組を取得するための取得ユニットと、
前記対象の関心ボリュームの画像を再構成するための画像処理装置と、
を有する検査装置であって、
前記画像処理装置が、フーリエ空間における同心円グリッドに基づいて、前記投影の組をフィルタリングして、フィルタリングされた投影データセットを生成し、前記フィルタリングされた投影データセットから、前記関心ボリュームを再構成して、前記関心ボリュームの再構成された画像を生成するように適応される、検査装置。 - 手荷物検査装置、医用診断装置、材料試験装置及び材料科学解析装置を含むグループのうちの1つとして構成される、請求項5に記載の検査装置。
- 検査装置によって、対象のコーンビーム投影の組から、該対象の関心ボリュームの3D画像を再構成するためのコンピュータプログラムが記憶されたコンピュータ可読媒体であって、前記プログラムは、プロセッサによって実行されるとき、
フーリエ空間における同心円グリッドに基づいて、前記投影の組をフィルタリングして、フィルタリングされた投影データセットを生成するステップと、
前記フィルタリングされた投影データセットから、前記関心ボリュームを再構成して、前記関心ボリュームの再構成された画像を生成するステップと、
を実行するように適応される、コンピュータ可読媒体。 - 対象のコーンビーム投影の組から、該対象の関心ボリュームの3D画像を再構成するためのプログラム要素であって、プロセッサによって実行されるとき、
フーリエ空間における同心円グリッドに基づいて、前記投影の組をフィルタリングして、フィルタリングされた投影データセットを生成するステップと、
前記フィルタリングされた投影データセットから、前記関心ボリュームを再構成して、前記関心ボリュームの再構成された画像を生成するステップと、
を実行するように適応されるプログラム要素。
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
EP05104916 | 2005-06-07 | ||
PCT/IB2006/051780 WO2006131872A2 (en) | 2005-06-07 | 2006-06-02 | Fast reconstruction algorithm for cone-beam ct |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2008541982A true JP2008541982A (ja) | 2008-11-27 |
Family
ID=37027052
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2008515352A Withdrawn JP2008541982A (ja) | 2005-06-07 | 2006-06-02 | コーンビームct用高速再構成アルゴリズム |
Country Status (4)
Country | Link |
---|---|
US (1) | US20080212860A1 (ja) |
EP (1) | EP1891603A2 (ja) |
JP (1) | JP2008541982A (ja) |
WO (1) | WO2006131872A2 (ja) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7792238B2 (en) | 2008-02-18 | 2010-09-07 | General Electric Company | Method and system for reconstructing cone-beam projection data with reduced artifacts |
US8708561B2 (en) | 2009-03-20 | 2014-04-29 | Orthoscan, Inc. | Mobile imaging apparatus |
CN101520899B (zh) * | 2009-04-08 | 2011-11-16 | 西北工业大学 | 一种锥束ct三维图像的并行重建方法 |
DE102010031740B4 (de) * | 2010-07-21 | 2019-05-02 | Siemens Healthcare Gmbh | Vorrichtung und Verfahren für ein Mammographiegerät |
WO2012082799A1 (en) | 2010-12-13 | 2012-06-21 | Orthoscan, Inc. | Mobile fluoroscopic imaging system |
US9600910B2 (en) * | 2014-01-08 | 2017-03-21 | Rensselaer Polytechnic Institute | Attenuation map reconstruction from TOF PET data |
CN105844690A (zh) * | 2016-03-29 | 2016-08-10 | 中北大学 | 基于gpu的快速三维ct迭代重建系统 |
CN105931280A (zh) * | 2016-03-29 | 2016-09-07 | 中北大学 | 基于gpu的快速三维ct迭代重建方法 |
US11436765B2 (en) * | 2018-11-15 | 2022-09-06 | InstaRecon | Method and system for fast reprojection |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
FR2615619B1 (fr) * | 1987-05-21 | 1991-07-19 | Commissariat Energie Atomique | Procede et dispositif d'imagerie tridimentionnelle a partir de mesures bidimensionnelles de l'attenuation d'un rayonnement |
DE10063442A1 (de) * | 2000-12-20 | 2002-07-04 | Philips Corp Intellectual Pty | Verfahren und Röntgeneinrichtung zur Ermittlung eines Satzes von Projektionsabbildungen eines Untersuchungsobjektes |
DE10116682A1 (de) * | 2001-04-04 | 2002-10-10 | Philips Corp Intellectual Pty | Verfahren und Vorrichtung zur Rekonstruktion dreidimensionaler Bilder aus Kegelstrahl-Projektionsdaten |
US7444011B2 (en) * | 2004-02-10 | 2008-10-28 | University Of Chicago | Imaging system performing substantially exact reconstruction and using non-traditional trajectories |
-
2006
- 2006-06-02 EP EP06756055A patent/EP1891603A2/en not_active Withdrawn
- 2006-06-02 WO PCT/IB2006/051780 patent/WO2006131872A2/en active Application Filing
- 2006-06-02 JP JP2008515352A patent/JP2008541982A/ja not_active Withdrawn
- 2006-06-02 US US11/916,540 patent/US20080212860A1/en not_active Abandoned
Also Published As
Publication number | Publication date |
---|---|
EP1891603A2 (en) | 2008-02-27 |
US20080212860A1 (en) | 2008-09-04 |
WO2006131872A3 (en) | 2008-03-13 |
WO2006131872A2 (en) | 2006-12-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Zou et al. | Exact image reconstruction on PI-lines from minimum data in helical cone-beam CT | |
Noo et al. | Exact helical reconstruction using native cone-beam geometries | |
JP2008541982A (ja) | コーンビームct用高速再構成アルゴリズム | |
CN100528088C (zh) | 用圆-线快速扫描算法重构计算机断层分析图像的方法 | |
US6990167B2 (en) | Image reconstruction method for divergent beam scanner | |
JP4714316B2 (ja) | イメージ再構成 | |
US6961404B2 (en) | Method and system for reconstructing an image from projection data acquired by a cone beam computed tomography system | |
JP4611168B2 (ja) | 画像再構成方法、およびx線ct装置 | |
JP5161888B2 (ja) | 円状及び直線状軌道を用いた医用撮像における断層撮影法の再構成のための方法及びシステム | |
EP2315178A1 (en) | Reconstruction of 3D image datasets from X-ray cone-beam data | |
CN102973291B (zh) | C型臂半精确滤波反投影断层成像方法 | |
JPH09149902A (ja) | X線断層撮影方法および装置 | |
EP1663004A2 (en) | Computer tomography method using a cone-shaped bundle of rays | |
JP4342164B2 (ja) | コンピュータ断層撮影装置 | |
WO2005010824A1 (en) | A fourier tomographic image reconstruction method for fan-beam data | |
WO2018179905A1 (ja) | インテリアct画像生成方法 | |
US7272205B2 (en) | Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix | |
US7187747B2 (en) | Computerized tomography method with helical relative movement and conical beam | |
JP2007198866A (ja) | 広義サドルコーンビームct装置および3次元再構成法 | |
JP2006505337A (ja) | 正確な円錐ビームコンピュータ断層撮影のための方法及び装置 | |
Dennerlein et al. | Exact and efficient cone-beam reconstruction algorithm for a short-scan circle combined with various lines | |
CN110974278B (zh) | 一种dsa锥束精确滤波反投影断层成像系统及成像方法 | |
Grangeat et al. | Indirect cone-beam three-dimensional image reconstruction | |
CN117425918A (zh) | 用于ct重建的方法 | |
JP2005522249A (ja) | コンピュータ断層撮影方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A300 | Application deemed to be withdrawn because no request for examination was validly filed |
Free format text: JAPANESE INTERMEDIATE CODE: A300 Effective date: 20090804 |