JP4555524B2 - イメージの高速な再投影のためのマルチレベル領域分解方法 - Google Patents
イメージの高速な再投影のためのマルチレベル領域分解方法 Download PDFInfo
- Publication number
- JP4555524B2 JP4555524B2 JP2001504966A JP2001504966A JP4555524B2 JP 4555524 B2 JP4555524 B2 JP 4555524B2 JP 2001504966 A JP2001504966 A JP 2001504966A JP 2001504966 A JP2001504966 A JP 2001504966A JP 4555524 B2 JP4555524 B2 JP 4555524B2
- Authority
- JP
- Japan
- Prior art keywords
- image
- sub
- sinogram
- reprojection
- images
- 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
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/001—Texturing; Colouring; Generation of texture or colour
Description
本発明は、ナショナル・サイエンス・ファウンデーションによって許可された、契約番号ＮＳＦ９５０２１３８及びＮＳＦ９１５７３７７の下に、政府の協力を得てなされた。政府は本発明に関し、一定の権利を有している。
【０００１】
これは、シリアル番号０９／３３８，０９２、出願日１９９９年６月２３日の一部継続出願である。
（技術分野）
本発明は、一般にイメージングに関する。特には、本発明は、断層面イメージの高速再投影の方法に関する。
【０００２】
（背景技術）
診断イメージは、細胞より小さい構成から環境への全てを観測する重要なツールであり、断層面再構成は、ほぼ全ての重要な診断イメージング法性の基板となるイメージング技術である。そのような断層面モードは、Ｘ線コンピュータトモグラフィ（ＣＴ）、陽子放出トモグラフィ（ＰＥＴ）及び単一光子放出数トモグラフィ（ＳＰＥＣＴ）を含んでいる。更に、磁気共鳴イメージング（ＭＲＩ）のためのある取得方法は、トモグラフィを使用し、電気インピーダンストモグラフィ（ＥＩＴ）及び光トモグラフィなどの新しい技術が出現している。断層面再構成は、電子顕微鏡で細胞より小さい構成を決定するためや、製造、地理学的研究、環境監視、及びリモートセンシングにおける非破壊評価（ＮＤＥ）を含む科学やエンジニアリングにおいて多くの他の応用の基礎ともなっている。
【０００３】
様々な断層面図の再構成処理において重要な要素は、前方投影あるいは再投影と呼ばれる、与えられたイメージの投影の集合を計算する処理である。再投影は、イメージから投影を作成し、投影がフィルタリングされ、逆投影されたときに、オリジナルのイメージを得るようにするための処理である。再投影は、様々な目的で広く用いられているが、再投影処理を利用する既知のアルゴリズムは、計算コストが高いという状況を起こしていた。
【０００４】
再投影は、幾つかの応用に於いて関心が持たれている。再投影は、Ｘ線ＣＴに用いられており、C.R.Crawford,J.G.Colsher,N.J.Pelec,and A.H.R.Lonn" High Speed Reprojection and its Application"、Proc. SPIE-Int.Soc.Opt.Eng.Conf.Medical Imaging II, Newport Beach, CA,Vol 914,pt A, pp.311-18, 1988を参照されたい。また、反復ビーム硬化修正アルゴリズムについては、P.M.Joseph and R.D.Spital,”A Method for Correcting Bone Induced Artifacts in Computed Tomography Scanners," JCAT vol 2 No.1, pp.100-108,1978 及びNaparsted への米国特許番号４，２１７，６４１及びCrawford等への５，４３８，６０２）を参照されたい。ストリーク（streak）抑圧アルゴリズムについては、G.Henrich," A Simple Computational Method for Reducing Streak Artifacts in CT Images, ”Computerized Tomography, vol.4, pp.67-71,1980 及びMattson 等への米国特許番号５，２２９，９３４を参照されたい。患者内に移植された金属の存在によって引き起こされるアーチファクトの除去のためのアルゴリズムについては、G.Glover and N.J.Pelc, "An Algorithm for the Reduction of Metal Clip Artifacts in CT Reconstructions."Medical Physics, vol.8 No.6, pp.799-807,1981及びTuy への米国特許番号５，２４３，６６４号あるいは、他の高密度の物体については、Crawfordへの米国特許番号４，７０９，３３３号、失われたデータの補正については、Hu等への米国特許番号５，３９６５２８号及び、部分体積効果については、Hsieh への米国特許番号５，７２７，０４１号を参照されたい。
【０００５】
ＰＥＴ及びＳＰＥＣＴイメージングにおいては、再投影は、減衰を補償するために使用されており、J.M.Ollinger, "Reconstruction-Reprojection Processing of Transmission Scans and the Variance of PET Images," IEEE Trans.Nuc.Sci. Vol.39 No.4, p.1122 Aug.1992 を参照されたい。再投影は、様々な取得エラーの検出と補償にも使用され、S.C.Huang and D.C.Yu,"Capability Evaluation of a Sinogram Error Detection and Correction Method in Computed Tomography," IEEE Trans. Nuc. Sci. Vol.39 No.4, pp.1106-1110 Aug. 1992及び、E.C.Frey, Z.W.Ju, and B.M.W.Tsui,"A Fast Projector-Backprojector Pair Modeling the Asymmetric, Spatially Varying Scatter Response Function for Scatter Compensation in SPECT Imaging," IEEE Trans. Nic. Sci, Vol.40, No.4,pp.1192-7,Aug.1993, 更に、患者の動きによる取得エラーについては、Nowak への米国特許番号４，８５８，１２８号、Arata への５，５５２，６０５、Lin への５，５７９，３５８、Kawai 等への５，８４８，１１４を参照されたい。ポアソンノイズ統計を考慮したものは、J.M.Ollinger," Iterative Reconstruction-Reprojection and Expectation-Maximization Algorithm," IEEE Trans. Med. Imag. Vol.9 No.1,p.94 Mar.1990 を参照されたい。再投影は、反復断層面再構成アルゴリズムの重要な要素であり、これらは、ＰＥＴ、ＳＰＥＣＴ、及び非破壊テストなどのイメージ法性における好適な方法である。反復断層面再構成アルゴリズムは、ＣＴ画像の中に金属片が存在することによって生じる様々な前述のアーチファクトを扱うためのＸ線ＣＴアルゴリズムにおいて重要である。
【０００６】
例えば、Ｘ線ＣＴを使ったマルチラインスパイラルＣＴ、カーディアックイメージングや、次のＣＴ蛍光顕微鏡法（fluoroscopy ）などのリアルタイムでデータを大量に集める能力を有する新技術の到来においては特に、そのような断増面イメージの高速な再構成が必要である。更に、例えば、手術を監視し、ガイドするために、リアルタイムの干渉（interventional）イメージングのために、このような要求が増えつつある。既知の技術においても非常にしばしば、そのようなデータからのイメージの再構成がボトルネックとなっている。ＣＴが導入されて以来、産業界や学界の継続的な努力が払われてきたが、今までの断増面再構成をスピードアップするための努力は成功していない。更に、断増面再構成のスピードアップの方法を開発することは、最近になって、National Institute of Health（NIH ）などのグループによって主要な推進分野として特定されている。
【０００７】
現在の反復法は、従来の直接フィルタド逆投影（ＦＢＰ）再構成よりもずっと高価であり、一般に、何オーダもの多くの計算を必要とする。これは、反復再構成を用いていた、１９７１年頃の最初のＣＴスキャナにまでさかのぼる反復法を高速化するために動機づけられた強力な仕事である。より最近になって、高速の反復再構成法を開発する必要性は、緊急の必要であると認識されており、Imaging Sciences Working Group,"Matching Clinical and Biological Needs with Emerging Imaging Technologies," tech.rep.,Diagnostic Imaging Program, National Cancer Institute,1998を参照されたい。様々な再構成アルゴリズムにおける再投影の処理の果たす重要な役割のために、これを高速化する重要性は、１９７８年には認識されていた。
【０００８】
再投影のための以前の方法は問題があった。直接法は、直接データを操作し、重み付け加算と補間の直接的な計算が含まれている。直接法は、計算量を増やせば幾らでも必要なだけ正確にすることができ、精度に関して基準となるものである。直接法の一般的な例は、P.M.Joseph,"An Improved Algorithm for Reprojecting Rays Through Pixel Images,"IEEE Trans. Med. Imag. Vol.1, No.3,pp.192-198 .1982, D.C.Yu and S.C.Huang," Study of Reprojection Methods in Terms of Their Resolution Loss and Sampling Errors," IEEE Trans. Nuc. Sci.,Vol.40,No.4,p1174,Aug.1993及びG.L.Zeng, Y.L.Hsieh, and G.T.Gullberg," A Rotating and Warping Projector/Backprojector for Fan-beam and Cone-beam Iterative Algorithm," IEEE Trans. Nuc. Sci., Vol.41,No.6,pp.2807-11,Dec.1994に記載されている。しかし、Zeng等への米国特許番号５，５５９，３３５号及び、Crandallへの５，６２５，１９０号などの様々な直接法の改良にも関わらず、直接法は、Ｎ×ＮピクセルのイメージのＮ個の投影を生成するのに、Ｎ3 に比例した（Ｏ（Ｎ3 ）と表す）高い計算コストが必要となる。例えば、一般的な４０９６×４０９６ピクセルのイメージの再投影を計算するのに、これらの方法の一つを使うと、２５６×２５６ピクセルのイメージに必要とされる計算の１６3 ＝４０９６倍が必要である。従って、様々な努力が学界や産業界に於いて、高速な再投影法の開発のために増えてきている。
【０００９】
製品における再投影処理を高速化する他の典型的なアプローチは、再投影を行う速度を大きくするための特別の目的を持ったハードウェアを開発することである。そのような方法は、E.B.Hinkle, J.L.C.Sanz, A.K.Jain, and D.Petkovic, "P/sup3/E:New Life for Projection-based Image Processing," J.Parallel&Distrib. Comput., Vol.4,No.1,pp.45-78, Feb.1987, T.M.Peters," Algorithms for Fast Back-and Re-Projection in Computed Tomography," IEEE Trans Nuc. Sci., Vol.NS-28, No.4,pp.3641-3647, Aug.1981, J.L.C.Sanz and E.B.Hinkle, "Computing Projections of Digital Images in Image Processing Pipeline Architectures, " IEEE Trans. Acoust., Speech&Signal Process., Vol.ASSP-35, No.2, pp.198-207, Feb.1987及びMeckley への米国特許番号４，９３０，０７６号、Brunnett等への５，００８，８２２号、Flickner等への５，１３６，６６０号、Jones 等への５，２２４，０３７号及び、Zeng等への５，５５９，３３５号に記載されている。C.R.Crawford,"Reprojection Using a Parallel Backprojector," Medical Physics, vol.13, No.4, pp.480-483, 1986 及びCrawford等への米国特許番号４，６２６，９９１号に記載された他のアプローチは、再構成システムに含まれることのあるハードウェアのバックプロジェクタを、バックプロジェクタが再投影を行えるように使用するものである。
【００１０】
特別の目的を持ったハードウェアを使ったそのような方法は、汎用コンピュータを使ったものよりも速い再構成を約束するが、特別な目的を持ったハードウェアは、高価であり、僅かなエキゾチックな実装でもほんの僅かなファクターの高速化しか得られない。通常、複数のプロセシングエレメント間で作業を分割する、これらの方法を実装するコストは、これらが提供するスピードアップに略比例している。更に、汎用コンピュータの性能の向上の驚くべき速度によれば、要求される特別な目的を持ったアーキテクチャは、すぐに時代遅れとなる。例えば、多くの現在のシステムは、もはやハードウェアのバックプロジェクタをもはやもっておらず、前述の方法は、これらのシステムでは時代遅れなものとしている。
【００１１】
T.Lu, et al., "Projection Iterative Reconstruction Technique and Optoelectronic Implementation." Proc. IEEE ISCAS '92, San Diego, May 10-13, pp.2469-2472, 1992 及びLu等への米国特許番号５，６５４，８２０号に記載された他の方法、アナログ光プロセシング方法は、他の機能の内、断層面再投影を実装するのに、アナログの光信号処理を使っている。この方法によれば、大きなスピードアップが得られようが、この方法は、他のハードウェア高速化スキームに対して、前に列挙した同じ欠点を持っている。更に、アナログ処理技術であるので、限定されたダイナミックレンジ、非線形性、不十分な精度、ロングタームドリフトなどの光アナログ信号処理の公知の限界を含んでいる。これらの限界は、大部分現在では、信号処理の現実的に全ての分野におけるデジタル実装に置き換わっているアナログ光信号処理の使用を制限する。
【００１２】
C.R.Crawford, J.G.Colsher, N.J.Pelc, and A. H. R. Lonn, "High Speed Reprojection and its Application," Proc. SPIE-Int. Soc. Opt. Eng. Conf. Medical Imaging II, Newport Beach, CA, vol 914, pt A, PP. 311-18,1988 及びCrawfordへの米国特許番号４，６１６，３１８号で提案されているような高速フーリエ変換（ＦＦＴ）を用いた方法などの他の方法は、理論的にはＯ（Ｎ2 ｌｏｇＮ）という許容可能な計算条件を示している。しかし、これらの方法は、フーリエ再構成アルゴリズム（ＦＲＡｓ）の基本である、フーリエスライスプロジェクション定理に基づいている。従って、これらは、現在３０年以上に渡って知られ、何回も再発見されてきているが、ＦＲＡｓが商業的に生き残ることを妨げている同様の問題を抱えている。
【００１３】
フーリエスライスプロジェクション定理は、角度θの投影の動径座標に沿ってのフーリエ変換は、θの方向に沿ったイメージの２次元フーリエ変換のスライスに等しいと述べている。このアルゴリズムの主要な問題は、ＦＦＴによってイメージの変換が計算されるフーリエ空間内の直交グリッドと、投影のフーリエ変換が評価されるべき極グリッド間の補間に必要とされるステップである。劣悪な理論的性質により、再構成において、シビアで、許容できないアーチファクトを顕在化させる。これらのアーチファクトを克服するために必要な多くの努力は、直接技術に対しての、僅かな、良くても少しの、例えば、２〜１０のファクタのスピードアップの結果しか生まない。
【００１４】
他の既知の方法、ＦＦＴを用いたリプロジェクタ−バックプロジェクタは、全ての場合ではないが、ある場合に有効な反復再構成アルゴリズムである結合再投影−逆投影（backprojection）のための処理である。この方法は、A.H. Delaney and Y. Bresler,"A Fast and Accurate Iterative Reconstruction Algorithm for Parallel-beam Tomography," IEEE Trans. Image Process., pp.740-753, May 1996 に記載されている。この方法は、ＦＲＡとは、異なった原理に基づいているが、やはりＦＦＴを使用している。これは、精度が良く、要求される計算量はＯ（Ｎ2 ｌｏｇＮ）である。しかし、このアルゴリズムは、平行ビーム幾何構成に本質的に限定されており、商業上の検査システムに一般に用いられるファンビームあるいは他の幾何構成に一般化することはできない。更に、この方法は、単独で再投影動作を行うことができず、従って、再投影アルゴリズムの重要な機能の多くを行うことができない。
【００１５】
高速離散ラドン変換として知られる更に他の方法は、再投影に似ているが、同一ではない処理を行う技術である。この方法は、M.L. Brady,"A Fast Discrete Approximation Algorithm for the Radon Transform," SIAM J. Comput., vol.27, No.1, pp.107-19, Feb.1998に記載されている。この方法は、離散ラドン変換を計算する。この方法は、イメージに対して、正方形のピクセル基底を仮定し、不均一に間隔の開いた角度に固定された集合に対して再投影を計算するのみである。この方法で計算される変換は、一般のトモグラフィにおける応用において考えられているラドン変換ではないので、問題である。むしろ、中心が所定の幅のストリップ内にあるピクセルの部分和について定義されている。それ故に、この方法は、２値イメージのハフ変換のような所定のイメージ処理アルゴリズムにおける応用が可能であり、トモグラフィにおける再構成には応用できない。
【００１６】
他の方法、Crawford等への米国特許番号４，７１４，９９７号において記載されている減縮データ処理は、再投影処理そのものをスピードアップすることは目的としていない。むしろ、このアプローチは、リプロジェクタの動作を改善するために加えられた前及び後処理サブシステムを使って、既存の再投影システムを使った処理時間を減縮するものである。特に、この方法は、減縮された空間分解能、すなわち、少ないピクセルのイメージについて動作し、たぶんに、オリジナルのデータ集合よりも少ないパスにおける再投影をするために、減縮されたサイズのイメージについて動作する。この方法は、再投影を利用する特定のアーチファクトを削減する応用に於いて十分な性能を示すが、精度について譲歩しないでこの方法が得るスピードアップは小さく、あるいは、最先端のシステムが要求する高精度に対してむしろ無視しても良いぐらいである。
【００１７】
図１を参照すると、既知の再投影技術において、スピードアップの必要性を示すために、辺の長さがＤの正方形領域が示されている。この正方形領域は、２次元イメージＦc （ｘ，ｙ）が正方形領域の外部で消えるように、すなわち、｜ｘ｜≦Ｄ／２、｜ｙ｜≦Ｄ／２であるように、定義される。イメージは、イメージのＮ×Ｎピクセル離散バージョンｆ（ｍ，ｎ）の補間によって十分な精度が得られるように近似される、すなわち、
【００１８】
【数１】
【００１９】
であり、ここで、ｐは、ピクセルの基底関数であり、Δ＝Ｄ／Ｎは、離散化サイズである。
ｘとｙの両方に於いて、ｐの台は、Δのサイズに比例する（Ｏ（Δ）と記載する）ように、Ｐの基本的帯域幅Ｂは、Ｏ（Δ-1）であると仮定する。ｐは空間的に限定されているので、ｐは、本質的に帯域制限されており、すなわち、ｐは、帯域の外では、無視できるエネルギーを有していることを意味している。例えば、サイズΔ正方ピクセル上のインジケータ関数のフーリエ変換は、各周波数軸に沿って、Δ-1のメインローブ半値幅を有しており、略帯域幅がＢ＝Δ-1である。同様に、以下に論じる数値実験において使用される３次元Ｂスプラインのテンソル積は、台が４Δ×４Δであり、帯域幅がΔ-1である。トモグラフィのための本質的に帯域が制限されたピクセル基底関数の別の選び方は、R.M.Lewitt,"Multidimentional digital image representations using generalized Kaiser-Bessel window functions," J. Opt. Soc. Am. A, vol 7, No. 10, pp. 1834-1846, Oct. 1990 に議論されている。実際は、標準のピクセル基底でしばしば充分である。というのも、トモグラフィにおける再構成の帯域幅は、ＦＢＰのフィルタリングあるいは畳み込みステップによって限定されるからである。
【００２０】
図１に、平行線投影の状況が記載されている。イメージのラドン変換
【００２１】
【数２】
【００２２】
は、角度の連続範囲、０≦θ≦πに渡って、イメージの平行線投影を生成する。投影のこの集合は、しばしばシノグラムと呼ばれる。
多くの現実的な状況に於いては、ｆc の投影は、時々制限された角度範囲｜θ｜≦Θの内部でのＰ個の角度｛θi ｝P i=1 とＯ（Ｎ）個の均一に隔てられた動径方向の点｛ｒk ｝の集合においてのみ測定される。与えられたイメージの要求される再投影動作は、全ての｛θi ｝P i=1 と｛ｒk ｝について、式（２）によって特定される計算に対応する。
【００２３】
【数３】
【００２４】
をピクセル基底関数の投影を示すとする。明らかに、これもサイズがＯ（Δ）の台を有している。式（１）の表示を使用すると、ある視角で、ある動径方向の位置における
【００２５】
【数４】
【００２６】
の計算は、Ｏ（Ｎ）の処理を必要とすることが分かる。従って、Ｐ個の視角とＮ個の動径方向の位置において投影を計算する処理数は、Ｏ（ＰＮ2 ）であり、Ｐ＝Ｏ（Ｎ）の典型的な場合には、計算コストは、Ｏ（Ｎ3 ）である。この好ましくないコストの大きさは、多くの応用に於いて、アルゴリズムの使用を制限し、反復イメージ再構成、例えば、データが収集された角度｛θi ｝に対してイメージを連続して近似する反復投影を要求する限定されたデータからの再構成において受け入れることができない。
【００２７】
以上を鑑み、本発明の目的は、既知の直接法の精度を得つつ、より少ない計算処理コスト、すなわち、Ｎ2 ｌｏｇＮあるいはそれ以上に良いスケーリングしか要求しない改良された方法を提供することである。典型的なイメージサイズにおいては、標準的なシリアル処理に対して、１００かそれ以上のファクターのスピードアップをする事ができる。
【００２８】
本発明の他の目的は、特別な目的のハードウェアを使うことができるが、要求はしないコストパフォーマンスの良い方法であって、標準のシリアルあるいはパラレルアーキテクチャに容易に実装できる方法を提供することである。
【００２９】
本発明の更に他の目的は、ＦＦＴを使った方法のような手法に起因したアーチファクトなしにスピードアップをすることができる方法を提供することである。
本発明の更なる他の目的は、柔軟で診察システムに一般に用いられるファンビームや他の幾何構成を適用できる方法を提供することである。
【００３０】
（発明の開示）
上記目的は、与えられた２次元イメージに対するラドン変換の投影あるいは角度サンプルの集合を計算することによって、再投影を行うために高速の計算方法を提供する本発明の方法によって達成できる。本発明の一実施形態に於いては、イメージのデータ表示を生成するために、再投影され、再結合された幾つかのサブイメージにイメージを分解する。イメージがサイズにおいて１ピクセルよりも大きい場合は、イメージをサブイメージに分解し、サブイメージは、サイズが１ピクセルになるまで、更に連続してより小さなサブイメージに分解される。サブイメージがサイズにおいて１ピクセルである場合には、これらは再投影される。従って、結果のイメージは、集められてトモグラフィの再投影イメージを生成し、表示される。収集は、全ての動径方向にシフトされた投影サブイメージが結合された後終わる。この処理は、反復処理、すなわち、有限数のステップでシステマチックに計算した処理値を許容する処理である。この方法の計算及びハードウェアのコストは、この方法がイメージを再構成するスピードに比べて低いものである。
【００３１】
略号テーブル
本特許は、幾つかの略号を使用する。以下のテーブルは、幾つかの略号の意味を読者が決定するための助けとして設けている。
【００３２】
ＣＰＵ＝中央処理ユニット
ＣＴ＝コンピュータトモグラフィ
ＥＩＴ＝電気インピーダンストモグラフィ
ＦＦＴ＝高速フーリエ変換
ＦＲＡ＝フーリエ再構成アルゴリズム
ＭＲＩ＝磁気共鳴イメージング
ＮＤＥ＝非破壊評価
ＮＩＨ＝National Institute of Health
ＰＥＴ＝陽子放出トモグラフィ
ｒｍｓ＝二乗平均平方根
ＳＰＥＣＴ＝単一光子放出数トモグラフィ
（発明の詳しい説明）
本発明は、コンピュータの処理能力及び装置のコストをかなり引き下げて、所定の角度範囲に渡って再投影を許容する処理である。この処理を理解するために、ラドン変換
【００３３】
【数５】
【００３４】
をその変数ｒとθの関数として見なしていくつかの知見を得る必要がある。
最初の２つの特徴は、ラドン変換の定義からすぐに得られるものである。特徴１は、ラドン変換は、線形変換であることである。
【００３５】
特徴２は、ラドン変換のシフト特性を定義するものである。ｆc,0 （ｘ、ｙ）＝ｆc （ｘ−ｘ0 、ｙ−ｙ0 ）が、ベクトル（ｘ0 、ｙ0 ）によってシフトされるイメージを示すとする。そのラドン変換は、外８ の動径方向にシフトされたものである。すなわち、
【００３６】
【外８】
【００３７】
【数６】
【００３８】
ここで、Ｓ（ｘ0 、ｙ0 ）は、元のものから、ベクトル（ｘ0 、ｙ0 ）によって移動された新しい原点に、投影を変換するシフトオペレータを示す。
特徴３は、投影の集合であるシノグラムの角度帯域制限である。第３の特徴は、「ボウネクタイ」性からくるもので、P.A. Rattey and A.G. Lindgren,"Sampling the 2-D Radon Transform," IEEE Trans. Acoust. Speech,Signal Proc., vol. 29, no. 5, pp.994-1002, 1981を参照されたい。これには、シノグラム 外９ の２次元フーリエ変換 外１０ の２つの変数ｒとθに対する本質的な台について書いている。
【００３９】
【外９】
【００４０】
【外１０】
【００４１】
なお、外１１ は、θについて周期的であるので、外１２は、ｕθ の直線スペクトルである。
【００４２】
【外１１】
【００４３】
【外１２】
【００４４】
あるいは、θに関するフーリエ級数を使うこともできる。本質的台外で 外１３ が０になる速度の厳密な評価は、F. Natterer,"The Mathematics of Computerized Tomography," Chichester, UK: Wiley, 1986 に与えられている。
【００４５】
【外１３】
【００４６】
外１４ のサポートが｜ｒ｜＜Ｒ（Ｒは、測定対象が収まるであろう最小円盤の半径である）で、外１５ が｜ｕr ｜＞Ｂr で実質的に０になるならば、
【００４７】
【外１４】
【００４８】
【外１５】
【００４９】
外１６ は、｜ｕθ ｜＞Ｒｂr で実質的に０になる。
【００５０】
【外１６】
【００５１】
従って、ｆc が｜ｘ｜≦Ｄ／２、｜ｙ｜≦Ｄ／２の領域で０でないならば、外１７ の台は、｜ｒ｜≦Ｄ／√２に限定される。
【００５２】
【外１７】
【００５３】
更に、Ｂ×Ｂの正方形に帯域制限された補間関数を使用するイメージｆc （ｘ、ｙ）の式（１）における近似は、ｆc が実質的にその変数について帯域制限されていることを示す。スライス−プロジェクション定理によれば、外１８ のｒについての帯域幅は、√２・Ｂで制限され、
【００５４】
【外１８】
【００５５】
特徴３から、外１９ のθに対する帯域幅は、√２・ＢＤ／√２＝ＤＢによって制限される。
【００５６】
【外１９】
【００５７】
前述の議論と、Ｂ＝Ｏ（Δ-1）であることに基づくと、特徴４は、重複のない投影の最小数は、Ｎ＝Ｄ／Δとして、イメージサイズＮに比例することである。これは、｜θ｜≦Θの範囲内の任意の値｛θi ｝に対して、正確に 外２０ を補間するために要求される、θの最低のサンプリングレートである。
【００５８】
【外２０】
【００５９】
本発明の一実施形態は、図２に見られるように、イメージの分解の単一レベルを行うものである。上記したように、既存の再投影は、ｃを定数として、ｃＮ3 の処理を必要とする。しかし、上記特徴によれば、本発明は、イメージの領域を小さく分割して、小さなサブイメージの集合とし、それらの投影を計算し、次に、結果を集めて、全体イメージの望まれる投影を得ることによって、処理数を減少している。この単一レベル領域分割処理では、イメージｆは、サイズが（Ｎ／Ｍ）×（Ｎ／Ｍ）のＭ2 個のサブイメージに細かく分解される。
【００６０】
図面、特に、図２を参照すると、正方形領域がＭ＝２の場合について示されている。各サブイメージに対しては、特徴４から、投影は、ｃ（Ｎ／Ｍ）3 のコストで、Ｏ（Ｎ／Ｍ）個の視角で、粗い角グリッド上の局所座標系で計算される。サブイメージの投影の集まりは、後に動径方向ｒにおけるシフトを伴う、角度座標θの補間を要求する。そして、線形性特徴１は、全体のイメージの投影を得るために、サブイメージの投影の加算に使用される。θ補間は、特徴４に適合するため、角度方向のサンプリングレートをＭ倍に増やす。というのも、イメージは全体として、サブイメージよりＭ倍大きいからである。θ補間による、Ｏ（Ｎ／Ｍ）個のサンプルをそれぞれ有するＯ（Ｎ）個の投影の計算コストは、ｃθ を定数として、ｃθ Ｎ2 ／Ｍで与えられる。上記特徴２によれば、動径方向のシフトは、サブイメージ座標系から全体イメージのそれへの平行移動を考慮する。この動径方向のシフトは、一般に、整数でない数のピクセル分である。しかし、それは、ｃr Ｎ2 ／Ｍ回の処理のコストを伴う補間となる。投影を計算する全コストは、
【００６１】
【数７】
【００６２】
となる。
上記の導出において、局所補間ができるように、θとｒの両方に適当なオーバサンプリングを適用することを仮定している。事実、そのようなオーバサンプリングによれば、補間誤差は、適切に選択された補間器の長さと共に、指数関数的に減少する。J.J. Knab."Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及び、O.M.Bucci, C. Gennarelli, and C. Savarese,"Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol 39, p. 1633, 1991 を参照されたい。Ｎが大きい場合には、Ｍ＝［ｃ／（ｃθ ＋ｃr ）］1/2 Ｎ1/2 という漸近的に最適な選択は、元のコストＯ（Ｎ3 ）と比較した場合、Ｏ（Ｎ2.5 ）のようにしか、Ｃ1 はスケーリングしないことを示す。精度と計算量は、投影収集ステップにおいて使用される特定の補間方法に依存する。
【００６３】
本発明の本実施形態は、図３と図４を比較することによってより良く理解される。既知の再投影処理では、完全なイメージ１２を、直接、投影の集合である単一のシノグラム１３に再投影している。図４においては、しかし、イメージ１２は、まず、独立にサブシノグラム２２、２４、２６及び２８にそれぞれ再投影される、４つのサブイメージ１４、１６、１８、及び２０に分割される。サブシノグラム２２、２４，２６、及び２８は、角度方向に補間され、動径方向にシフトされて、それぞれサブシノグラム３０、３２、３４、３６にそれぞれなる。サブシノグラムは、サブシノグラム２２、２４、２６、及び２８のθ’とｒ’のスケールをサブシノグラム３０、３２、３４及び３６のθ及びｒに効率的に調整することにより、補間とシフトを受け、サブシノグラム３０、３２、３４、及び３６は、直接単一のシノグラム１３を生成するように、足し合わされる。
【００６４】
他の実施形態においては、Ｏ（Ｎ2.5 ）の計算コストは、単一レベル領域分割方法ではなく、本発明のマルチレベル方法を採用することによって、更に減少される。本発明の説明を簡単にするために、イメージのサイズＮを２のべき乗とし、イメージｆは、図２に示されるように、同じ大きさの４つのサブイメージに細かく分割されるとする。イメージを小さなサブイメージに細かく分割する処理は、サブイメージが単一のピクセルの大きさのレベルになるまで、繰り返し、反復的に行うことができることが理解されよう。単一のピクセルレベルにおいては、投影の計算は、ピクセル基底補間関数 外２１ の投影は、解析的に知られていると仮定
【００６５】
【外２１】
【００６６】
しており、Ｏ（１）だけの投影が必要であるとされているので、自明なものとなる。従って、本発明の細分化の処理は、投影の計算をサブイメージ投影のマルチレベル集積に帰する。
【００６７】
集積処理は、座標系の原点、すなわち、回転の中心をサブイメージの中心から、複合イメージのそれへの平行移動を必要とする。これは、式（３）で定義されるシフトオペレータＳを用いることによって達成される。本発明のマルチレベル計算シーケンスは、反復処理として、最も簡単に定式化できる。以下の疑似コードは、イメージｆの投影 外２２ の計算を記述する。
【００６８】
【外２２】
【００６９】
外２３ ＝Ｆａｓｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆ）
【００７０】
【外２３】
【００７１】
ｉｆ ｓｉｚｅ（ｆ）＝１×１（１ピクセル）あるいは、ある最小のサイズ
外２４ ＝Ｄｉｒｅｃｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆ）
【００７２】
【外２４】
【００７３】
ｅｌｓｅ
外２５ ＝０
【００７４】
【外２５】
【００７５】
ｆｏｒ ｉ＝１ ｔｏ サブイメージの数
外２６ ＝Ｆａｓｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆS i ）
【００７６】
【外２６】
【００７７】
外２７＝ 外２８＋Ｓ[ｘi 、ｙi ］Ｉθ 外２９
【００７８】
【外２７】
【００７９】
【外２８】
【００８０】
【外２９】
【００８１】
ｅｎｄ ｆｏｒ
ｅｎｄ ｉｆ
ここで、Ｉθ は、θ補間オペレータであり、ｆS i は、イメージｆのｉ番目のサブイメージである。また、（ｘi 、ｙi ）は中心を、外３０ は、ｆS i の投影を示す。
【００８２】
【外３０】
【００８３】
オペレータＩθ は、補間によって、投影の数を２倍にする。反復の各レベルのコストは、Ｏ（Ｎ2 ）である。というのも、サブイメージの数は、各レベルで４倍になり、サブイメージのピクセル数に比例するサブイメージ毎の計算量は、４のファクターで減少する。従って、ｌｏｇＮレベルでは、マルチレベル領域分割処理の計算の複雑さは、Ｏ（Ｎ2 ｌｏｇＮ）である。
【００８４】
前節で述べた処理の数値的実装の前に、補間、基底関数の選択、及びオーバサンプリングの量に関する選択をしなくてはならない。これらのパラメータは、互いに関連し、処理の精度と複雑さに影響する。
【００８５】
反復処理におけるイベントの構成とシーケンスが、図５に示されている。図４のサブイメージ２０は、それ自身サブイメージ４０、４２、４４及び４６に細かく分割される。これらの細分イメージは、反復的に１６個のサブイメージに細かく分割され、各サブイメージが１ピクセルの大きさになるまで細分化処理は続けられる。その点に於いて、各ピクセルは、要求されるように補間され、シフトされるシノグラムに再投影され、サブイメージのシノグラムが合わせられると、シノグラム１３が生成される。
【００８６】
２つの補間ステップ、細かいレベルから粗いレベルに行く角度補間と動径方向のシフトを考慮する動径方向補間を、より詳しく述べる。先に述べたように、適切なオーバサンプリングによれば、補間器の長さに伴って補間エラーが漸近的に指数関数で減少すること約束する補間器が使用される。J.J. Knab."Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及びO.M. Bucci, C.Gennarelli, and C. Savarese, " Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol 39, p.1633, 1991を参照されたい。簡単のため、動径及び角度変数の両方について、線形補間のみを実装する。しかし、ちょうど良いイメージのサイズ及び近似的にのみ帯域制限された基底関数に対しては、漸近的振る舞いは得られず、角度及び動径補間器長、角度及び動径オーバサンプリングの量、全計算コスト、及び、エラーの性質のトレードオフを考えなければならない。
【００８７】
線形補間を正確にするために、基底関数は、十分滑らかでなくてはならず、オーバサンプリングも必要である。Crawfordへの米国特許番号４，７０９，３３３号にある、３次元Ｂスプラインのテンソル積からなる基底関数
【００８８】
【数８】
【００８９】
が、良い近似及び滑らかさの特性のために選択される。３次元Ｂスプラインは、単位パルスの４次畳み込み、特に、｜ｔ｜≦１／２に対して、ｂ0 ＝１で、その他は０とする場合、ｂ＝ｂ0 ＊ｂ0 ＊ｂ0 ＊ｂ0 で定義されることを思い出されたい。
【００９０】
ここで、図６を参照すると、既知のシェップ−ローガンファントムが示されている。L.A. Shepp and B.F. Logan," The Fourier Reconstruction of a Head Section," IEEE Trans. Nucl. Sci., vol. 21, no. 1, pp. 21-43, 1974を参照されたい。シェップ−ローガンファントムは、本発明の有効性を示すために選択された。既知の直接投影と本発明の処理を使った、Ｎ＝２５６の場合に計算されたシェップ−ローガンファントムの投影が、図７ａ及び図７ｂのそれぞれに等高線プロットとして示されている。図７ａと７ｂを見て比較して明らかなように、本発明の方法は、正確な投影を実現している。更に、数値実験は、再投影の精度は、角度及び動径方向のオーバサンプリングの量に依存していることが示されている。図７ｂは、動径方向及び角度方向の４倍のオーバサンプリングの結果を示している。
【００９１】
ここで、図８ａ及び８ｂを参照すると、図７ａ及び７ｂに示されるシノグラムの２つの断面、特に、外３１ 及び 外３２ が示されている。
【００９２】
【外３１】
【００９３】
【外３２】
【００９４】
実線は、本発明の高速処理を示しており、ダイヤ印は、既知の直接再投影を示している。本発明のマルチレベル処理によって得られる結果は、視覚的に、既知の直接再投影に一致している。従って、本発明の精度が数値的に確かめられた。
【００９５】
再投影の精度を測定するために、規格化二乗平均平方根（ｒｍｓ）エラーを定義する。
【００９６】
【数９】
【００９７】
ここで、外３３ は、本発明の高速処理及び既知の直接再投影によって、計算され
【００９８】
【外３３】
【００９９】
た投影 外３４ を示している。
【０１００】
【外３４】
【０１０１】
ここで、インデックスｉ及びｊは、それぞれサンプルされた角度及び動径方向に対応する。
図９を参照すると、Δrms 対ｒのオーバサンプリング量のプロットが、θについての４倍角度オーバサンプリングについて示されている。最適な補間方法を使用することによって、よりよい性能のトレードオフが、少ないオーバサンプリングファクタで可能である。J.J. Knab, " Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及び、O.M. Bucci, C.Gennarelli, and C. Savarese, " Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol. 39, p.1633, 1991 を参照されたい。
【０１０２】
本発明の計算上の利点は、様々なサイズのイメージを考慮することによって、確かめられる。図１０において、ＣＰＵ時間対イメージサイズＮが既知の直接再投影と本発明の高速処理による投影の計算について示されている。直接再投影のカーブの傾きは、予測されるＯ（Ｎ3 ）の振る舞いに対応している。また、高速処理のカーブの観測される傾向は、本発明の処理の予想されるＯ（Ｎ2 ｌｏｇＮ）の複雑さであることを支持している。従って、このプロットから、本発明は、既知の直接再投影法よりも高速な結果を達成することが明らかである。
【０１０３】
本発明の原理に従って作られたイメージング装置５０が、図１１に示されている。イメージング装置５０は、ＣＴスキャナあるいは、様々な他のイメージング装置である。イメージング装置５０は、頭などの対象物からの生のデータを生成するスキャナ５２を含んでいる。データは、レシーバ５４に送られ、それから、後処理装置（あるいはステップ）５６に送られる。リビニング（re-binning）などの処理は、後処理装置５６で行うことが可能である。後処理装置５６の出力は、装置（あるいはステップ）５８において再構成され、結果のイメージは、ディスプレイ装置６０に表示される。しかし、イメージがアーチファクト（例えば、頭の中の金属片などによって、）を有する場合は、先に述べたように、イメージをエラー修正装置（あるいは、ステップ）６４と再投影装置（あるいは、ステップ）６２に供給することによって、結果としてのエラーを取り除くことが可能である。再投影後のシノグラム出力は、再構成装置（あるいは、ステップ）５８の入力に供給される。再投影とエラー修正は、アーチファクトによって起きたエラーが修正されるまで、繰り返される。
【０１０４】
前述の記載から、イメージの再投影あるいはトモグラフィの投影を計算する改良された方法は、多くの望ましい性質及び利点を有することが示され記載されたことが理解される。これは、既存の技術と同じ品質の投影を１００倍速く生成することができる。本発明は、キーとなる機能、従って、トモグラフィのスキャナ、例えば、ＣＴ、ＰＥＴ、ＳＰＥＣＴなどの全体の再構成を、高価な特別の目的を持ったコンピュータハードウェアを要求することなく、スピードアップすることができる。にもかかわらず、これは、同じ品質の投影を生成する特別目的のハードウェアに実装することが可能であるが、同様のコストの特別目的ハードウェアシステムよりずっと速く動作する。
【０１０５】
本発明の方法の特別な実施形態について示し、記載したが、当業者によれば、発明の最も広い側面及び請求項に記載されるものから離れることなく、これに変更などを施すことができることが理解されるであろう。
【０１０６】
本発明のこれら及び他の特徴や利点は、詳細な説明及び図面を参照すれば、当業者によれば明確であろう。
【図面の簡単な説明】
【図１】 既知の平行線投影。
【図２】 本発明の例示的投影。
【図３】 既知の再投影処理の図である。
【図４】 本発明の原理に基づいて行われる階層的分解と再投影の図である。
【図５】 図４のサブイメージに行われる反復的サブイメージ分解の図である。
【図６】 本発明の動作の図示に使用されたシェップ−ローガンファントム（Shepp-Logan phantom ）の幾何構成の表示である。
【図７】 図７ａ及び図７ｂは本発明の高速処理と既知の直接再投影それぞれによって計算された図６のシェップ−ローガンファントムのシノグラム投影である。
【図８】 図８ａ及び図８ｂは、
【数１０】
それぞれに対する図７ａと図７ｂのシノグラムの断面図である。
【図９】 図６のシェップ−ローガンファントムに対する、二乗平均平方根（ｒｍｓ）再投影エラーΔrms 対動径方向のオーバサンプリングの量のプロットである。
【図１０】 既知の直接再投影と本発明の方法による投影を計算する際の中央処理ユニット（ＣＰＵ）時間対イメージサイズＮを比較する図である。
【図１１】 本発明の原理に基づいたイメージ装置のブロック図である。
【０００１】
これは、シリアル番号０９／３３８，０９２、出願日１９９９年６月２３日の一部継続出願である。
（技術分野）
本発明は、一般にイメージングに関する。特には、本発明は、断層面イメージの高速再投影の方法に関する。
【０００２】
（背景技術）
診断イメージは、細胞より小さい構成から環境への全てを観測する重要なツールであり、断層面再構成は、ほぼ全ての重要な診断イメージング法性の基板となるイメージング技術である。そのような断層面モードは、Ｘ線コンピュータトモグラフィ（ＣＴ）、陽子放出トモグラフィ（ＰＥＴ）及び単一光子放出数トモグラフィ（ＳＰＥＣＴ）を含んでいる。更に、磁気共鳴イメージング（ＭＲＩ）のためのある取得方法は、トモグラフィを使用し、電気インピーダンストモグラフィ（ＥＩＴ）及び光トモグラフィなどの新しい技術が出現している。断層面再構成は、電子顕微鏡で細胞より小さい構成を決定するためや、製造、地理学的研究、環境監視、及びリモートセンシングにおける非破壊評価（ＮＤＥ）を含む科学やエンジニアリングにおいて多くの他の応用の基礎ともなっている。
【０００３】
様々な断層面図の再構成処理において重要な要素は、前方投影あるいは再投影と呼ばれる、与えられたイメージの投影の集合を計算する処理である。再投影は、イメージから投影を作成し、投影がフィルタリングされ、逆投影されたときに、オリジナルのイメージを得るようにするための処理である。再投影は、様々な目的で広く用いられているが、再投影処理を利用する既知のアルゴリズムは、計算コストが高いという状況を起こしていた。
【０００４】
再投影は、幾つかの応用に於いて関心が持たれている。再投影は、Ｘ線ＣＴに用いられており、C.R.Crawford,J.G.Colsher,N.J.Pelec,and A.H.R.Lonn" High Speed Reprojection and its Application"、Proc. SPIE-Int.Soc.Opt.Eng.Conf.Medical Imaging II, Newport Beach, CA,Vol 914,pt A, pp.311-18, 1988を参照されたい。また、反復ビーム硬化修正アルゴリズムについては、P.M.Joseph and R.D.Spital,”A Method for Correcting Bone Induced Artifacts in Computed Tomography Scanners," JCAT vol 2 No.1, pp.100-108,1978 及びNaparsted への米国特許番号４，２１７，６４１及びCrawford等への５，４３８，６０２）を参照されたい。ストリーク（streak）抑圧アルゴリズムについては、G.Henrich," A Simple Computational Method for Reducing Streak Artifacts in CT Images, ”Computerized Tomography, vol.4, pp.67-71,1980 及びMattson 等への米国特許番号５，２２９，９３４を参照されたい。患者内に移植された金属の存在によって引き起こされるアーチファクトの除去のためのアルゴリズムについては、G.Glover and N.J.Pelc, "An Algorithm for the Reduction of Metal Clip Artifacts in CT Reconstructions."Medical Physics, vol.8 No.6, pp.799-807,1981及びTuy への米国特許番号５，２４３，６６４号あるいは、他の高密度の物体については、Crawfordへの米国特許番号４，７０９，３３３号、失われたデータの補正については、Hu等への米国特許番号５，３９６５２８号及び、部分体積効果については、Hsieh への米国特許番号５，７２７，０４１号を参照されたい。
【０００５】
ＰＥＴ及びＳＰＥＣＴイメージングにおいては、再投影は、減衰を補償するために使用されており、J.M.Ollinger, "Reconstruction-Reprojection Processing of Transmission Scans and the Variance of PET Images," IEEE Trans.Nuc.Sci. Vol.39 No.4, p.1122 Aug.1992 を参照されたい。再投影は、様々な取得エラーの検出と補償にも使用され、S.C.Huang and D.C.Yu,"Capability Evaluation of a Sinogram Error Detection and Correction Method in Computed Tomography," IEEE Trans. Nuc. Sci. Vol.39 No.4, pp.1106-1110 Aug. 1992及び、E.C.Frey, Z.W.Ju, and B.M.W.Tsui,"A Fast Projector-Backprojector Pair Modeling the Asymmetric, Spatially Varying Scatter Response Function for Scatter Compensation in SPECT Imaging," IEEE Trans. Nic. Sci, Vol.40, No.4,pp.1192-7,Aug.1993, 更に、患者の動きによる取得エラーについては、Nowak への米国特許番号４，８５８，１２８号、Arata への５，５５２，６０５、Lin への５，５７９，３５８、Kawai 等への５，８４８，１１４を参照されたい。ポアソンノイズ統計を考慮したものは、J.M.Ollinger," Iterative Reconstruction-Reprojection and Expectation-Maximization Algorithm," IEEE Trans. Med. Imag. Vol.9 No.1,p.94 Mar.1990 を参照されたい。再投影は、反復断層面再構成アルゴリズムの重要な要素であり、これらは、ＰＥＴ、ＳＰＥＣＴ、及び非破壊テストなどのイメージ法性における好適な方法である。反復断層面再構成アルゴリズムは、ＣＴ画像の中に金属片が存在することによって生じる様々な前述のアーチファクトを扱うためのＸ線ＣＴアルゴリズムにおいて重要である。
【０００６】
例えば、Ｘ線ＣＴを使ったマルチラインスパイラルＣＴ、カーディアックイメージングや、次のＣＴ蛍光顕微鏡法（fluoroscopy ）などのリアルタイムでデータを大量に集める能力を有する新技術の到来においては特に、そのような断増面イメージの高速な再構成が必要である。更に、例えば、手術を監視し、ガイドするために、リアルタイムの干渉（interventional）イメージングのために、このような要求が増えつつある。既知の技術においても非常にしばしば、そのようなデータからのイメージの再構成がボトルネックとなっている。ＣＴが導入されて以来、産業界や学界の継続的な努力が払われてきたが、今までの断増面再構成をスピードアップするための努力は成功していない。更に、断増面再構成のスピードアップの方法を開発することは、最近になって、National Institute of Health（NIH ）などのグループによって主要な推進分野として特定されている。
【０００７】
現在の反復法は、従来の直接フィルタド逆投影（ＦＢＰ）再構成よりもずっと高価であり、一般に、何オーダもの多くの計算を必要とする。これは、反復再構成を用いていた、１９７１年頃の最初のＣＴスキャナにまでさかのぼる反復法を高速化するために動機づけられた強力な仕事である。より最近になって、高速の反復再構成法を開発する必要性は、緊急の必要であると認識されており、Imaging Sciences Working Group,"Matching Clinical and Biological Needs with Emerging Imaging Technologies," tech.rep.,Diagnostic Imaging Program, National Cancer Institute,1998を参照されたい。様々な再構成アルゴリズムにおける再投影の処理の果たす重要な役割のために、これを高速化する重要性は、１９７８年には認識されていた。
【０００８】
再投影のための以前の方法は問題があった。直接法は、直接データを操作し、重み付け加算と補間の直接的な計算が含まれている。直接法は、計算量を増やせば幾らでも必要なだけ正確にすることができ、精度に関して基準となるものである。直接法の一般的な例は、P.M.Joseph,"An Improved Algorithm for Reprojecting Rays Through Pixel Images,"IEEE Trans. Med. Imag. Vol.1, No.3,pp.192-198 .1982, D.C.Yu and S.C.Huang," Study of Reprojection Methods in Terms of Their Resolution Loss and Sampling Errors," IEEE Trans. Nuc. Sci.,Vol.40,No.4,p1174,Aug.1993及びG.L.Zeng, Y.L.Hsieh, and G.T.Gullberg," A Rotating and Warping Projector/Backprojector for Fan-beam and Cone-beam Iterative Algorithm," IEEE Trans. Nuc. Sci., Vol.41,No.6,pp.2807-11,Dec.1994に記載されている。しかし、Zeng等への米国特許番号５，５５９，３３５号及び、Crandallへの５，６２５，１９０号などの様々な直接法の改良にも関わらず、直接法は、Ｎ×ＮピクセルのイメージのＮ個の投影を生成するのに、Ｎ3 に比例した（Ｏ（Ｎ3 ）と表す）高い計算コストが必要となる。例えば、一般的な４０９６×４０９６ピクセルのイメージの再投影を計算するのに、これらの方法の一つを使うと、２５６×２５６ピクセルのイメージに必要とされる計算の１６3 ＝４０９６倍が必要である。従って、様々な努力が学界や産業界に於いて、高速な再投影法の開発のために増えてきている。
【０００９】
製品における再投影処理を高速化する他の典型的なアプローチは、再投影を行う速度を大きくするための特別の目的を持ったハードウェアを開発することである。そのような方法は、E.B.Hinkle, J.L.C.Sanz, A.K.Jain, and D.Petkovic, "P/sup3/E:New Life for Projection-based Image Processing," J.Parallel&Distrib. Comput., Vol.4,No.1,pp.45-78, Feb.1987, T.M.Peters," Algorithms for Fast Back-and Re-Projection in Computed Tomography," IEEE Trans Nuc. Sci., Vol.NS-28, No.4,pp.3641-3647, Aug.1981, J.L.C.Sanz and E.B.Hinkle, "Computing Projections of Digital Images in Image Processing Pipeline Architectures, " IEEE Trans. Acoust., Speech&Signal Process., Vol.ASSP-35, No.2, pp.198-207, Feb.1987及びMeckley への米国特許番号４，９３０，０７６号、Brunnett等への５，００８，８２２号、Flickner等への５，１３６，６６０号、Jones 等への５，２２４，０３７号及び、Zeng等への５，５５９，３３５号に記載されている。C.R.Crawford,"Reprojection Using a Parallel Backprojector," Medical Physics, vol.13, No.4, pp.480-483, 1986 及びCrawford等への米国特許番号４，６２６，９９１号に記載された他のアプローチは、再構成システムに含まれることのあるハードウェアのバックプロジェクタを、バックプロジェクタが再投影を行えるように使用するものである。
【００１０】
特別の目的を持ったハードウェアを使ったそのような方法は、汎用コンピュータを使ったものよりも速い再構成を約束するが、特別な目的を持ったハードウェアは、高価であり、僅かなエキゾチックな実装でもほんの僅かなファクターの高速化しか得られない。通常、複数のプロセシングエレメント間で作業を分割する、これらの方法を実装するコストは、これらが提供するスピードアップに略比例している。更に、汎用コンピュータの性能の向上の驚くべき速度によれば、要求される特別な目的を持ったアーキテクチャは、すぐに時代遅れとなる。例えば、多くの現在のシステムは、もはやハードウェアのバックプロジェクタをもはやもっておらず、前述の方法は、これらのシステムでは時代遅れなものとしている。
【００１１】
T.Lu, et al., "Projection Iterative Reconstruction Technique and Optoelectronic Implementation." Proc. IEEE ISCAS '92, San Diego, May 10-13, pp.2469-2472, 1992 及びLu等への米国特許番号５，６５４，８２０号に記載された他の方法、アナログ光プロセシング方法は、他の機能の内、断層面再投影を実装するのに、アナログの光信号処理を使っている。この方法によれば、大きなスピードアップが得られようが、この方法は、他のハードウェア高速化スキームに対して、前に列挙した同じ欠点を持っている。更に、アナログ処理技術であるので、限定されたダイナミックレンジ、非線形性、不十分な精度、ロングタームドリフトなどの光アナログ信号処理の公知の限界を含んでいる。これらの限界は、大部分現在では、信号処理の現実的に全ての分野におけるデジタル実装に置き換わっているアナログ光信号処理の使用を制限する。
【００１２】
C.R.Crawford, J.G.Colsher, N.J.Pelc, and A. H. R. Lonn, "High Speed Reprojection and its Application," Proc. SPIE-Int. Soc. Opt. Eng. Conf. Medical Imaging II, Newport Beach, CA, vol 914, pt A, PP. 311-18,1988 及びCrawfordへの米国特許番号４，６１６，３１８号で提案されているような高速フーリエ変換（ＦＦＴ）を用いた方法などの他の方法は、理論的にはＯ（Ｎ2 ｌｏｇＮ）という許容可能な計算条件を示している。しかし、これらの方法は、フーリエ再構成アルゴリズム（ＦＲＡｓ）の基本である、フーリエスライスプロジェクション定理に基づいている。従って、これらは、現在３０年以上に渡って知られ、何回も再発見されてきているが、ＦＲＡｓが商業的に生き残ることを妨げている同様の問題を抱えている。
【００１３】
フーリエスライスプロジェクション定理は、角度θの投影の動径座標に沿ってのフーリエ変換は、θの方向に沿ったイメージの２次元フーリエ変換のスライスに等しいと述べている。このアルゴリズムの主要な問題は、ＦＦＴによってイメージの変換が計算されるフーリエ空間内の直交グリッドと、投影のフーリエ変換が評価されるべき極グリッド間の補間に必要とされるステップである。劣悪な理論的性質により、再構成において、シビアで、許容できないアーチファクトを顕在化させる。これらのアーチファクトを克服するために必要な多くの努力は、直接技術に対しての、僅かな、良くても少しの、例えば、２〜１０のファクタのスピードアップの結果しか生まない。
【００１４】
他の既知の方法、ＦＦＴを用いたリプロジェクタ−バックプロジェクタは、全ての場合ではないが、ある場合に有効な反復再構成アルゴリズムである結合再投影−逆投影（backprojection）のための処理である。この方法は、A.H. Delaney and Y. Bresler,"A Fast and Accurate Iterative Reconstruction Algorithm for Parallel-beam Tomography," IEEE Trans. Image Process., pp.740-753, May 1996 に記載されている。この方法は、ＦＲＡとは、異なった原理に基づいているが、やはりＦＦＴを使用している。これは、精度が良く、要求される計算量はＯ（Ｎ2 ｌｏｇＮ）である。しかし、このアルゴリズムは、平行ビーム幾何構成に本質的に限定されており、商業上の検査システムに一般に用いられるファンビームあるいは他の幾何構成に一般化することはできない。更に、この方法は、単独で再投影動作を行うことができず、従って、再投影アルゴリズムの重要な機能の多くを行うことができない。
【００１５】
高速離散ラドン変換として知られる更に他の方法は、再投影に似ているが、同一ではない処理を行う技術である。この方法は、M.L. Brady,"A Fast Discrete Approximation Algorithm for the Radon Transform," SIAM J. Comput., vol.27, No.1, pp.107-19, Feb.1998に記載されている。この方法は、離散ラドン変換を計算する。この方法は、イメージに対して、正方形のピクセル基底を仮定し、不均一に間隔の開いた角度に固定された集合に対して再投影を計算するのみである。この方法で計算される変換は、一般のトモグラフィにおける応用において考えられているラドン変換ではないので、問題である。むしろ、中心が所定の幅のストリップ内にあるピクセルの部分和について定義されている。それ故に、この方法は、２値イメージのハフ変換のような所定のイメージ処理アルゴリズムにおける応用が可能であり、トモグラフィにおける再構成には応用できない。
【００１６】
他の方法、Crawford等への米国特許番号４，７１４，９９７号において記載されている減縮データ処理は、再投影処理そのものをスピードアップすることは目的としていない。むしろ、このアプローチは、リプロジェクタの動作を改善するために加えられた前及び後処理サブシステムを使って、既存の再投影システムを使った処理時間を減縮するものである。特に、この方法は、減縮された空間分解能、すなわち、少ないピクセルのイメージについて動作し、たぶんに、オリジナルのデータ集合よりも少ないパスにおける再投影をするために、減縮されたサイズのイメージについて動作する。この方法は、再投影を利用する特定のアーチファクトを削減する応用に於いて十分な性能を示すが、精度について譲歩しないでこの方法が得るスピードアップは小さく、あるいは、最先端のシステムが要求する高精度に対してむしろ無視しても良いぐらいである。
【００１７】
図１を参照すると、既知の再投影技術において、スピードアップの必要性を示すために、辺の長さがＤの正方形領域が示されている。この正方形領域は、２次元イメージＦc （ｘ，ｙ）が正方形領域の外部で消えるように、すなわち、｜ｘ｜≦Ｄ／２、｜ｙ｜≦Ｄ／２であるように、定義される。イメージは、イメージのＮ×Ｎピクセル離散バージョンｆ（ｍ，ｎ）の補間によって十分な精度が得られるように近似される、すなわち、
【００１８】
【数１】
であり、ここで、ｐは、ピクセルの基底関数であり、Δ＝Ｄ／Ｎは、離散化サイズである。
ｘとｙの両方に於いて、ｐの台は、Δのサイズに比例する（Ｏ（Δ）と記載する）ように、Ｐの基本的帯域幅Ｂは、Ｏ（Δ-1）であると仮定する。ｐは空間的に限定されているので、ｐは、本質的に帯域制限されており、すなわち、ｐは、帯域の外では、無視できるエネルギーを有していることを意味している。例えば、サイズΔ正方ピクセル上のインジケータ関数のフーリエ変換は、各周波数軸に沿って、Δ-1のメインローブ半値幅を有しており、略帯域幅がＢ＝Δ-1である。同様に、以下に論じる数値実験において使用される３次元Ｂスプラインのテンソル積は、台が４Δ×４Δであり、帯域幅がΔ-1である。トモグラフィのための本質的に帯域が制限されたピクセル基底関数の別の選び方は、R.M.Lewitt,"Multidimentional digital image representations using generalized Kaiser-Bessel window functions," J. Opt. Soc. Am. A, vol 7, No. 10, pp. 1834-1846, Oct. 1990 に議論されている。実際は、標準のピクセル基底でしばしば充分である。というのも、トモグラフィにおける再構成の帯域幅は、ＦＢＰのフィルタリングあるいは畳み込みステップによって限定されるからである。
【００２０】
図１に、平行線投影の状況が記載されている。イメージのラドン変換
【００２１】
【数２】
は、角度の連続範囲、０≦θ≦πに渡って、イメージの平行線投影を生成する。投影のこの集合は、しばしばシノグラムと呼ばれる。
多くの現実的な状況に於いては、ｆc の投影は、時々制限された角度範囲｜θ｜≦Θの内部でのＰ個の角度｛θi ｝P i=1 とＯ（Ｎ）個の均一に隔てられた動径方向の点｛ｒk ｝の集合においてのみ測定される。与えられたイメージの要求される再投影動作は、全ての｛θi ｝P i=1 と｛ｒk ｝について、式（２）によって特定される計算に対応する。
【００２３】
【数３】
をピクセル基底関数の投影を示すとする。明らかに、これもサイズがＯ（Δ）の台を有している。式（１）の表示を使用すると、ある視角で、ある動径方向の位置における
【００２５】
【数４】
の計算は、Ｏ（Ｎ）の処理を必要とすることが分かる。従って、Ｐ個の視角とＮ個の動径方向の位置において投影を計算する処理数は、Ｏ（ＰＮ2 ）であり、Ｐ＝Ｏ（Ｎ）の典型的な場合には、計算コストは、Ｏ（Ｎ3 ）である。この好ましくないコストの大きさは、多くの応用に於いて、アルゴリズムの使用を制限し、反復イメージ再構成、例えば、データが収集された角度｛θi ｝に対してイメージを連続して近似する反復投影を要求する限定されたデータからの再構成において受け入れることができない。
【００２７】
以上を鑑み、本発明の目的は、既知の直接法の精度を得つつ、より少ない計算処理コスト、すなわち、Ｎ2 ｌｏｇＮあるいはそれ以上に良いスケーリングしか要求しない改良された方法を提供することである。典型的なイメージサイズにおいては、標準的なシリアル処理に対して、１００かそれ以上のファクターのスピードアップをする事ができる。
【００２８】
本発明の他の目的は、特別な目的のハードウェアを使うことができるが、要求はしないコストパフォーマンスの良い方法であって、標準のシリアルあるいはパラレルアーキテクチャに容易に実装できる方法を提供することである。
【００２９】
本発明の更に他の目的は、ＦＦＴを使った方法のような手法に起因したアーチファクトなしにスピードアップをすることができる方法を提供することである。
本発明の更なる他の目的は、柔軟で診察システムに一般に用いられるファンビームや他の幾何構成を適用できる方法を提供することである。
【００３０】
（発明の開示）
上記目的は、与えられた２次元イメージに対するラドン変換の投影あるいは角度サンプルの集合を計算することによって、再投影を行うために高速の計算方法を提供する本発明の方法によって達成できる。本発明の一実施形態に於いては、イメージのデータ表示を生成するために、再投影され、再結合された幾つかのサブイメージにイメージを分解する。イメージがサイズにおいて１ピクセルよりも大きい場合は、イメージをサブイメージに分解し、サブイメージは、サイズが１ピクセルになるまで、更に連続してより小さなサブイメージに分解される。サブイメージがサイズにおいて１ピクセルである場合には、これらは再投影される。従って、結果のイメージは、集められてトモグラフィの再投影イメージを生成し、表示される。収集は、全ての動径方向にシフトされた投影サブイメージが結合された後終わる。この処理は、反復処理、すなわち、有限数のステップでシステマチックに計算した処理値を許容する処理である。この方法の計算及びハードウェアのコストは、この方法がイメージを再構成するスピードに比べて低いものである。
【００３１】
略号テーブル
本特許は、幾つかの略号を使用する。以下のテーブルは、幾つかの略号の意味を読者が決定するための助けとして設けている。
【００３２】
ＣＰＵ＝中央処理ユニット
ＣＴ＝コンピュータトモグラフィ
ＥＩＴ＝電気インピーダンストモグラフィ
ＦＦＴ＝高速フーリエ変換
ＦＲＡ＝フーリエ再構成アルゴリズム
ＭＲＩ＝磁気共鳴イメージング
ＮＤＥ＝非破壊評価
ＮＩＨ＝National Institute of Health
ＰＥＴ＝陽子放出トモグラフィ
ｒｍｓ＝二乗平均平方根
ＳＰＥＣＴ＝単一光子放出数トモグラフィ
（発明の詳しい説明）
本発明は、コンピュータの処理能力及び装置のコストをかなり引き下げて、所定の角度範囲に渡って再投影を許容する処理である。この処理を理解するために、ラドン変換
【００３３】
【数５】
をその変数ｒとθの関数として見なしていくつかの知見を得る必要がある。
最初の２つの特徴は、ラドン変換の定義からすぐに得られるものである。特徴１は、ラドン変換は、線形変換であることである。
【００３５】
特徴２は、ラドン変換のシフト特性を定義するものである。ｆc,0 （ｘ、ｙ）＝ｆc （ｘ−ｘ0 、ｙ−ｙ0 ）が、ベクトル（ｘ0 、ｙ0 ）によってシフトされるイメージを示すとする。そのラドン変換は、外８ の動径方向にシフトされたものである。すなわち、
【００３６】
【外８】
【数６】
ここで、Ｓ（ｘ0 、ｙ0 ）は、元のものから、ベクトル（ｘ0 、ｙ0 ）によって移動された新しい原点に、投影を変換するシフトオペレータを示す。
特徴３は、投影の集合であるシノグラムの角度帯域制限である。第３の特徴は、「ボウネクタイ」性からくるもので、P.A. Rattey and A.G. Lindgren,"Sampling the 2-D Radon Transform," IEEE Trans. Acoust. Speech,Signal Proc., vol. 29, no. 5, pp.994-1002, 1981を参照されたい。これには、シノグラム 外９ の２次元フーリエ変換 外１０ の２つの変数ｒとθに対する本質的な台について書いている。
【００３９】
【外９】
【外１０】
なお、外１１ は、θについて周期的であるので、外１２は、ｕθ の直線スペクトルである。
【００４２】
【外１１】
【外１２】
あるいは、θに関するフーリエ級数を使うこともできる。本質的台外で 外１３ が０になる速度の厳密な評価は、F. Natterer,"The Mathematics of Computerized Tomography," Chichester, UK: Wiley, 1986 に与えられている。
【００４５】
【外１３】
外１４ のサポートが｜ｒ｜＜Ｒ（Ｒは、測定対象が収まるであろう最小円盤の半径である）で、外１５ が｜ｕr ｜＞Ｂr で実質的に０になるならば、
【００４７】
【外１４】
【外１５】
外１６ は、｜ｕθ ｜＞Ｒｂr で実質的に０になる。
【００５０】
【外１６】
従って、ｆc が｜ｘ｜≦Ｄ／２、｜ｙ｜≦Ｄ／２の領域で０でないならば、外１７ の台は、｜ｒ｜≦Ｄ／√２に限定される。
【００５２】
【外１７】
更に、Ｂ×Ｂの正方形に帯域制限された補間関数を使用するイメージｆc （ｘ、ｙ）の式（１）における近似は、ｆc が実質的にその変数について帯域制限されていることを示す。スライス−プロジェクション定理によれば、外１８ のｒについての帯域幅は、√２・Ｂで制限され、
【００５４】
【外１８】
特徴３から、外１９ のθに対する帯域幅は、√２・ＢＤ／√２＝ＤＢによって制限される。
【００５６】
【外１９】
前述の議論と、Ｂ＝Ｏ（Δ-1）であることに基づくと、特徴４は、重複のない投影の最小数は、Ｎ＝Ｄ／Δとして、イメージサイズＮに比例することである。これは、｜θ｜≦Θの範囲内の任意の値｛θi ｝に対して、正確に 外２０ を補間するために要求される、θの最低のサンプリングレートである。
【００５８】
【外２０】
本発明の一実施形態は、図２に見られるように、イメージの分解の単一レベルを行うものである。上記したように、既存の再投影は、ｃを定数として、ｃＮ3 の処理を必要とする。しかし、上記特徴によれば、本発明は、イメージの領域を小さく分割して、小さなサブイメージの集合とし、それらの投影を計算し、次に、結果を集めて、全体イメージの望まれる投影を得ることによって、処理数を減少している。この単一レベル領域分割処理では、イメージｆは、サイズが（Ｎ／Ｍ）×（Ｎ／Ｍ）のＭ2 個のサブイメージに細かく分解される。
【００６０】
図面、特に、図２を参照すると、正方形領域がＭ＝２の場合について示されている。各サブイメージに対しては、特徴４から、投影は、ｃ（Ｎ／Ｍ）3 のコストで、Ｏ（Ｎ／Ｍ）個の視角で、粗い角グリッド上の局所座標系で計算される。サブイメージの投影の集まりは、後に動径方向ｒにおけるシフトを伴う、角度座標θの補間を要求する。そして、線形性特徴１は、全体のイメージの投影を得るために、サブイメージの投影の加算に使用される。θ補間は、特徴４に適合するため、角度方向のサンプリングレートをＭ倍に増やす。というのも、イメージは全体として、サブイメージよりＭ倍大きいからである。θ補間による、Ｏ（Ｎ／Ｍ）個のサンプルをそれぞれ有するＯ（Ｎ）個の投影の計算コストは、ｃθ を定数として、ｃθ Ｎ2 ／Ｍで与えられる。上記特徴２によれば、動径方向のシフトは、サブイメージ座標系から全体イメージのそれへの平行移動を考慮する。この動径方向のシフトは、一般に、整数でない数のピクセル分である。しかし、それは、ｃr Ｎ2 ／Ｍ回の処理のコストを伴う補間となる。投影を計算する全コストは、
【００６１】
【数７】
となる。
上記の導出において、局所補間ができるように、θとｒの両方に適当なオーバサンプリングを適用することを仮定している。事実、そのようなオーバサンプリングによれば、補間誤差は、適切に選択された補間器の長さと共に、指数関数的に減少する。J.J. Knab."Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及び、O.M.Bucci, C. Gennarelli, and C. Savarese,"Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol 39, p. 1633, 1991 を参照されたい。Ｎが大きい場合には、Ｍ＝［ｃ／（ｃθ ＋ｃr ）］1/2 Ｎ1/2 という漸近的に最適な選択は、元のコストＯ（Ｎ3 ）と比較した場合、Ｏ（Ｎ2.5 ）のようにしか、Ｃ1 はスケーリングしないことを示す。精度と計算量は、投影収集ステップにおいて使用される特定の補間方法に依存する。
【００６３】
本発明の本実施形態は、図３と図４を比較することによってより良く理解される。既知の再投影処理では、完全なイメージ１２を、直接、投影の集合である単一のシノグラム１３に再投影している。図４においては、しかし、イメージ１２は、まず、独立にサブシノグラム２２、２４、２６及び２８にそれぞれ再投影される、４つのサブイメージ１４、１６、１８、及び２０に分割される。サブシノグラム２２、２４，２６、及び２８は、角度方向に補間され、動径方向にシフトされて、それぞれサブシノグラム３０、３２、３４、３６にそれぞれなる。サブシノグラムは、サブシノグラム２２、２４、２６、及び２８のθ’とｒ’のスケールをサブシノグラム３０、３２、３４及び３６のθ及びｒに効率的に調整することにより、補間とシフトを受け、サブシノグラム３０、３２、３４、及び３６は、直接単一のシノグラム１３を生成するように、足し合わされる。
【００６４】
他の実施形態においては、Ｏ（Ｎ2.5 ）の計算コストは、単一レベル領域分割方法ではなく、本発明のマルチレベル方法を採用することによって、更に減少される。本発明の説明を簡単にするために、イメージのサイズＮを２のべき乗とし、イメージｆは、図２に示されるように、同じ大きさの４つのサブイメージに細かく分割されるとする。イメージを小さなサブイメージに細かく分割する処理は、サブイメージが単一のピクセルの大きさのレベルになるまで、繰り返し、反復的に行うことができることが理解されよう。単一のピクセルレベルにおいては、投影の計算は、ピクセル基底補間関数 外２１ の投影は、解析的に知られていると仮定
【００６５】
【外２１】
しており、Ｏ（１）だけの投影が必要であるとされているので、自明なものとなる。従って、本発明の細分化の処理は、投影の計算をサブイメージ投影のマルチレベル集積に帰する。
【００６７】
集積処理は、座標系の原点、すなわち、回転の中心をサブイメージの中心から、複合イメージのそれへの平行移動を必要とする。これは、式（３）で定義されるシフトオペレータＳを用いることによって達成される。本発明のマルチレベル計算シーケンスは、反復処理として、最も簡単に定式化できる。以下の疑似コードは、イメージｆの投影 外２２ の計算を記述する。
【００６８】
【外２２】
外２３ ＝Ｆａｓｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆ）
【００７０】
【外２３】
ｉｆ ｓｉｚｅ（ｆ）＝１×１（１ピクセル）あるいは、ある最小のサイズ
外２４ ＝Ｄｉｒｅｃｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆ）
【００７２】
【外２４】
ｅｌｓｅ
外２５ ＝０
【００７４】
【外２５】
ｆｏｒ ｉ＝１ ｔｏ サブイメージの数
外２６ ＝Ｆａｓｔ Ｒｅｐｒｏｊｅｃｔｉｏｎ（ｆS i ）
【００７６】
【外２６】
外２７＝ 外２８＋Ｓ[ｘi 、ｙi ］Ｉθ 外２９
【００７８】
【外２７】
【外２８】
【外２９】
ｅｎｄ ｆｏｒ
ｅｎｄ ｉｆ
ここで、Ｉθ は、θ補間オペレータであり、ｆS i は、イメージｆのｉ番目のサブイメージである。また、（ｘi 、ｙi ）は中心を、外３０ は、ｆS i の投影を示す。
【００８２】
【外３０】
オペレータＩθ は、補間によって、投影の数を２倍にする。反復の各レベルのコストは、Ｏ（Ｎ2 ）である。というのも、サブイメージの数は、各レベルで４倍になり、サブイメージのピクセル数に比例するサブイメージ毎の計算量は、４のファクターで減少する。従って、ｌｏｇＮレベルでは、マルチレベル領域分割処理の計算の複雑さは、Ｏ（Ｎ2 ｌｏｇＮ）である。
【００８４】
前節で述べた処理の数値的実装の前に、補間、基底関数の選択、及びオーバサンプリングの量に関する選択をしなくてはならない。これらのパラメータは、互いに関連し、処理の精度と複雑さに影響する。
【００８５】
反復処理におけるイベントの構成とシーケンスが、図５に示されている。図４のサブイメージ２０は、それ自身サブイメージ４０、４２、４４及び４６に細かく分割される。これらの細分イメージは、反復的に１６個のサブイメージに細かく分割され、各サブイメージが１ピクセルの大きさになるまで細分化処理は続けられる。その点に於いて、各ピクセルは、要求されるように補間され、シフトされるシノグラムに再投影され、サブイメージのシノグラムが合わせられると、シノグラム１３が生成される。
【００８６】
２つの補間ステップ、細かいレベルから粗いレベルに行く角度補間と動径方向のシフトを考慮する動径方向補間を、より詳しく述べる。先に述べたように、適切なオーバサンプリングによれば、補間器の長さに伴って補間エラーが漸近的に指数関数で減少すること約束する補間器が使用される。J.J. Knab."Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及びO.M. Bucci, C.Gennarelli, and C. Savarese, " Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol 39, p.1633, 1991を参照されたい。簡単のため、動径及び角度変数の両方について、線形補間のみを実装する。しかし、ちょうど良いイメージのサイズ及び近似的にのみ帯域制限された基底関数に対しては、漸近的振る舞いは得られず、角度及び動径補間器長、角度及び動径オーバサンプリングの量、全計算コスト、及び、エラーの性質のトレードオフを考えなければならない。
【００８７】
線形補間を正確にするために、基底関数は、十分滑らかでなくてはならず、オーバサンプリングも必要である。Crawfordへの米国特許番号４，７０９，３３３号にある、３次元Ｂスプラインのテンソル積からなる基底関数
【００８８】
【数８】
が、良い近似及び滑らかさの特性のために選択される。３次元Ｂスプラインは、単位パルスの４次畳み込み、特に、｜ｔ｜≦１／２に対して、ｂ0 ＝１で、その他は０とする場合、ｂ＝ｂ0 ＊ｂ0 ＊ｂ0 ＊ｂ0 で定義されることを思い出されたい。
【００９０】
ここで、図６を参照すると、既知のシェップ−ローガンファントムが示されている。L.A. Shepp and B.F. Logan," The Fourier Reconstruction of a Head Section," IEEE Trans. Nucl. Sci., vol. 21, no. 1, pp. 21-43, 1974を参照されたい。シェップ−ローガンファントムは、本発明の有効性を示すために選択された。既知の直接投影と本発明の処理を使った、Ｎ＝２５６の場合に計算されたシェップ−ローガンファントムの投影が、図７ａ及び図７ｂのそれぞれに等高線プロットとして示されている。図７ａと７ｂを見て比較して明らかなように、本発明の方法は、正確な投影を実現している。更に、数値実験は、再投影の精度は、角度及び動径方向のオーバサンプリングの量に依存していることが示されている。図７ｂは、動径方向及び角度方向の４倍のオーバサンプリングの結果を示している。
【００９１】
ここで、図８ａ及び８ｂを参照すると、図７ａ及び７ｂに示されるシノグラムの２つの断面、特に、外３１ 及び 外３２ が示されている。
【００９２】
【外３１】
【外３２】
実線は、本発明の高速処理を示しており、ダイヤ印は、既知の直接再投影を示している。本発明のマルチレベル処理によって得られる結果は、視覚的に、既知の直接再投影に一致している。従って、本発明の精度が数値的に確かめられた。
【００９５】
再投影の精度を測定するために、規格化二乗平均平方根（ｒｍｓ）エラーを定義する。
【００９６】
【数９】
ここで、外３３ は、本発明の高速処理及び既知の直接再投影によって、計算され
【００９８】
【外３３】
た投影 外３４ を示している。
【０１００】
【外３４】
ここで、インデックスｉ及びｊは、それぞれサンプルされた角度及び動径方向に対応する。
図９を参照すると、Δrms 対ｒのオーバサンプリング量のプロットが、θについての４倍角度オーバサンプリングについて示されている。最適な補間方法を使用することによって、よりよい性能のトレードオフが、少ないオーバサンプリングファクタで可能である。J.J. Knab, " Interpolation of Band-Limited Functions Using the Approximate Prolate Series," IEEE Trans. Inf. Theory, vol. 25, No. 6, pp. 717-20, 1979 及び、O.M. Bucci, C.Gennarelli, and C. Savarese, " Optimal Interpolation of Radiated Fields Over the Sphere," IEEE Trans. Antennas and Propag., vol. 39, p.1633, 1991 を参照されたい。
【０１０２】
本発明の計算上の利点は、様々なサイズのイメージを考慮することによって、確かめられる。図１０において、ＣＰＵ時間対イメージサイズＮが既知の直接再投影と本発明の高速処理による投影の計算について示されている。直接再投影のカーブの傾きは、予測されるＯ（Ｎ3 ）の振る舞いに対応している。また、高速処理のカーブの観測される傾向は、本発明の処理の予想されるＯ（Ｎ2 ｌｏｇＮ）の複雑さであることを支持している。従って、このプロットから、本発明は、既知の直接再投影法よりも高速な結果を達成することが明らかである。
【０１０３】
本発明の原理に従って作られたイメージング装置５０が、図１１に示されている。イメージング装置５０は、ＣＴスキャナあるいは、様々な他のイメージング装置である。イメージング装置５０は、頭などの対象物からの生のデータを生成するスキャナ５２を含んでいる。データは、レシーバ５４に送られ、それから、後処理装置（あるいはステップ）５６に送られる。リビニング（re-binning）などの処理は、後処理装置５６で行うことが可能である。後処理装置５６の出力は、装置（あるいはステップ）５８において再構成され、結果のイメージは、ディスプレイ装置６０に表示される。しかし、イメージがアーチファクト（例えば、頭の中の金属片などによって、）を有する場合は、先に述べたように、イメージをエラー修正装置（あるいは、ステップ）６４と再投影装置（あるいは、ステップ）６２に供給することによって、結果としてのエラーを取り除くことが可能である。再投影後のシノグラム出力は、再構成装置（あるいは、ステップ）５８の入力に供給される。再投影とエラー修正は、アーチファクトによって起きたエラーが修正されるまで、繰り返される。
【０１０４】
前述の記載から、イメージの再投影あるいはトモグラフィの投影を計算する改良された方法は、多くの望ましい性質及び利点を有することが示され記載されたことが理解される。これは、既存の技術と同じ品質の投影を１００倍速く生成することができる。本発明は、キーとなる機能、従って、トモグラフィのスキャナ、例えば、ＣＴ、ＰＥＴ、ＳＰＥＣＴなどの全体の再構成を、高価な特別の目的を持ったコンピュータハードウェアを要求することなく、スピードアップすることができる。にもかかわらず、これは、同じ品質の投影を生成する特別目的のハードウェアに実装することが可能であるが、同様のコストの特別目的ハードウェアシステムよりずっと速く動作する。
【０１０５】
本発明の方法の特別な実施形態について示し、記載したが、当業者によれば、発明の最も広い側面及び請求項に記載されるものから離れることなく、これに変更などを施すことができることが理解されるであろう。
【０１０６】
本発明のこれら及び他の特徴や利点は、詳細な説明及び図面を参照すれば、当業者によれば明確であろう。
【図面の簡単な説明】
【図１】 既知の平行線投影。
【図２】 本発明の例示的投影。
【図３】 既知の再投影処理の図である。
【図４】 本発明の原理に基づいて行われる階層的分解と再投影の図である。
【図５】 図４のサブイメージに行われる反復的サブイメージ分解の図である。
【図６】 本発明の動作の図示に使用されたシェップ−ローガンファントム（Shepp-Logan phantom ）の幾何構成の表示である。
【図７】 図７ａ及び図７ｂは本発明の高速処理と既知の直接再投影それぞれによって計算された図６のシェップ−ローガンファントムのシノグラム投影である。
【図８】 図８ａ及び図８ｂは、
【数１０】
【図９】 図６のシェップ−ローガンファントムに対する、二乗平均平方根（ｒｍｓ）再投影エラーΔrms 対動径方向のオーバサンプリングの量のプロットである。
【図１０】 既知の直接再投影と本発明の方法による投影を計算する際の中央処理ユニット（ＣＰＵ）時間対イメージサイズＮを比較する図である。
【図１１】 本発明の原理に基づいたイメージ装置のブロック図である。
Claims (9)
- 電子形式のイメージ（１２）をイメージシノグラム（１３）に再投影する方法であって、
該イメージ（１２）をサブイメージ（１４、１６、１８、２０）に分割し、
該サブイメージ（１４、１６、１８、２０）を、関連したサブイメージシノグラム（２２、２４、２６、２８）を生成するために再投影し、
角度方向のサンプリングレートを増加するように、角度方向に、該サブイメージシノグラム（２２、２４、２６、２８）を補間し、各角度方向に補間されたサブイメージシノグラムの回転の中心が、該イメージシノグラム（１３）の回転の中心にシフトするように、該角度方向に補間されたサブイメージシノグラムを動径方向にシフトし、該シフトされ、補間されたサブイメージシノグラム（３０、３２、３４、３６）を合わせて、該イメージシノグラム（１３）を得ることによって、該サブイメージ再投影ステップの結果を集めて、該イメージシノグラム（１３）を得ることを特徴とする方法。 - 前記再投影は、前記サブイメージ（１４、１６、１８、２０）のラドン変換を行なうことを特徴とする請求項１に記載の方法。
- 前記サブイメージ（１４、１６、１８、２０）が、１ピクセルの大きさになるまで、反復して前記サブイメージ（１４、１６、１８、２０）を細かく分割し、前記再投影及び収集ステップを行うステップを更に備えることを特徴とする請求項１に記載の方法。
- 前記収集ステップが、反復して行われることを特徴とする請求項３に記載の方法。
- 前記ステップの少なくとも１つは、特別な目的を持ったハードウェアで行われることを特徴とする請求項３に記載の方法。
- 対象物のイメージを生成する装置であって、
該対象物からデータを生成するスキャナ（５２）と、
該スキャナ（５２）データからの該イメージの少なくとも１つの投影を生成するプロセッサ（５６）と、
該少なくとも１つの投影から該イメージを再構成する手段（５８）と、
該再構成手段によって生成された該イメージのエラーを修正する手段（６４）と、
角度方向のサンプリングレートを増加するように、角度方向に、該サブイメージシノグラム（２２、２４、２６、２８）を補間し、各角度方向に補間されたサブイメージシノグラムの回転の中心が、イメージシノグラム（１３）の回転の中心にシフトするように、該角度方向に補間されたサブイメージシノグラムを動径方向にシフトし、該シフトされ、補間されたサブイメージシノグラム（３０、３２、３４、３６）を合わせて、該イメージシノグラム（１３）を得ることによって、エラー修正の後、該イメージをサブイメージに分割し、該サブイメージをサブイメージシノグラムに再投影し、該サブイメージシノグラムを集積し、該イメージシノグラムを得ることによって、エラー修正後のイメージに対し該イメージを再投影する手段（６２）と、
該イメージシノグラムを該再構成手段に供給する手段と、
エラーが修正された後、該再構成手段によって生成された該イメージを表示する手段（６０）と、
を備えることを特徴とする装置。 - 前記動径方向へのシフトステップは、更に、前記角度方向に補間されたサブイメージシノグラムを動径方向に補間することを含むことを特徴とする請求項１に記載の方法。
- 前記角度方向の補間は、異なる精度を持つことを特徴とする請求項１に記載の方法。
- 前記動径方向の補間は、異なる精度を持つことを特徴とする請求項７に記載の方法。
Applications Claiming Priority (5)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US33809299A | 1999-06-23 | 1999-06-23 | |
| US09/338,092 | 1999-06-23 | ||
| US09/419,415 US6263096B1 (en) | 1999-06-23 | 1999-10-15 | Multilevel domain decomposition method for fast reprojection of images |
| US09/419,415 | 1999-10-15 | ||
| PCT/US2000/017125 WO2000079473A1 (en) | 1999-06-23 | 2000-06-22 | A multilevel domain decomposition method for fast reprojection of images |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2003502766A JP2003502766A (ja) | 2003-01-21 |
| JP4555524B2 true JP4555524B2 (ja) | 2010-10-06 |
Family
ID=26991031
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2001504966A Expired - Lifetime JP4555524B2 (ja) | 1999-06-23 | 2000-06-22 | イメージの高速な再投影のためのマルチレベル領域分解方法 |
Country Status (7)
| Country | Link |
|---|---|
| US (1) | US6263096B1 (ja) |
| EP (1) | EP1221129B1 (ja) |
| JP (1) | JP4555524B2 (ja) |
| AT (1) | ATE373286T1 (ja) |
| CA (1) | CA2375318A1 (ja) |
| DE (1) | DE60036393T2 (ja) |
| WO (1) | WO2000079473A1 (ja) |
Families Citing this family (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR19990007929A (ko) | 1995-04-26 | 1999-01-25 | 데이빗로스클리블랜드 | 다면 반복 노광 방법 및 장치 |
| DE10135994A1 (de) * | 2001-07-24 | 2003-02-20 | Siemens Ag | Verfahren zur Reduzierung von Artefakten in CT-Bildern, die durch Strukturen hoher Dichte hervorgerufen werden |
| US6771732B2 (en) * | 2002-02-28 | 2004-08-03 | The Board Of Trustees Of The University Of Illinois | Methods and apparatus for fast divergent beam tomography |
| JP2006523891A (ja) * | 2003-03-26 | 2006-10-19 | ブリガム・ヤング・ユニバーシティ | 局所洗練を用いてｔ−スプライン及びｔ−ｎｕｒｃｃ表面を定義するためのシステム及び方法 |
| CN100407995C (zh) * | 2003-09-09 | 2008-08-06 | 伊利诺斯州立大学托管会 | 快速分层x线断层造影方法和装置 |
| WO2005065090A2 (en) * | 2003-12-30 | 2005-07-21 | The Mitre Corporation | Techniques for building-scale electrostatic tomography |
| US8774355B2 (en) * | 2004-06-30 | 2014-07-08 | General Electric Company | Method and apparatus for direct reconstruction in tomosynthesis imaging |
| US7602879B2 (en) * | 2005-11-02 | 2009-10-13 | Wisconsin Alumni Research Foundation | Method for increasing the resolution of a CT image during image reconstruction |
| DE102006004618A1 (de) * | 2006-02-01 | 2007-08-09 | Siemens Ag | Arbeitsablauf-basiertes Management von medizinischen Bilddaten |
| US20070291277A1 (en) | 2006-06-20 | 2007-12-20 | Everett Matthew J | Spectral domain optical coherence tomography system |
| DE102006044189A1 (de) * | 2006-09-20 | 2008-04-17 | Siemens Ag | Verfahren zum Bereitstellen von vielfältig verarbeiteten Bilddaten, Verfahren zum Verarbeiten vielfältiger Bilddaten und Röntgenbildsystem |
| FR2923152A1 (fr) * | 2007-11-06 | 2009-05-08 | Gen Electric | Procede d'acquisition d'une image radiologique tridimensionnelle d'un organe en mouvement |
| US8135186B2 (en) * | 2008-01-25 | 2012-03-13 | Purdue Research Foundation | Method and system for image reconstruction |
| US9354455B2 (en) | 2009-01-23 | 2016-05-31 | Rodenstock Gmbh | Controlling designs using a polygonal design |
| JP5455210B2 (ja) * | 2009-10-29 | 2014-03-26 | ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー | Ｘ線ｃｔ装置およびｘ線ｃｔ画像再構成方法 |
| WO2011100575A2 (en) * | 2010-02-11 | 2011-08-18 | Emory University | Systems, methods and computer readable storage mediums storing instructions for applying multiscale bilateral filtering to magnetic resonance (mr) images |
| US9569891B2 (en) * | 2010-07-16 | 2017-02-14 | Tyoterveyslaitos | Method, an apparatus and an arrangement for visualizing information |
| CN102456228B (zh) * | 2010-10-29 | 2015-11-25 | Ge医疗系统环球技术有限公司 | 图像重建方法和装置及ct机 |
| RU2605519C2 (ru) * | 2011-10-24 | 2016-12-20 | Конинклейке Филипс Н.В. | Двухпроходная коррекция металлического артефакта с компенсацией движения для изображений срезов компьютерной томографии |
| KR102220226B1 (ko) * | 2014-02-12 | 2021-02-25 | 삼성전자주식회사 | 컴퓨터 단층 촬영 장치 및 그에 따른 ct 영상 복원 방법 |
| US11210804B2 (en) * | 2018-02-23 | 2021-12-28 | Sony Group Corporation | Methods, devices and computer program products for global bundle adjustment of 3D images |
Family Cites Families (45)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4042811A (en) | 1975-11-28 | 1977-08-16 | Picker Corporation | Tomography system having an ultrahigh-speed processing unit |
| US4149247A (en) | 1975-12-23 | 1979-04-10 | Varian Associates, Inc. | Tomographic apparatus and method for reconstructing planar slices from non-absorbed and non-scattered radiation |
| US4217641A (en) | 1978-04-28 | 1980-08-12 | U.S. Philips Corporation | Correction for polychromatic X-ray distortion in CT images |
| US4491932A (en) | 1981-10-01 | 1985-01-01 | Yeda Research & Development Co. Ltd. | Associative processor particularly useful for tomographic image reconstruction |
| US4590558A (en) * | 1981-12-30 | 1986-05-20 | General Electric Company | Method and apparatus for removing objects from CT images |
| US4626991A (en) | 1983-04-21 | 1986-12-02 | Elscint Incorporated | System for reprojecting images acquired by backprojection |
| US4714997A (en) | 1983-06-02 | 1987-12-22 | Elscint Incorporated | Data reduction in reprojection systems |
| US4616318A (en) | 1983-06-07 | 1986-10-07 | Elscint, Inc. | System for reprojecting images using transform techniques |
| GB2169180B (en) * | 1984-12-28 | 1988-06-15 | Toshiba Kk | Ct apparatus and operating method therefor |
| US4991093A (en) | 1985-08-21 | 1991-02-05 | Exxon Research And Engineering Company | Method for producing tomographic images using direct Fourier inversion |
| US4709333A (en) | 1986-01-03 | 1987-11-24 | General Electric Company | Method and apparatus for imaging in the presence of multiple high density objects |
| US4858128A (en) | 1986-08-11 | 1989-08-15 | General Electric Company | View-to-view image correction for object motion |
| US4930076A (en) | 1988-11-09 | 1990-05-29 | General Electric Company | Systolic radon transform processor |
| US5008822A (en) | 1988-11-25 | 1991-04-16 | Picker International, Inc. | Combined high speed backprojection and forward projection processor for CT systems |
| US5253308A (en) * | 1989-06-21 | 1993-10-12 | Amber Engineering, Inc. | Massively parallel digital image data processor using pixel-mapped input/output and relative indexed addressing |
| US5136660A (en) | 1989-10-13 | 1992-08-04 | International Business Machines Corporation | Apparatus and method for computing the radon transform of digital images |
| US5229934A (en) | 1990-06-18 | 1993-07-20 | Picker International, Inc. | Post-processing technique for cleaning up streaks and artifacts in diagnostic images |
| US5253171A (en) * | 1990-09-21 | 1993-10-12 | General Electric Company | Parallel processing method and apparatus based on the algebra reconstruction technique for reconstructing a three-dimensional computerized tomography (CT) image from cone beam projection data |
| US5224037A (en) | 1991-03-15 | 1993-06-29 | Cti, Inc. | Design of super-fast three-dimensional projection system for Positron Emission Tomography |
| US5396528A (en) | 1991-06-28 | 1995-03-07 | General Electric Company | Tomographic image reconstruction using cross-plane rays |
| US5243664A (en) | 1991-09-16 | 1993-09-07 | Picker International, Inc. | Post-processing technique for reducing metallic clip artifacts in CT images |
| US5375156A (en) | 1992-03-31 | 1994-12-20 | Siemens Medical Systems, Inc. | Method and apparatus for 3-D computer tomography |
| US5414623A (en) | 1992-05-08 | 1995-05-09 | Iowa State University Research Foundation | Optoelectronic system for implementation of iterative computer tomography algorithms |
| US5300782A (en) | 1992-06-26 | 1994-04-05 | General Electric Company | Gamma ray detector for pet scanner |
| JPH06259541A (ja) * | 1992-10-30 | 1994-09-16 | Toshiba Corp | 画像歪み補正方法およびそのシステム |
| JPH06186341A (ja) * | 1992-12-21 | 1994-07-08 | Hamamatsu Photonics Kk | ポジトロンコンピュータ断層撮影装置 |
| FR2701135B1 (fr) | 1993-01-29 | 1995-03-10 | Commissariat Energie Atomique | Procédé de reconstruction d'images tridimensionnelles d'un objet évoluant. |
| US5406479A (en) * | 1993-12-20 | 1995-04-11 | Imatron, Inc. | Method for rebinning and for correcting cone beam error in a fan beam computed tomographic scanner system |
| US5438602A (en) | 1993-12-23 | 1995-08-01 | General Electric Company | Correction of CT attenuation data using fan beam reprojections |
| US5552605A (en) | 1994-11-18 | 1996-09-03 | Picker International, Inc. | Motion correction based on reprojection data |
| US5559335A (en) | 1994-12-28 | 1996-09-24 | The University Of Utah | Rotating and warping projector/backprojector for converging-beam geometries |
| US5579358A (en) | 1995-05-26 | 1996-11-26 | General Electric Company | Compensation for movement in computed tomography equipment |
| US5625190A (en) | 1995-07-11 | 1997-04-29 | General Electric Company | Forward projection algorithm |
| JP3373720B2 (ja) | 1996-03-25 | 2003-02-04 | 株式会社日立メディコ | Ｘ線断層撮影装置 |
| JP4163767B2 (ja) | 1996-05-02 | 2008-10-08 | シーメンス アクチエンゲゼルシヤフト | コンピュータ断層撮影装置の画像再構成方法 |
| US5778038A (en) * | 1996-06-06 | 1998-07-07 | Yeda Research And Development Co., Ltd. | Computerized tomography scanner and method of performing computerized tomography |
| US5825031A (en) * | 1996-10-11 | 1998-10-20 | Board Of Regents The University Of Texas System | Tomographic pet camera with adjustable diameter detector ring |
| US5878102A (en) * | 1996-10-15 | 1999-03-02 | International Business Machines Corporation | Use of computer tomography scout images as an alternative to sinogram data |
| US5727041A (en) | 1996-11-13 | 1998-03-10 | General Electric Company | Methods and apparatus for reducing partial volume image artifacts |
| US6008493A (en) * | 1997-05-30 | 1999-12-28 | Adac Laboratories | Method and apparatus for performing correction of emission contamination and deadtime loss in a medical imaging system |
| US5862198A (en) | 1997-09-30 | 1999-01-19 | Siemens Corporate Research, Inc. | Pre-calculated hitlist for reducing run-time processing of an exact cone beam reconstruction algorithm |
| US5901196A (en) * | 1997-09-30 | 1999-05-04 | Siemens Corporate Research, Inc. | Reduction of hitlist size in spiral cone beam CT by use of local radon origins |
| US6108007A (en) | 1997-10-09 | 2000-08-22 | Silicon Graphics, Inc. | Method, system, and computer program product for increasing interpolation precision using multi-channel texture mapping |
| US6026142A (en) * | 1998-04-02 | 2000-02-15 | International Business Machines Corporation | System and method for boundary detection in tomographic images by geometry-constrained edge detection of projected data |
| US6028907A (en) * | 1998-05-15 | 2000-02-22 | International Business Machines Corporation | System and method for three-dimensional geometric modeling by extracting and merging two-dimensional contours from CT slice data and CT scout data |
-
1999
- 1999-10-15 US US09/419,415 patent/US6263096B1/en not_active Expired - Lifetime
-
2000
- 2000-06-22 DE DE60036393T patent/DE60036393T2/de not_active Expired - Lifetime
- 2000-06-22 AT AT00943031T patent/ATE373286T1/de not_active IP Right Cessation
- 2000-06-22 CA CA002375318A patent/CA2375318A1/en not_active Abandoned
- 2000-06-22 EP EP00943031A patent/EP1221129B1/en not_active Expired - Lifetime
- 2000-06-22 WO PCT/US2000/017125 patent/WO2000079473A1/en active IP Right Grant
- 2000-06-22 JP JP2001504966A patent/JP4555524B2/ja not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JP2003502766A (ja) | 2003-01-21 |
| CA2375318A1 (en) | 2000-12-28 |
| EP1221129B1 (en) | 2007-09-12 |
| US6263096B1 (en) | 2001-07-17 |
| ATE373286T1 (de) | 2007-09-15 |
| EP1221129A1 (en) | 2002-07-10 |
| DE60036393D1 (de) | 2007-10-25 |
| EP1221129A4 (en) | 2006-06-21 |
| DE60036393T2 (de) | 2008-06-12 |
| WO2000079473A1 (en) | 2000-12-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4555524B2 (ja) | イメージの高速な再投影のためのマルチレベル領域分解方法 | |
| JP4500312B2 (ja) | 高速階層型断層撮影法および装置 | |
| Natterer | The mathematics of computerized tomography | |
| JP4714316B2 (ja) | イメージ再構成 | |
| US4991093A (en) | Method for producing tomographic images using direct Fourier inversion | |
| O'Connor et al. | Fourier-based forward and back-projectors in iterative fan-beam tomographic image reconstruction | |
| CN107016672B (zh) | 医学扫描图像的重建方法和装置以及医学成像系统 | |
| EP1644897B1 (en) | A fourier tomographic image reconstruction method for fan-beam data | |
| EP0928458B1 (en) | A computer tomographic method and a computer tomograph | |
| Boag et al. | A multilevel domain decomposition algorithm for fast O (N/sup 2/logN) reprojection of tomographic images | |
| US6351548B1 (en) | Fast hierarchical reprojection algorithm for tomography | |
| Wen et al. | An investigation on the property and fast implementation of a ray-driven method for inversion of the attenuated Radon transform with variable focusing fan-beam collimators | |
| Yan et al. | Derivation and analysis of a filtered backprojection algorithm for cone beam projection data | |
| Bleuet et al. | An adapted fan volume sampling scheme for 3-D algebraic reconstruction in linear tomosynthesis | |
| Wang et al. | Exact and approximate conebeam x-ray microtomography | |
| Fessler | Iterative tomographic image reconstruction using nonuniform fast Fourier transforms | |
| Taguchi et al. | Cone-beam image reconstruction using spherical harmonics | |
| US6307911B1 (en) | Fast hierarchical backprojection for 3D Radon transform | |
| Hamill et al. | Iterative reconstruction methods for high-throughput PET tomographs | |
| Rajgopal et al. | Image reconstruction from incomplete projection data: A linear prediction approach | |
| You et al. | A harmonic decomposition reconstruction algorithm for spatially varying focal length collimators | |
| Wen et al. | Speed up of an analytical algorithm for nonuniform attenuation correction by using PC video/graphics card architecture | |
| Zeng et al. | Cone-beam and fan-beam image reconstruction algorithms based on spherical and circular harmonics | |
| Varga et al. | Iterative high resolution tomography from combined high-low resolution sinogram pairs | |
| Chen et al. | Accurate perspective projection calculation using a pixel-pyramid model for iterative cone beam reconstruction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20070613 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20100223 |
|
| A601 | Written request for extension of time |
Free format text: JAPANESE INTERMEDIATE CODE: A601 Effective date: 20100519 |
|
| A602 | Written permission of extension of time |
Free format text: JAPANESE INTERMEDIATE CODE: A602 Effective date: 20100526 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20100623 |
|
| 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: 20100713 |
|
| 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: 20100716 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130723 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4555524 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| 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 |