JP5543448B2 - 高効率コンピュータ断層撮影 - Google Patents
高効率コンピュータ断層撮影 Download PDFInfo
- Publication number
- JP5543448B2 JP5543448B2 JP2011516784A JP2011516784A JP5543448B2 JP 5543448 B2 JP5543448 B2 JP 5543448B2 JP 2011516784 A JP2011516784 A JP 2011516784A JP 2011516784 A JP2011516784 A JP 2011516784A JP 5543448 B2 JP5543448 B2 JP 5543448B2
- Authority
- JP
- Japan
- Prior art keywords
- data
- projection
- object space
- space
- voxels
- 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.)
- Active
Links
- 238000002591 computed tomography Methods 0.000 title description 3
- 238000000034 method Methods 0.000 claims description 64
- 230000004044 response Effects 0.000 claims description 49
- 238000013500 data storage Methods 0.000 claims description 31
- 230000005284 excitation Effects 0.000 claims description 14
- 238000009499 grossing Methods 0.000 claims description 12
- 238000001514 detection method Methods 0.000 claims description 11
- 238000007639 printing Methods 0.000 claims description 10
- 238000005259 measurement Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 4
- 239000002245 particle Substances 0.000 claims description 2
- 230000002093 peripheral effect Effects 0.000 claims description 2
- 238000000605 extraction Methods 0.000 claims 1
- 230000006870 function Effects 0.000 description 33
- 239000011159 matrix material Substances 0.000 description 17
- 238000004364 calculation method Methods 0.000 description 15
- 230000004048 modification Effects 0.000 description 7
- 238000012986 modification Methods 0.000 description 7
- 238000012545 processing Methods 0.000 description 7
- 238000012546 transfer Methods 0.000 description 7
- 210000000078 claw Anatomy 0.000 description 6
- 238000012937 correction Methods 0.000 description 6
- 238000003384 imaging method Methods 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000005316 response function Methods 0.000 description 5
- 230000015572 biosynthetic process Effects 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 238000003786 synthesis reaction Methods 0.000 description 4
- 238000003325 tomography Methods 0.000 description 4
- 238000002583 angiography Methods 0.000 description 3
- 238000004891 communication Methods 0.000 description 3
- 230000007812 deficiency Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000012800 visualization Methods 0.000 description 3
- 238000010009 beating Methods 0.000 description 2
- 230000015556 catabolic process Effects 0.000 description 2
- 210000004027 cell Anatomy 0.000 description 2
- 210000000170 cell membrane Anatomy 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 238000006731 degradation reaction Methods 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000009472 formulation Methods 0.000 description 2
- 239000004973 liquid crystal related substance Substances 0.000 description 2
- 230000015654 memory Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 238000012014 optical coherence tomography Methods 0.000 description 2
- 238000002600 positron emission tomography Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000002603 single-photon emission computed tomography Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000004627 transmission electron microscopy Methods 0.000 description 2
- 238000009966 trimming Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- ZCYVEMRRCGMTRW-UHFFFAOYSA-N 7553-56-2 Chemical compound [I] ZCYVEMRRCGMTRW-UHFFFAOYSA-N 0.000 description 1
- 241000202917 Spiroplasma Species 0.000 description 1
- 241000116871 Spiroplasma melliferum BC-3 Species 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000002872 contrast media Substances 0.000 description 1
- 238000002586 coronary angiography Methods 0.000 description 1
- 210000004292 cytoskeleton Anatomy 0.000 description 1
- 238000013501 data transformation Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001493 electron microscopy Methods 0.000 description 1
- 210000001723 extracellular space Anatomy 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 210000002458 fetal heart Anatomy 0.000 description 1
- 230000008571 general function Effects 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 238000003706 image smoothing Methods 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000003834 intracellular effect Effects 0.000 description 1
- 239000000193 iodinated contrast media Substances 0.000 description 1
- 229910052740 iodine Inorganic materials 0.000 description 1
- 239000011630 iodine Substances 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000000386 microscopy Methods 0.000 description 1
- 238000000465 moulding Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000010349 pulsation Effects 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 210000003705 ribosome Anatomy 0.000 description 1
- 235000002020 sage Nutrition 0.000 description 1
- 239000000243 solution Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000013179 statistical model Methods 0.000 description 1
- 238000000859 sublimation Methods 0.000 description 1
- 230000008022 sublimation Effects 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 239000002023 wood Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/424—Iterative
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Biochemistry (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Health & Medical Sciences (AREA)
- Algebra (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Mathematical Optimization (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Image Processing (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Description
本発明は、画像再構成（image reconstruction）の分野に関する。
断層撮影（tomography）、すなわちｎ次元空間内の領域のｍ次元投影（ピクセルとして表され、通常０＜ｍ＜ｎである）に基づく、その領域中の密度の計算（推定）は、２つの主要カテゴリ、すなわちフィルタ補正逆投影（ＦＢＰ：filtered back projection）および代数的再構成法（ＡＲＴ：algebraic reconstruction technique）の変形としても知られる反復再構成（ＩＲ）に分けられる。ＦＢＰおよび従来のＩＲ法には両方とも、特に、１つの投影図が２０００ピクセル（１次元）から最大２０００×２０００ピクセル（２次元）までの範囲になることが多い場合に、１０６〜１０１０個の画像要素（すなわち、ボクセル）を計算することがしばしば求められる最新の設定において、特有の欠陥がある。
ＦＢＰに関連した重大な欠陥は、多数の投影が必要とされるのに、達成される定量的精度は限られることである。投影数は一般に、数百にもなるが、投影は、可能な限り効率的には使われない。例えば、物理的強度が負になることはないと分かっているのに、計算された密度推定値が負になる場合がある。限られた定量的精度は、後から改良を繰り返すことによって向上する場合もある。ＦＢＰに関連した他の欠陥は、例えば、非特許文献１で論じられているように、ボクセルの場所に依存してデータの重みを変更できないこと、および様々な制約を効果的に考慮できないことである。
ＦＢＰ法の強みは、その計算速度であるが、これは、ＦＢＰ法が、一般に高速フーリエ変換、中央断面定理、および場合によっては適切な事前計算フィルタの使用に依拠するからである。
ＦＢＰ法および従来のＩＲ法に対する代替技法は、行列演算を使用し、ボリューム要素またはピクチャ要素（ボクセルまたはピクセル）の数が数千の単位である「小規模」問題に適用可能であり得る。ただし、典型的なトモグラフィ再構成の設定用途では、計算負荷は、当面はこうした行列に基づく技法によっては対処され得ない。というのは、コンピュータ演算の数がＮの累乗の規模だからである。例えば、３次元のケースでは、演算カウントは、非特許文献１および特許文献１で論じられているように、Ｎ５より速く増え得る（ここで、Ｎは、対象領域を表す再構成立方体の一辺のピクセル数である）。
ＩＲ法の利点は、制約に対処できること、特に密度推定値が負にならないように保証できることである。負でない密度推定値を保証することは、大幅な画像コントラスト向上、およびより優れた定量的結果につながり得る。
従来のＩＲ法に関連した欠陥は、後に続く反復のための補正項を取得するために、反転演算および順投影演算を繰り返し解く必要があることである。これに伴い、こうした演算を何度も再帰的にペアリングする必要が生じることにより、通常、ＩＲ法は遅くなり、特定の実装形態によっては不安定になる。
こうした従来のＩＲ法は、様々な（反復）最適化方法を用いるので、解への収束率は、ある特定の技法がどの程度ヘシアン行列（Hessian matrix）に対処するかに依存し、ヘシアン行列は、最適化されるべき目的関数（すなわち、画像要素に関連した性能基準）の二次導関数行列である。最も普及している技法はＩＲの変形である。ただし、ヘシアン行列のサイズは画像再構成にとって大きいので、行列の構造（および行列の逆の構造）は一般に、無視され、または不十分な近似をなす。さらに、ヘシアンの固有値は広く分布するので、現在の最適化技法は、一定数の反復（一般に、数十から数百のカウントである）を越えると改善を示さない傾向があり、大きい固有値のみを扱い得る。
こうしたアルゴリズムのマルチグリッド変形体は有用であり得るが、結局、細かいグリッドに伴うヘシアンのサイズのせいで、やはりうまくいかない。本明細書でのマルチグリッド解像度は、反復が実施されると次第に細かくなっていく解像度の使用を指す。
Wood, S.L., Morf, M., A fast implementation of a minimum variance estimator for computerized tomography image reconstruction, IEEE, Trans. On Biomed. Eng., Vol. BME-28, No. 2, pp. 56-68, 1981
Box, G.E.P., Jenkins, G.M., Time Series Analysis: Forecasting and Control, San Francisco, CA:Holden-Day, 1976
Sage, A.P., J.L. Melsa, Estimation Theory with Applications to Communications and Control, New York:McGraw-Hill, 1971, page 262
Wiener, Norbert, Extrapolation, Interpolation, and Smoothing of Stationary Time Series. New York: Wiley, 1949
Trachtenberg,S., Dorward,L.M., Speransky,V.V., Jaffe,H., Andrews,S.B., Leapman,R.D., "Structure of the Cytoskeleton of Spiroplasma melliferum BC3 and Its Interactions with the Cell Membrane", J. Mol. Biol., 378, pp.776-787, 2008
Sage, Jarisch W. R. & Detwiler J. S., "Statistical Modeling of Fetal Heart Rate Variability", IEEE Trans. On Biomed. Eng., Vol. BME-27, No. 10, pp. 582-589, Oct. 1980
Jarisch W., Determination of Material Properties by Limited Scan X-ray Tomography, Technical Report USAF AFWAL-TR-81-4050, NTIS, HC A07/MF A01, Defense Technical Information Center, Cameron Station, Alexandria, Virginia 22314, 1981
Kolehmainen,V., Brandt,S.S., Structure-From-Motion Without Correspondence From Tomography Projections by Bayesian Inversion Theory, IEEE Trans. Medical Imaging, ITMID4, Vol.26, No. 2, pp. 238-248, Feb. 2007
本発明は、観察物体内に励起エネルギーを送信するための１つまたは複数の送信機と、観察物体への送信励起エネルギーに応答して１つまたは複数の検出装置によって受信されるエネルギーをエンコードする投影空間データを生成するための１つまたは複数の検出装置と、励起エネルギーを送信するための１つまたは複数の送信機および投影空間データを生成するための１つまたは複数の検出装置を制御するためのコントローラと、投影空間データと予測投影データの間の差を特徴づける第１の量を計算し、各インパルス応答が単位強度の参照ボクセルに対応する少なくとも１つのインパルス応答に対応する第２の量を計算し、第１の量および第２の量の反転関数を使ってアップデート値を計算し、アップデート値を使って観察物体を表す物体空間画像を再構成することによって、投影空間データを受信し、投影空間データを処理するための画像再構成装置とを備えるシステムを提供する。
本発明は、１つまたは複数のプロセッサと、１つまたは複数のプロセッサによって実行されるソフトウェアを格納する１つまたは複数のデータ記憶装置とを備え、ソフトウェアは、物体空間から投影空間に投影される物体を表す入力データを受信するためのソフトウェアコードと、入力データと予測投影データの間の差を特徴づける第１の量を計算するためのソフトウェアコードと、少なくとも１つのインパルス応答に対応する第２の量を計算するためのソフトウェアコードであって、各インパルス応答は、物体空間内の場所での単位強度の参照ボクセルに対応するソフトウェアコードと、第１の量および第２の量の反転関数を使ってアップデート値を計算するためのソフトウェアコードと、アップデート値を使って観察物体を表す物体空間画像を再構成するためのソフトウェアコードとを含む、ワークステーションを提供する。
本発明は、１つまたは複数のデータ記憶装置に格納されたソフトウェアコードを実行する１つまたは複数のプロセッサによって実行される方法を提供することができ、この方法は、投影画像取得システム、１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから、物体空間から投影空間に投影される物体を表す入力データを受信すること、解像度グリッド上の入力データと解像度グリッド上の予測投影データの間の差を特徴づける第１の量を計算すること、少なくとも１つのインパルス応答に対応する第２の量を計算することであって、各インパルス応答は、物体空間内の場所での単位強度の参照ボクセルに対応すること、第１の量および第２の量の反転関数を使ってアップデート値を計算すること、並びにアップデート値を使って観察物体を表す画像を再構成すること、再構成画像を出力装置に出力することであって、出力装置は、１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つであることを含み、入力データを受信すること、第１の量を計算すること、第２の量を計算すること、アップデート値を計算すること、画像を再構成すること、および再構成画像を出力することは、１つまたは複数のデータ記憶装置に格納されたソフトウェアコードに従って１つまたは複数のプロセッサによって実施される。
例示的な実施形態を関連図面とともに記載する。
本発明の例示的な実施形態による体系図を示す図である。
本発明の例示的な実施形態による撮像システムを示す図である。
本発明の例示的な実施形態によるフローチャートである。
従来のフィルタ補正逆投影（ＦＢＰ）再構成に基づく細胞の再構成画像を示す図である。
本発明の例示的な実施形態による方法に基づく細胞の再構成画像を示す図である。
模擬ファントムの正面像サンプル投影を示す図である。
同時反復再構成法（ＳＩＲＴ）に基づく再構成誤差を示す図である。
本発明の例示的な実施形態による方法に基づく再構成誤差を示す図である。
拍動液圧式の冠動脈ファントムにヨード造影剤を注入している間の１枚のＸ線投影図を示す図である。
本発明の例示的な実施形態による方法に基づいて取得されるファントムの再構成画像を示す図である。
本発明の例示的な実施形態を以下で詳述する。例示的な実施形態を記載する際、明確にするために特殊な用語が利用される。ただし、本発明は、選択された特殊な用語に限定されることは意図していない。本発明の広範な概念から逸脱することなく、他の等価な構成要素を利用され得ること、および他の方法が開発され得ることが当業者には理解されよう。本明細書で引用する全ての参考文献は、それぞれが個々に組み込まれているものとして参照によって組み込まれる。
図１は、本発明のある実施形態による体系図を示す。物体の内部構造をエンコードする情報を含む投影ピクセルデータ１０４は、画像再構成装置１０５を使って再構成画像１０６に変換されて、出力装置１０７上で可視化され得る。投影ピクセルデータ１０４は、実験による取得１０１、またはコンピュータ上でのシミュレーション１０３、またはより以前の実験もしくはシミュレーションにより記録された投影ピクセルデータ１０４を含むデータ記憶装置１０２からのものでよい。実験による取得１０１、データ記憶装置１０２、およびシミュレーション１０３は、直接またはネットワーク、例えば、ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）、インターネットなどを介して画像再構成装置１０５に投影ピクセルデータ１０４をリモートに与えることができる。投影ピクセルデータ１０４のみを図示しているが、概してデータは、観察物体の中にどの形の励起エネルギーを伝えた結果でもよく、従って、ＣＴ（コンピュータ断層撮影）に加え、例えば、電子顕微鏡法、ＭＲ（磁気共鳴）、ＰＥＴ（陽電子放出断層撮影）、ＳＰＥＣＴ（単光子放出コンピュータ断層撮影）、超音波、蛍光、ＭＰＭ（多光子顕微鏡法）、ＯＣＴ（光干渉断層法）などからのデータを含んでもよい。データ記憶装置１０２は、例えば、コンピュータメモリ、ＣＤ−ＲＯＭ、ＤＶＤ、光磁気（ＭＯ）ディスク、ハードディスク、フロッピー（登録商標）ディスク、ジップディスク、フラッシュドライブなどでよい。
図２は、本発明のある実施形態による撮像システムを示す。撮像システムは、コントローラ２０５と、調査者２１１（例えば、人間、コンピュータなど）との信号通信およびコントローラ２０５との同期のためのインタフェース２０６（例えば、グラフィックユーザインタフェース、キーボード、ジョイスティックなど）と、コントローラ２０５からの制御信号に応答して励起エネルギー２０２を生成し放出する１つまたは複数の送信機２０１と、投影ピクセルデータ１０４を生成するように構成された１つまたは複数の検出装置２０７とを含むことができる。投影ピクセルデータ１０４は、記憶装置に格納されるデジタル形式とすることができる。撮像システムはまた、その上にある物体２０３を受信するように構成され、送信機２０１および検出装置２０７と関係して翻訳し（translate）、または循環させるように動作可能な翻訳または循環テーブル２０４と、電気的に、またはインターネット、他のネットワーク、もしくはデータバスなどの任意選択のデータ輸送装置２０８を介して検出装置２０７およびコントローラ２０５に結合された画像再構成装置１０５、２０９とを含むことができる。データバスは、例えば、コンピュータの内側のコンピュータ構成要素の間またはコンピュータの間でデータを転送する電気線のサブシステムでよい。接続構成要素間のデータフローは単方向だが、接続構成要素間の通信は双方向でもよい。
励起エネルギー２０２は、例えば、Ｘ線エネルギー、電磁（ＥＭ）エネルギー、光エネルギー、赤外線（ＩＲ）エネルギー、粒子エネルギー（例えば、電子、中性子、原子ビーム）、超音波などの振動エネルギーなど、放射エネルギーの形とすることができる。励起エネルギー２０２は、例えば、ファントム、人間の患者、標本、またはそのどの組合せでもよい物体２０３の上に照射され得る。励起エネルギー２０２は、対応するエネルギータイプの送信機２０１によって放出され得る。励起エネルギー２０２は、物体２０３の中を伝播することができ、一部分は、適切な検出装置２０７によって受信され得る。検出装置２０７は、受信エネルギーを測定可能な電気信号にコンバートすることができ、測定された電気信号を、デジタル形式の投影ピクセルデータ１０４にさらにコンバートすることができる。
コントローラ２０５は、送信機２０１および検出装置２０７に接続された回路でよく、励起エネルギー２０２の送信と検出装置２０７の動作とを同期させるための制御信号を送信機２０１および検出装置２０７に送ることができる。回路は、アナログ、デジタル、または混合信号回路でよい。コントローラ２０５は、送信機２０１および検出装置２０７に制御信号を送るための１つまたは複数のプロセッサおよび１つまたは複数のメモリを有するコンピュータとすることもできる。
画像再構成装置１０５、２０９は、本発明のある実施形態による方法により物体２０３の画像を再構成するために、コントローラ２０５に応答し、投影ピクセルデータ１０４を受信して、高い計算効率で物体の高忠実度画像を生成することができる。画像再構成装置１０５、２０９は、本発明の例示的な実施形態に従ってコンピュータを動作させるためのソフトウェアコードを格納する１つまたは複数のデータ記憶装置をもつコンピュータとすることができる。コンピュータは、１つまたは複数のデータ記憶装置に格納されたソフトウェアコードを読み取り、実行するための１つまたは複数のプロセッサおよび１つまたは複数のメモリを含むことができる。別の例示的実施形態では、画像再構成装置１０５、２０９は、本発明の例示的な実施形態に従って動作可能な１つまたは複数のプログラム格納／実行装置を含むことができる。画像再構成装置１０５、２０９は、再構成画像１０６を生成することができる。
出力装置１０７、２１０は、再構成画像１０６を受信することができる。出力装置１０７、２１０は、可視化装置でもデータ記憶装置でもよい。可視化装置は、例えば、表示装置でも印刷装置でもよい。表示装置は、例えば、ＣＲＴ（陰極線チューブ）、ＬＥＤ（発光ダイオード）ディスプレイ、ＬＣＤ（液晶ディスプレイ）、ＤＬＰ（デジタルライトプロジェクション）モニタ、ＶＦＤ（真空蛍光ディスプレイ）、ＳＥＤ（表面伝導型電子放出素子ディスプレイ）、ＦＥＤ（電界放出ディスプレイ）、ＬＣＯＳ（反射型液晶素子）ディスプレイなどを含むことができる。印刷装置は、例えば、トナー式プリンタ、液体インクジェットプリンタ、固体インクプリンタ、染料昇華型プリンタ、並びに感熱式プリンタおよび紫外線（ＵＶ）プリンタなどのインクレスプリンタなどを含むことができる。印刷装置は、３次元立体成型装置など、３次元（３−Ｄ）で印刷することができる。
撮像システムは調査者２１１をさらに含み得る。調査者２１１は、出力装置２１０または画像再構成装置２０９から再構成画像１０６を受信し、アルゴリズム（例えば、予めプログラムされたルーチン、人工知能、機械学習など）を適用して物体２０３についての診断情報を抽出し、または送信機２０１、検出装置２０７、テーブル２０４、画像再構成装置２０９などのための制御パラメータを微調整するための人間またはプログラムされたコンピュータとすることができる。ある実施形態では、インタフェース２０６または出力装置２１０は必要ない場合がある。ある実施形態では、調査者２１１、コントローラ２０５、および画像再構成装置１０５、２０９は、同じコンピュータ上にあってもよいし、または別個のコンピュータ上にあってもよい。
本発明のある実施形態は、画像再構成装置１０５、２０９と同様にして画像を再構成するように構成された１つまたは複数のプロセッサを備えるワークステーションを提供することができる。ワークステーションは、撮像システム、データ記憶装置、またはコンピュータの少なくとも１つから入力データを受信することができる。入力データは、データバス、ケーブル、有線ネットワーク、無線ネットワークなどを介して受信され得る。ワークステーションは、再構成画像を受信するための出力装置１０７、２１０をさらに備えることができる。出力装置は、データ記憶装置、表示装置、印刷装置などでよい。データ記憶装置、表示装置、および印刷装置例は、上述した通りである。
図３は、本発明の例示的な実施形態による方法のフローチャートを示す。
ボックス３０１で、例えば、
−ａ＜ｆ（ｘ）＜ｂ 式（１）
など、適切な限定（トリミング）関数ｆ（ｘ）を選ぶことができ、上式で、−ａ＜０＜ｂであり、ｘ＝０のある程度の近傍に対して、関数ｆ（ｘ）は、
適当な小さいε＞０に対して、｜ｆ（ｘ）−ｘ｜＜ε 式（２）
に従ってほぼｘである。
−ａ＜ｆ（ｘ）＜ｂ 式（１）
など、適切な限定（トリミング）関数ｆ（ｘ）を選ぶことができ、上式で、−ａ＜０＜ｂであり、ｘ＝０のある程度の近傍に対して、関数ｆ（ｘ）は、
適当な小さいε＞０に対して、｜ｆ（ｘ）−ｘ｜＜ε 式（２）
に従ってほぼｘである。
関数ｆは、特に初期反復において過度の値をトリミングすることができる。例えば、このような関数の１つがａｔａｎ（ｘ）である。限定関数は、軽減してよく、厳密なカットオフ値ａ、ｂがなくてよい。別の限定関数例は、ｈ（ｘ）＝ｃ・ｘ＋（１−ｃ）・ａｔａｎ（ｘ））、０＜ｃ＜１とすることができる。例えば、非特許文献２は、関数ｆに対して使うことができる様々な非線形データ変換を含む。
ボックス３０２で、ｍ次元（ｍ−Ｄ）ピクセルの入力投影データ１０４を、所与のすべての投影方向に対して、可変解像度の連続する解像度グリッド（例えば、各後続解像度グリッドは、画像範囲を維持したまま、直前グリッドよりピクセルが少なくなる）上に配置することができる。例えば、連続グリッド中の解像度グリッドは、連続グリッド中の直前グリッドの解像度の半分の解像度をもつ。最も粗いサイズは、例えば、単一ｍ次元ピクセルであり得る。
通常、１＜ｍ＜ｎであり、ｍは全投影に対して等しくてよいことに留意されたい。ただし、概してｍは、全ての投影とともに変化してよく、ｍ＞＝０（例えば、特定のノイズあり推定問題の場合）である。
またボックス３０２で、ｎ次元（ｎ−Ｄ）ボクセルの初期セットを、その投影が連続解像度グリッド上の利用可能ｍ次元投影ピクセルをカバーするように、ｎ次元空間内に配置することができる。例えば、こうした初期ボクセル用の適切な値セットは、平均投影密度に対応する値でよく、各ボクセルｊに対して、ｖＬｊ＝ｌｏｇ（ｖｊ）を設定することができる。対数以外の単調関数を使ってもよい。例えば、初期単一ｎ次元ボクセル密度ｖｉを平均投影密度に等しく設定してよい。
ボックス３０３で、現在のグリッド解像度に対するｎ次元ボクセルｖおよびｍ次元ピクセルｐは両方とも、任意で適切なカーネル（kernel）で平滑化することができ、ここで、ｍ次元カーネルはほぼｎ次元カーネルの投影である。この演算は、ノイズあり観察を伴うカルマンフィルタの予測子アルゴリズム（アップデートなし）に対応し得る（ｖは、ｖＬから取得される）。関数は、例えば、ｖ＝ｅｘｐ（ｖＬ）でよい。不連続を含む他の関数を、追加密度制約に対処するのに利用してよい。例示的な平滑カーネル（smoothing kernel）は、例えば、０．４ボクセル間隔と２ボクセル間隔との間のシグマ幅のガウス形状をもち得る。ただし、ボックス３０７を参照して後述するように修正が可能である。
必要に応じて、ボクセルｖは、一定の値へのマップにより制限することができる。例えば、一部の領域では、ボクセル密度がゼロであることが事前に分かっている場合がある。また、非常に小さいボクセル密度に対する数値アンダーフローを避けるために、ｖＬｊは無視できる最小値に限定してよい。例示的な限定値は、最大ボクセル密度値ｖｋ ｍａｘの約１０−１５の密度を表し得る。
ボックス３０４で、予測投影ｍ次元ピクセルｐｐを、ｎ次元空間内のｎ次元ボクセルｖのセットに基づいて計算することができる。予測投影ｍ次元ピクセルｐｐは、連続解像度グリッドの中の１つの解像度グリッド上にあってよい。グリッド解像度が、直前の解像度グリッドから増大している場合、予測投影およびボクセル推定値ｖは、より粗い／先行解像度グリッドからの補間により取得することができる。
ボックス３０５で、解像度グリッド上の予測投影データと入力投影データ１０４との間の差を特徴づける第１の量を計算することができる。
入力投影データ１０４の投影ピクセル値ｐｉがゼロより大きい場合、第１の量は、例えば、
ｕｉ＝ｆ（ｌｏｇ（Ｌｉ）） 式（３）
となり、上式で、Ｌｉは、入力投影データ１０４の測定されたｐｉと、予測投影データの予測ピクセル値ｐｐｉとの間の誤差比であり、例えば、
Ｌｉ＝ｐｉ／ｐｐｉ 式（４）
で表される。
ｕｉ＝ｆ（ｌｏｇ（Ｌｉ）） 式（３）
となり、上式で、Ｌｉは、入力投影データ１０４の測定されたｐｉと、予測投影データの予測ピクセル値ｐｐｉとの間の誤差比であり、例えば、
Ｌｉ＝ｐｉ／ｐｐｉ 式（４）
で表される。
例えば、Ｌｉがｐｉとともに増大し、増大するｐｐｉとともに減少する関数など、Ｌｉに対して他の適切な関係も用いることができる。こうした関数は、入力投影データ１０４と予測投影データとの間のある程度の誤差を取り込むことができる。
上記のログ関数は、誤差比の対数でよいが、入力対予測投影値のより一般的な関数（例えば、透過型電子顕微鏡法における情報欠落（missing wedge）問題において使われる関数）でもよい。
入力投影データ１０４の投影ピクセル値ｐｉがゼロ以下の場合、第１の量は、例えば、
ｕｉ＝−ａ 式（５）
でよく、上式でａは限界を示す。
ｕｉ＝−ａ 式（５）
でよく、上式でａは限界を示す。
ボックス３０６で、ｎ次元空間内の単位強度の参照ボクセルに対する少なくとも１つのインパルス応答Ｒｐに対応する量を計算することができる。インパルス応答Ｒｐは、現在の解像度グリッドに依存して、ｐ個すべての投影図および逆投影図用の所望の各参照ｎ次元ボクセルに関して計算することができる。インパルス応答Ｒｐは、参照ｎ次元ボクセルの場所に依存し得る。インパルス応答Ｒｐおよびその逆は、好ましい一連の演算に依存して、ｍ次元投影それぞれに対して、またはｎ次元ボクセルに対して計算することができ、一連の演算は所望の精度に依存する。このことはＦＢＰ法に類似している。参照ｎ次元ボクセルの限定サブセットが、インパルス応答Ｒｐの評価のために選択され得る。一部のアプリケーションでは、再構成空間の中心にあるわずか１つのｎ次元参照ボクセルを選択すればよい。
インパルス応答は、対応する参照ｎ次元ボクセルに関していくつかのやり方で計算することができ、こうしたインパルス応答の逆は、ボックス３０７を参照して後述するように、第１の量をもつ後続の組合せに対して使われ得る。
例示的な実施形態では、ｎ次元再構成／物体空間の領域中のある特定の場所でのインパルス応答Ｒｐは、単位強度の参照ボクセルが考慮される場合、線形システムのインパルス応答でよい。
ある特定の場所でのこの単一参照ｎ次元ボクセルの中心が、ｐ個すべての投影方向に投影されると、ｐ個の単一ピクセル投影図（必ずしも同じ次元ｍのすべてではない）を作成することができる。こうしたｐ個の投影図が、再構成のために逆投影されると、元のボクセルの劣化版が作成され得る。劣化版は、元のボクセルと線形的に関係するが、ｎ次元空間全体に分布し得る。この劣化関数Ｒｐは、この特定の場所に対する再構成システムのインパルス応答（または点像分布関数ＰＳＦ）でよい。
有限再構成ボリュームの場合、このインパルス応答Ｒｐは、参照ｎ次元ボクセルの場所に依存し得る。例えば、投影図の数が小さく、グリッドが十分に細かいとき、インパルス応答Ｒｐは蜘蛛状でよい（例えば、「蜘蛛の脚」状の伸長または「蜘蛛」の中心場所から発する筋を示す）。隣接し合うボクセルは、同様の形状の「蜘蛛の脚」状の伸長を有し得るが、有限再構成領域が境界の所で伸長を切り得るので、異なる長さにもなり得る。さらに、ｐ個の円錐ビーム投影図を使う円錐ビームなどの幾何構造では、「蜘蛛の脚」がｐ個の円錐の固定先端を指し得るので、「蜘蛛の脚」の方向および相互角（mutual angle）は場所によって変わり得る。
ただし、近似目的のために、こうしたインパルス応答関数の緩慢な空間的変化を利用して、ｎ次元空間内の十分に狭いボクセル近傍に渡って同じインパルス応答関数を共有することができる。「蜘蛛の脚」の類似度を幾何学的に考察することにより、同じ近似インパルス応答を共有する十分に一致したボクセル領域のサイズを判定することができる。例えば、その内容に相対して、「蜘蛛の脚」のエンドポイントでのいくつか（例えば、数十）パーセントの積分誤差は、特に後続反復がこうした誤差を訂正することができる場合、全く許容可能であり得る。対照的に、ＦＢＰ法では通常、この切捨て効果を考慮しないが、非特許文献１に示すような行列定式化では考慮する。
またボックス３０６で、所与のボクセルに対するインパルス応答Ｒｐは、ノイズ分散からの寄与率と合成されてＲを取得することができ、Ｒは第２の量となる。入力投影データを測定してノイズ分散を得ることができる。この合成により、「蜘蛛」の中心での値が増大し得る。Ｒおよびその逆は、好ましい一連の演算に依存して、ｍ次元投影図それぞれに関しても、ｎ次元ボクセルに関しても計算することができ、演算は、所望の精度に依存し得る。想定ノイズレベルが低いときなどのように増大が小さい方が、（例えば、カルマン／ウィーナーフィルタ利得に対応する）Ｒの逆を精密に計算することが強く求められ得る。ただし、好都合には、多重解像度グリッドにおける計算により、残留投影誤差における信号部分が適度／低レベルに保たれ得るので、Ｒ−１の効率的計算が支援される。信号部分は、入力投影データから測定されるノイズ分散によって対処されない、残留投影誤差中の構造的情報でよい。
ボックス３０７で、計算された第２の量は逆にされて、例えば、（逆ヘシアンを使う代わりに）Ｐ＝Ｒ−１として表されるＰを得ることができる。このステップはフーリエ領域において実施することができる。本発明の例示的な実施形態は、後述するように、反転関数をフーリエ領域におけるカルマンまたはウィーナー利得として特徴づけることができる。
測定された入力投影データＹに対する観察モデルは、例えば、Ｙ＝ＨＸ＋Ｗでよく、ここでＸは物体密度を表し、Ｈはシステム変換関数を表し、Ｙは入力投影データ１０４を表し得る。
行列表記において、Ｘは、入手可能データＹから推定することができる。データＹは、測定ノイズから得た寄与率を含み得る。Ｘは、以下の統計的プロパティを有し得る。
Ｘ０＝Ｅ［Ｘ］ 式（６）
ＶＸ＝Ｅ［（Ｘ−Ｘ０）・（Ｘ−Ｘ０）ｔ］ 式（７）
上式で、Ｘ０はＸの期待値Ｅ（すなわち、平均）であり、ＶＸは物体密度Ｘの共分散行列である。
Ｘ０＝Ｅ［Ｘ］ 式（６）
ＶＸ＝Ｅ［（Ｘ−Ｘ０）・（Ｘ−Ｘ０）ｔ］ 式（７）
上式で、Ｘ０はＸの期待値Ｅ（すなわち、平均）であり、ＶＸは物体密度Ｘの共分散行列である。
投影ノイズＷは、以下のように表すことができる。
Ｅ［Ｗ］＝０、 式（８）
Ｖ＝Ｅ［Ｗ・Ｗｔ］ 式（９）
Ｅ［Ｗ・Ｘｔ］＝０ 式（１０）
上式で、Ｖは投影ノイズＷの共分散行列である。
Ｅ［Ｗ］＝０、 式（８）
Ｖ＝Ｅ［Ｗ・Ｗｔ］ 式（９）
Ｅ［Ｗ・Ｘｔ］＝０ 式（１０）
上式で、Ｖは投影ノイズＷの共分散行列である。
例えば、非特許文献３で論じられる最小分散（カルマン）利得、即ちＫは、ＶＸＨｔが可逆であるとすると、
Ｋ＝ＶＸＨｔ［ＨＶＸＨｔ＋Ｖ］−１＝ＶＸＨｔＲＹ −１ 式（１１）
によって、または
Ｋ＝｛Ｈｔ［Ｉ＋Ｖ（ＨＶ Ｘ Ｈｔ）−１］｝−１ 式（１２）
によって与えることができ、上式で、Ｉは識別行列である。この条件は、入力投影データ中のピクセルの数が推定画像ボクセルの数を超え得る推奨計算の初期の低解像度反復において満たされ得る。断層撮影のコンテキストでは、ＲＹ −１は、反転には大き過ぎる場合があるが、式（１２）は、実地のＫを得るためのＨの有用な修正を推奨し得る。式（１１）により、最大階数Ｖを与えると、行列ＲＹ −１から可能なゼロ固有値を除くことができる。ＨＶＸＨｔとＶとの間の関係により、画像Ｘの推定がＸの事前統計にどの程度まで依拠し得るか判定することができる。
Ｋ＝ＶＸＨｔ［ＨＶＸＨｔ＋Ｖ］−１＝ＶＸＨｔＲＹ −１ 式（１１）
によって、または
Ｋ＝｛Ｈｔ［Ｉ＋Ｖ（ＨＶ Ｘ Ｈｔ）−１］｝−１ 式（１２）
によって与えることができ、上式で、Ｉは識別行列である。この条件は、入力投影データ中のピクセルの数が推定画像ボクセルの数を超え得る推奨計算の初期の低解像度反復において満たされ得る。断層撮影のコンテキストでは、ＲＹ −１は、反転には大き過ぎる場合があるが、式（１２）は、実地のＫを得るためのＨの有用な修正を推奨し得る。式（１１）により、最大階数Ｖを与えると、行列ＲＹ −１から可能なゼロ固有値を除くことができる。ＨＶＸＨｔとＶとの間の関係により、画像Ｘの推定がＸの事前統計にどの程度まで依拠し得るか判定することができる。
複合スペクトル表記を用いると、物体のパワースペクトル密度はＳＸＸになり、観察転送フィルタはＨになり、投影ノイズスペクトルパワー密度はＳＶＶになり、ウィーナー（例えば、非特許文献４を参照）利得は、
Ｗ＝Ｈ＊ＳＸＸ／（Ｈ＊ＳＸＸＨ＋ＳＶＶ） 式（１３）
＝［Ｈ（Ｕ＋ＳＶＶ／（Ｈ＊ＳＸＸＨ）］−１ 式（１４）
として表すことができ、上式で、＊は共役関数であり、Ｕは単位伝達関数である。
Ｗ＝Ｈ＊ＳＸＸ／（Ｈ＊ＳＸＸＨ＋ＳＶＶ） 式（１３）
＝［Ｈ（Ｕ＋ＳＶＶ／（Ｈ＊ＳＸＸＨ）］−１ 式（１４）
として表すことができ、上式で、＊は共役関数であり、Ｕは単位伝達関数である。
スペクトル定式化は、断層アプリケーションにとって計算可能であり得る。トモグラフィ再構成において、物体密度の不確実性（例えば、ＳＸＸまたはＶＸは比較的大きくてよい）が、投影ノイズ（例えば、ＳＶＶまたはＶは比較的小さくてよい）に優先し、ウィーナー利得Ｗ（またはカルマン利得Ｋ）を見つける主要タスクは、式（１３）（または等価には式（１１））への有用な近似を見つけることであり得る。
従って、入力投影データ１０４を測定して、上述したように、ＶおよびＶＸを得ることができる。ＨＶＸＨｔおよびＨ＊ＳＸＸＨは、例えば、測定されたＶＸに基づいて計算することができる。ボックス３０６で言及したノイズ分散から得た寄与率は、例えば、Ｖ、ＨＶＸＨｔ、またはＨ＊ＳＸＸＨに起因する寄与率を指し得る。カルマン利得のための式（１２）およびウィーナーフィルタ利得のための式（１４）は直接計算することができ、さらに単純化することが可能であり得る。例えば、測定ノイズＷ（例えば、Ｖは対角優位でよい）の小さい相互相関、および投影される物体密度（例えば、ＨＶＸＨｔは、特に反復プロセス中は対角優位でよい）の適度な相互相関に対して、式（１２）または（１４）の複合表現は、個々の投影図（またはそのグループ）に対するブロック単位の利得を計算することによって単純化することができる。例えば、ブロック単位処理は、例えば、ボクセル（ｎ次元処理）またはピクセル（ｍ次元近似処理）単位で修正ローカルインパルス応答の逆を使うことによってさらに単純化することができる。ローカルインパルス応答の修正は、例えば、式（１２）に見られるノイズ対信号構造Ｖ（ＨＶＸＨｔ）−１の要素によって、または同様に式（１４）に見られるＳＶＶ／（Ｈ＊ＳＸＸＨ）によって与えられ得る。この修正はボクセル単位で適用することができる。修正は、「蜘蛛」中心領域（例えば、式（１２）の行列Ｖ（ＨＶＸＨｔ）−１または式（１４）のＳＶＶ／（Ｈ＊ＳＸＸＨ）の（近）対角線要素に対応する、インパルス応答関数のピーク付近）の値を増大させることによって近似することができる。Ｗが白色ノイズであり、投影される物体の不確定性が白色であり、推定二乗平均平方根（ＲＭＳ）ノイズ対信号比がｒのケースでは、蜘蛛の中心は、例えば、（１＋ｒ２）倍だけ増大し得る。この修正の後、カルマン利得およびウィーナーフィルタ利得を、例えば、式（１２）、（１４）による数値反転によって計算することができる。
Ｒｐは、上述したように、式（１１）、（１３）に表した上記行列（または転送フィルタ）Ｈと関係することに留意されたい。処理シーケンスは、カルマン利得Ｋまたはウィーナー利得Ｗ中に表されるシングルステップを、２つのステップのベクトル演算に分割することによって変更することができる。２つのステップは、最初は、例えば、投影データの逆投影であり、次いで、参照ボクセルＫのｎ次元フィルタＰＫ＝ＲＫ −１でのフィルタリングでよい。ある程度の近傍周辺ボクセルＫは、結果として生じたフィルタ出力を使うことができる。２つのステップは、最初は、例えば、ｐ個のｍ次元フィルタＰｉＫ（０＜＝ｉ＜ｐ）での投影データの事前フィルタリングであり、ここで各ＰｉＫは、ノイズ対信号比だけ修正されたｎ次元インパルス応答（ＰＳＦ）ＲＰＫの投影の反転から得ることができ、次いで、フィルタ出力の逆投影でよい。やはり、ある程度の近傍周辺ボクセルＫ（またはその投影図）は、その結果生じたフィルタ出力（アップデート）を使うことができる。ｎ次元空間全体をカバーするために、いくつかのボクセル場所およびその近傍を使ってよい。従って、式（１１）、（１３）に示す極めて大きい行列の反転は、はるかに小さい数字セットＲＫからなる適度なセットの反転に分解することができ、ここで各ＲＫは、物体空間中で極めて緩慢に変化することができ、繰り返し使うことが可能である。さらに、参照ボクセルに対する各ＲＫ −１は、フーリエ技法によって非常に効率的に計算することができる。
全体として、利得Ｐを計算するのに使われる近似の２つの側面は、高い統計的効率に寄与し得る。
第１に、低い周波数の場合、カルマン／ウィーナーフィルタ利得（ここでは利得Ｐで置き換えられる）は、十分に広い周波数範囲に渡るノイズ対信号比に対処し得るインパルス応答関数の原点の所で小さい定数を加えることによって推定することができる。インパルス応答関数の中心近くの小さい近傍は、原点の所で単なる定数ではなく「狭い」ガウスベルを加えることによって、例えば、「ピンク」測定ノイズに対処するように修正することができる。
第２に、より高い周波数の場合、利得Ｐに対応するフィルタは、高周波信号対ノイズ比によって特徴づけられるものとしてテーパリング／平滑化することができる。数値の精度を向上させるために、このテーパリング／平滑化関数は、現在の解像度グリッドにおける画像密度および／または推定画像密度のアップデートに対して実施することができる。一連の解像度グリッドにおける反復計算中にすべての段階で推定画像密度を低域フィルタリングすることにより、高周波画像成分を漸減するのを助けることができ、計算負荷が削減される。
ボックス３０８で、第１の量およびＰを使ってアップデート値を計算することができる。例えば、Ｐを第１の量と合成して、所与のボクセル（またはピクセル）画像要素に対する生（raw）アップデート値ｉを得ることができる。例えば、第１の量でＰを畳み込んで、生アップデート値を得ることができる。Ｐがｎ次元再構成空間内でゆっくり変化するのに従って、適切だが十分に小さい所与のボクセル近傍にＰを近似として印加することによって、計算時間を節約することができる。上述した式（１）、（２）のｆと同様に構築される限定関数ｇを、生アップデート値ｉに適用して、修正アップデート値ｉｇ＝ｇ（ｉ）を得ることができる。この修正は過度の訂正を防ぎ得る。本発明の例示的な実施形態を開発する間、生アップデート値を直接使用する多くの実験が失敗した。本発明者は、過度のアップデート値を制限するための限定関数およびその賢明な使用を思いつき、開発した。修正アップデートは、後続処理のためのアップデート値となり得る。
ボックス３０４〜３０８で、単純化をさらに検討することが、円錐ビーム幾何構造の再構成（例えば、血管造影法で使われる）にとって有用な場合があり、ここで投影の数ｐは、小さくてよい（ｐ＜＜Ｎであり、Ｎは、再構成ボリュームを通る代表的投影パスに沿ったボクセルの数を表し得る）。この場合、個々の「蜘蛛の脚」は、特に多重再構成グリッドの高解像度最終ステージにおいては相異なってよい。この高解像度ｎ次元空間内で多くの小領域に対するインパルス応答Ｒｐを計算するのは、時間がかかり得る。ｎ次元インパルス応答使用の近似として、次元ｍｉの比較的低次元で投影されたｐ個のインパルス応答（０＜−ｉ＜ｐ）を使うことができる。対応する修正された逆またはカルマン／ウィーナー利得を、投影の個々の誤差画像に繰り返し適用してよい。こうしたｍ次元インパルス応答は、蜘蛛状でよい上記のｎ次元インパルス応答の、それぞれのｍ次元投影空間への投影でよい。こうしたインパルス応答のプロパティは上述した通りでよく、すなわち、投影空間内でゆっくり変化する。ただし、ｐ個の投影のインパルス応答、こうした計算されたインパルス応答の逆、および後続畳込み（ここではｍ次元投影に適用される）を計算するための計算負荷は、特にｐ＜＜Ｎのとき、ｎ次元空間内での計算よりはるかに小さくなり得る。計算されたインパルス応答は、Ｒの場合と同様にして修正することができる。計算されたインパルス応答は、トモグラフィ再構成におけるノイズを削減するように帯域通過フィルタリングすることもできる。この断層反転技法は、投影自体が再構成に先立ってフィルタリングされるフィルタ補正逆投影と同様でよい。この技法は、円錐ビーム幾何構造に限定されなくてよく、例えば、並列ビーム幾何構造など、他の幾何構造にも適用することができる。
ボックス３０９で、修正アップデート値ｉｇを、ｎ次元ボクセルセットｖ中の対応する画像密度の値に加えて、ｎ次元ボクセルセットｖを更新することができる。加算は、上述した式（１）、（２）と同様に、限定関数によって境界を定められ得る。加算は、ｎ次元ボクセルセット全体の一部分に対して実施することができる。精度のために、ボクセル値ｖは、その投影が生画像データと最も合致するようにスケーリングしてよい。
ボックス３１０で、解像度グリッドが最も細かい場合、また、計算された予測と入力投影データとの間の差（例えば、残差）が（統計的に）十分に小さい場合、密度推定値を表す更新ｎ次元ボクセルセットｖ、ｖＬ、または両方（並びに必要に応じて推定予測および誤差）を、ボックス３１１に示すように出力装置に出力することができる。出力に先立って平滑化を適用することができる。それ以外の場合、解像度グリッドは、現在の解像度グリッドより大きいサイズの解像度グリッドに変更することができ（または可能性としては、ｎ次元平滑カーネルの重大度／範囲を低下させることができ）、ボックス３１２で、フローはボックス３０３に進み得る。
ボックス３１１で終了すると、物体１０６の高忠実度の再構成画像を取得することができる。この再構成画像１０６は、物体２０３の内部情報を非侵襲的および非破壊的に明らかにすることができ、診断用医用撮像および非破壊的評価において広範に適用される。再構成画像１０６は、出力装置１０７に関して上述したように、可視化装置上で可視化することも、コンピュータメモリまたは記憶媒体に格納することもできる。
例えば、非特許文献５にある投影ピクセルデータ１０４を使って、本発明の例示的な実施形態により達成可能な改良を証明した。図４Ａは、投影ピクセルデータ１０４のフィルタ補正逆投影（ＦＢＰ）再構成に基づく再構成画像１０６を示す。図４Ｂは、本発明の例示的な実施形態による再構成画像１０６を示す。図４Ａと４Ｂとの間の比較は、例えば、細胞内対細胞外空間、細胞膜、リボソーム、スピロプラズマの細部など、様々な密度のより明確な対比によって測定される、図４Ｂにおける優れた画像鮮明度を証明している。図４Ｂの画像はユーザにとってより望ましい。両方の再構成は、ほぼ＋／−７０度の範囲に渡るほぼ１６６個の二軸傾斜式の透過型電子顕微鏡法（ＴＥＭ）投影を用いた。
図５は、模擬ファントムの正面像サンプル投影を示す。模擬ファントムは、サイズが２５６×２５６×６４ボクセルの３次元ボックス内のランダムの場所に、サイズが６×６×６ボクセルの１００個程度の立方体を含む。こうしたボクセルは、０．８５単位の周辺密度に組み込まれた５単位の密度をもつ。白色ノイズが加算された合計で４９個の投影（単軸傾斜、＋／−７２度）が、コンピュータシミュレーション１０３によって取得されて、投影ピクセルデータ１０４を生成した。
統計モデル検証における公知のゴールドスタンダードによると、再構成品質は、残余検査、すなわち、実像（例えば、予め分かっている図５の模擬ファントム）と再構成画像１０６との間の不一致によって最も厳しく評価され得る。高忠実再構成は、ランダムであるとともに３次元物体（例えば、図５の模擬ファントム内の立方体）の構造からの非常に少ない寄与率を含む残余につながり得る。図６Ａは、同時反復再構成法（ＳＩＲＴ）の残余を示し、図６Ｂは、本発明の例示的な実施形態による方法から取得されるものを示す。図６Ａのかなり著しい残差平方パターン（１５回の反復の後）および図６Ｂの顕著な白色ノイズ出現（１度の高解像度反復の後）は、画像忠実度および計算効率における本発明の優越性をさらに証明する。ＳＩＲＴのための反復カウントは、ＳＩＲＴのためのそれより多いまたはより少ない反復が画像品質の劣化につながっただけであったので、この例のために最適化してある。
本発明は、投影の数が小さい状況、例えば、血管造影法における円錐ビーム幾何構造にも当てはまり得る。図７Ａは、プラスチックチューブで作られた鼓動式水圧冠動脈ファントムにヨード造影剤（通例、例えば、介入的心臓サブトラクション血管造影法で利用される）を注入している間の１枚のＸ線投影図を示す。用いられた投影ピクセルデータ１０４は、注入の一過性のため、基底ファントムのわずか６度のデジタルサブトラクション投影によるものであった。図７Ｂは、本発明の例示的な実施形態による再構成画像１０６を示し、鼓動式水圧ファントムの幾何構造を正確に表している。再構成画像は、フィルタバック投影（ＦＢＰ）によっても同時反復再構成法（ＳＩＲＴ）によっても、最大エントロピー（ＭＥＮＴ）法によってさえも取得することはできない。というのは、投影回数が異常に小さく、かつ／または反復回数が現在の計算機資源では実際的でなくなるからである。
上述したように、本発明の例示的ないくつかの実施形態では、計算／処理技法を合成し、結果として、対象領域内で所望のｎ次元ボクセルの高速で正確な計算／推定をもたらすとともに、任意のｍ次元投影（および方向）を使用できるようにする。ｍは固定しなくてよく、投影指標に依存し得ることに留意されたい。合成技法の効率は、広範な投影回数に当てはまり得る。入手可能な投影（測定）情報を効率的に使用することにより、投影回数が非常に大きいとき、および冠動脈血管造影法などにおいて制限される（小さいまたは不完全）ときの両方で本方法は効果的になり得る。
本発明の非線形反転技法は、（ｉ）再構成の領域内部の２つの非線形関連ドメインへのｎ次元画像推定値の分解による、非線形カルマンフィルタの一定部分の線形化（ベイズの推定の近似）、（ｉｉ）画像平滑化の低下度を用いることによって解像度をさらに細分した「階段」に匹敵する（「ステップ」に匹敵する）マルチグリッド計算、および（ｉｉｉ）必要とされるとき、逆ヘシアンを使う代わりに、画像ノイズの推定値によって修正された線形システムの逆インパルス応答（または伝達関数）の効率的計算であって、物体空間における選択領域に対する近似カルマン／ウィーナーフィルタ利得を生じる計算を合成している。
具体的には、非線形カルマンフィルタは、例えば、非特許文献６で、画像再構成のためのカルマンフィルタリングの適用は、例えば、非特許文献７で論じられている。例えば、ＦＢＰ、ＦＦＴ、およびウィーナーフィルタなどの代替方法を、反転の線形部分の代わりに使うことができる。再構成の領域内部の２つの非線形関連ドメインのどちらも、カルマンシステム状態を表し得る。
上記合成によって、最適化におけるヘシアンの効果的使用と同様の効率的な誤差削減を、すべての反復で達成することができる。達成された効率は、二次収束のものに近くなり得る。この手法により、例えば、非特許文献８にある手法など、他の手法によってとられる手法で見られる高次元であり乏しい収束率の重大な問題が避けられ得る。
あるいは、本発明の例示的な実施形態は、ＩＲと、本明細書ではＳ−ＦＢＰと呼ばれる（空間可変／非定常／分割）ＦＢＰとの合成と見なすことができる。近似の要件に依存して、Ｓ−ＦＢＰに関連づけられた１つまたはいくつかのセグメントがあり得る。概して、ＦＢＰは、再構成のフィールドにおいて定常ではない。ただし、フィールドプロパティをゆっくりと変更することにより、再構成のフィールドを、領域の間の遷移が重みづけされた重なり領域／セグメントに分けることによって近似が可能になり得る。Ｓ−ＦＢＰの速度、効率、および柔軟性は、後述するいくつかの実施形態の適切な選択および合成により達成することができる。
例えば、Ｓ−ＦＢＰ部分は、生投影データではなく、非線形変換されたデータに適用することができる。
例えば、マルチグリッド計算は、線形化をサポートし、各反復での演算カウントを非常に低く保つことができる。例えば、最高解像度グリッドでの最終計算に先立つ合計演算カウントは、最高解像度グリッドでのＳ−ＦＢＰのものよりも低くなり得る。
例えば、ボクセル密度は、２つの非線形関連ドメイン（マップにより接続される）、すなわちＳ−ＦＢＰによるアップデート用ドメイン（所望の近似技法に依存して、ｎ次元空間内において可変でよい）、および投影予測の計算用ドメインで表すことができる。Ｓ−ＦＢＰ部分のフィルタ重みは、ボクセルの場所に依存し得るが、同じ重みをもって使用するためのボクセルのグループ化により、計算効率が向上し得る。
例えば、線形化システムインパルス応答（または伝達関数）の逆から、Ｓ−ＦＢＰ部分のフィルタの重みを計算することによって、任意の投影次元および方向を用いることができる。線形化システムインパルス応答（または伝達関数）は、入力投影データの測定値から判定されたノイズ分散と合成済みでよい。
例えば、適切に選ばれた限定（トリミング）関数は、計算中の数値特異性および不安定性を防ぎ得る。
例えば、各反復においてボクセルおよびピクセル値に慎重に適用される、適切に選ばれた平滑化関数は、数値計算の安定性および精度に寄与し得る。このボクセルおよびピクセル値の平滑化は、アップデート値の生成に先立つカルマン予測子アルゴリズムに機能的に匹敵し得る。平滑化の反復適用は、カルマンフィルタにおける状態予測子アルゴリズムと同様でよい。
最高解像度グリッド上での計算に先行する反復の精度が十分であり得るので、最高解像度グリッドを通る反復は計算上必要なくてよく、従って、誤差訂正項の計算が節約され得る。こうした事例では、再構成速度は、Ｓ−ＦＢＰのものに匹敵するようになり得る。というのは、最高解像度グリッド上での計算により、ＦＢＰのものに匹敵する演算カウント（１より大きい因数の範囲内）が達成され、より粗い解像度グリッド上での先行するより低い解像度の反復すべてに対する合成演算カウントは、特により高次元の画像再構成に対しては、最高解像度グリッド上での計算の場合より少なくなり得るからであり得る。
図３に示す方法を参照して論じた本発明の例示的な実施形態は、例えば、ＣＤ−ＲＯＭ、ＤＶＤ、ＭＯ（光磁気）ディスク、ハードディスク、フロッピー（登録商標）ディスク、ジップディスク、フラッシュドライブなどのデータ記憶装置に格納されるソフトウェアコードとして提供され得る。格納されたソフトウェアコードは、図３を参照して上述した例示的な方法を実施するための、１つまたは複数のプロセッサと、例えば、ＲＡＭ（ランダムアクセスメモリ）装置、ＤＲＡＭ（ダイナミックＲＡＭ）装置、フラッシュメモリ装置、およびＳＲＡＭ（スタティックＲＡＭ）装置などの、１つまたは複数のメモリ装置とを有するコンピュータによって可読であり実行可能でよい。
本発明の例示的な実施形態は、図３を参照して上述した例示的な方法を格納し実行する１つまたは複数のプログラム記憶および実行装置、例えば、ＡＳＩＣ（特定用途向け集積回路）、ＦＰＧＡ（フィールドプログラム可能ゲートアレイ）、ＣＰＬＤ（複合プログラム可能論理素子）、ＡＳＩＰ（特定用途命令セットプロセッサ）などを提供し得る。
本明細書に記載した例および実施形態は、非限定的な例である。本発明を、例示的な実施形態を参照して詳しく記載したので、上記内容から、本発明からより広範な態様において逸脱することなく変更および修正を行ってよく、したがって、請求項において定義される本発明は、このような本発明の真の精神の範囲内であるあらゆる変更および修正をカバーすることを意図していることが当業者には明らかであろう。
Claims (32)
- 観察物体内に励起エネルギーを送信する１つまたは複数の送信機と、
前記観察物体への前記送信励起エネルギーに応答して１つまたは複数の検出装置によって受信されるエネルギーをエンコードする投影空間データを生成する前記１つまたは複数の検出装置と、
前記励起エネルギーを送信する前記１つまたは複数の送信機および前記投影空間データを生成する１つまたは複数の検出装置を制御するコントローラと、
前記投影空間データを受信する画像再構成装置であって、
複数の物体空間ボクセルのそれぞれについて第１の量を計算し、前記第１の量は、前記投影空間データの値および予測投影データの対応する値の間の比率の対数を示し、
複数の物体空間ボクセルのそれぞれについて第２の量を計算し、前記第２の量は、ノイズ分散からの寄与率と結合される少なくとも１つの反転インパルス応答に対応し、各インパルス応答は、単位強度の参照ボクセルに対応し、前記第２の量は、ノイズ分散からの寄与率、およびカルマンフィルタ利得、近似カルマンフィルタ利得、ウィーナーフィルタ利得または近似ウィーナーフィルタ利得の少なくとも一つを使用して計算され、
複数の物体空間ボクセルのそれぞれについてアップデート値を計算し、それぞれのアップデート値は、対応する第１の量および対応する第２の量を使用して計算され、
前記アップデート値を使用して前記観察物体を表す物体空間画像を再構成する
ことによって、
前記投影空間データを処理する画像再構成装置と
を備えることを特徴とするシステム。 - 前記励起エネルギーは、電磁（ＥＭ）エネルギー、Ｘ線エネルギー、粒子線、赤外線エネルギー、光エネルギー、または振動エネルギーの少なくとも１つであることを特徴とする請求項１に記載のシステム。
- 前記予測投影データは、前記物体空間ボクセルのセットに基づいて計算され、前記物体空間ボクセルのセットは、前記物体空間ボクセルのセットが投影空間に投影されると、可変サイズの複数の解像度グリッドをカバーすることを特徴とする請求項１に記載のシステム。
- 前記第１の量および前記第２の量を使用してアップデート値を計算することは、前記第１の量の重みを前記第２の量で畳み込むことを含むことを特徴とする請求項３に記載のシステム。
- 前記アップデート値を使用して前記観察物体を表す前記物体空間画像を再構成することは、前記物体空間ボクセルのセットのデータを対応するアップデート値で更新することを含むことを特徴とする請求項１に記載のシステム。
- 前記観察物体を表す前記物体空間画像を受信する出力装置をさらに備え、
前記出力装置は、データ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項１に記載のシステム。 - 前記物体空間画像を受信する調査用コンピュータであって、前記物体空間画像からの診断情報の抽出を実施するように、あるいは前記１つもしくは複数の送信機、前記１つもしくは複数の検出装置、または前記画像再構成装置の少なくとも１つのためのパラメータを微調整するようにプログラムされた、調査用コンピュータをさらに備えることを特徴とする請求項１に記載のシステム。
- １つまたは複数のプロセッサと、
前記１つまたは複数のプロセッサによって実行されるソフトウェアを格納する１つまたは複数のデータ記憶装置と
を備えた、１つまたは複数のプログラム記憶および実行装置であって、前記プロセッサによって実行される時、前記ソフトウェアに含まれるソフトウェアコードは、前記プロセッサに
物体空間から投影空間に投影される物体を表す入力データを受信すること、
複数の物体空間ボクセルのそれぞれについて第１の量を計算することであって、前記第１の量は、前記入力データの値および予測投影データの対応する値の間の比率の対数を示す、こと、
複数の物体空間ボクセルのそれぞれについて第２の量を計算することであって、前記第２の量は、ノイズ分散からの寄与率と結合される少なくとも１つの反転インパルス応答に対応し、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応し、前記第２の量は、ノイズ分散からの寄与率、およびカルマンフィルタ利得、近似カルマンフィルタ利得、ウィーナーフィルタ利得または近似ウィーナーフィルタ利得の少なくとも一つを使用して計算される、こと、
複数の物体空間ボクセルのそれぞれについてアップデート値を計算することであって、それぞれのアップデート値は、対応する第１の量および対応する第２の量を使用して計算される、こと、
前記アップデート値を使用して前記観察物体を表す物体空間画像を再構成すること
を実行させることを特徴とするプログラム記憶および実行装置。 - 前記入力データは、投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから受信されることを特徴とする請求項８に記載のプログラム記憶および実行装置。
- 前記入力データは、データバス、ケーブル、有線ネットワーク、または無線ネットワークを介して受信されることを特徴とする請求項８に記載のプログラム記憶および実行装置。
- 前記物体空間画像を受信するための出力装置をさらに備え、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項８に記載のプログラム記憶および実行装置。
- 前記ソフトウェアは、前記第１の量の重みを前記第２の量で畳み込むためのソフトウェアコードをさらに含むことを特徴とする請求項８に記載のプログラム記憶および実行装置。
- １つまたは複数のデータ記憶装置に格納されたソフトウェアコードを実行する１つまたは複数のプロセッサによって実行される方法であって、
前記方法は、
投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから、物体空間から投影空間に投影される物体を表す入力データを受信するステップと、
複数の物体空間ボクセルのそれぞれについて第１の量を計算するステップであって、前記第１の量は、解像度グリッド上の前記入力データの値および前記解像度グリッド上の予測投影データの対応する値の間の比率の対数を示す、ステップと、
複数の物体空間ボクセルのそれぞれについて第２の量を計算するステップであって、前記第２の量は、ノイズ分散からの寄与率と結合される少なくとも１つの反転インパルス応答に対応し、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応し、前記第２の量は、ノイズ分散からの寄与率、およびカルマンフィルタ利得、近似カルマンフィルタ利得、ウィーナーフィルタ利得または近似ウィーナーフィルタ利得の少なくとも一つを使用して計算される、ステップと、
複数の物体空間ボクセルのそれぞれについてアップデート値を計算するステップであって、それぞれのアップデート値は、対応する第１の量および対応する第２の量を使用して計算され、前記第２の量は、前記第１の量と結合され、および対数的に変換された負でない予備物体の密度に寄与する、ステップと、
前記アップデート値を使用して前記観察物体を表す画像を再構成するステップと、
前記再構成画像を出力装置に出力するステップであって、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つである、ステップと
を含み、
前記入力データを受信するステップ、前記第１の量を計算するステップ、前記第２の量を計算するステップ、前記アップデート値を計算するステップ、前記画像を再構成するステップ、および前記再構成画像を出力するステップは、前記１つまたは複数のデータ記憶装置に格納された前記ソフトウェアコードに従って前記１つまたは複数のプロセッサによって実施されることを特徴とする方法。 - 前記解像度グリッドは、可変サイズの複数の解像度グリッドから取得されることを特徴とする請求項１３に記載の方法。
- 前記予測投影データは、前記物体空間ボクセルのセットに基づいて計算され、前記物体空間ボクセルのセットは、前記物体空間ボクセルのセットが投影空間に投影されると、可変サイズの前記複数の解像度グリッドをカバーすることを特徴とする請求項１４に記載の方法。
- 前記解像度グリッド上の前記入力データまたは前記物体空間ボクセルのセットの少なくとも１つを平滑カーネルで平滑化するステップをさらに含むことを特徴とする請求項１５に記載の方法。
- 前記物体空間ボクセルのセットの１つまたは複数の特徴は、値の１つまたは複数の範囲内に抑制されることを特徴とする請求項１５に記載の方法。
- 前記第１の量は、前記入力データと前記予測投影データとの間の誤差の測度として特徴づけられることを特徴とする請求項１３に記載の方法。
- 前記第２の量は、前記入力データの測定値から取得されるノイズ分散からの寄与率を含むことを特徴とする請求項１３に記載の方法。
- 前記第２の量は、前記物体空間内の単位強度の前記参照ボクセルの前記場所をカバーする周辺の量、または前記参照ボクセルの周辺の量を含む、前記物体空間内の単位強度の前記参照ボクセルの投影に関して計算されることを特徴とする請求項１３に記載の方法。
- 前記画像を再構成するステップは、
前記第１の量の重みを前記第２の量で畳み込んで、前記アップデート値を取得するステップを含むことを特徴とする請求項１３に記載の方法。 - 前記アップデート値は、限定関数によってさらに限定されることを特徴とする請求項１３に記載の方法。
- 対応するアップデート値を使用して前記物体空間ボクセルのセットの密度を更新するステップをさらに含むことを特徴とする請求項１３に記載の方法。
- 前記物体空間ボクセルのセットが更新された後、前記解像度グリッドより大きいサイズをもつ異なる解像度グリッドに変更するステップをさらに含むことを特徴とする請求項２３に記載の方法。
- １つまたは複数のプロセッサによって実行されるソフトウェアを格納する記憶装置であって、
前記ソフトウェアは、実行されると、前記１つまたは複数のプロセッサに、
投影画像取得システム、前記記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから、物体空間から投影空間に投影される物体を表す入力データを受信すること、
複数の物体空間ボクセルのそれぞれについて第１の量を計算することであって、前記第１の量は、解像度グリッド上の前記入力データの値および前記解像度グリッド上の予測投影データの対応する値の間の比率の対数を示す、こと、
複数の物体空間ボクセルのそれぞれについて第２の量を計算することであって、前記第２の量は、ノイズ分散からの寄与率と結合される少なくとも１つのインパルス応答に対応し、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応し、前記第２の量は、ノイズ分散からの寄与率、およびカルマンフィルタ利得、近似カルマンフィルタ利得、ウィーナーフィルタ利得または近似ウィーナーフィルタ利得の少なくとも一つを使用して計算され、前記第２の量を計算することは、
前記参照ボクセルの近傍について前記投影空間データを逆投影して１つまたは複数の逆投影された出力を生成し、および前記ノイズ分散からの寄与率と結合される前記１つまたは複数の逆投影された出力を反転すること、
信号対雑音比によって修正されるインパルス応答の投影を反転すること、
の少なくとも１つをさらに含む、こと、
複数の物体空間ボクセルのそれぞれについてアップデート値を計算することであって、それぞれのアップデート値は、対応する第１の量および対応する第２の量を使用して計算される、こと、
前記アップデート値を使用して前記観察物体を表す画像を再構成すること、
前記再構成画像を出力装置に出力することであって、前記出力装置は、前記記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つである、こと
を実行させることを特徴とする記憶装置。 - 前記投影空間データを受信し、および前記投影空間データを処理する画像再構成装置は、予測空間データおよび物体空間ボクセルデータのセットの少なくとも１つを平滑化することによって前記投影空間データを処理することをさらに含むことを特徴とする請求項１に記載のシステム。
- 前記ソフトウェアは、予測空間データおよび物体空間ボクセルデータのセットの少なくとも１つを平滑化するソフトウェアコードをさらに含むことを特徴とする請求項８に記載のプログラム記憶および実行装置。
- 予測空間データおよび物体空間ボクセルデータのセットの少なくとも１つ平滑化するステップをさらに含むことを特徴とする請求項１３に記載の方法。
- 前記ソフトウェアは、実行されると、前記１つまたは複数のプロセッサに、
予測空間データおよび物体空間ボクセルデータのセットの少なくとも１つ平滑化することをさらに実行させることを特徴とする請求項２５に記載の記憶装置。 - 前記第２の量を計算することは、
前記参照ボクセルの近傍について前記投影空間データを逆投影して１つまたは複数の逆投影された出力を生成し、および前記ノイズ分散からの寄与率と結合される前記１つまたは複数の逆投影された出力を反転すること、または、
信号対雑音比によって修正されるインパルス応答の投影を反転すること、
の少なくとも１つを含むことを特徴とする請求項１に記載のシステム。 - 前記第２の量を計算することは、
前記参照ボクセルの近傍について前記投影空間データを逆投影して１つまたは複数の逆投影された出力を生成し、および前記ノイズ分散からの寄与率と結合される前記１つまたは複数の逆投影された出力を反転すること、または、
信号対雑音比によって修正されるインパルス応答の投影を反転すること、
の少なくとも１つを含むことを特徴とする請求項８に記載のプログラム記憶および実行装置。 - 前記第２の量を計算することは、
前記参照ボクセルの近傍について前記投影空間データを逆投影して１つまたは複数の逆投影された出力を生成し、および前記ノイズ分散からの寄与率と結合される前記１つまたは複数の逆投影された出力を反転すること、または、
信号対雑音比によって修正されるインパルス応答の投影を反転すること、
の少なくとも１つを含むことを特徴とする請求項１３に記載の方法。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US12945908P | 2008-06-27 | 2008-06-27 | |
| US61/129,459 | 2008-06-27 | ||
| PCT/US2009/049102 WO2009158718A1 (en) | 2008-06-27 | 2009-06-29 | High efficiency computed tomography |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2011526367A JP2011526367A (ja) | 2011-10-06 |
| JP5543448B2 true JP5543448B2 (ja) | 2014-07-09 |
Family
ID=41444997
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2011516784A Active JP5543448B2 (ja) | 2008-06-27 | 2009-06-29 | 高効率コンピュータ断層撮影 |
Country Status (9)
| Country | Link |
|---|---|
| US (2) | US8660328B2 (ja) |
| EP (1) | EP2310840B1 (ja) |
| JP (1) | JP5543448B2 (ja) |
| KR (1) | KR20120138256A (ja) |
| CN (1) | CN102132149B (ja) |
| AU (1) | AU2009261945A1 (ja) |
| CA (1) | CA2729166A1 (ja) |
| IN (1) | IN2010DN09038A (ja) |
| WO (1) | WO2009158718A1 (ja) |
Families Citing this family (23)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050273009A1 (en) * | 2004-06-02 | 2005-12-08 | Harald Deischinger | Method and apparatus for co-display of inverse mode ultrasound images and histogram information |
| US8660328B2 (en) | 2008-06-27 | 2014-02-25 | Wolfram R. JARISCH | High efficiency computer tomography |
| US8660330B2 (en) * | 2008-06-27 | 2014-02-25 | Wolfram Jarisch | High efficiency computed tomography with optimized recursions |
| EP2385494A1 (en) * | 2010-05-07 | 2011-11-09 | IBBT vzw | A method and device for estimating noise in a reconstructed image |
| DE102010040041B3 (de) * | 2010-08-31 | 2012-01-26 | Siemens Aktiengesellschaft | Verfahren zur Korrektur von durch zeitliche Veränderungen von Schwächungswerten auftretenden Artefakten |
| CN102540276B (zh) | 2010-12-29 | 2014-11-26 | 通用电气公司 | 软场层析成像系统和方法 |
| JP2012189456A (ja) * | 2011-03-10 | 2012-10-04 | Ihi Corp | 潤滑剤分布取得装置及び潤滑剤分布取得方法 |
| US10186056B2 (en) | 2011-03-21 | 2019-01-22 | General Electric Company | System and method for estimating vascular flow using CT imaging |
| JP5818588B2 (ja) * | 2011-09-05 | 2015-11-18 | 株式会社東芝 | 放射線検出データ処理装置及び方法 |
| DE102011117654B4 (de) * | 2011-11-04 | 2013-09-05 | Eizo Gmbh | Verfahren zum Betreiben einer Bildverarbeitungseinrichtung sowie entsprechende Bildverarbeitungseinrichtung |
| US9677869B2 (en) | 2012-12-05 | 2017-06-13 | Perimeter Medical Imaging, Inc. | System and method for generating a wide-field OCT image of a portion of a sample |
| WO2014097065A1 (en) * | 2012-12-21 | 2014-06-26 | Koninklijke Philips N.V. | Image processing apparatus and method for filtering an image |
| US9415546B2 (en) * | 2014-01-29 | 2016-08-16 | Xerox Corporation | System and method for controlling material drop volume in three dimensional object printing |
| RU2720161C2 (ru) * | 2015-10-16 | 2020-04-24 | Эмтензор Гмбх | Электромагнитная томография с распознаванием картин интерференции |
| TWI550556B (zh) * | 2015-12-28 | 2016-09-21 | A Method of Establishing Stereo Image Model by Kalman Filtering | |
| US10650621B1 (en) | 2016-09-13 | 2020-05-12 | Iocurrents, Inc. | Interfacing with a vehicular controller area network |
| US10607378B2 (en) | 2016-11-18 | 2020-03-31 | Wolfram R. JARISCH | Extended high efficiency computed tomography with optimized recursions and applications |
| CA3044844A1 (en) | 2016-11-23 | 2018-05-31 | Emtensor Gmbh | Use of electromagnetic field for tomographic imaging of head |
| EP3413033B1 (en) * | 2017-06-09 | 2020-09-23 | Roche Diagnostics GmbH | Method and apparatus for determining properties of a laboratory sample contained in a laboratory sample container |
| WO2019014767A1 (en) | 2017-07-18 | 2019-01-24 | Perimeter Medical Imaging, Inc. | SAMPLE CONTAINER FOR STABILIZING AND ALIGNING EXCISED ORGANIC TISSUE SAMPLES FOR EX VIVO ANALYSIS |
| EP3814758A1 (en) | 2018-06-29 | 2021-05-05 | Universiteit Antwerpen | Item inspection by dynamic selection of projection angle |
| CN109091109B (zh) * | 2018-07-02 | 2021-04-20 | 南京大学 | 基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法 |
| CN110032758B (zh) * | 2019-02-28 | 2023-06-06 | 苏州瑞派宁科技有限公司 | 计算电信号的能量的方法、装置和计算机存储介质 |
Family Cites Families (23)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4593355A (en) * | 1983-11-21 | 1986-06-03 | American Science And Engineering, Inc. | Method of quick back projection for computed tomography and improved CT machine employing the method |
| US5566341A (en) * | 1992-10-05 | 1996-10-15 | The Regents Of The University Of California | Image matrix processor for fast multi-dimensional computations |
| US6018562A (en) * | 1995-11-13 | 2000-01-25 | The United States Of America As Represented By The Secretary Of The Army | Apparatus and method for automatic recognition of concealed objects using multiple energy computed tomography |
| US6110113A (en) | 1998-12-15 | 2000-08-29 | Siemens Medical Systems, Inc. | Method and apparatus for removing transients and gaps from ultrasound echo signals |
| US6345113B1 (en) | 1999-01-12 | 2002-02-05 | Analogic Corporation | Apparatus and method for processing object data in computed tomography data using object projections |
| 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 |
| US6859564B2 (en) * | 2001-02-15 | 2005-02-22 | James N. Caron | Signal processing using the self-deconvolving data reconstruction algorithm |
| US7825667B2 (en) * | 2003-04-04 | 2010-11-02 | Microwave Imaging Systems Technologies, Inc. | Microwave imaging system and processes, and associated software products |
| BRPI0511772A (pt) * | 2004-06-01 | 2008-01-08 | Exxonmobil Upstream Res Co | método para uso de uma abordagem de filtro de kalman para o processamento de dados eletromagnéticos |
| US7386088B2 (en) * | 2004-09-24 | 2008-06-10 | General Electric Company | Method and system for iterative image reconstruction |
| US7215732B2 (en) | 2004-09-30 | 2007-05-08 | General Electric Company | Method and system for CT reconstruction with pre-correction |
| US7251306B2 (en) * | 2004-11-17 | 2007-07-31 | General Electric Company | Methods, apparatus, and software to facilitate iterative reconstruction of images |
| US8687869B2 (en) | 2005-11-30 | 2014-04-01 | The Research Foundation Of State Of University Of New York | System and method for acceleration of image reconstruction |
| DE102006003609B4 (de) | 2006-01-25 | 2014-09-04 | Siemens Aktiengesellschaft | Tomographie-System und Verfahren zur Visualisierung einer tomographischen Darstellung |
| CN100591269C (zh) | 2006-02-17 | 2010-02-24 | 株式会社东芝 | 数据校正装置、数据校正方法、磁共振成像装置和x射线ct装置 |
| JP5248010B2 (ja) * | 2006-02-17 | 2013-07-31 | 株式会社東芝 | データ補正装置、データ補正方法、磁気共鳴イメージング装置およびｘ線ｃｔ装置 |
| CN101405619B (zh) | 2006-03-16 | 2013-01-02 | 皇家飞利浦电子股份有限公司 | 计算机断层造影数据采集装置和方法 |
| US8090182B2 (en) | 2006-11-13 | 2012-01-03 | National University Corporation Kyoto Institute Of Technology | Image reconstruction device, image reconstruction method, image reconstruction program, and CT apparatus |
| US8175115B2 (en) * | 2006-11-17 | 2012-05-08 | General Electric Company | Method and system for iterative reconstruction |
| US8194961B2 (en) | 2008-04-21 | 2012-06-05 | Kabushiki Kaisha Toshiba | Method, apparatus, and computer-readable medium for pre-reconstruction decomposition and calibration in dual energy computed tomography |
| US8660328B2 (en) | 2008-06-27 | 2014-02-25 | Wolfram R. JARISCH | High efficiency computer tomography |
| US8234824B2 (en) | 2008-06-27 | 2012-08-07 | Sunpower Corporation | Photovoltaic module with removable wind deflector |
| US7852977B2 (en) * | 2008-09-11 | 2010-12-14 | Samplify Systems, Inc. | Adaptive compression of computed tomography projection data |
-
2009
- 2009-06-29 US US12/493,989 patent/US8660328B2/en active Active
- 2009-06-29 US US13/001,592 patent/US20110158492A1/en not_active Abandoned
- 2009-06-29 AU AU2009261945A patent/AU2009261945A1/en not_active Abandoned
- 2009-06-29 WO PCT/US2009/049102 patent/WO2009158718A1/en active Application Filing
- 2009-06-29 KR KR1020107028604A patent/KR20120138256A/ko not_active Application Discontinuation
- 2009-06-29 EP EP09771249.1A patent/EP2310840B1/en active Active
- 2009-06-29 CA CA2729166A patent/CA2729166A1/en not_active Abandoned
- 2009-06-29 JP JP2011516784A patent/JP5543448B2/ja active Active
- 2009-06-29 CN CN200980123522.8A patent/CN102132149B/zh active Active
-
2010
- 2010-12-17 IN IN9038DEN2010 patent/IN2010DN09038A/en unknown
Also Published As
| Publication number | Publication date |
|---|---|
| AU2009261945A2 (en) | 2011-04-21 |
| US20110158492A1 (en) | 2011-06-30 |
| US20090324047A1 (en) | 2009-12-31 |
| US8660328B2 (en) | 2014-02-25 |
| EP2310840B1 (en) | 2014-03-19 |
| WO2009158718A1 (en) | 2009-12-30 |
| AU2009261945A1 (en) | 2009-12-30 |
| EP2310840A1 (en) | 2011-04-20 |
| EP2310840A4 (en) | 2012-10-03 |
| KR20120138256A (ko) | 2012-12-26 |
| CA2729166A1 (en) | 2009-12-30 |
| CN102132149A (zh) | 2011-07-20 |
| CN102132149B (zh) | 2014-05-28 |
| JP2011526367A (ja) | 2011-10-06 |
| IN2010DN09038A (ja) | 2015-07-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP5543448B2 (ja) | 高効率コンピュータ断層撮影 | |
| TWI555514B (zh) | 用於實施具遞迴最佳化之高效率電腦斷層掃描之系統、工作站及方法 | |
| JP6280700B2 (ja) | 反復式再構成の方法、非一時的なコンピュータ可読の媒体及びイメージング・システム | |
| Sidky et al. | Constrained ${\rm T} p {\rm V} $ Minimization for Enhanced Exploitation of Gradient Sparsity: Application to CT Image Reconstruction | |
| JP5530637B2 (ja) | 画像再構成の方法及びシステム | |
| US8571287B2 (en) | System and method for iterative image reconstruction | |
| US8971599B2 (en) | Tomographic iterative reconstruction | |
| JP6212294B2 (ja) | 反復式再構成の方法及び装置 | |
| US8897528B2 (en) | System and method for iterative image reconstruction | |
| US9020230B2 (en) | Method and apparatus for iterative reconstruction | |
| US20170340287A1 (en) | Method And Apparatus For Motion Correction In CT Imaging | |
| Wang et al. | Sparse-view cone-beam CT reconstruction by bar-by-bar neural FDK algorithm | |
| Sidky et al. | Accurate image reconstruction in circular cone-beam computed tomography by total variation minimization: a preliminary investigation | |
| Pohlmann et al. | Estimation of missing fan-beam projections using frequency consistency conditions | |
| Karimi et al. | A hybrid stochastic-deterministic gradient descent algorithm for image reconstruction in cone-beam computed tomography | |
| Karimi et al. | Sparse-view image reconstruction in cone-beam computed tomography with variance-reduced stochastic gradient descent and locally-adaptive proximal operation | |
| Gao et al. | Distributed computation of linear inverse problems with application to computed tomography | |
| Hassan et al. | Study of compressed sensing application to low dose computed tomography data collection | |
| Cho | Improving statistical image reconstruction for cardiac X-ray computed tomography | |
| Colmeiro et al. | Reconstruction of positron emission tomography images using adaptive sliced remeshing strategy | |
| CN114581545A (zh) | 医学图像重建方法、装置、计算机设备及存储介质 | |
| BOONE et al. | Constrained TpV Minimization for Enhanced Exploitation of Gradient Sparsity: Application to CT Image Reconstruction | |
| Niinimaki et al. | Multiresolution local tomography in dental radiology using wavelets | |
| Wang et al. | System modeling studies in iterative X-ray CT reconstruction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20120216 |
|
| A601 | Written request for extension of time |
Free format text: JAPANESE INTERMEDIATE CODE: A601 Effective date: 20130911 |
|
| A602 | Written permission of extension of time |
Free format text: JAPANESE INTERMEDIATE CODE: A602 Effective date: 20130919 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20131011 |
|
| 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: 20140408 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20140508 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 5543448 Country of ref document: JP 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 |