JP4550426B2 - 高速発散ビーム断層撮影のための方法及び装置 - Google Patents

高速発散ビーム断層撮影のための方法及び装置 Download PDF

Info

Publication number
JP4550426B2
JP4550426B2 JP2003573443A JP2003573443A JP4550426B2 JP 4550426 B2 JP4550426 B2 JP 4550426B2 JP 2003573443 A JP2003573443 A JP 2003573443A JP 2003573443 A JP2003573443 A JP 2003573443A JP 4550426 B2 JP4550426 B2 JP 4550426B2
Authority
JP
Japan
Prior art keywords
image
partial
sinogram
projection
desired number
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.)
Expired - Lifetime
Application number
JP2003573443A
Other languages
English (en)
Other versions
JP2005518892A (ja
Inventor
シャオ,シュー
ブレスラー,ヨラム
デイビッド,シー.,ジュニアー. マンソン
Original Assignee
ザ、ボード、オブ、トラスティーズ、オブ、ザ、ユニバシティー、オブ、イリノイ
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by ザ、ボード、オブ、トラスティーズ、オブ、ザ、ユニバシティー、オブ、イリノイ filed Critical ザ、ボード、オブ、トラスティーズ、オブ、ザ、ユニバシティー、オブ、イリノイ
Publication of JP2005518892A publication Critical patent/JP2005518892A/ja
Application granted granted Critical
Publication of JP4550426B2 publication Critical patent/JP4550426B2/ja
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/027Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis characterised by the use of a particular data acquisition trajectory, e.g. helical or spiral
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10TECHNICAL SUBJECTS COVERED BY FORMER USPC
    • Y10STECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10S378/00X-ray or gamma ray systems or devices
    • Y10S378/901Computer tomography program or processor

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Description

本発明は、発散ビーム断層撮影のための方法及び装置に係り、特には、視認可能なほどの劣化や重大な数値上の誤りなしに計算上の複雑さを低減する、発散ビーム断層撮影のための方法及び装置に関する。
断層撮影の再構成(tomographic reconstruction)は、X線コンピュータ断層撮影法(CT)、ポジトロン放射断層撮影法(PET)、単一フォトン放射断層撮影法(SPECT)、及び磁気共鳴イメージング法(MRI)用の或る捕捉(acquisition)法を含む診断用イメージング法(diagnostic imaging modalities)のほとんど全ての基礎となる周知の技術である。それはまた、非破壊評価(NDE)や、より最近では空港での荷物検査にも広く利用されている。現在のCTスキャナにおける主データ捕捉モードは、ファンビームジオメトリであり、しばしばヘリカルスキャン又はスパイラルスキャンと組み合わされる。
より高速の画像取得のために、CTスキャナ技術は、面検出器(area detectors)を用いたコーンビーム取得へと移行しつつある。生産されたCTスキャナの多くは、すでに、多列(multirow)検出器(典型的には16列)を有することで、(狭角の)コーンビーム捕捉が得られている。広角のコーンビーム捕捉用の面検出器は、今後数年に渡って医用CTスキャナに導入されるであろう。現在、隠された武器や爆破装置を検出するための高速の3D断層撮影画像法への必要性が存在する。PETやSPECTのようなその他の方法においては、コーンビーム捕捉がすでに主フォーマットとなっている。
ファンビームCT及びコーンビームCTにおける再構成上の問題は、特定の捕捉軌跡又は捕捉軌道上を進むソース(source)からそれぞれ扇形又は円錐として発散する線に沿った1組の線積分投影から、オブジェクトの2D又は3D画像を復元することである。従って、ファンビーム及びコーンビームのジオメトリは、発散ビームジオメトリの一例である。
断層撮影の再構成のための選択方法は、フィルタされた(filtered)逆投影(FBP)又は回旋(convolution)逆投影(CBP)であり、このどちらも重み付け(weighted)逆投影ステップを使用する。このステップは、2DにおけるN×Nピクセル画像にとってのNのような尺度化(scaling)や、3DにおけるN×N×Nボクセル画像にとっての少なくともNのような尺度化を行う計算上の必要性を伴う技術において、計算上のネックとなっている。従って、画像分解能を2倍にしてNから2Nにしたとすると、計算がおおよそ8倍(3Dでは16倍)増加する。これまで以上の大量データをリアルタイムで収集可能な新技術(例えば、多列検出器を備えた心臓画像法、干渉画像法)の出現や、3Dコーンビームジオメトリへの移行により、コンピュータがより一層高速になる一方で、高速再構成技術の必要性がまだ高まりつつある。高速再構成は、画像形成プロセスをスピードアップし、特殊目的の画像再構成コンピュータ(以下を参照)のコストを低減し、或いはそれら両方を可能にする。
逆投影の双対(dual)演算は、再投影(reprojection)であり、これは電子的に格納された画像の投影を計算するプロセスである。このプロセスも、断層撮影の再構成において基本的な役割を演じる。逆投影と再投影の組み合わせは、被験体の実際の3D画像法へのキーとなる、ヘリカルコーンビームジオメトリにおける長いオブジェクトの問題のための高速再構成アルゴリズムを構築するのにも利用可能である。更に、各種応用において、逆投影ステップと再投影ステップの両方が1つの画像の再構成のために数回行われる、反復の再構成アルゴリズムを使用することは、有利であり、或いは必要でさえもある。逆投影ステップと再投影ステップをスピードアップすることは、そのような反復方法の経済的な実行可能性を決定することになろう。
高速方法について議論する場合、断層撮影データ用の2つの主フォーマット、すなわち(i)平行ビームと(ii)発散ビーム及びその特別な場合、を互いに区別することが重要である。平行ビームフォーマットは理論上の操作と計算上の操作の両方にとって最も便利であるのに対し、発散ビームは市販の医用CTスキャナや産業用CTスキャナにおいて最も一般的に見出されるフォーマットである。元々、2D法のファンビーム再構成は、ヘリカル及びマルチスライスヘリカル3D(体積測定の)再構成のための技術状態法(stage-of-the-art)におけるキー要素であるが、新たな設計及び将来の設計は、コーンビームジオメトリへと移行することが予測される。従って、2Dにせよ3Dにせよ、発散ビームジオメトリは主たる画像フォーマットであり、おそらく近い将来においてもそのままであろう。
オブジェクトからの大きなソース距離の限界では、発散ビームジオメトリは平行ビームジオメトリにまで減少する。従って、(本発明の方法を含む)発散ビームジオメトリのための全てのデータ処理方法は、平行ビームジオメトリにも応用可能である。しかし、その逆は真ではない。すなわち、実際に見出された小さな又は適度のソース距離では、2つのジオメトリは十分に相違しており、平行ビーム処理方法は、これを発散ビームデータに直接応用した場合は劣った結果又は不満足な結果を生じる。従って、発散ビームジオメトリは、平行ビームジオメトリとは異なる処理を必要とする。
本発明は、発散ビーム断層撮影法を扱う。再構成方法の幾つかは、再ビニング(rebinning)と呼ばれるプロセスを利用することで、発散ビームデータを平行ビームデータに再配置(又は補間)し、これは次に平行ビームアルゴリズムによって処理される。3D発散ビームジオメトリのための他の方法は、コーンビームデータを3Dランダム変換データ(又はそれの派生データ)に変換することを利用する。その代わりに、限定による元々の発散ビーム法(本発明はその高速バージョンである)が、先立って平行ビームデータを再ビニングすることや、そのデータを3Dランダム変換領域に変換することなしに、発散ビームデータを直接処理する。元々の発散ビーム法における再構成プロセスは、投影の前処理(例えば重み付けやフィルタリング)、それに続く重み付けされた発散ビームの逆投影演算、及びおそらくは発散ビームの再投影演算、又は一連のそのような演算によって構成されている。従って、そのような方法は、発散ビームのフィルタされた逆投影(DB−FBP)アルゴリズムと呼ばれる。それらの計算は、平行ビームの場合には、通常、逆投影演算及び再投影演算によって支配される。
発散ビームの再投影は、まず平行ビームの再投影を行い、次にその結果を再ビニングすることによってなし得るが、ここでも、再ビニングを使用せずに元々の発散ビーム再投影アルゴリズムを使用するのが有利である。
再ビニングの欠点は、スピードアップを制限する、アーチファクト(artifacts)及び追加の計算を含むことである。再ビニングはまた、処理が開始可能になる前にデータの多くの部分の捕捉及び操作を必要とし、これはスピードを再び制限することになる。3Dランダム領域への変換を利用する方法は、同様な欠点を有している。従って、発散ビーム断層撮影のための方法及び装置にとって、これらの欠点を解消し、高度に柔軟性があり、従来の発散ビーム逆投影を使用する方法と比べて大きな性能向上を提供する必要がある。
特殊目的のハードウェアが、逆投影プロセスをスピードアップするための伝統的な手段である。カスタムチップ、マルチプルプロセッサ、又はそれらの組み合わせを使用する特別なタイプのコンピュータが、必要な計算を行うために設計されている。
このアプローチの欠点は、マルチプルプロセッサやカスタムハードウェアのコストを含み、また、一般目的のコンピュータの性能が急速に向上しているので、要求される特殊目的の構造がすぐに時代後れになってしまうという事実を含んでいることである。従って、発散ビームジオメトリにおける高速処理は、特殊目的のハードウェアを必要とせず、かつ、標準のシリアル又はパラレル構造において容易に実行されることで、より高いコスト効率にする必要がある。
従って、本発明の目的は、発散ビーム断層撮影のための新規で改良された方法及び装置を提供することにある。
他の目的は、視認可能なほどの劣化や重大な数値上の誤りなしに計算上の複雑さを低減する、発散ビーム断層撮影のための新規で改良された方法及び装置を提供することにある。
本発明の1つのアスペクトによれば、逆投影を施せる、前処理された発散ビームシノグラムから電子画像を生成する方法が提供され、上記シノグラムは発散ビーム投影の集合である。この方法は、上記シノグラムを複数の部分シノグラムに細分するステップと、撮像されるオブジェクト上に中心を置く全体座標系内で上記部分シノグラムの重み付けされた逆投影を行って、全体座標系内の正しい位置に中心を置く複数の対応する部分画像を生成するステップと、上記部分画像を集合させて上記電子画像を生成するステップと、を含んでいる。上記シノグラムの部分シノグラムへの細分は、所望回数の正確な細分と、所望回数の近似の(概算的な)細分とを行うことにより、なされる
ここで、上記所望回数の正確な細分は、シノグラムにおける投影のトランケーションを行うこと、によってなされる。また、上記所望回数の近似の(概算的な)細分は、シノグラムにおける投影のトランケーションを行い、このトランケーションの行われた投影を全体座標系の中心に対して、選択された距離だけシフトさせ、このシフトされた投影のデシメーションを、選択された係数によって行い、このデシメーションの行われた投影をその元の位置へ、上記選択された距離を用いてシフトさせること、によってなされる。
本発明の他のアスペクトによれば、電子画像を再投影する方法、すなわち該画像から発散ビームシノグラムを生成する方法が提供される。この方法は、上記電子画像のトランケーションを行うことにより、上記電子画像を、上記電子画像上に中心を置く全体座標系内の正しい位置に中心を置く複数の部分画像に分割するステップと、全体座標系内の各部分画像の部分シノグラムを計算するステップと、上記部分シノグラムを集合させて上記シノグラムを生成するステップと、を含んでいる。上記画像の細分は、所望回数の正確な再投影と、所望回数の近似の(概算的な)再投影とを行うことにより、なされる。
ここで、上記所望回数の正確な再投影は、全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成すること、によって行われる。また、上記所望回数の近似の(概算的な)再投影は、全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成し、この部分シノグラムにおける投影を、全体座標系の中心に対して、選択された距離だけシフトさせ、上記部分シノグラムにおける追加の投影を、補間を行うことにより生成し、上記部分シノグラムをその元の位置へ、上記選択された距離を用いてシフトさせること、によって行われる。
これらの方法は、従来の逆投影アルゴリズムや再投影アルゴリズムと比べて、FFTと同じ程のスピードアップを提供する。認識可能なほどの劣化や重大な数値上の誤りなしに、実質的な計算上の節約が得られる。
添付の図面を参照しつつ、本発明の実施形態についての以下の記載を参照することにより、本発明における上述した構成及びその他の構成と、それらを得るための方法が、一層明らかになり、そして、本発明が最も良く理解されることになろう。
本発明は、CTスキャナを含む様々なイメージング装置への応用を有している。典型的なイメージング装置10(図1)は、頭のようなオブジェクトからデータを取得して、発散ビーム投影(DB投影)14に相当する生(raw)データを投影プリプロセッサ(前処理装置)16へ送るスキャナ12を含んでいる。投影プリプロセッサ16は、重み付けや、シフト変更であるフィルタリングだけでなく、データに対して様々な変換、正規化、及び訂正を加える。投影プリプロセッサ16の出力は、発散ビームシノグラム(DBシノグラム1)18であり、これはシノグラム更新プロセッサ20に与えられる。このシノグラム更新プロセッサ20は、発散ビームシノグラム(DBシノグラム2)34からの情報を用いて、入力されたシノグラム18を変更可能であり、これは例えば、ビーム硬化(beam-hardening)を含む様々なアーチファクトのための訂正を行いつつ、又はマルチステップの中の一部として、又は反復する再構成手順として行われる。
シノグラム更新プロセッサ20の出力は、発散ビームシノグラム(DBシノグラム3)22であり、これは高速逆投影プロセッサ24への入力となる。この高速逆投影プロセッサ24は、一般に、コンピュータ又は特殊目的のハードウェア、又はそれらの何らかの組み合わせであり、この中で記述されている正確な又は近似の逆投影アルゴリズムを行うようプログラムされた何らかの適切なタイプのものである。
高速逆投影プロセッサ24の出力は、電子画像(電子画像1)26であり、これは画像コンディショニングプロセッサ28に入力される。この画像コンディショニングプロセッサ28は、電子画像の必要な後処理を行うものであり、アーチファクト画像や、マルチステップ又は反復再構成プロセスにおける更なる処理のための画像を確認及び抽出することを含んでいてもよい。
必要ならば、画像コンディショニングプロセッサ28は、高速再投影プロセッサ32に与えられる電子画像(電子画像2)30を生成可能である。高速再投影プロセッサ32は、一般に、コンピュータ又は特殊目的のハードウェア、又はそれらの何らかの組み合わせであり、この中で記述されている正確な又は近似の再投影アルゴリズムを行うようプログラムされた何らかの適切なタイプのものである。必要ならば、このプロセッサは、逆投影プロセッサ24に採用されたのと同じコンピュータ及びハードウェアの使用を共有することも可能である。
高速再投影プロセッサ32の出力は、発散ビームシノグラム(DBシノグラム2)34であり、これはシノグラム更新プロセッサ20にフィードバックされる。逆投影/再投影のプロセスは、適切な結果が得られるまで継続可能である。再投影は必ずしも必要とされないが、多くの状況において役立つものである。
電子画像26が適切である場合、画像コンディショニングプロセッサ28が電子画像(電子画像3)36を生成し、これは格納/解析/表示装置38に与えられる。これにより、電子画像36は、コンピュータメモリに格納可能であり、及び/又は、例えば異常や危険物質のために電子的に解析可能であり、及び/又は、表示可能であり、及び/又は、幾つかの視覚可能な形態で印刷可能である。
本発明の高速発散ビーム逆投影アルゴリズムは、米国特許第6,287,257号に記載の平行ビーム断層撮影用の高速階層逆投影アルゴリズム(fast hierarchical backprojection algorithm;FHBP)に基づいている。同様に、本発明の高速発散ビーム再投影アルゴリズムは、米国特許第6,263,096号及び同第6,351,548号に記載の平行ビーム断層撮影用の高速階層再投影アルゴリズム(fast hierarchical reprojection algorithm;FHRP)に基づいている。逆投影と再投影は、密接に関連しており、それらの高速アルゴリズムは同一原理から得られる。より一般的に使用される高速逆投影は、より詳細に記述することとする。
平行ビームFHBPを理解するために、画像の中心に原点を有する全体座標系の付随したN×N画像について考える。この大きなN×N画像を、部分画像と呼ばれるN/2×N/2のサイズの小さな4つの画像に分解することを考える。部分画像の各々は、全体座標系の別々の象限に配置される。ラドン変換(Radon transform)の2つの特性が、FHBPアルゴリズムの派生に使用される。(i)ボウ−タイ特性(bow-tie property)については、全体座標系の原点に中心を置いた半分サイズの画像にとって、角度変数に関する投影のスペクトル支持体(spectral support)も半分だけ減少される、と言われている。これは、中心に置かれた半分サイズの画像は投影の数の半分を備えるシノグラムから再構成されることが可能であることを意味している。(ii)シフト特性(shift property)については、シフトされた画像の平行投影が原画像のシフトされた平行投影に相当する、と言われている。
今、P個のフィルタされた投影を備えるシノグラムが、全体画像fを再構成するのに利用可能であるとする。部分画像の1つfを再構成するために、P個の投影がまず部分画像の投影の支持体にまでトランケートされ(以下、「トランケートする」とは、「トランケーションを行う」の意味である)、次に部分画像の中心に置かれたバージョンに相当するようシフトされ(shifted)、最後に角度方向に関してP/2個の投影にまでデシメートされる(以下、「デシメートする」とは、「デシメーションを行う」の意味である)。部分画像の中心に置かれたバージョンfが、続いてこれら新たなP/2個の投影から再構成され、この部分画像が全体座標系中の正しい位置に戻るようシフトされる。これは、4つの部分画像のそれぞれについて行われ、この4つの部分画像が一緒に集合されると全体画像fの再構成が提供される。合計で4c(P/2)(N/2)=cPN/2の演算が、4つの部分画像を再構成するのに必要であり、ここでcは定数である。これは、元々の再構成コストを1/2だけ減少させることになる。反復して分解することで、その各ステップで画像サイズを1/2だけ減少させ、1/2の数の投影を用いて半分のサイズの部分画像の形成を許容するので、計算コストの合計をO(NlogN)にまで低減させることが可能である。
この分解及び集合のアプローチの双対(dual)は、米国特許第6,263,096号及び第6,351,548号に記載の高速階層平行ビーム再投影アルゴリズム(FHRP)で使用されている。
発散ビームの逆投影及び再投影のための本発明は、画像を階層的に分解し、次に、より少ない発散ビーム投影を備える発散ビーム部分シノグラムを用いて部分画像を再構成する(又は再投影する)、という考えに基づいている。しかし、平行ビームFHBP又はFHRPの考えを発散ビームに適用しようとする場合に生じる直接の困難性は、上述のシフト特性(ii)がもはや当てはまらないということである。すなわち、シフトされた部分画像の発散ビーム投影は、元の部分画像のシフトされた発散ビーム投影に相当しないということである。平行ビームFHBP及びFHRPはこの特性を使用するので、平行ビームアプローチは発散ビーム用に変更しなければならない。
本発明のアルゴリズムは、分解用の異なる手順を使用するものであり、この手順は、オブジェクトに固定された、すなわちフルサイズ画像に固定された全体座標系の中心に対する部分画像のシフティングを含まない。平行ビームFHBPと同様、新規な発散ビーム逆投影アルゴリズムは、2つのタイプの分解、すなわち1つは正確だが低速のもの、もう1つは近似だが高速のものを有している。この新規なアルゴリズムの正確な発散ビーム分解は、データのシフティングを全く必要とせず、部分画像に相当する部分シノグラムへのシノグラムのトランケーションを必要とするのみである。これらの部分シノグラムは、次に、全体座標系内に逆投影されることで、対応する部分画像を全体座標系内の正しい位置に生成する。近似の発散ビーム分解は、平行ビームFHBPに使用されるものとは異なる手順を使用することで、より小さなサイズの部分画像を再構成するのに使用される投影の数を減少させる。この新規な手順は、投影のみのトランケーション及びシフティングを含み、部分画像のものは含まない。この新規な手順は、FHBPには使用されていなかった重み付けを含んでいてもよい。同様なことが、発散ビーム再投影のための新規な手順に適用される。
半分サイズの部分画像を再構成するために、適切にトランケート(重み付け及びフィルタ)された発散ビーム投影が、全体座標系における部分画像の中心の投影位置によって決定される距離だけシフトされ、次に、ソース位置パラメータにおいて2だけデシメートされ、そして、同一量だけシフトバックされる。この部分画像は、続いて、この減少された投影群から良好な近似へと再構成されることが可能である。部分画像に対する正確な重み付け逆投影と近似の重み付け逆投影の両方において、重み付けを含む、画像座標を必要とする全ての計算は、元々のフルサイズ画像を意味する絶対画像座標又は全体画像座標の中でなされる。正確な分解と近似の分解とを反復形式で組み合わせることで、発散ビームジオメトリのための高速階層逆投影アルゴリズムが得られる。その平行ビームFHBPにおける場合と同様、これらのアルゴリズムは、近似の分解に対する正確な分解の数を選択し、かつ、各種のアルゴリズムパラメータを調整することにより、計算と正確度(精度)との間のトレードオフを提供する。

〔発散ビームイメージングジオメトリ〕
3D断層撮影のための一般的な発散ビーム捕捉ジオメトリは、画像化されるオブジェクトの周囲の3D空間内の軌跡上を移動する発散放射線の頂点(通常はソース)と、検出器面とを含んでおり、この面上において、画像化されるオブジェクトを通るソース放射線に沿った線積分が計測される。発散ビームジオメトリではなく収束ビームジオメトリを有する、放射断層撮影法(SPECT及びPET)に見られるような他のイメージング法は、このジオメトリに変換可能である。
上記の捕捉ジオメトリは、通常、検出器の形状、それに関連する放射線サンプリングスキーム、及びソースの軌跡又は軌道によって分類される。2つの最も一般的な検出器の種類は、(i)等間隔の検出器を有する平面状検出器面(planar detector surface)と、(ii)方位角が等間隔にサンプルされ、かつ高さが直線的等間隔にサンプルされる放射線を有する円柱状検出器面(cylindrical detector surface)である。面検出器を有する3Dイメージングにおいて、発散ビームジオメトリはコーンビームジオメトリとしても知られている。
発散ビーム2D捕捉の特別な場合において、検出器面は、単一線状又は弧状の検出器に減少され、それらはそれぞれ、いわゆる同一直線等間隔ファンビームジオメトリと等角ファンビームジオメトリ(colinear equispaced, and equiangular fan-beam geometries)となる。これらのファンビームジオメトリは、現在市販されているシステムにおいて広く配備されている。
提案されている各種ソース軌道の中で最も一般的なのは、単一円形軌道と、ヘリカル(又はスパイラル)軌道である。しかし、2つの円や、円と線からなるような他の軌道も、実際的な関心を有するものであり、これら各種の軌道を任意に組み合わせたバージョンも実際には生じ得る。
本発明のプロセス、アルゴリズム、及び装置は、一般的な発散ビームジオメトリに適用可能であり、また、先に述べた特別なジオメトリのうちのいずれかや、発散ビームジオメトリの他の例に専門に使用されることも可能である。同様に、本発明の方法は、収束ビームジオメトリを含む(これには限定されない)何らかの数学的に同等なデータ群にも適用可能である。前記の目的のために、この方法は、以下の場合について少々詳細に説明する。

1.単一円のソース軌道
(a)同一直線等間隔ファンビームジオメトリ
(b)等角ファンビームジオメトリ
(c)平面コーンビームジオメトリ
2.平面検出器を有する一般のコーンビームジオメトリ
3.ヘリカルソース軌道

一般のプロセスがこれらの場合に対してどのように専門化されるのかを、実証する。このことから、或る1つの検出器形状及びソース軌道のためのアルゴリズムはそれとは異なる検出器形状及びソース軌道のために容易に変更されるものであるということが明らかになる。

〔任意の軌道を有する発散ビームイメージングジオメトリ〕
平面等間隔コーンビームジオメトリは、検出器が1つの平表面上に均一に配置されるものであり、これは図2に示されている。この場合、検出器平面は、各ソース位置毎に任意の位置及び方向を有することが可能となっている。
発散放射線のソースSは、3Dオブジェクトの周囲の軌道上を進行し、オブジェクトに固定された全体座標系の中心Oに関するソース位置ベクトル
によって定義される。このベクトルは、パラメータθによって特定される軌道上のソース位置を有しており、θmin <θ<θmaxである。原点Oθが検出器平面上への軌道点
の正投影(orthogonal projection)となるよう、各θ毎に、検出器に付随した3D座標系が定義される。正規直交ベクトル
が検出器平面内で選択され、ここで、
は検出器平面に直交する単位ベクトルであり、Oθからソースを向いている。検出器平面内の点
が、その座標(t,t)を用いて定義される。この系において、ソース軌道点は
と表現され、ここでD(θ)はソースと検出器との間の距離である。(SPECTでは、D(θ)はコーンビームコリメータの焦点距離である。)
は、3D画像(すなわちオブジェクト)内の位置を示すものとし、
は、点
(図2を参照)を通過するソース放射線が検出器平面と交差するt,t位置を示すものとする。ソース軌道位置θ及び検出器位置(t,t)でのオブジェクトfの発散ビーム投影(Pf)(θ,t,t)は、(θ,t,t)によってパラメータ化されたソース放射線に沿った線積分であり、
によって与えられる。ここで、
はディラックのデルタ関数である。ソース放射線がソースから発散するのではなくソースへ向けて収束するような状況を記述するのにも、同じ捕捉ジオメトリを使用可能である。従って、本発明のプロセスは、収束ビームジオメトリにも等しく適用される。
投影は、ソース軌道に沿ったパラメータによって定義されるP個の別々のソース位置θp=pΔθ,p=0,・・・P−1、で得られるが、通常、必ずしも均一な間隔Δθ=(θmax−θmin)/Pである必要はない。(Pf)(θ,・)は、
の全ての値において、ソース位置θでの発散ビーム投影と呼ばれる。検出器平面は、通常、t及びt軸上で異なる間隔T及びTを有する、均一な直角グリッド上でサンプルされる。
異なるソース位置での発散ビーム投影の集合は、発散ビームシノグラムと呼ばれる。これと同一の用語が、一組の前処理された発散ビーム投影や、その他の数学的に同等なデータ群においても使用される。シノグラム中の投影の集合は、不完全な場合もあり、単一の投影を備えていてもよい。
一般の発散ビームジオメトリの幾つかの特別な場合について、これから記述する。

〔平面検出器を有するコーンビームヘリカル断層撮影〕
この場合のソース軌道は、図3に示されているように、軸回りの螺旋又は渦巻きであり、
で定義される。ここで、半径Dは一定であり、hは螺旋のピッチとして知られている。
この場合のパラメータθは、x軸に対する(x,y)平面内のソース角度を示している。検出器平面は、一般性を損なうことなしに、x軸を含み、ソースの動きと共に回転すると考えることで、それは各θ毎に(x,y)平面上へのソースベクトル
の投影に垂直となる。検出器平面内の原点Oθは、この場合、z軸上に落ちる。t軸はz軸と一致し、t軸は(x,y)平面に平行である。(その他の一般的な選択、例えば、
が位置θで軌道の接線に平行であり、
がそれに直角であるといったことは、ここで座標回転によって記述される場合まで減少可能である。)
実際には、検出器平面は異って配置されてもよいが、単純な直線座標変換が、実際の検出器平面上で得られる投影を、(t,t)座標に沿って均一な間隔を有する上述の仮想の検出器上で得られる投影に変換する。

〔単一円軌道及び平面検出器を有するコーンビームイメージング〕
この場合のソース軌道は、図4に示されているように、ヘリカル軌道の特殊な場合であり、ピッチがゼロ h=0、すなわち、軌道がx,y平面内において原点Oでの半径Dの単一円
の上にある。検出器平面は、一般性を損なうことなしに、z軸を含み、ソースから中心への線
に垂直であると考える。
オブジェクト内の点
について考える。検出器角度θでの、検出器上へのそのコーンビーム投影の位置
は、
によって与えられる。
これらの明白な表現は、後に、高速アルゴリズムの説明において有用である。

〔円形軌道及び同一直線等間隔検出器を有するファンビーム断層撮影〕
この場合は、図5に示されるように、オブジェクトの単一断面がイメージされ、2D画像が再構成される。ソース軌道はまた、(x,y)平面内の半径Dの円であり、原点Oに中心が置かれ、ソース角度θで示されるソース位置を有している。しかし、検出器は(x,y)平面内の単一直線(t軸)に限定され、これは一般性を損なうことなしに、ソースの回転中心Oを通過し、ソースから中心への線
に垂直である。検出器要素はt軸上に均一に配置されている。
は2D画像中の位置を示すものとし、
はt軸と
(図5参照)を通過するソース放射線との交点のt位置を示すものとする。ソース角度θ及び検出器位置tでのオブジェクトfのファンビーム投影(Pf)(θ,t)は、(θ,t)によってパラメータ化されたソース放射線に沿った線積分であり、
によって与えられる。
投影は、均一な間隔Δθ=θmax/Pを有するP個の別々のソース角度θp=pΔθ,p=0,・・・P−1、で得られる。(Pf)(θ,・)は、tの全ての値において、ソース角度θでの投影であり、異なるPでの投影の集合がシノグラムである。最大のソース角度θmaxは、フルスキャンの場合には2πに等しいが、ショートスキャンやオーバースキャンの場合には他の値をとることも可能である。部分的なデータの集合は、角度の何らかの部分集合での投影からなっていてもよい。
発散ビームジオメトリのそれぞれの場合において、平面ではなく曲面の検出器面が、変更された放射線サンプリングスキームと共に、しばしば使用されてもよい。例えば、ヘリカル又は単一円形軌道のコーンビームイメージングでは、検出器面はしばしば、ソース上に中心の置かれた円柱の表面となるよう選択され、放射線サンプリングは方位が等角で、高さzが等間隔となる。これは、図6に単一円形ソース軌道で記載されており、ここでは、曲面検出器がz軸を含むように選択されている。
円柱検出器形状と平面検出器形状との関係は、一般に、ファンビームジオメトリを研究することによって得られる。図7に示されているような、等角放射線センプリングを有するファンビームジオメトリにおいて、検出器dは、ソースSに中心の置かれた弧上で、角度がΔtの間隔で等間隔に配置されている。従って、このイメージングジオメトリにおいては、同一直線等間隔ファンビームジオメトリに見られる検出器位置tの代わりに、ファン(扇形)角度tが使用される。今、図7に示されるように、
が点
を通過する放射線のファン角度位置を示すものとする。その時、投影は再び式(4)によって与えられる。
同様に、本発明の適用される、発散ビームジオメトリの特別な場合であるファンビームジオメトリの他の例が存在するということを考えてみる。

〔発散ビームのフィルタされた逆投影再構成及び再投影〕
発散ビームデータの本来の(直接の)逆転は、近似であっても正確であっても、重み付けされフィルタされた逆投影(weighted filtered backprojection)として明確化可能である。先ず、発散ビームシノグラムが前処理される。例えば、シノグラムを構成する投影が重み付け及びフィルタされて、ソース位置θに相当する投影
を構成する、変更された発散ビームシノグラムを生成する。このフィルタリングは、ファンビーム又は近似のコーンビーム再構成の場合のように、単純なシフト不変ランプフィルタリング(shift-invariant ramp filtering)であってもよく、又は、正確なコーンビーム再構成の場合のように、もっと複雑なシフト変化フィルタリング(shift varying filtering)であってもよい。重み付け及びフィルタリングは、部分スキャンされた長いオブジェクトを有するヘリカルコーンビームデータからの再構成のための正確又はほぼ正確な方法の場合のように、2以上の投影群を生成してそれらの分担が加算される、もっと一般的なものであってもよい。簡潔性及び一般性のために、逆投影されるデータは、それを元々の発散ビーム投影データから得るのに使用される前処理がどんなものであろうと、「発散ビーム投影データ」又は「発散ビームシノグラム」と呼ぶこととする。固定されたpでは、それは(ソース位置θでの)「発散ビーム投影」と呼ばれ、関数g(p,・)で示される。
その場合、3D画像は、従来の発散ビームの離散θでの逆投影(divergent-beam discrete-θ backprojection)
によって再構成される。ここで、
は適当な重み関数である。この離散(discrete)逆投影式は、全てのθで計測される投影を使用する逆投影のための積分表現を近似する。離散検出器を使用するために、実際には、
はまた
でサンプルされる。これは、
が通常は利用可能なサンプル位置に相当しないために、逆投影を実行するには、変数
でgの補間を必要とする。2Dファンビームの場合における逆投影式は、τがスカラーで
が2次元であることを除き、同一の形式を有する。
N×N×N画像での3Dコーンビーム逆投影の計算コストは、cNPである。ここで、定数cは補間の複雑さのような実行詳細に依存する。これに対して、重み付け及びランプフィルタリングの計算コストは、フィルタリングがシフト不変であって、コンボルーション(convolution)がFFTを用いて行われる場合、ほんのO(PNlogN)にすぎない。フィルタリングのもっと手の込んだ形式も、しばしば同様なコストで行うことが可能である。従って、逆投影のコストは、多くの場合のようにP=O(N)である時においてO(PN)又はO(PN)のコストを有する従来のコーンビーム再構成のコストを卓越している。この状況は、フィルタリング及び逆投影のステップの複雑さがそれぞれO(NlogN)及びO(N)である2Dファンビーム再構成においても、同様である。
電子的に格納された画像fの本来の(直接の)発散ビーム再投影は、式(1)の実行である。実際には、fは離散形式でのみ利用可能であり、また、

でのサンプルのみが計算される。例えば、fは次の級数
によって表すことが可能である。ここで、係数βijkはfの離散表現であり、
は、スプラインのテンソル積のような局所化された基本関数、又は標準ピクセル基本関数、すなわち第ijk番目(ijk-th)のボクセルでの指示関数であり、これらは、2Dについては例えば米国特許第6,263,096号及び第6,287,257号に記載されており、3Dについては米国特許第6,332,035号に記載されている。その場合の投影は、
によって計算される。ここで、
は基本関数
の発散ビーム投影であり、式(1)によって定義される。これら基本関数の投影は、通常、解析表現を用いて計算されるか、又はルックアップテーブルに格納されている。逆投影の場合と同様、本来の再投影のコストは、3DではO(N)であり、2DではO(N)である。
,tが検出器上の適当な座標である曲面検出器の場合、今、
が点
を通過するソース放射線の検出器平面上の(t,t)位置を示すものとする。投影は、再び式(1)によって与えられ、曲面検出器上の基本関数
の発散ビーム投影を示す
と共に、式(7)を使って計算可能である。同様に、再構成は、適当に変更された投影
の重み付け逆投影として再び表すことが可能であり、ここで、t,tは曲面検出器上の適当な座標である。逆投影は前記と同一の式(5)によって与えられるが、異なる重み付け関数
を伴っている。例えば、2D等角ファンビームジオメトリにおいては、等間隔同一直線ジオメトリの線検出器上の変位tの代わりに、ファン角度tが使用される。
の離散性に由来する考察は、平面検出器アルゴリズムにおけるものと同じであり、計算の複雑さについての考察である。
既に論じたように、2Dであろうが3Dであろうが、逆投影ステップ(及び、もし使用されるのであれば再投影)は、通常、計算上は最も高くつくものであり、従って、それは本発明の目標でもある。あらゆる発散ビームジオメトリに応用可能な、高速逆投影の新規な一般プロセスについて記述する。例を単純化するために、まず、2Dファンビームジオメトリのための高速逆投影及び再投影について述べる。

〔同一直線等間隔ファンビームジオメトリのための高速の本来の逆投影及び再投影アルゴリズム〕
図8に示されるような、fの部分画像
のための逆投影演算について考える。
が、
によって定義される画像トランケーション演算子であるとする。ここで、
である。従って、
が、fの部分画像であり、そのサイズはM×Mで、
に中心が置かれている。ソース角度pΔθ,p=0・・・Q−1、でのQ個の投影を用いた、全体座標系内での(全体座標系におけるその正しい位置に配置された)部分画像f′上への逆投影は、
によって与えられる。
は、関係する逆投影演算子を示している。従って、特に、
である。
逆投影の局所性のために、投影g(p,・)の一部分のみがf′上への逆投影に寄与する。この部分が
によって示される。ここで、
は、各θ毎に変数tでのg(p,t)を部分画像
(図8参照)の支持体(support)の投影の支持体
にトランケートする演算子である。
を決定するトランケーション間隔は、全てのソース角度θと、対象の部分画像のサイズM及び位置
とで、前もって計算可能である。この代わりに、より単純で控えめな近似を使用可能であるが、シフティングとデシメーションを必要とする投影の必要な支持体を過剰評価することになるため、階層アルゴリズムにおいて計算コストが幾分上昇する。
逆投影の局所化から、もし
である場合には
となる。すなわち、f′上への逆投影は、近似的に適切にトランケートされた投影のサイズ(M,P)の逆投影によって得ることができる。
今、画像fを、それぞれのサイズがM×MのオーバーラップなしのJ=(N/M)個の部分画像へと区切ることを考える。
ここで、ベクトル
は、図9に示されるように、全体座標系における部分画像の中心である。式(10)を適用すると、以下の、逆投影のための画像上への、逆投影への正確な分解が得られる。
この分解は、N/M=2におけるものが図10に示されている。式(12)中のgによっても示される発散ビームシノグラム22は、今述べた部分画像に対応する複数の部分シノグラムに分割される。この部分シノグラムは、1ピクセル又は1ボクセルと同じ程度に小さな、何らかの所望のサイズであってもよい。トランケーション演算子1010、1012、1014及び1016によって生成される部分シノグラムは、個々に、それぞれ1018、1020、1022及び1024において全体座標系内に逆投影される。逆投影されたシノグラムは部分画像f、f、f及びfを生成し、これらは1026で集合されて、図1中の電子画像26を生成する。このプロセスは、ファンビーム投影の正確な分解を生成する。N/M=2では、図10に示されるように、4つの部分シノグラム及び部分画像へと分解される。
図10を再び参照すると、線22の下のPは、シノグラム又は部分シノグラム中の投影の数を示す。この分解は、それ自身では、単一ステップの逆投影BN,Pと比べてスピードアップを提供するものではない。事実、単一の部分画像逆投影BM,Pの計算コストcMPは、全画像での単一ステップBN,Pのものよりも(N/M)=J倍だけ小さいのに対し、必要とされるJ個の部分画像逆投影における合計コストは、依然として同じである(簿記のための僅かな追加の諸経費を伴う可能性もある)。正確な分解は、J個の並列プロセッサ間で逆投影を区切ることによるスピードアップを提供したり、或いは、J倍小さなメモリを必要とする単一プロセッサの繰り返し使用を許容するのに、まだ使用可能である。
高速逆投影アルゴリズムは、以下の追加の特性を伴う分解の考え方を使用する。固定された画像分解能(帯域幅)では、画像を再構成するのに必要とされる投影の数が、画像のサイズに比例する。すなわち、N×N画像を再構成するのにP個の投影が必要とされる場合、それと同一分解能でM×M画像を再構成するのにはP′=(M/N)P個の投影で十分である。
従って、高速逆投影アルゴリズムは、減少された数の投影からの部分画像
の再構成に基づいている。この減少プロセスは、演算子
によって行われ、この演算子自体は、今から定義する幾つかの演算子から構成されている。
投影シフト演算子
及び
は、それぞれt軸に沿って
だけ投影をシフトさせるものであり、
として定義される。
L倍(L-fold)角度デシメーション演算子
は、角度フィルタリングの後に角度サブサンプリングを行うことにより、1組のP個の投影における視界(view)の数をP/L個の投影に減少させる。このフィルタリングは、隣接する角度で投影を平均化するのと同じくらいに単純であってもよく、或いは、もっと手の込んだものであってもよく、そしてこれは、
と表される。ここで、「*」はコンボルーション(convolution)を示し、hはフィルタ零空間(filter kernel)である。このフィルタリングは、それぞれの固定されたtにおけるp(即ち、離散したθ)にのみ発生してもよく、或いは、θ及びtでの組み合わされたフィルタリングであってもよく、それは分離可能でもそうでなくてもよい。角度サブサンプリングは、全L個のフィルタされた投影
の中から1つを保持する。
演算子
は、トランケーション演算子、シフト演算子、及び角度デシメーション演算子の合成であり、
によって定義される。ここで、P′=P/Lである。演算子
は、デシメーションの前後に、重み付け関数W(θ,t),k=1,2 による乗法重み付け(multiplicative weighting)を含んでいる。1つの可能な選択は、
及び
である。もっと一般的に言えば、重み付は、アルゴリズムの正確度を最適にするよう選択可能であり、又は全く省略することも可能である。
図11は、部分シノグラムの処理され得る方法を一層詳細に示している。例えば、トランケーション
がまず初めにステップ1110にて起こり、その後にステップ1112でシフティング演算子
が続き、そして重み付けステップ1114、その後にデシメーションステップ1116、もう1つの重み付けステップ1118、そして第2のシフティングステップ1120となる。このプロセスは、シノグラム分解演算子
を記述している。
これらの定義により、部分画像f′上への逆投影のための正確な式(10)は、次の近似によって置き換えられる。
ここで、BM,P/Lは、式(9)によって定義されるP/L個の投影を用いた、M×M部分画像上への全体座標系内での逆投影である。これにより、式(12)に類似の区切られた画像のための逆投影演算子の近似分解が得られる。
この分解は図12に示されており、この図は、L=N/M=2での投影数の減少を伴うファンビーム逆投影のための近似分解を示している。この近似分解を実行するために、高速逆投影プロセッサ24は、図12中にgとして示された、P個の投影からなる発散ビームシノグラムを分割し、1210、1212、1214及び1216によって、g、g、g及びgで示されるそれぞれP/2個の投影からなる4つの部分シノグラムを生成するように、プログラムされる。部分シノグラムg、g、g及びgが、1218、1220、1222及び1224によって全体座標系内に逆投影されて、全体座標系内での正しい位置に部分画像f、f、f及びfが生成される。部分画像が1226で集合されて、図12中にfとしても示されている電子画像26が生成される。
再び図12を参照して、J=(N/M)の部分画像のそれぞれへの逆投影のための投影数はcMP/Lであり、全ての逆投影では合計cNP/Lである。これは、単一ステップの逆投影BN,PにおけるよりもL倍少ない。画像サイズでの投影の数の尺度化についての先の観測に基づけば、LはL=(N/M)と同じ大きさに選択可能である。
の計算コストが含まれる場合、部分画像の最適な数は、米国特許第6,263,096号に示されるような単一レベル分解のために決定可能である。
高速逆投影アルゴリズムの好ましい実施の形態においては、式(12)及び(18)中の分解は反復して適用され、各レベルで式(11)に従い部分画像をより小さな部分画像へと更に区切る。画像のダイアディック(N/M=2)2レベル区切りの一例が図13に示されている。異なるレベルでの幾つかの部分画像の中心位置も示されている。それらは全て、ゼロ番目のレベルの画像原点O、すなわち全体座標系に言及している。それに対応する、2レベル近似分解の一例と、1つの正確な分解及び1つの近似分解の段階からなる2レベル分解の一例が、それぞれ図14及び15に示されている。図14において、P個の投影を備える発散ビーム部分シノグラム22は、高速逆投影プロセッサ24で処理されて、それぞれがP/2個の投影を備える部分シノグラム1410、1412、1414及び1416に分割される。この近似の細分は、図11のフローチャートによって記述された近似分解演算子によって行われる。部分シノグラム1410の更なる近似分解が、1418、1420、1422及び1424で、図11のフローチャートによって記述される演算子を再び用いて生成されるが、今回は図13中のより小さな部分画像に対応する異なるパラメータで行われる。P/4個の投影を備える上記部分画像のそれぞれは、1426、1428、1430及び1432で全体座標系内で逆投影されて、図13に示された全体座標系内における正しい位置に複数の部分画像f1,1、f1,2、f1,3、f1,4が生成される。これら部分画像は1434で集合されて、やはり全体座標系内の正しい位置に部分画像fが生成される。
1416で表される部分シノグラムは、1436、1438、1440及び1442で同様に細分される。これらの部分シノグラムは、それぞれ1444、1446、1448及び1450で逆投影され、そして1452で集合されて部分画像fが生成される。点で示されたブロック1454及び1456は、それぞれf及びfのための対応するレベル2の分解を表している。図14において、分解の全てが、L=N/M=2におけるファンビーム投影の近似分解である。
図15は、L=N/M=2におけるファンビーム投影の2レベル混合(1つは正確な、もう1つは近似の)分解を示している。図15において、P個の投影を備えるシノグラム22は、やはりそれぞれがP個の投影を備える正確な部分シノグラム1510、1512、1514及び1516に分解される。これら正確な部分シノグラムは、それぞれg1、g2、g3及びg4で示され、(適当なパラメータと共に)図11に示された近似分解演算子を用いて、それぞれP/2個の投影を備える近似部分シノグラムに分解される。例えば、gは、1518、1520、1522及び1524で分解されて、近似部分シノグラムg′1,1、g′1,2、g′1,3及びg′1,4が生成される。それらの部分シノグラムは、1526、1528、1530及び1532で全体座標系内に逆投影されて、全体座標系内の正しい位置に部分画像f1,1、f1,2、f1,3及びf1,4が生成される。それらの部分画像は、1534で集合されて部分画像fが生成される。
1516で生成された部分シノグラムは、図15中でgで示されており、これはまた更に、1536、1538、1540及び1542で4つの部分シノグラムに分解されて、部分シノグラムg′4,1、g′4,2、g′4,3及びg′4,4が生成される。1544、1546、1548及び1550での逆投影は部分画像f4,1、f4,2、f4,3及びf4,4を生成し、これらが1552で集合されて部分画像fが生成される。点で示されたブロック1554及び1556は、それぞれg及びgのための対応する2レベル混合分解を表している。それらのプロセスは部分画像f及びf
を生成し、部分画像f、f、f及びfが1558で集合されて、図15中にfで示される画像26が生成される。
部分画像が或る所望の最小サイズMminとなるまで反復を継続可能であり、その後、逆投影BPmin,Mminが式(5)を用いて全体座標系内で行われる。最適なMminは、通常、インプリメンテーションに依存する(implementation-dependent)。それは、最小の部分画像が1ピクセル又は1ボクセルのサイズとなるよう、Mmin=1と同じくらい小さく選択することが可能である。Nは2の累乗であり、分解の全てのレベルでN/M=2であるとすると、logNのレベルが存在することになる。更にL=(N/M)=2であり、かつP=kNであり、kを小さな整数であるとすると、最後のレベルはN個の単一ピクセル逆投影BP/N,1を含むこととなり、合計コストがcPNになる。
の実行におけるデシメーション及びシフトの演算子のために固定長の補間が使用されるものとすると、各レベルでの計算コストがcNP(ここでcは定数)となることが示される。従って、logNのレベルの合計コストはO(PNlogN)であり、そしてP=O(N)なので、反復分解アルゴリズムの合計コストはO(NlogN)になる。
実際には、離散した検出器を使用するので、投影データは間隔Δtを有するtでサンプルされる。従って、投影シフト演算子
及び
は、離散形式で実行されることになるが、Δtの非整数倍によるシフトが必要となる。投影は逆投影の前にフィルタされているので、それらは帯域制限され、そして、例えば米国特許第6,287,257号に記載されているように、補間及び再サンプリングの組み合わせを用いてシフトされることが可能である。特に、Δtの非整数倍によるシフティングは、Δtの分数によるシフティングを伴って段階的に行われる、Δtの整数倍によるシフティングへと分離されることが可能である。分数シフトはt変数内の補間によって行われるが、整数シフトは単純な再インデクシングに相当する。更に、この補間は、式(14)内の一般化されたデシメーションフィルタリングステップにおけるフィルタh(p,t)を有する放射次元内でのコンボルーション(convolution)に対応する放射フィルタリングによって実行可能である。
重要なことには、ここで提案される特定の投影シフティングスキームによれば、投影のサンプリング間隔が均一かつ一定のままとなり、その結果、補間及び再サンプリングは離散対離散として行われ、単純なディジタルフィルタとして効率的に実行可能である。そのような分数tシフティングに含まれるエラーは、補間の長さで指数関数的に衰微させることが可能であり、その結果、通常は短い補間で十分となり、単純な線形補間でさえ十分となる。t補間の正確度はまた、tで投影をオーバーサンプリングすることによって向上し、これは必要ならば、投影のディジタル補間(アップサンプリングの後にフィルタリング)によって効率的に行われ得る。
近似分解を用いて得られる逆投影の正確度は、P,N,M,L、画像及び投影の帯域幅中の分解能、オブジェクトと交差するソース放射線の最大ファン角度、検出器のレスポンス、及びアルゴリズム中で使用される補間を含む、プロセス中の各種パラメータに依存する。これらのパラメータは、発散ビームアルゴリズム中で適当に最適化されてもよい。正確度に影響を与える重要なパラメータは、逆投影に使用される投影の数と、部分画像のサイズとの間の比(P/L)/Mである。角度オーバーサンプリング(angular oversampling)とも呼ばれる、この比を増大させることにより、近似分解の正確度が向上する。
以下の(i)〜(iii)を含む角度オーバーサンプリングを増大するための幾つかの方法がある:(i)分解の開始前のPを増加させるための、投影gのディジタル角度補間、(ii)近似分解の各段階でのL<(N/M)の選択、及び(iii)例えば図15中に例として示されているような正確な分解ステップの使用。方法(iii)を理解するために、式(12)における正確な分解は部分画像のサイズを減少させるが、部分シノグラム中の投影の数を減少させることはない、ということを思い出すこと。従って、これにより、因子(N/M)による角度オーバーサンプリングの増大と、改善された正確度とが得られる。
事実、上記方法の中の1つである(ii)又は(iii)に他と比べての小さな計算上の利点を与える実行詳細は別として、

に中心の置かれた部分画像の部分画像の中心である時にはいつでも成立する恒等式
のために、角度オーバーサンプリングのための方法(iii)は方法(ii)に等しい。従って、ダイアディック区切りを有する2段階の分解、-1つは正確、もう1つは近似- は、可能な最大値L=4の代わりにデシメーションL=2での、それぞれサイズがN/4の16個の部分画像への単一近似分解に等しい。従って、両方法は、同一因子2の角度オーバーサンプリングを提供する。一方、L=3を選択することにより、16個の部分画像への単一レベル近似分解は、因子4/3の角度オーバーサンプリングを提供可能であり、従って、動作点の選択に追加の柔軟性を提供する。一般には、角度オーバーサンプリングを増大させるための方法により、一定因子によるアルゴリズムの計算コストが増加する。従って、それらは計算と正確度との間のトレードオフを調整する方法を提供する。
本アルゴリズムの好ましい実施の形態は、幾つかの正確な細分、及び幾つかの近似の細分、又は細分因子よりも小さなデシメーション因子を有する近似の細分を含む。正確な細分及び/又はデシメーションの因子の数はオーバーサンプリングの因子を調整するものであり、また、計算と正確度との間の所望のトレードオフが得られるよう選択される。
最も単純な実施形態は、近似の細分レベルではL=2で、各レベルでN/M=2の一定因子による細分を使用するものであるが、もちろんその他の選択も可能である。
先に示したように、本アルゴリズムは、全てのレベルにおいてデシメーション因子Lの積によるP個への分割が必要である。これが、デシメーション因子の選択によって好都合又は効率的に達成され得ない場合は、角度補間及び投影の再サンプリングによってPを変更することが可能である。
本アルゴリズムの好ましい実施形態の反復構成を更に示すために、M/N=L=2における反復関数として書かれたアルゴリズムのための翻訳言語(pseudo-code)の一例が、図16に示されている。この翻訳言語は自明ではあるが、簡単に説明する。
関数FAN_FHBPは、全シノグラムG上の第1トップレベルの援助コール(invocation call)における反復関数演算子である。それ自身への次のコールにおいて、この関数は画像を部分画像に分解し、最小画像サイズNminに達するまでシノグラム分解を行う。これらの分解のQは正確であり、その残りは近似である。画像サイズNminに達すると、全体座標系内での正確な逆投影が、部分画像に対応する部分シノグラム上で行われる。Nminよりも大きな部分画像では、部分画像に相当するよう投影がトランケートされる。これは、正確なシノグラム細分を構成する。正確な分解が完了すると(Q<1の時)、近似分解のための追加ステップが行われる。図16において、投影はシフトされ、重みデシメートされ(weight-decimated)、そして再度シフトされる。この部分画像は
で集合される。
部分画像の中心はまた、
として表現される、
に中心の置かれたより大きな部分画像の部分画像iの中心
を用いて反復的に計算される。ここで、
は、
に対する部分画像fの近似の位置ベクトルである。フルサイズ画像の中心は、定義により、全体座標系の原点にある。従って、その後の全ての部分画像の中心は全体座標系において表現される。
高速逆投影アルゴリズムの好ましい実施の形態においては、反復分解が以下の観点により単純化される。tでの投影シフト、投影トランケーション、及び投影重み付けの演算が、(適当な座標変換を用いて)代わりに使用できる。従って、1つの段階の演算
がその後の段階の
と組み合わされて、1つのシフト及び1つの重み付けの演算になる。同様に、最後の部分画像の逆投影BMmin,Pminの直前に生じる、このシフティング及び重み付けの演算
は、逆投影自身における重み付け及び補間の演算と組み合わせ可能である。従って、反復実行において、
は、1つの投影シフトステップと、各段階毎の(多くとも)1つの重み付けのみが必要となる。
信号処理アルゴリズムの設計及び実行における当業者にとっては明らかなことであろうが、本発明のプロセスを定義する各種演算は、関数又はプロセスの原理を変化させなければ、添付の図によって記載された順番とは異なる各種順番で再配置することや、並列に実行することも可能である。幾つかのそのような再配置は、実行、演算、又はコストの観点から有利である場合がある。例えば、部分画像のための投影の数をデシメートする演算は、P個の全ての投影がトランケート及びシフトされるのを待つことなしに、2以上の投影がこの部分画像のためにトランケート及びシフトされるとすぐに開始されるようにすることが可能である。そのような再配置は、再構成アルゴリズムを配備したり、データの捕捉でオーバラップしたりするのに使用可能である。
更に、投影は適時に連続的に得られるので、本発明の方法は連続的に更新される断層撮影画像に適用可能である。これは、最も最近に得られた投影とそれよりも古い投影群との間の相違から構成されるシノグラムを逆投影するための提案された方法を用いることにより達成可能である。その逆投影の結果が、続いて現在の画像に加算されて、更新画像が生成される。そのような更新は、連続的手法で続行可能である。

〔高速再投影〕
逆投影と再投影は、双対(dual)の演算、又は互いの随伴(adjoint)である。よって、その実行のための高速の本来のファンビームアルゴリズムは密接に関連している。従って、高速の本来のファンビームアルゴリズムについては、高速逆投影アルゴリズムの全ての記述が再投影アルゴリズムにおいても同様であるとの理解のもとに、それほど詳しくは記述しない。
図8に示された、fのM×M部分画像
について考える。そのソース角度θq=qΔθ,q=0・・・Q−1でのファンビーム投影は、
である。
位置
に中心を有するM×M部分画像のための全体座標系内の、その関係する再投影演算子を
として示す。特に、
は、全体画像のファンビーム投影である。PM,Qの計算コストは、或る定数cにおいてcMQである。
次に、図9に示された式(11)の画像区切りについて考える。投影演算子の直線性のために、以下の、同一数Pのソース角度でより小さいM×M画像のJ個の再投影の合計への、再投影の正確な分解は、次のように得られる。
再投影のためのこの正確な分解は、N/M=2のものが図17に示されている。図17には、図17中にfでも示されている電子画像30は、トランケーション演算子1710、1712、1714及び1716によって、それぞれ部分画像f、f、f及びfに細分される。この部分画像が1718、1720、1722及び1724によって全体座標系内に再投影されて、部分シノグラムg、g、g及びgが生成される。この部分シノグラムが1726で集合されて、図17中にgでも示されているシノグラム34が生成される。
逆投影の式(12)の正確な分解の場合と同じ理由により、再投影の分解は単一ステップの逆投影
と比べてスピードアップを提供するものではないが、特に近似分解との組み合わせにより、その他の可能な適用を有している。
高速再投影アルゴリズムは、先に述べた以下の追加の特性を有する分解の考えを使用する。すなわち、固定された画像分解能(帯域幅)においては、画像を特徴づけるファンビーム投影の数、又は同様に、θ変数における投影の帯域幅は、画像のサイズに比例する、ということである。従って、もしP個の投影がN×N画像fを特徴づける場合は、P′=P/(N/M)個の投影が部分画像f′を特徴づけ、その投影の残りはこの小さい方の組から得られる。
従って、高速再投影アルゴリズムは、M×Mの部分画像f′の1組のP′個の投影を計算することと、より大きい組のP個の投影を得るのに「投影拡大(projection augmentation)プロセス」を使用することに基づいている。この投影拡大プロセスは演算子
によって行われ、これはそれ自体が幾つかの演算子から構成されている。投影シフト演算子
及び
は、先に式(13)で定義されている。L倍角度の補間演算子
は、L倍角度のアップサンプリング(upsampling)とその後の角度フィルタリングにより、一組のQ個の投影における視界(views)の数をQL個の投影に増加させる。角度アップサンプリングは、L−1個の、全ての2投影間のゼロ値(zero-valued)の投影を導入する。アップサンプルされた投影をg(p,t)で示すと、フィルタリングは式(14)及びそれに関係する議論によって記述される。その結果得られる角度補間演算子
は、高速逆投影アルゴリズムにおいて生じる角度デシメーション演算子
の随伴(adjoint)である。投影拡大演算子
は、シフト及び角度補間演算子の組み合わせであり、
によって定義される。ここで、P′=QLである。演算子
は、補間の前後に、重み付け関数W(θ,t),k=1,2 による乗法重み付けを含んでおり、これは同様な関数及び
におけるような選択基準を有している。例えば、図18に示されるように、第1のシフティング関数1810の後に、第1の重み付け関数1812、補間関数1814、第2の重み付け関数1816、及び第2のシフティング関数1818が続く。もしその方がよければ、投影拡大演算子の他の組み合わせも使用可能である。
これらの定義により、部分画像f′の近似再投影は、
である。ここで、PM,P/Lは、P/L個の投影のみを計算する、全体座標系内の再投影である。これにより、式(20)に類似の区切られた画像のための再投影演算の近似分解が得られる。
この近似分解は、L=N/M=2におけるものが図19に示されている。図19において、図1中の高速再投影プロセッサ32が、図19中にfでも示された電子画像30を受け入れて、1910、1912、1914及び1916においてそれぞれ部分画像f、f、f及びfに分割するようプログラムされている。全体座標系内での再投影は、1918、1920、1922及び1924で行われて、それぞれP/2個の投影を有する部分シノグラムg′、g′、g′及びg′が生成される。この部分シノグラムは、前述したように、1926、1928、1930及び1932で拡大され、それによって得られるP個の投影を有する部分シノグラムg、g、g及びgのそれぞれが1934で集合されて、図19中にgでも示されるシノグラム34が生成される。
逆投影の近似分解の場合のように、計算コストは単一ステップの再投影PN,PよりもL倍小さく、同じ議論が適用される。
高速再投影アルゴリズムの好ましい実施の形態においては、式(20)及び(23)における分解が反復して適用され、式(11)に従って、各レベルで部分画像がもっと小さな部分画像へと区切られる。これは図13に示されている。これに対応する、2レベル近似分解の例と、1つの正確な分解及び1つの近似分解の段階から構成された2レベル分解の例が、それぞれ図20及び21に示されている。図20においては、L=N/M=2でのファンビーム再投影の2レベル近似分解が、図20中にfでも示された電子画像30で開始される。この画像fは2010、2012、2014及び2016でそれぞれ部分画像f、f、f及びfに分割され、部分画像fは更に2018、2020、2022及び2024で部分画像f1,1、f1,2、f1,3及びf1,4に更に細分される。それらの部分画像はそれぞれ2026、2028、2030及び2032によって全体座標系内に再投影されて、それぞれP/4個の投影を備える部分シノグラムg″1,1、g″1,2、g″1,3及びg″1,4が生成される。これらの部分シノグラムは2034、2036、2038及び2040によって拡大されて、それぞれP/2個の投影を備える部分シノグラムg′1,1、g′1,2、g′1,3及びg′1,4が生成される。これらは2042で集合されて、P/2個の投影を備える部分シノグラムg′が生成される。2044での更なる拡大の後、P個の投影を備える部分シノグラムgが2046で部分シノグラムg、g及びgと集合される。
部分画像fは2048、2050、2052及び2054で更に分解され、それら部分画像が2056、2058、2060及び2062で同様に再投影されて、それぞれ部分シノグラムg″4,1、g″4,2、g″4,3及びg″4,4が生成される。2064、2066、2068及び2070での拡大の後、それによって得られた部分シノグラムg′4,1、g′4,2、g′4,3及びg′4,4が2072で集合され、そして、それによって得られた部分シノグラムが2074で更に拡大される。その結果得られたgが2046で他の部分シノグラムと集合されて、図20中ではgでも示された、図1中の部分シノグラム34が生成される。点で示されたブロック2076及び2078は、それぞれf及びfのための、対応するレベル2の分解を表している。
図21は、L=N/M=2でのファンビーム再投影の2レベル混合(1つは正確な、もう1つは近似の)分解を示している。図21において、図中にfでも示された、図1中の電子画像30が、2110、2112、2114及び2116によって、それぞれ部分画像f、f、f及びfに分解される。部分画像fは、2118、2120、2122及び2124によって部分画像f1,1、f1,2、f1,3及びf1,4に細分される。これらの部分画像は、2126、2128、2130及び2132によって全体座標系内に再投影されて、それぞれP/2個の投影を備える部分シノグラムg′1,1、g′1,2、g′1,3及びg′1,4が生成される。2134、2136、2138及び2140による拡大の後に、それぞれP個の投影を備える部分シノグラムg1,1、g1,2、g1,3及びg1,4が生成される。2142での集合により、P個の投影を備える部分シノグラムgが生成され、これが2144で部分シノグラムg、g及びgと集合されて、P個の投影を備えるシノグラムgが生成される。同様に、部分画像fは、2146、2148、2150及び2152で更に細分されて、部分画像f4,1、f4,2、f4,3及びf4,4が生成される。2154、2156、2158及び2160による再投影によって、部分シノグラムg′4,1、g′4,2、g′4,3及びg′4,4が生成され、これらが2162、2164、2166及び2168によって拡大されることで、部分シノグラムg4,1、g4,2、g4,3及びg4,4が生成される。これらの部分シノグラムは、2170で集合されて部分シノグラムgが生成される。点で示されたブロック2172及び2174は、それぞれg及びgのための、対応するレベル2の分解を表している。
部分画像が何らかの所望の最小サイズMminになるまで反復を続け、その後に、再投影PPmin,Mminを式(7)を用いて計算することが可能である。
高速再投影アルゴリズムのための以下の考察は、逆投影の場合に議論したものと類似している。

1.P、N、M、L及びMminを含む、高速再投影アルゴリズムにおける各種パラメータの選択。
2.tにおける離散の操作。
3.正確度とスピードとの間のトレードオフを制御するための、角度オーバーサンプリングと、正確及び近似の分解ステップの組み合わせの使用。
4.連続する段階間でのシフティング及び重み付けの演算を組み合わせることによる演算子
及び
の単純化を行い、その結果としての、各段階で単一のシフティング及び(多くとも)単一の重み付け演算。
5.演算の順番の可能な再配置。
6.O(N2logN)の、その結果としての階層再投影アルゴリズムの計算の複雑性。

〔等角ファンビームジオメトリのための高速の本来のアルゴリズム〕
等角検出器の場合にも、画像をもっと小さな部分画像に階層的に分解するために、同様なスキームを適用可能である。図22は、この場合における部分画像f′のための投影ジオメトリを示している。
式(8)と同様な、Kの定義を用いて、逆投影演算子
は式(9)によって与えられる。ここで、
は、等角ジオメトリ(図6参照)のために先に定義されている。f′上への逆投影に寄与するg(p,・)の一部分は、再び
である。ここで今、
は、部分画像
の支持体の投影の角度t及びtによって決定される角度支持体へ、ファン角度tでg(p,・)をトランケートする。同一直線等間隔ファンビームジオメトリにおけるのと同じ議論及び討議が適用され、これにより、式(10)−(12)と、図9に示された正確な分解とが得られる。
近似分解はまた、以下の手順を用いて、部分画像f′を再構成するために使用される投影の数を減少させる。tOO′が、f及び部分画像f′のそれぞれの中心O及びO′を通過する2つの放射線間の角度であるとする。同一直線等間隔検出器の場合におけるf′の中心O′の投影に対応する間隔
によるトランケートされた投影をtシフティングするための直接の類比において、部分画像のためのトランケートされた投影は、tOO′によってファン角度でtシフトされる。同一の角度間隔Δtで均一に配置されたサンプル点が得られ、これはt=0で中心にある。従って、等角の場合は、角度変数tが、等間隔ジオメトリに見られるtの役割を担い、
が図7のように定義される。或いは、前述したのと同じ反復分解プロセスが、ここで適用される。特に、部分画像に対応するファンビーム投影が一度トランケートされると、それらはソース角度θに関してLによりデシメートされ、その残りの手順に続く。各種演算子の定義は、既にリストされた変化を除き、ほとんど変更がない。近似分解は、図11によって再び記述される。
同様な議論及び変更が、高速の本来の再投影アルゴリズムに適用される。

〔平面等間隔検出器を有するコーンビームジオメトリのための高速の本来の逆投影及び再投影アルゴリズム〕
一般の3D発散ビームジオメトリのためのプロセスは、2Dの場合において述べたのと同じ要素を有しており、以下の主な相違を備えている。
1.画像は、全ての3座標に沿って3Dで再構成される。
2.投影が2次元であるため、個々の投影のトランケーション、シフティング及び補間が2次元である。
3.投影は、円上のソース角度ではなく、軌道に沿ったソース位置θでデシメート(又は補間)される。
4.ソース軌道はしばしば3Dにおいて等方的でないため、好ましい分解は等方的でなく、空間的に変化している場合がある。
図23に示された、fの部分画像
のための逆投影演算について考える。
が式(8)によって定義された画像トランケーション演算子であるとし、これにより、
は、サイズM×M×Mで、
に中心の置かれた、fの部分画像となるようにする。式(9)及び(10)がまだ適用され、この場合、τが2Dベクトル
によって置換され、また、部分画像
のM×M×Mの支持体の投影の2D支持体Ωへのオブジェクト全体の投影をトランケートするトランケーション演算子
を用いる(図23を参照)。
今、画像fを、それぞれのサイズがM×M×Mである、J=(N/M)個のオーバーラップの無い部分画像に区切ることについて考える。
ここで、ベクトル
は、図24に示された、M/N=2の場合における部分画像の中心である。
部分画像上への逆投影への逆投影のための正確な分解を提供する式(12)が適用される。それは、4つの「経路(channels)」の代わりに8つの経路を備えていることを除き、図9に類似した図によって記述されている。計算コストの議論は、BN,P及びBM,Pの計算コストがそれぞれcNP及びcMPであることを除き、2Dの場合におけるものと同様であり、2Dの場合とちょうど同じ比Jを有している。
2Dファンビームの場合もそうであるが、高速3Dコーンビーム逆投影アルゴリズムは、以下の更なる特性を有する分解の考えを使用する。固定された画像分解能(帯域幅)においては、画像を再構成するのに必要とされる投影の数は画像のサイズに比例する。すなわち、N×N×Nの画像を再構成するのにP個の投影が必要とされる場合、同一分解能のM×M×Mの画像を再構成するにはP′=(M/N)P個の投影で十分である。
減少を行うための演算子
を用いて、部分画像
が、減少された数の投影から再び再構成される。それは、以下の変更を伴う、以前のものと同様な演算子の構成によって定義される。
1.今、投影トランケーション演算子
は、先に説明したように、2D領域Ωにトランケートする。
2.投影シフト演算子
は、ファンビームジオメトリの場合における1Dシフトの代わりに、t、t座標内でベクトル
だけ2Dシフトを行う。このシフトは、例えばまずtでシフティングし、次にtで、といったような分離可能な方法で行われてもよい。
3.演算子
は、角度ではなく、ソース位置でデシメートする。軌道が、直角をなす2つの円のような、2つ以上の別々の部分の結合から構成されている場合、ソース位置でのデシメーションは、各部分毎に分離して行われる。
4.式(14)当たりのデシメーションに含まれる
変数におけるフィルタリングは、1次元ではなく2次元である。それも、
フィルタリングはθでのフィルタリングから分離可能であってもなくてもよいが、t及びtでの分離可能な形式で行われる場合がある。
近似分解は、再び式(18)によって与えられる。今ここで、BM,P/Lは、式(9)によって定義されたP/L個の投影を用いた、M×M×Mの部分画像上への全体座標系内での発散ビーム逆投影である。L=N/M=2において、この分解は、4つの「経路(channels)」の代わりに8つの経路を備えていることを除き、図10に似た図によって記載される。計算コストについての議論は、BN,P及びBM,P/Lの計算コストがそれぞれcNP及びcMP/Lであって、J=(N/M)であることを除き、2Dの場合のものと同様である。J個の部分画像の逆投影の合計コストは、従来の全体画像の逆投影BN,Pにおけるものよりも再びL倍小さい。
高速逆投影アルゴリズムの好ましい実施の形態においては、式(12)及び(18)の分解が反復して適用され、各レベルで式(24)に従って部分画像をもっと小さな部分画像へと更に区切る。ファンビームジオメトリにおいて先に示された、プロセス、各種の特徴、及び改善についての記載は、正確又は近似の分解それ自身においてリストされたのと同様な変更を伴って、同様に一般のコーンビームジオメトリに適用される。特に、反復分解アルゴリズムを用いたコーンビーム逆投影の合計コストはO(NlogN)となり、これは従来のコーンビーム逆投影におけるO(N)と比較される。
ここまで、我々は、画像が3軸の全てに等しく分解されるという、3D画像分解の最も単純な実施の形態を記述してきた。しかし、コーンビームジオメトリの幾つかの例においては、ソース軌道が3Dに等方的でないために、好ましい分解も等方的ではなく、更には空間的に変化している場合もある。重要な例は、単一円及びヘリカルコーンビーム軌道である(図4、6及び3)。これらの場合は、z軸回りのソース軌道の円柱対称(又はほぼ対称)を示し、z軸で選択された画像区切りはx軸及びy軸でのものと異なる場合がある。軌道の全体が対称軸を持たない時でさえ、それは時々、直角をなす2つの円から構成されるソース軌道の場合のように、部分の結合へと分解される場合がある。それで、逆投影演算は、それぞれそれ自身の分離した軸対称軌道を有する2つの逆投影の加算として行われ得る。
そのような例においては、画像の部分をもっと少なく区切るか、又はz軸においてもっと大きなサイズの部分画像へと区切るだけで十分である。極端な例では、図25に示されるように、区切りは全体的に(x,y)方向に限定されてもよい。もっと一般的な、非等方的で空間的に不均一な区切りが、図26に示されている。これは、t軸に沿った投影への必要とされるシフト、補間又はフィルタリングの演算の数を減少させることにより、計算の節約へとつながる。
そのように区切り作業を減少させる理由は、単一円形のソース軌道の場合における式(2)及び(3)について検討することにより明らかになる。これらの式は、検出器上への点の投影の位置がθに依存することを示している。小さなz及び/又は小さな
を有する点を考える。すなわち、ソース軌道の平面に近い点、及び/又はz軸に近い点である。なお、そのような点では、
が、
よりも、θでもっとゆっくりと変化する。これは、たとえそのz次元がその(x,y)次元と同じ比率だけ減少されたとしても、中心
がそのように配置された部分画像が、減少された数の視界(views)から正確に再構成され得る、ということを示すために使用可能である。結果として、(x,y)平面内におけるもっと細かい区切りが、減少された数の投影からの正確な再構成を確実にするために十分である。もし、tでのhθ/(2π)だけの全体的なシフトが投影にまず適用されて、θの変化に伴う、z軸に沿った検出器上での座標の中心の動きが補償されるのであれば(図3を参照)、同様な議論がヘリカル軌道の場合にも適用される。
そのような非等方的で不均一な区切りでは、それに応じて反復の分解が変更される。もし図25の区切りが反復的に使用されるのであれば、それに対応する分解は、反復の各レベルで4つの経路を有し、それはt変数でのみシフトを含むことになり、また、これは図10−15によって記載されることになる。一方、もし図26の区切りが2レベルの近似分解に使用されるのであれば、分解の第1レベルは、(立方体の粗い区切りからなる8つの八分体(octants)に対応する)8つの経路を有することになるが、次のレベルでのそれぞれの分解は、それぞれの粗い画像八分体を区切って得られた7つの画像に対応する7つの経路を有するのみとなる。
先に議論した発散ビームジオメトリの特定の例を含む、一般の発散ビームジオメトリのための高速再投影の設計は、それに対応する高速逆投影アルゴリズムの双対(dual)である。それは、高速ファンビーム再投影が高速ファンビーム逆投影に対してもたらす関係と同様な関係を、高速発散ビームジオメトリ逆投影に対してもたらす。後者の関係は既に少々詳細に記載してあるので、ここではそれを繰り返し述べない。

〔実行及び実験〕
新規なアルゴリズムは、MATLABTM及びCプログラミング言語で実行され、2Dの場合にはサイズN×N=512×512で、3Dの場合にはサイズN×N×N=128×128×128のシェップローガンヘッドファントム(Shepp-Logan head phantom:断層撮影アルゴリズムの数値評価に使用される標準テスト画像)の再構成を含む幾つかの数値実験でのテストに成功した。我々は、頭蓋骨と脳との間のコントラストが明確な(及び高い)、いわゆる未修正ファントム(unmodified phantom)を使用した。これは、このアルゴリズムのための魅力的なテスト(challenging test)を提示するものであり、それは頭蓋骨の再構成によって生成された僅かなアーティファクトでさえ脳の背景で目立つためである。
まず、2Dファンビーム再構成の結果について記述する。等間隔と等角度のファンビームイメージングジオメトリの両方において、ソースと画像中心間の距離を画像サイズの1.25倍に、すなわちD=1.25Nに設定し、ファン上に1025個の検出器を有する場合における合計のファン角度を2α=1.17ラジアンに設定した。P=1536個の均一間隔で配置されたソース角度、10楕円(10-ellipse)のファントムにおけるθ∈[0,2π]、を有するフルスキャンのファンビーム投影データを、楕円を通る線積分のための分析表現を用いて生成した。
我々のシミュレーションにおいては、各レベルで画像を半分サイズの4つの部分画像に分解した。近似分解では、デシメーション因子L=2を用いた。画像は、単一ピクセルの部分画像、すなわちMmin=1、に繰り返し分解された。角度オーバーサンプリングを増大させるために、分解の最初のQレベルは、角度デシメーションなしの正確な分解とした。Q=2及びQ=3を用いた。ここで、Qは「ホールドオフパラメータ(holdoff parameter)」と言うことにする。投影トランケーション演算子
のためのパラメータ(「足跡(footprint)」)は、計算が少々増えるコストで、単純な控えめな近似を用いて大急ぎで計算された。
図27及び28は、正確な逆投影アルゴリズムと、2つのファンビームイメージングジオメトリのための我々の高速ファンビーム逆投影アルゴリズムとを比較している。参考及び目安として、図27(a)及び28(a)は元々のファントムを示し、図27(b)及び28(b)は、等間隔及び等角度の検出器の両方のための正確なファンビーム逆投影を用いて再構成された画像である。図27(c)、28(c)、27(d)及び28(d)は、それぞれホールドオフパラメータQ=3及びQ=2を有する我々の高速アルゴリズムを用いた再構成である。高速アルゴリズムのスピードアップは、ホールドオフがQ=3及びQ=2においてそれぞれ30×及び60×である。正確かつ高速のアルゴリズムを用いた再構成は、知覚可能な品質相違について、Q=2においてさえもほとんど示していない。図28中の再構成された画像の列及びコラム(縦行)を通してのプロットが、図29(a)及び29(b)に示されている。図29(a)は第410番目の列を通るスライスを比較しており、図29(b)は図28の第356番目のコラムを通るスライスを比較している。これらのプロットにおいて、実線は元々のファントムを表し、破線は正確なファンビーム逆投影を表し、鎖線はホールドオフ3の高速アルゴリズムを表し、点線はホールドオフ2の高速アルゴリズムを表している。
次に、等間隔検出器ジオメトリのための正確及び高速のファンビーム逆投影アルゴリズムの点分布関数(PSF’s:point-spread functions)を比較する。位置(10,10)、(128,128)、(49,37)及び(201,93)に配置された4つの点ターゲットを有するサイズ256×256の画像上で、正確及び高速のアルゴリズムを実行した。先に使用したものと同様な以下のパラメータを用いた。すなわち、ソースと中心間の距離D=1.25N、ソース角度の数P=2N、及び737個の等間隔検出器での合計ファン角度が1.20ラジアンである。点応答(point responses)は、4つの位置の全てで非常に類似しており、本質的にシフト不変の応答を示唆した。(Q=3の)正確及び高速のアルゴリズムで得られた(128,128)での点ターゲットにおける点分布関数(PSF’s)が、図30(a)及び30(b)において比較されている。図30(a)は正確なアルゴリズムにおける中心ピクセルPSFを示し、図30(b)はホールドオフが2での高速アルゴリズムにおける中心ピクセルPSFを示している。これらPSF’sの2つの軸を通るスライスが図31(a)及び31(b)に示されている。図31(a)は、正確なアルゴリズム(実線)及びQ=3の高速アルゴリズム(破線)における第128番目の列を通るスライスを表している。図31(b)は、正確なアルゴリズム(実線)及びQ=3の高速アルゴリズム(破線)における第128番目のコラムを通るスライスを表している。2つのアルゴリズムのPSF’sはほとんど同一であり、これは高速アルゴリズムの正確さを裏付けている。
次に、(x,y)平面内で単一円形ソース軌道を有するコーンビーム再構成を用いた実験の結果について議論する。ソースからボリュームの中心までの距離が、側面部の長さの1.25倍である。コーンビーム投影は、[0,2π]内に均一に配置された512個のソース角度で行われた。(方位角における)ファン角度は1.20ラジアンであり、コーン角度は1.5ラジアンであった。それぞれ図4及び図6に示された平面等間隔及び円柱状の検出器ジオメトリの両方について考える。それぞれの場合において、検出器面は、t座標において等間隔又は等角度のいずれかの間隔で、垂直軸(t)に沿って等間隔配置された、375×375の検出器から構成された。
平面検出器を有する単一円コーンビーム軌道のための従来の再構成アルゴリズムは、よく知られたフェルドカンプ、デービス及びクレス(FDK:Feldkamp, Davis and Kress)アルゴリズムであり、これは、投影に重み付けをし、t1座標に沿ってフィルタリングし、そして次に、重み付けされたコーンビーム逆投影を行うことから構成されている。この重み付けとフィルタリングは、平面検出器を有するファンビームジオメトリにおけるものと同じであり、相違するのは、ファンビーム逆投影ではなくコーンビーム逆投影を用いる点である。これは近似アルゴリズムであり、大きなコーン角度では幾つかの固有のアーティファクトを受ける。それでも、それはコーンビームジオメトリのための標準的な目安である。円柱状検出器ジオメトリにおいては、重み付けされたコーンビーム逆投影を円柱状検出器に組み込むことにより、等角検出器ジオメトリのためのファンビームアルゴリズムをコーンビームジオメトリに拡張することによって得られた、類似のアルゴリズムを用いた。
デモンストレーション目的のために、対応する高速アルゴリズムは、前述した反復分解を有する逆投影を用いており、これは図25に示された単純な区切りスキームで行われる。我々のシミュレーションでは、アルゴリズムの特別な実行中に含まれるオーバーヘッドが小ボリュームのための逆投影演算子のスピードアップを相殺するので、部分画像サイズが4×4×128にまで減少された時に分解プロセスを終了した。高速アルゴリズムは、従来のFDKアルゴリズムと同じフィルタリング及び重み付けを用いた。
従来のFDK再構成の結果と、本発明のプロセスを用いる高速再構成の結果が、それぞれ平面及び円柱状検出器ジオメトリにおいて図32及び図33で比較されている。図32は、単一円形軌道及び平面等間隔検出器ジオメトリにおける3Dコーンビーム再構成を比較している。最上位の列はz=−0.25でのxyスライスであり、中央の列はx=−0.13でのyzスライスであり、最下位の列はy=−0.04でのxzスライスである。左のコラムは従来のFDKアルゴリズムを用いた再構成の結果であり、中央及び右のコラムはそれぞれホールドオフ2及び1での本発明の高速アルゴリズムを用いている。図33は、単一ソース軌道及び円柱状検出器ジオメトリにおける3Dコーンビーム再構成の比較である。最上位の列はz=−0.25でのxyスライスであり、中央の列はx=−0.13でのyzスライスであり、最下位の列はy=−0.04でのxzスライスである。左のコラムは従来のFDKアルゴリズムを用いた再構成結果であり、中央及び右のコラムはそれぞれホールドオフ2及び1での本発明の高速アルゴリズムを用いている。図34は、図32中の対応する列における画像を通るスライスを表している。最上位のプロットは列102を示し、中央のプロットは列80を示し、最下位のプロットはコラム81を示している。実線は従来のFDK再構成のものである。ホールドオフ2での高速アルゴリズムが鎖線で示されており、ホールドオフ1での高速アルゴリズムが点線で示されている。図35は、図33中の対応する列における画像を通るスライスを示している。最上位のプロットは列102を示し、中央のプロットは列80を示し、最下位のプロットはコラム81を示している。実線は従来のFDK再構成のものであり、鎖線はホールドオフ2での高速アルゴリズムを用いた結果を表し、点線はホールドオフ1での高速アルゴリズムを用いた結果を表している。更なる詳細が、図34及び35によって提供され、これらは図32及び33の画像を通る或る線に沿った、再構成された密度のプロットを示している。ホールドオフ=1及びホールドオフ=2における結果が与えられている。
ホールドオフ因子Q=2での本発明を用いたFDKアルゴリズムの高速バージョンは、従来のフェルドカンプ(Feldkamp)アルゴリズムよりも約3倍速く、画像品質における劣化はほとんど見られない。ホールドオフ因子Q=1でさえも、画像品質は幾つかの応用において容認できるものであり、7倍のスピードアップが図れる。512×512×512の画像ボリューム、図24又は26におけるように十分になされた3D区切り、及び対応する分解においては、スピードアップは10倍の向上が見込める。(事実、円形又はヘリカル軌道のような、ほぼ円柱対称を有するコーンビームジオメトリにおいては、何らかのzサイズHにおける512×512×Hでは10倍のスピードアップが得られる。)

〔実験の結論〕
2Dにおける我々の実験において、新規なファンビーム逆投影アルゴリズムは、512×512の画像において30×−60×の実際のスピードアップを提供し、正確度における認識可能な損失はなかった。3Dにおいて、単一円軌道及び平面検出器のコーンビームジオメトリで、分解の最も単純な形式を用いた場合、スピードアップは128×128×128において3×−7×であった。本発明で提案された分解を十分に実行することで、512×512×512の画像において、10倍のスピードアップ因子が可能である。
この実験により、異なる検出器ジオメトリにおける本発明のプロセスの有効性が実証された。単純な単一円形軌道ジオメトリにおけるこれらの結果は、ヘリカルコーンビームのような、もっと一般的な発散ビームジオメトリに拡張される。再投影アルゴリズムと逆投影アルゴリズムとの間の、原理及び構造の類似性のために、本発明の高速発散ビーム再投影アルゴリズムにおいても同様な結果が維持される。再投影アルゴリズムと逆投影アルゴリズムの両方が、個々にも組み合わせても有用であり、それらは、例えば米国特許第6,263,096号及び第6,351,548号において説明したように、断層撮影再構成への各種の反復アプローチや、断層撮影におけるアーティファクト修正に使用可能である。逆投影と再投影の組み合わせは、ヘリカルコーンビームジオメトリにおける長いオブジェクトの問題のための高速再構成アルゴリズムを構築するためにも使用可能である。
以上のように、本発明の原理は特定の装置及び応用と関係付けて記述されているが、この記述はほんの一例であって、本発明の範囲を限定するものではないことを理解すべきである。
本発明で使用する装置のブロック図である。 平面等間隔コーンビームジオメトリ(a planar eqispaced cone-beam geometry)の図である。 平面等間隔ヘリカルコーンビームジオメトリ(a planar eqispaced helical cone-beam geometry)の図である。 円形ソース軌道及び平面等間隔検出器を有するコーンビーム断層撮影を示す図である。 同一直線等間隔ファンビームジオメトリ(colinear eqispaced fan-beam geometry)の図である。 円形ソース軌道及び円柱状検出器を有するコーンビーム断層撮影を示す図である。 等角コーンビームジオメトリ(an equiangular fan-beam geometry)の図である。 MxM部分画像における同一直線等間隔ファンビーム投影ジオメトリの図である。 画像の部分画像への区切りの図である。 ファンビーム逆投影の正確な分解を示す図である。 シノグラム分解演算子O[L,M,δ]によって行われる処理のフローチャートである。 ファンビーム逆投影の近似分解の図である。 2レベルのダイアディック(dyadic)画像の区切りの図である。 ファンビーム逆投影の2レベル近似分解の図である。 ファンビーム逆投影の2レベル混合(1つは正確な、1つは近似の)分解の図である。 高速ファンビーム階層逆投影のための翻訳言語である。 ファンビーム再投影の正確な分解を示す図である。 投影拡大演算子(projection augmentation operator)V[L,M,δ]によって行われる処理のフローチャートである。 ファンビーム再投影の近似分解の図である。 ファンビーム再投影の2レベル近似分解の図である。 ファンビーム再投影の2レベル混合(1つは正確な、1つは近似の)分解の図である。 部分画像のための等角ファンビーム投影ジオメトリの図である。 部分画像のための平面等間隔コーンビーム投影ジオメトリの図である。 3D画像の部分画像への区切りの図である。 3D画像の部分画像への非等方区切りの図である。 3D画像の2レベル非一様反復区切りの図である。 同一直線等間隔検出器における正確及び高速のファンビーム逆投影を比較する写真である。 等角検出器における正確及び高速のファンビーム逆投影を比較する写真である。 図28(a)〜28(d)の逆投影のスライスを通って得られたプロットである。 点分布関数(point spread function)の比較である。 点分布関数を通るスライスのプロットの比較である。 3Dコーンビーム再構成の写真の集まりである。 3Dコーンビーム再構成の写真の他の集まりである。 図32中の画像を通るスライスのプロットの集まりである。 図33中の画像を通るスライスのプロットの集まりである。

Claims (29)

  1. 逆投影を施せる発散ビームシノグラム(22)から電子画像を生成する方法であって、
    所望回数の正確な細分と所望回数の概算的な細分を行うことにより、前記シノグラム(22)を複数の部分シノグラム(図10中のg、g、g、g)に細分するステップ(1010、1012、1014、1016)と、
    撮像されるオブジェクト上に中心を置く全体座標系内で、前記部分シノグラムの重み付けされた逆投影を行って(1018、1020、1022、1024)、前記全体座標系内の正しい位置に中心を置く複数の対応する部分画像(図10中のf、f、f、f)を生成するステップと、
    前記部分画像(図10中のf、f、f、f)を集合させて(1026)、前記電子画像(図10中のf、26)を生成するステップと、
    を備え
    前記所望回数の正確な細分は、前記シノグラムにおける投影のトランケーションを行うこと、によってなされ、
    前記所望回数の概算的な細分は、前記シノグラムにおける投影のトランケーションを行い、該トランケーションの行われた投影を前記全体座標系の中心に対して、選択された距離だけシフトさせ、該シフトされた投影のデシメーションを、選択された係数によって行い、該デシメーションの行われた投影をその元の位置へ、前記選択された距離を用いてシフトさせること、によってなされる、方法。
  2. 前記シノグラムは、各部分シノグラムが所望サイズの部分画像を示すまで、反復手法で複数の部分シノグラムに細分され請求項1記載の方法。
  3. 前記部分シノグラムは、サイズが1ピクセル又は1ボクセルと同程度に小さい部分画像に対応する請求項記載の方法。
  4. 前記集合のステップは反復手法で行われる請求項1記載の方法。
  5. 前記電子画像は断層撮影画像である請求項1記載の方法。
  6. データ座標内でのオーバーサンプリングを使用して前記電子画像の精度を向上させる前処理ステップ(16)を更に備える請求項1記載の方法。
  7. 前記シノグラムは前処理された発散ビーム投影を含む請求項1記載の方法。
  8. 前記所望回数の概算的な細分は、前記シフトされた投影をデシメーションの前に重み付けすること、及び、前記デシメーションの行われた投影をシフトの前に重み付けすること(1114、1118)を更に含んでいる請求項記載の方法。
  9. 前記所望回数の正確な細分及び前記所望回数の概算的な細分は反復的に行われ、前記所望回数の概算的な細分中に前記反復して連続的に行なわれる前記シフの演算が互いに結合されることで計算コストを低減する請求項記載の方法。
  10. 前記所望回数の正確な細分及び前記所望回数の概算的な細分が空間において非等方的かつ非均一に行われることで、反復の或るレベルでの部分画像が同一形状又は同一サイズであることを必要としない請求項1記載の方法。
  11. 前記所望回数の正確な細分及び前記所望回数の概算的な細分がソース軌道の対称性に従って選択されることで、計算コストと画像精度との間で所望のトレードオフを達成する請求項10記載の方法。
  12. 前記方法が、複数のソース軌道から得られた複数のシノグラムに対して個々に適用されることで、最終画像を生成するのに結合される画像を生成する請求項1記載の方法。
  13. 前記所望回数の正確な細分は、部分画像用のシノグラム又は部分シノグラムにおける投影のトランケーション演算が、1以上の投影が該部分画像に利用可能になった後に開始されるようになされている請求項1記載の方法。
  14. 前記所望回数の正確な細分及び前記所望回数の概算的な細分は、部分画像用のシノグラムの細分が、1以上の投影が該部分画像に利用可能になった後に開始されるようになされている請求項1記載の方法。
  15. 複数の発散ビームシノグラムが利用可能になるにつれて、該複数の発散ビームシノグラムに対し時間的に連続して適用される請求項1記載の方法。
  16. 前記発散ビームシノグラムは収束ビームジオメトリ内で得られる請求項1記載の方法。
  17. オブジェクトの電子画像(36)を生成する装置であって、
    前記オブジェクトを走査して前記オブジェクトの画像を表すデータを生成する手段(12)と、
    前記データを操作して、処理されたシノグラム(22)を生成する手段(16、20)と、
    所望回数の正確な細分と所望回数の概算的な細分を行うことにより、前記シノグラム(22)を複数の部分シノグラム(1010、1012、1014、1016)に細分する手段(24)と、
    前記オブジェクト上に中心を置く全体座標系内で、前記部分シノグラムのそれぞれを逆投影して(1018、1020、1022、1024)、前記全体座標系内の正しい位置に中心を置く複数の対応する部分画像(図10中のf、f、f、f)を生成する手段(24)と、
    前記部分画像を集合させて前記電子画像を生成する手段(24)と、
    前記電子画像を格納、及び/又は表示、及び/又は解析する手段(38)と、
    を備え
    前記所望回数の正確な細分は、前記シノグラムにおける投影のトランケーションを行うこと、によってなされ、
    前記所望回数の概算的な細分は、前記シノグラムにおける投影のトランケーションを行い、該トランケーションの行われた投影を前記全体座標系の中心に対して、選択された距離だけシフトさせ、該シフトされた投影のデシメーションを、選択された係数によって行い、該デシメーションの行われた投影をその元の位置へ、前記選択された距離を用いてシフトさせること、によってなされる、装置。
  18. 前記シノグラムは、各部分シノグラムが所望サイズの部分画像に相当するまで、反復手法で複数の部分シノグラムに細分され請求項17記載の装置。
  19. 前記部分シノグラムは、サイズが1ピクセル又は1ボクセルと同程度に小さい部分画像に対応する請求項18記載の装置。
  20. 前記集合させる手段は反復手法で演算する請求項17記載の装置。
  21. 前記電子画像は断層撮影画像である請求項17記載の装置。
  22. 前記操作する手段は、データ座標内でオーバーサンプリングを行って(16)、前記電子画像の精度を向上させる請求項17記載の装置。
  23. 前記所望回数の概算的な細分は、前記シフトされた投影をデシメーションの前に重み付けすること、及び、前記デシメーションの行われた投影をシフトの前に重み付けすること(1114、1118)を更に含んでいる請求項17記載の装置。
  24. 電子画像から発散ビームシノグラム(34)を生成する方法であって、
    前記電子画像のトランケーションを行うことにより、前記電子画像を、前記電子画像上に中心を置く全体座標系内の正しい位置に中心を置く複数の部分画像(図17及び図9中のf、f、f、f)に分割するステップ(1710、1712、1714、1716)と、
    所望回数の正確な再投影と所望回数の概算的な再投影を行うことにより、前記全体座標系内で前記部分画像それぞれの部分シノグラム(g、g、g、g)を計算するステップ(1718、1720、1722、1724)と、
    前記部分シノグラムを集合させて前記シノグラム(34)を生成するステップ(1726)と、
    を備え
    前記所望回数の正確な再投影は、前記全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成すること、によって行われ、
    前記所望回数の概算的な再投影は、前記全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成し、該部分シノグラムにおける投影を、前記全体座標系の中心に対して、選択された距離だけシフトさせ、前記部分シノグラムにおける追加の投影を、補間を行うことにより生成し、前記部分シノグラムをその元の位置へ、前記選択された距離を用いてシフトさせること、によって行われる、方法。
  25. 前記分割は、前記部分画像のそれぞれが所望のサイズを有するまで、反復手法で行われ請求項24記載の方法。
  26. 前記集合は反復手法で行われる請求項24記載の方法。
  27. 前記所望回数の概算的な再投影、前記シフトされた投影を、補間を行う前に重み付けすること、及び、前記補間された投影をシフトの前に重み付けすること(1812、1816)を更に含んでいる請求項24記載の方法。
  28. オブジェクトの画像(38)を生成する装置であって、
    前記オブジェクトからデータを生成するスキャナ(12)と、
    前記画像の少なくとも1つの発散ビーム投影を生成するプロセッサ(16)と、
    少なくとも1つの投影から画像を再構成する、或いは反復再構成装置においては少なくとも1つの投影から画像を逆投影する、手段(24)と、
    前記再構成する手段によって生成された画像中のエラーを検出する手段(28)と、
    エラーの訂正後に画像を再投影してシノグラムを生成し、前記再構成する手段によって生成されたエラーを更に訂正するために前記シノグラムに更なる訂正を加え、該訂正されたシノグラムを前記再構成又は逆投影する手段(24)に与える(20)再投影手段(32)と、
    前記エラーが訂正された後に前記再構成する手段によって生成された画像を表示する手段(38)と、を備え、
    前記再投影手段(32)は、前記画像のトランケーションを行うことにより前記画像を複数の部分画像(図17中のf、f、f、f)に分割し、所望回数の正確な再投影及び所望回数の概算的な再投影を行うことにより前記オブジェクト上に中心を置く全体座標系内で前記部分画像の部分シノグラム(図17中のg、g、g、g)を計算し、該部分シノグラムを集合させて(1726)前記画像シノグラム(図17中のg)を生成し、
    前記所望回数の正確な再投影は、前記全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成すること、によって行われ、
    前記所望回数の概算的な再投影は、前記全体座標系内の部分画像を再投影して、対応する部分シノグラムを生成し、該部分シノグラムにおける投影を、前記全体座標系の中心に対して、選択された距離だけシフトさせ、前記部分シノグラムにおける追加の投影を、補間を行うことにより生成し、前記部分シノグラムをその元の位置へ、前記選択された距離を用いてシフトさせること、によって行われる、装置。
  29. 前記分割は、前記部分画像のそれぞれが所望のサイズを有するまで、反復手法で行われ請求項28記載の装置。
JP2003573443A 2002-02-28 2003-02-28 高速発散ビーム断層撮影のための方法及び装置 Expired - Lifetime JP4550426B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US36039002P 2002-02-28 2002-02-28
US10/190,295 US6771732B2 (en) 2002-02-28 2002-07-05 Methods and apparatus for fast divergent beam tomography
PCT/US2003/006040 WO2003075036A2 (en) 2002-02-28 2003-02-28 Methods and apparatus for divergent fast beam tomography

Publications (2)

Publication Number Publication Date
JP2005518892A JP2005518892A (ja) 2005-06-30
JP4550426B2 true JP4550426B2 (ja) 2010-09-22

Family

ID=27760117

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2003573443A Expired - Lifetime JP4550426B2 (ja) 2002-02-28 2003-02-28 高速発散ビーム断層撮影のための方法及び装置

Country Status (8)

Country Link
US (1) US6771732B2 (ja)
EP (1) EP1478274B1 (ja)
JP (1) JP4550426B2 (ja)
CN (1) CN100475146C (ja)
AT (1) ATE362152T1 (ja)
AU (1) AU2003230578A1 (ja)
DE (1) DE60313742T2 (ja)
WO (1) WO2003075036A2 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US12014557B2 (en) 2021-12-02 2024-06-18 V5Med Inc. High-speed automatic scanning system for interpreting images with AI assistance and method using the same

Families Citing this family (51)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6775264B1 (en) 1997-03-03 2004-08-10 Webley Systems, Inc. Computer, internet and telecommunications based network
US7516190B2 (en) 2000-02-04 2009-04-07 Parus Holdings, Inc. Personal voice-based information retrieval system
US6721705B2 (en) 2000-02-04 2004-04-13 Webley Systems, Inc. Robust voice browser system and voice activated device controller
JP2005519688A (ja) * 2002-03-13 2005-07-07 ブレークアウェイ・イメージング・エルエルシー 擬似同時多平面x線画像化システムおよび方法
CN1643371B (zh) 2002-03-19 2011-07-06 麦德特尼克航空公司 带有跟随数轴x射线源移动的探测器的计算机x光断层摄影装置
US20030190065A1 (en) * 2002-03-26 2003-10-09 Cti Pet Systems, Inc. Fast iterative image reconstruction from linograms
WO2003103496A1 (en) * 2002-06-11 2003-12-18 Breakaway Imaging, Llc Cantilevered gantry apparatus for x-ray imaging
WO2004019279A2 (en) * 2002-08-21 2004-03-04 Breakaway Imaging, Llc Apparatus and method for reconstruction of volumetric images in a divergent scanning computed tomography system
US7889835B2 (en) * 2003-08-07 2011-02-15 Morpho Detection, Inc. System and method for detecting an object by dynamically adjusting computational load
US7729526B2 (en) * 2003-09-09 2010-06-01 The Board Of Trustees Of The University Of Illinois Fast hierarchical tomography methods and apparatus
WO2005059840A1 (en) * 2003-12-16 2005-06-30 Philips Intellectual Property & Standards Gmbh Imaging method with filtered back projection
US7376255B2 (en) * 2004-06-23 2008-05-20 General Electric Company System and method for image reconstruction
US8774355B2 (en) * 2004-06-30 2014-07-08 General Electric Company Method and apparatus for direct reconstruction in tomosynthesis imaging
US7203267B2 (en) * 2004-06-30 2007-04-10 General Electric Company System and method for boundary estimation using CT metrology
US7215734B2 (en) * 2004-06-30 2007-05-08 General Electric Company Method and system for three-dimensional reconstruction of images
US7372937B2 (en) * 2004-07-16 2008-05-13 University Of Iowa Research Foundation Systems and methods of non-standard spiral cone-beam computed tomograpy (CT)
US7583777B2 (en) * 2004-07-21 2009-09-01 General Electric Company Method and apparatus for 3D reconstruction of images
US7362843B2 (en) * 2004-09-23 2008-04-22 General Electric Company System and method for reconstruction of cone beam tomographic projections with missing data
US20060072800A1 (en) * 2004-09-24 2006-04-06 General Electric Company Method and system for distributed iterative image reconstruction
US7385200B2 (en) * 2004-09-27 2008-06-10 Siemens Medical Solutions Usa, Inc. Re-binning method for nuclear medicine imaging devices
US7251306B2 (en) * 2004-11-17 2007-07-31 General Electric Company Methods, apparatus, and software to facilitate iterative reconstruction of images
US7840249B2 (en) * 2004-11-24 2010-11-23 University Of Iowa Research Foundation Clinical micro-CT (CMCT) methods, techniques and apparatus
US20060198491A1 (en) * 2005-03-04 2006-09-07 Kabushiki Kaisha Toshiba Volumetric computed tomography system for imaging
DE112006001584T5 (de) * 2005-06-16 2008-05-29 Ii-Vi Inc. Energie unterscheidendes Streuabbildungssystem
US7187750B1 (en) * 2005-09-20 2007-03-06 General Electric Company Method and apparatus for compensating non-uniform detector collimator plates
US7646842B2 (en) * 2005-09-23 2010-01-12 General Electric Company Methods and apparatus for reconstructing thick image slices
US7672421B2 (en) * 2005-10-12 2010-03-02 Siemens Medical Solutions Usa, Inc. Reduction of streak artifacts in low dose CT imaging through multi image compounding
US7215731B1 (en) * 2005-11-30 2007-05-08 General Electric Company Fast backprojection/reprojection with hexagonal segmentation of image
US7613275B2 (en) * 2005-12-19 2009-11-03 General Electric Company Method and apparatus for reducing cone beam artifacts using spatially varying weighting functions
US8120358B2 (en) * 2006-04-13 2012-02-21 The Regents Of The University Of California Magnetic resonance imaging with high spatial and temporal resolution
EP2011085A1 (en) * 2006-04-25 2009-01-07 Wisconsin Alumni Research Foundation System and method for estimating data missing from ct imaging projections
DE102007020879A1 (de) * 2006-05-10 2009-04-02 Gachon University Of Medicine & Science Industry-Academic Cooperation Foundation Verfahren und Vorrichtung für die äußerst schnelle Symmetrie- und SIMD- gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion
EP1860458A1 (en) * 2006-05-22 2007-11-28 Interuniversitair Microelektronica Centrum Detection of resonant tags by UWB radar
US20080075342A1 (en) * 2006-09-27 2008-03-27 Lazuka David M Pet scanner with digital trigger
US8270559B2 (en) * 2006-11-24 2012-09-18 Kabushiki Kaisha Toshiba Method and system for tomographic reconstruction in medical imaging using the circle and line trajectory
US8217937B2 (en) * 2007-03-28 2012-07-10 The Aerospace Corporation Isosurfacial three-dimensional imaging system and method
WO2008146186A2 (en) * 2007-05-30 2008-12-04 Koninklijke Philips Electronics, N.V. Pet local tomography
US8184887B2 (en) * 2008-08-29 2012-05-22 General Electric Company System and method for image reconstruction
US8755586B2 (en) * 2009-07-14 2014-06-17 Koninklijke Philips N.V. Image reconstruction including shift-variant blur compensation
US8315428B2 (en) * 2009-10-22 2012-11-20 Siemens Medical Solutions Usa, Inc. Automatic line identification and pairing for nuclear imaging collimator vector map characterization
US8340386B2 (en) * 2011-01-10 2012-12-25 Siemens Medical Solutions Usa, Inc. System and method for measuring hole orientation for SPECT collimators
GB201112359D0 (en) * 2011-07-19 2011-08-31 Univ Antwerpen Filter for tomographic reconstruction
US9224216B2 (en) 2013-07-31 2015-12-29 Kabushiki Kaisha Toshiba High density forward projector for spatial resolution improvement for medical imaging systems including computed tomography
US10653339B2 (en) 2014-04-29 2020-05-19 Nxp B.V. Time and frequency domain based activity tracking system
US9898840B2 (en) * 2014-05-15 2018-02-20 General Electric Company Systems and methods for continuous motion breast tomosynthesis
US9715744B2 (en) * 2014-12-11 2017-07-25 General Electric Company Method and device of obtaining beam hardening correction coefficient for carrying out beam hardening correction on computed tomography data
EP3136345A1 (en) * 2015-08-24 2017-03-01 FEI Company Positional error correction in a tomographic imaging apparatus
US10628973B2 (en) 2017-01-06 2020-04-21 General Electric Company Hierarchical tomographic reconstruction
CN111177497B (zh) * 2019-12-15 2023-11-24 腾讯科技(深圳)有限公司 层次数据的关联关系可视化处理方法、服务器及存储介质
CN113011121B (zh) * 2021-03-22 2022-04-22 华南理工大学 一种超高频开关变换器的变步长精细离散映射建模方法
CN114494952B (zh) * 2022-01-19 2024-04-23 杭州电子科技大学 一种基于感知损失的乳腺mri影像时间序列生成方法

Family Cites Families (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS60194939A (ja) * 1984-03-16 1985-10-03 横河メディカルシステム株式会社 計算機トモグラフイ装置
JPS628740A (ja) * 1985-07-04 1987-01-16 株式会社東芝 断層検査装置
JPS63243852A (ja) * 1987-03-31 1988-10-11 Toshiba Corp Ctスキヤナの画像再構成方式
JP2732595B2 (ja) * 1988-07-14 1998-03-30 株式会社東芝 バックプロジェクション処理用定数発生装置
JP3018194B2 (ja) * 1989-10-18 2000-03-13 ジーイー横河メディカルシステム株式会社 X線ctスキャナ
US5778038A (en) * 1996-06-06 1998-07-07 Yeda Research And Development Co., Ltd. Computerized tomography scanner and method of performing computerized tomography
SE9602594D0 (sv) * 1996-07-01 1996-07-01 Stefan Nilsson Förfarande och anordning vid datortomografi
US5805098A (en) * 1996-11-01 1998-09-08 The United States Of America As Represented By The Secretary Of The Army Method and system for forming image by backprojection
JP3672399B2 (ja) * 1996-11-21 2005-07-20 株式会社東芝 Ct画像再構成方法
US5848117A (en) * 1996-11-27 1998-12-08 Analogic Corporation Apparatus and method for computed tomography scanning using halfscan reconstruction with asymmetric detector system
US6351548B1 (en) * 1999-06-23 2002-02-26 The Board Of Trustees Of The University Of Illinois Fast hierarchical reprojection algorithm for tomography
US6282257B1 (en) * 1999-06-23 2001-08-28 The Board Of Trustees Of The University Of Illinois Fast hierarchical backprojection method for imaging
US6263096B1 (en) * 1999-06-23 2001-07-17 The Board Of Trustees Of The University Of Illinois Multilevel domain decomposition method for fast reprojection of images
US6332035B1 (en) * 1999-06-23 2001-12-18 The Board Of Trustees Of The University Of Illinois Fast hierarchical reprojection algorithms for 3D radon transforms
US6307911B1 (en) * 1999-06-23 2001-10-23 The Board Of Trustees Of The University Of Illinois Fast hierarchical backprojection for 3D Radon transform
US6287257B1 (en) * 1999-06-29 2001-09-11 Acuson Corporation Method and system for configuring a medical diagnostic ultrasound imaging system

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US12014557B2 (en) 2021-12-02 2024-06-18 V5Med Inc. High-speed automatic scanning system for interpreting images with AI assistance and method using the same

Also Published As

Publication number Publication date
EP1478274A4 (en) 2006-08-16
EP1478274B1 (en) 2007-05-09
CN1658795A (zh) 2005-08-24
AU2003230578A8 (en) 2003-09-16
EP1478274A2 (en) 2004-11-24
DE60313742D1 (de) 2007-06-21
DE60313742T2 (de) 2008-01-24
WO2003075036A3 (en) 2003-12-11
JP2005518892A (ja) 2005-06-30
AU2003230578A1 (en) 2003-09-16
US20030161443A1 (en) 2003-08-28
ATE362152T1 (de) 2007-06-15
CN100475146C (zh) 2009-04-08
WO2003075036A2 (en) 2003-09-12
US6771732B2 (en) 2004-08-03

Similar Documents

Publication Publication Date Title
JP4550426B2 (ja) 高速発散ビーム断層撮影のための方法及び装置
US7729526B2 (en) Fast hierarchical tomography methods and apparatus
US7215731B1 (en) Fast backprojection/reprojection with hexagonal segmentation of image
EP1953700B1 (en) System and method for reconstructing an image by rectilinear trajectory scanning
US6907102B1 (en) Iterative reconstruction methods for multi-slice computed tomography
JP4714316B2 (ja) イメージ再構成
JPH0910202A (ja) 診断撮像方法と装置およびヘリカル部分コーンビームデータからの画像再構成
EP1644897B1 (en) A fourier tomographic image reconstruction method for fan-beam data
WO1992005507A1 (en) Parallel processing method and apparatus based on the algebra reconstruction technique for reconstructing a three-dimensional computerized tomography
US5341460A (en) Method and apparatus for producing a three-dimensional computerized tomography image of an object with improved conversion of cone beam data to radon data
EP1749281A2 (en) Image reconstruction method for divergent beam scanner
JP3787374B2 (ja) 断層撮影画像内のアーチファクトを低減するための方法および装置
JP2001057976A (ja) 立体画像再構成方法及び装置並びにctスキャナー
US7602879B2 (en) Method for increasing the resolution of a CT image during image reconstruction
JP2000126173A (ja) 円錐ビ―ムデ―タのための画像再構築
JP4106251B2 (ja) 3次元逆投影方法およびx線ct装置
US7272205B2 (en) Methods, apparatus, and software to facilitate computing the elements of a forward projection matrix
US8379948B2 (en) Methods and systems for fast iterative reconstruction using separable system models
Grangeat et al. Indirect cone-beam three-dimensional image reconstruction
US6307911B1 (en) Fast hierarchical backprojection for 3D Radon transform
Sourbelle et al. Performance evaluation of local ROI algorithms for exact ROI reconstruction in spiral cone-beam computed tomography
JP2006527618A (ja) 余剰な測定値を使用するコンピュータ断層撮影法

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20060216

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20090127

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090423

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20090630

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090925

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20091124

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20100219

A602 Written permission of extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A602

Effective date: 20100226

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20100426

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20100622

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20100708

R150 Certificate of patent or registration of utility model

Ref document number: 4550426

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130716

Year of fee payment: 3

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

EXPY Cancellation because of completion of term