高効率コンピュータ断層撮影
システムは、観察物体内に励起エネルギーを送信するための送信機と、物体への送信励起エネルギーに応答して受信されるエネルギーをエンコードする投影空間データを生成するための検出装置と、投影空間データを生成するための送信機および受信機を制御するためのコントローラと、投影空間データを受信するための画像再構成装置であって、投影空間データと予測投影データの間の差を特徴づける第１の量を計算し、各インパルス応答が単位強度の参照ボクセルに対応する少なくとも１つのインパルス応答に対応する第２の量を計算し、第１の量および第２の量の反転関数を使ってアップデート値を計算し、アップデート値を使って観察物体を表す物体空間画像を再構成することによって投影空間データを処理するための画像再構成装置とを含む。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像再構成（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）アップデート値ｉを得ることができる。例えば、第１の量でＰを畳み込んで、生アップデート値を得ることができる。Ｐがｎ次元再構成空間内でゆっくり変化するのに従って、適切だが十分に小さい所与のボクセル近傍にＰを近似として印加することによって、計算時間を節約することができる。上述した式（１）、（２）のｆと同様に構築される限定関数ｇを、生アップデート値ｉに適用して、修正アップデート値ｉｇ＝ｇ（ｉ）を得ることができる。この修正は過度の訂正を防ぎ得る。本発明の例示的な実施形態を開発する間、生アップデート値を直接使用する多くの実験が失敗した。本発明者は、過度のアップデート値を制限するための限定関数およびその賢明な使用を思いつき、開発した。修正アップデートは、後続処理のためのアップデート値となり得る。
【００５９】
ボックス３０４〜３０８で、単純化をさらに検討することが、円錐ビーム幾何構造の再構成（例えば、血管造影法で使われる）にとって有用な場合があり、ここで投影の数ｐは、小さくてよい（ｐ＜＜Ｎであり、Ｎは、再構成ボリュームを通る代表的投影パスに沿ったボクセルの数を表し得る）。この場合、個々の「蜘蛛の脚」は、特に多重再構成グリッドの高解像度最終ステージにおいては相異なってよい。この高解像度ｎ次元空間内で多くの小領域に対するインパルス応答Ｒｐを計算するのは、時間がかかり得る。ｎ次元インパルス応答使用の近似として、次元ｍｉの比較的低次元で投影されたｐ個のインパルス応答（０＜−ｉ＜ｐ）を使うことができる。対応する修正された逆またはカルマン／ウィーナー利得を、投影の個々の誤差画像に繰り返し適用してよい。こうしたｍ次元インパルス応答は、蜘蛛状でよい上記のｎ次元インパルス応答の、それぞれのｍ次元投影空間への投影でよい。こうしたインパルス応答のプロパティは上述した通りでよく、すなわち、投影空間内でゆっくり変化する。ただし、ｐ個の投影のインパルス応答、こうした計算されたインパルス応答の逆、および後続畳込み（ここではｍ次元投影に適用される）を計算するための計算負荷は、特にｐ＜＜Ｎのとき、ｎ次元空間内での計算よりはるかに小さくなり得る。計算されたインパルス応答は、Ｒの場合と同様にして修正することができる。計算されたインパルス応答は、トモグラフィ再構成におけるノイズを削減するように帯域通過フィルタリングすることもできる。この断層反転技法は、投影自体が再構成に先立ってフィルタリングされるフィルタ補正逆投影と同様でよい。この技法は、円錐ビーム幾何構造に限定されなくてよく、例えば、並列ビーム幾何構造など、他の幾何構造にも適用することができる。
【００６０】
ボックス３０９で、修正アップデート値ｉｇを、ｎ次元ボクセルセットｖ中の対応する画像密度の値に加えて、ｎ次元ボクセルセットｖを更新することができる。加算は、上述した式（１）、（２）と同様に、限定関数によって境界を定められ得る。加算は、ｎ次元ボクセルセット全体の一部分に対して実施することができる。精度のために、ボクセル値ｖは、その投影が生画像データと最も合致するようにスケーリングしてよい。
【００６１】
ボックス３１０で、解像度グリッドが最も細かい場合、また、計算された予測と入力投影データとの間の差（例えば、残差）が（統計的に）十分に小さい場合、密度推定値を表す更新ｎ次元ボクセルセットｖ、ｖＬ、または両方（並びに必要に応じて推定予測および誤差）を、ボックス３１１に示すように出力装置に出力することができる。出力に先立って平滑化を適用することができる。それ以外の場合、解像度グリッドは、現在の解像度グリッドより大きいサイズの解像度グリッドに変更することができ（または可能性としては、ｎ次元平滑カーネルの重大度／範囲を低下させることができ）、ボックス３１２で、フローはボックス３０３に進み得る。
【００６２】
ボックス３１１で終了すると、物体１０６の高忠実度の再構成画像を取得することができる。この再構成画像１０６は、物体２０３の内部情報を非侵襲的および非破壊的に明らかにすることができ、診断用医用撮像および非破壊的評価において広範に適用される。再構成画像１０６は、出力装置１０７に関して上述したように、可視化装置上で可視化することも、コンピュータメモリまたは記憶媒体に格納することもできる。
【００６３】
例えば、非特許文献５にある投影ピクセルデータ１０４を使って、本発明の例示的な実施形態により達成可能な改良を証明した。図４Ａは、投影ピクセルデータ１０４のフィルタ補正逆投影（ＦＢＰ）再構成に基づく再構成画像１０６を示す。図４Ｂは、本発明の例示的な実施形態による再構成画像１０６を示す。図４Ａと４Ｂとの間の比較は、例えば、細胞内対細胞外空間、細胞膜、リボソーム、スピロプラズマの細部など、様々な密度のより明確な対比によって測定される、図４Ｂにおける優れた画像鮮明度を証明している。図４Ｂの画像はユーザにとってより望ましい。両方の再構成は、ほぼ＋／−７０度の範囲に渡るほぼ１６６個の二軸傾斜式の透過型電子顕微鏡法（ＴＥＭ）投影を用いた。
【００６４】
図５は、模擬ファントムの正面像サンプル投影を示す。模擬ファントムは、サイズが２５６×２５６×６４ボクセルの３次元ボックス内のランダムの場所に、サイズが６×６×６ボクセルの１００個程度の立方体を含む。こうしたボクセルは、０．８５単位の周辺密度に組み込まれた５単位の密度をもつ。白色ノイズが加算された合計で４９個の投影（単軸傾斜、＋／−７２度）が、コンピュータシミュレーション１０３によって取得されて、投影ピクセルデータ１０４を生成した。
【００６５】
統計モデル検証における公知のゴールドスタンダードによると、再構成品質は、残余検査、すなわち、実像（例えば、予め分かっている図５の模擬ファントム）と再構成画像１０６との間の不一致によって最も厳しく評価され得る。高忠実再構成は、ランダムであるとともに３次元物体（例えば、図５の模擬ファントム内の立方体）の構造からの非常に少ない寄与率を含む残余につながり得る。図６Ａは、同時反復再構成法（ＳＩＲＴ）の残余を示し、図６Ｂは、本発明の例示的な実施形態による方法から取得されるものを示す。図６Ａのかなり著しい残差平方パターン（１５回の反復の後）および図６Ｂの顕著な白色ノイズ出現（１度の高解像度反復の後）は、画像忠実度および計算効率における本発明の優越性をさらに証明する。ＳＩＲＴのための反復カウントは、ＳＩＲＴのためのそれより多いまたはより少ない反復が画像品質の劣化につながっただけであったので、この例のために最適化してある。
【００６６】
本発明は、投影の数が小さい状況、例えば、血管造影法における円錐ビーム幾何構造にも当てはまり得る。図７Ａは、プラスチックチューブで作られた鼓動式水圧冠動脈ファントムにヨード造影剤（通例、例えば、介入的心臓サブトラクション血管造影法で利用される）を注入している間の１枚のＸ線投影図を示す。用いられた投影ピクセルデータ１０４は、注入の一過性のため、基底ファントムのわずか６度のデジタルサブトラクション投影によるものであった。図７Ｂは、本発明の例示的な実施形態による再構成画像１０６を示し、鼓動式水圧ファントムの幾何構造を正確に表している。再構成画像は、フィルタバック投影（ＦＢＰ）によっても同時反復再構成法（ＳＩＲＴ）によっても、最大エントロピー（ＭＥＮＴ）法によってさえも取得することはできない。というのは、投影回数が異常に小さく、かつ／または反復回数が現在の計算機資源では実際的でなくなるからである。
【００６７】
上述したように、本発明の例示的ないくつかの実施形態では、計算／処理技法を合成し、結果として、対象領域内で所望のｎ次元ボクセルの高速で正確な計算／推定をもたらすとともに、任意のｍ次元投影（および方向）を使用できるようにする。ｍは固定しなくてよく、投影指標に依存し得ることに留意されたい。合成技法の効率は、広範な投影回数に当てはまり得る。入手可能な投影（測定）情報を効率的に使用することにより、投影回数が非常に大きいとき、および冠動脈血管造影法などにおいて制限される（小さいまたは不完全）ときの両方で本方法は効果的になり得る。
【００６８】
本発明の非線形反転技法は、（ｉ）再構成の領域内部の２つの非線形関連ドメインへのｎ次元画像推定値の分解による、非線形カルマンフィルタの一定部分の線形化（ベイズの推定の近似）、（ｉｉ）画像平滑化の低下度を用いることによって解像度をさらに細分した「階段」に匹敵する（「ステップ」に匹敵する）マルチグリッド計算、および（ｉｉｉ）必要とされるとき、逆ヘシアンを使う代わりに、画像ノイズの推定値によって修正された線形システムの逆インパルス応答（または伝達関数）の効率的計算であって、物体空間における選択領域に対する近似カルマン／ウィーナーフィルタ利得を生じる計算を合成している。
【００６９】
具体的には、非線形カルマンフィルタは、例えば、非特許文献６で、画像再構成のためのカルマンフィルタリングの適用は、例えば、非特許文献７で論じられている。例えば、ＦＢＰ、ＦＦＴ、およびウィーナーフィルタなどの代替方法を、反転の線形部分の代わりに使うことができる。再構成の領域内部の２つの非線形関連ドメインのどちらも、カルマンシステム状態を表し得る。
【００７０】
上記合成によって、最適化におけるヘシアンの効果的使用と同様の効率的な誤差削減を、すべての反復で達成することができる。達成された効率は、二次収束のものに近くなり得る。この手法により、例えば、非特許文献８にある手法など、他の手法によってとられる手法で見られる高次元であり乏しい収束率の重大な問題が避けられ得る。
【００７１】
あるいは、本発明の例示的な実施形態は、ＩＲと、本明細書ではＳ−ＦＢＰと呼ばれる（空間可変／非定常／分割）ＦＢＰとの合成と見なすことができる。近似の要件に依存して、Ｓ−ＦＢＰに関連づけられた１つまたはいくつかのセグメントがあり得る。概して、ＦＢＰは、再構成のフィールドにおいて定常ではない。ただし、フィールドプロパティをゆっくりと変更することにより、再構成のフィールドを、領域の間の遷移が重みづけされた重なり領域／セグメントに分けることによって近似が可能になり得る。Ｓ−ＦＢＰの速度、効率、および柔軟性は、後述するいくつかの実施形態の適切な選択および合成により達成することができる。
【００７２】
例えば、Ｓ−ＦＢＰ部分は、生投影データではなく、非線形変換されたデータに適用することができる。
【００７３】
例えば、マルチグリッド計算は、線形化をサポートし、各反復での演算カウントを非常に低く保つことができる。例えば、最高解像度グリッドでの最終計算に先立つ合計演算カウントは、最高解像度グリッドでのＳ−ＦＢＰのものよりも低くなり得る。
【００７４】
例えば、ボクセル密度は、２つの非線形関連ドメイン（マップにより接続される）、すなわちＳ−ＦＢＰによるアップデート用ドメイン（所望の近似技法に依存して、ｎ次元空間内において可変でよい）、および投影予測の計算用ドメインで表すことができる。Ｓ−ＦＢＰ部分のフィルタ重みは、ボクセルの場所に依存し得るが、同じ重みをもって使用するためのボクセルのグループ化により、計算効率が向上し得る。
【００７５】
例えば、線形化システムインパルス応答（または伝達関数）の逆から、Ｓ−ＦＢＰ部分のフィルタの重みを計算することによって、任意の投影次元および方向を用いることができる。線形化システムインパルス応答（または伝達関数）は、入力投影データの測定値から判定されたノイズ分散と合成済みでよい。
【００７６】
例えば、適切に選ばれた限定（トリミング）関数は、計算中の数値特異性および不安定性を防ぎ得る。
【００７７】
例えば、各反復においてボクセルおよびピクセル値に慎重に適用される、適切に選ばれた平滑化関数は、数値計算の安定性および精度に寄与し得る。このボクセルおよびピクセル値の平滑化は、アップデート値の生成に先立つカルマン予測子アルゴリズムに機能的に匹敵し得る。平滑化の反復適用は、カルマンフィルタにおける状態予測子アルゴリズムと同様でよい。
【００７８】
最高解像度グリッド上での計算に先行する反復の精度が十分であり得るので、最高解像度グリッドを通る反復は計算上必要なくてよく、従って、誤差訂正項の計算が節約され得る。こうした事例では、再構成速度は、Ｓ−ＦＢＰのものに匹敵するようになり得る。というのは、最高解像度グリッド上での計算により、ＦＢＰのものに匹敵する演算カウント（１より大きい因数の範囲内）が達成され、より粗い解像度グリッド上での先行するより低い解像度の反復すべてに対する合成演算カウントは、特により高次元の画像再構成に対しては、最高解像度グリッド上での計算の場合より少なくなり得るからであり得る。
【００７９】
図３に示す方法を参照して論じた本発明の例示的な実施形態は、例えば、ＣＤ−ＲＯＭ、ＤＶＤ、ＭＯ（光磁気）ディスク、ハードディスク、フロッピー（登録商標）ディスク、ジップディスク、フラッシュドライブなどのデータ記憶装置に格納されるソフトウェアコードとして提供され得る。格納されたソフトウェアコードは、図３を参照して上述した例示的な方法を実施するための、１つまたは複数のプロセッサと、例えば、ＲＡＭ（ランダムアクセスメモリ）装置、ＤＲＡＭ（ダイナミックＲＡＭ）装置、フラッシュメモリ装置、およびＳＲＡＭ（スタティックＲＡＭ）装置などの、１つまたは複数のメモリ装置とを有するコンピュータによって可読であり実行可能でよい。
【００８０】
本発明の例示的な実施形態は、図３を参照して上述した例示的な方法を格納し実行する１つまたは複数のプログラム記憶および実行装置、例えば、ＡＳＩＣ（特定用途向け集積回路）、ＦＰＧＡ（フィールドプログラム可能ゲートアレイ）、ＣＰＬＤ（複合プログラム可能論理素子）、ＡＳＩＰ（特定用途命令セットプロセッサ）などを提供し得る。
【００８１】
本明細書に記載した例および実施形態は、非限定的な例である。本発明を、例示的な実施形態を参照して詳しく記載したので、上記内容から、本発明からより広範な態様において逸脱することなく変更および修正を行ってよく、したがって、請求項において定義される本発明は、このような本発明の真の精神の範囲内であるあらゆる変更および修正をカバーすることを意図していることが当業者には明らかであろう。
【技術分野】
【０００１】
本発明は、画像再構成（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）アップデート値ｉを得ることができる。例えば、第１の量でＰを畳み込んで、生アップデート値を得ることができる。Ｐがｎ次元再構成空間内でゆっくり変化するのに従って、適切だが十分に小さい所与のボクセル近傍にＰを近似として印加することによって、計算時間を節約することができる。上述した式（１）、（２）のｆと同様に構築される限定関数ｇを、生アップデート値ｉに適用して、修正アップデート値ｉｇ＝ｇ（ｉ）を得ることができる。この修正は過度の訂正を防ぎ得る。本発明の例示的な実施形態を開発する間、生アップデート値を直接使用する多くの実験が失敗した。本発明者は、過度のアップデート値を制限するための限定関数およびその賢明な使用を思いつき、開発した。修正アップデートは、後続処理のためのアップデート値となり得る。
【００５９】
ボックス３０４〜３０８で、単純化をさらに検討することが、円錐ビーム幾何構造の再構成（例えば、血管造影法で使われる）にとって有用な場合があり、ここで投影の数ｐは、小さくてよい（ｐ＜＜Ｎであり、Ｎは、再構成ボリュームを通る代表的投影パスに沿ったボクセルの数を表し得る）。この場合、個々の「蜘蛛の脚」は、特に多重再構成グリッドの高解像度最終ステージにおいては相異なってよい。この高解像度ｎ次元空間内で多くの小領域に対するインパルス応答Ｒｐを計算するのは、時間がかかり得る。ｎ次元インパルス応答使用の近似として、次元ｍｉの比較的低次元で投影されたｐ個のインパルス応答（０＜−ｉ＜ｐ）を使うことができる。対応する修正された逆またはカルマン／ウィーナー利得を、投影の個々の誤差画像に繰り返し適用してよい。こうしたｍ次元インパルス応答は、蜘蛛状でよい上記のｎ次元インパルス応答の、それぞれのｍ次元投影空間への投影でよい。こうしたインパルス応答のプロパティは上述した通りでよく、すなわち、投影空間内でゆっくり変化する。ただし、ｐ個の投影のインパルス応答、こうした計算されたインパルス応答の逆、および後続畳込み（ここではｍ次元投影に適用される）を計算するための計算負荷は、特にｐ＜＜Ｎのとき、ｎ次元空間内での計算よりはるかに小さくなり得る。計算されたインパルス応答は、Ｒの場合と同様にして修正することができる。計算されたインパルス応答は、トモグラフィ再構成におけるノイズを削減するように帯域通過フィルタリングすることもできる。この断層反転技法は、投影自体が再構成に先立ってフィルタリングされるフィルタ補正逆投影と同様でよい。この技法は、円錐ビーム幾何構造に限定されなくてよく、例えば、並列ビーム幾何構造など、他の幾何構造にも適用することができる。
【００６０】
ボックス３０９で、修正アップデート値ｉｇを、ｎ次元ボクセルセットｖ中の対応する画像密度の値に加えて、ｎ次元ボクセルセットｖを更新することができる。加算は、上述した式（１）、（２）と同様に、限定関数によって境界を定められ得る。加算は、ｎ次元ボクセルセット全体の一部分に対して実施することができる。精度のために、ボクセル値ｖは、その投影が生画像データと最も合致するようにスケーリングしてよい。
【００６１】
ボックス３１０で、解像度グリッドが最も細かい場合、また、計算された予測と入力投影データとの間の差（例えば、残差）が（統計的に）十分に小さい場合、密度推定値を表す更新ｎ次元ボクセルセットｖ、ｖＬ、または両方（並びに必要に応じて推定予測および誤差）を、ボックス３１１に示すように出力装置に出力することができる。出力に先立って平滑化を適用することができる。それ以外の場合、解像度グリッドは、現在の解像度グリッドより大きいサイズの解像度グリッドに変更することができ（または可能性としては、ｎ次元平滑カーネルの重大度／範囲を低下させることができ）、ボックス３１２で、フローはボックス３０３に進み得る。
【００６２】
ボックス３１１で終了すると、物体１０６の高忠実度の再構成画像を取得することができる。この再構成画像１０６は、物体２０３の内部情報を非侵襲的および非破壊的に明らかにすることができ、診断用医用撮像および非破壊的評価において広範に適用される。再構成画像１０６は、出力装置１０７に関して上述したように、可視化装置上で可視化することも、コンピュータメモリまたは記憶媒体に格納することもできる。
【００６３】
例えば、非特許文献５にある投影ピクセルデータ１０４を使って、本発明の例示的な実施形態により達成可能な改良を証明した。図４Ａは、投影ピクセルデータ１０４のフィルタ補正逆投影（ＦＢＰ）再構成に基づく再構成画像１０６を示す。図４Ｂは、本発明の例示的な実施形態による再構成画像１０６を示す。図４Ａと４Ｂとの間の比較は、例えば、細胞内対細胞外空間、細胞膜、リボソーム、スピロプラズマの細部など、様々な密度のより明確な対比によって測定される、図４Ｂにおける優れた画像鮮明度を証明している。図４Ｂの画像はユーザにとってより望ましい。両方の再構成は、ほぼ＋／−７０度の範囲に渡るほぼ１６６個の二軸傾斜式の透過型電子顕微鏡法（ＴＥＭ）投影を用いた。
【００６４】
図５は、模擬ファントムの正面像サンプル投影を示す。模擬ファントムは、サイズが２５６×２５６×６４ボクセルの３次元ボックス内のランダムの場所に、サイズが６×６×６ボクセルの１００個程度の立方体を含む。こうしたボクセルは、０．８５単位の周辺密度に組み込まれた５単位の密度をもつ。白色ノイズが加算された合計で４９個の投影（単軸傾斜、＋／−７２度）が、コンピュータシミュレーション１０３によって取得されて、投影ピクセルデータ１０４を生成した。
【００６５】
統計モデル検証における公知のゴールドスタンダードによると、再構成品質は、残余検査、すなわち、実像（例えば、予め分かっている図５の模擬ファントム）と再構成画像１０６との間の不一致によって最も厳しく評価され得る。高忠実再構成は、ランダムであるとともに３次元物体（例えば、図５の模擬ファントム内の立方体）の構造からの非常に少ない寄与率を含む残余につながり得る。図６Ａは、同時反復再構成法（ＳＩＲＴ）の残余を示し、図６Ｂは、本発明の例示的な実施形態による方法から取得されるものを示す。図６Ａのかなり著しい残差平方パターン（１５回の反復の後）および図６Ｂの顕著な白色ノイズ出現（１度の高解像度反復の後）は、画像忠実度および計算効率における本発明の優越性をさらに証明する。ＳＩＲＴのための反復カウントは、ＳＩＲＴのためのそれより多いまたはより少ない反復が画像品質の劣化につながっただけであったので、この例のために最適化してある。
【００６６】
本発明は、投影の数が小さい状況、例えば、血管造影法における円錐ビーム幾何構造にも当てはまり得る。図７Ａは、プラスチックチューブで作られた鼓動式水圧冠動脈ファントムにヨード造影剤（通例、例えば、介入的心臓サブトラクション血管造影法で利用される）を注入している間の１枚のＸ線投影図を示す。用いられた投影ピクセルデータ１０４は、注入の一過性のため、基底ファントムのわずか６度のデジタルサブトラクション投影によるものであった。図７Ｂは、本発明の例示的な実施形態による再構成画像１０６を示し、鼓動式水圧ファントムの幾何構造を正確に表している。再構成画像は、フィルタバック投影（ＦＢＰ）によっても同時反復再構成法（ＳＩＲＴ）によっても、最大エントロピー（ＭＥＮＴ）法によってさえも取得することはできない。というのは、投影回数が異常に小さく、かつ／または反復回数が現在の計算機資源では実際的でなくなるからである。
【００６７】
上述したように、本発明の例示的ないくつかの実施形態では、計算／処理技法を合成し、結果として、対象領域内で所望のｎ次元ボクセルの高速で正確な計算／推定をもたらすとともに、任意のｍ次元投影（および方向）を使用できるようにする。ｍは固定しなくてよく、投影指標に依存し得ることに留意されたい。合成技法の効率は、広範な投影回数に当てはまり得る。入手可能な投影（測定）情報を効率的に使用することにより、投影回数が非常に大きいとき、および冠動脈血管造影法などにおいて制限される（小さいまたは不完全）ときの両方で本方法は効果的になり得る。
【００６８】
本発明の非線形反転技法は、（ｉ）再構成の領域内部の２つの非線形関連ドメインへのｎ次元画像推定値の分解による、非線形カルマンフィルタの一定部分の線形化（ベイズの推定の近似）、（ｉｉ）画像平滑化の低下度を用いることによって解像度をさらに細分した「階段」に匹敵する（「ステップ」に匹敵する）マルチグリッド計算、および（ｉｉｉ）必要とされるとき、逆ヘシアンを使う代わりに、画像ノイズの推定値によって修正された線形システムの逆インパルス応答（または伝達関数）の効率的計算であって、物体空間における選択領域に対する近似カルマン／ウィーナーフィルタ利得を生じる計算を合成している。
【００６９】
具体的には、非線形カルマンフィルタは、例えば、非特許文献６で、画像再構成のためのカルマンフィルタリングの適用は、例えば、非特許文献７で論じられている。例えば、ＦＢＰ、ＦＦＴ、およびウィーナーフィルタなどの代替方法を、反転の線形部分の代わりに使うことができる。再構成の領域内部の２つの非線形関連ドメインのどちらも、カルマンシステム状態を表し得る。
【００７０】
上記合成によって、最適化におけるヘシアンの効果的使用と同様の効率的な誤差削減を、すべての反復で達成することができる。達成された効率は、二次収束のものに近くなり得る。この手法により、例えば、非特許文献８にある手法など、他の手法によってとられる手法で見られる高次元であり乏しい収束率の重大な問題が避けられ得る。
【００７１】
あるいは、本発明の例示的な実施形態は、ＩＲと、本明細書ではＳ−ＦＢＰと呼ばれる（空間可変／非定常／分割）ＦＢＰとの合成と見なすことができる。近似の要件に依存して、Ｓ−ＦＢＰに関連づけられた１つまたはいくつかのセグメントがあり得る。概して、ＦＢＰは、再構成のフィールドにおいて定常ではない。ただし、フィールドプロパティをゆっくりと変更することにより、再構成のフィールドを、領域の間の遷移が重みづけされた重なり領域／セグメントに分けることによって近似が可能になり得る。Ｓ−ＦＢＰの速度、効率、および柔軟性は、後述するいくつかの実施形態の適切な選択および合成により達成することができる。
【００７２】
例えば、Ｓ−ＦＢＰ部分は、生投影データではなく、非線形変換されたデータに適用することができる。
【００７３】
例えば、マルチグリッド計算は、線形化をサポートし、各反復での演算カウントを非常に低く保つことができる。例えば、最高解像度グリッドでの最終計算に先立つ合計演算カウントは、最高解像度グリッドでのＳ−ＦＢＰのものよりも低くなり得る。
【００７４】
例えば、ボクセル密度は、２つの非線形関連ドメイン（マップにより接続される）、すなわちＳ−ＦＢＰによるアップデート用ドメイン（所望の近似技法に依存して、ｎ次元空間内において可変でよい）、および投影予測の計算用ドメインで表すことができる。Ｓ−ＦＢＰ部分のフィルタ重みは、ボクセルの場所に依存し得るが、同じ重みをもって使用するためのボクセルのグループ化により、計算効率が向上し得る。
【００７５】
例えば、線形化システムインパルス応答（または伝達関数）の逆から、Ｓ−ＦＢＰ部分のフィルタの重みを計算することによって、任意の投影次元および方向を用いることができる。線形化システムインパルス応答（または伝達関数）は、入力投影データの測定値から判定されたノイズ分散と合成済みでよい。
【００７６】
例えば、適切に選ばれた限定（トリミング）関数は、計算中の数値特異性および不安定性を防ぎ得る。
【００７７】
例えば、各反復においてボクセルおよびピクセル値に慎重に適用される、適切に選ばれた平滑化関数は、数値計算の安定性および精度に寄与し得る。このボクセルおよびピクセル値の平滑化は、アップデート値の生成に先立つカルマン予測子アルゴリズムに機能的に匹敵し得る。平滑化の反復適用は、カルマンフィルタにおける状態予測子アルゴリズムと同様でよい。
【００７８】
最高解像度グリッド上での計算に先行する反復の精度が十分であり得るので、最高解像度グリッドを通る反復は計算上必要なくてよく、従って、誤差訂正項の計算が節約され得る。こうした事例では、再構成速度は、Ｓ−ＦＢＰのものに匹敵するようになり得る。というのは、最高解像度グリッド上での計算により、ＦＢＰのものに匹敵する演算カウント（１より大きい因数の範囲内）が達成され、より粗い解像度グリッド上での先行するより低い解像度の反復すべてに対する合成演算カウントは、特により高次元の画像再構成に対しては、最高解像度グリッド上での計算の場合より少なくなり得るからであり得る。
【００７９】
図３に示す方法を参照して論じた本発明の例示的な実施形態は、例えば、ＣＤ−ＲＯＭ、ＤＶＤ、ＭＯ（光磁気）ディスク、ハードディスク、フロッピー（登録商標）ディスク、ジップディスク、フラッシュドライブなどのデータ記憶装置に格納されるソフトウェアコードとして提供され得る。格納されたソフトウェアコードは、図３を参照して上述した例示的な方法を実施するための、１つまたは複数のプロセッサと、例えば、ＲＡＭ（ランダムアクセスメモリ）装置、ＤＲＡＭ（ダイナミックＲＡＭ）装置、フラッシュメモリ装置、およびＳＲＡＭ（スタティックＲＡＭ）装置などの、１つまたは複数のメモリ装置とを有するコンピュータによって可読であり実行可能でよい。
【００８０】
本発明の例示的な実施形態は、図３を参照して上述した例示的な方法を格納し実行する１つまたは複数のプログラム記憶および実行装置、例えば、ＡＳＩＣ（特定用途向け集積回路）、ＦＰＧＡ（フィールドプログラム可能ゲートアレイ）、ＣＰＬＤ（複合プログラム可能論理素子）、ＡＳＩＰ（特定用途命令セットプロセッサ）などを提供し得る。
【００８１】
本明細書に記載した例および実施形態は、非限定的な例である。本発明を、例示的な実施形態を参照して詳しく記載したので、上記内容から、本発明からより広範な態様において逸脱することなく変更および修正を行ってよく、したがって、請求項において定義される本発明は、このような本発明の真の精神の範囲内であるあらゆる変更および修正をカバーすることを意図していることが当業者には明らかであろう。
【特許請求の範囲】
【請求項１】
観察物体内に励起エネルギーを送信するための１つまたは複数の送信機と、
前記観察物体への前記送信励起エネルギーに応答して１つまたは複数の検出装置によって受信されるエネルギーをエンコードする投影空間データを生成するための前記１つまたは複数の検出装置と、
前記励起エネルギーを送信するための前記１つまたは複数の送信機および前記投影空間データを生成するための１つまたは複数の受信機を制御するためのコントローラと、
前記投影空間データと予測投影データの間の差を特徴づける第１の量を計算し、
各インパルス応答が単位強度の参照ボクセルに対応する少なくとも１つのインパルス応答に対応する第２の量を計算し、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算し、
前記アップデート値を使って前記観察物体を表す物体空間画像を再構成する
ことによって、
前記投影空間データを受信し、前記投影空間データを処理するための画像再構成装置と
を備えることを特徴とするシステム。
【請求項２】
前記励起エネルギーは、電磁（ＥＭ）エネルギー、Ｘ線エネルギー、粒子線、赤外線エネルギー、光エネルギー、または振動エネルギーの少なくとも１つであることを特徴とする請求項１に記載のシステム。
【請求項３】
前記予測投影データは、物体空間ボクセルセットに基づいて計算され、前記物体空間ボクセルセットは、前記物体空間ボクセルセットが投影空間に投影されると、可変サイズの複数の所定の解像度グリッドをカバーすることを特徴とする請求項１に記載のシステム。
【請求項４】
前記第１の量および前記第２の量の逆関数を使ってアップデート値を計算することは、
前記第１の量を前記第２の量の前記反転関数で畳み込むことを含むことを特徴とする請求項３に記載のシステム。
【請求項５】
前記アップデート値を使って前記観察物体を表す前記物体空間画像を再構成することは、
前記物体空間ボクセルセットを前記アップデート値で更新することを含むことを特徴とする請求項１に記載のシステム。
【請求項６】
前記観察物体を表す前記物体空間画像を受信するための出力装置をさらに備え、前記出力装置は、データ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項１に記載のシステム。
【請求項７】
前記物体空間画像を受信するための調査用コンピュータであって、前記物体空間画像からの診断情報の抽出を実施するように、あるいは前記１つもしくは複数の送信機、前記１つもしくは複数の検出装置、または前記画像再構成装置の少なくとも１つのためのパラメータを微調整するようにプログラムされた、調査用コンピュータをさらに備えることを特徴とする請求項１に記載のシステム。
【請求項８】
１つまたは複数のプロセッサと、
前記１つまたは複数のプロセッサによって実行されるソフトウェアを格納する１つまたは複数のデータ記憶装置と
を備え、前記ソフトウェアは、
物体空間から投影空間に投影される物体を表す入力データを受信するためのソフトウェアコードと、
前記入力データと予測投影データの間の差を特徴づける第１の量を計算するためのソフトウェアコードと、
少なくとも１つのインパルス応答に対応する第２の量を計算するためのソフトウェアコードであって、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応する、ソフトウェアコードと、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算するためのソフトウェアコードと、
前記アップデート値を使って前記観察物体を表す物体空間画像を再構成するためのソフトウェアコードと
を備えることを特徴とするワークステーション。
【請求項９】
前記入力データは、投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから受信されることを特徴とする請求項８に記載のワークステーション。
【請求項１０】
前記入力データは、データバス、ケーブル、有線ネットワーク、または無線ネットワークを介して受信されることを特徴とする請求項８に記載のワークステーション。
【請求項１１】
前記画像を受信するための出力装置をさらに備え、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項８に記載のワークステーション。
【請求項１２】
前記ソフトウェアは、前記第１の量を前記第２の量の前記反転関数で畳み込むためのソフトウェアコードをさらに備えることを特徴とする請求項８に記載のワークステーション。
【請求項１３】
１つまたは複数のデータ記憶装置に格納されたソフトウェアコードを実行する１つまたは複数のプロセッサによって実装される方法であって、
投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから、物体空間から投影空間に投影される物体を表す入力データを受信するステップと、
解像度グリッド上の前記入力データと前記解像度グリッド上の予測投影データの間の差を特徴づける第１の量を計算するステップと、
少なくとも１つのインパルス応答に対応する第２の量を計算するステップであって、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応する、ステップと、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算するステップと、
前記アップデート値を使って前記観察物体を表す画像を再構成するステップと、
前記再構成画像を出力装置に出力するステップであって、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つである、ステップと
を含み、
前記入力データを受信するステップ、前記第１の量を計算するステップ、前記第２の量を計算するステップ、前記アップデート値を計算するステップ、前記画像を再構成するステップ、および前記再構成画像を出力するステップは、前記１つまたは複数のデータ記憶装置に格納された前記ソフトウェアコードに従って前記１つまたは複数のプロセッサによって実施されることを特徴とする方法。
【請求項１４】
前記解像度グリッドは、可変サイズの複数の所定の解像度グリッドから取得されることを特徴とする請求項１３に記載の方法。
【請求項１５】
前記予測投影データは、物体空間ボクセルセットに基づいて計算され、前記物体空間ボクセルセットは、前記物体空間ボクセルセットが投影空間に投影されると、可変サイズの前記複数の所定の解像度グリッドをカバーすることを特徴とする請求項１４に記載の方法。
【請求項１６】
前記解像度グリッド上の前記入力データまたは前記物体空間ボクセルセットの少なくとも１つを平滑カーネルで平滑化するステップをさらに含むことを特徴とする請求項１５に記載の方法。
【請求項１７】
前記物体空間ボクセルセットは、値の範囲内に抑制されることを特徴とする請求項１５に記載の方法。
【請求項１８】
前記第１の量は、前記入力データと前記予測投影データとの間の誤差の測度として特徴づけられることを特徴とする請求項１３に記載の方法。
【請求項１９】
前記第２の量は、前記入力データの測定値から取得されるノイズ分散からの寄与率を含むことを特徴とする請求項１３に記載の方法。
【請求項２０】
前記反転関数は、前記物体空間内の単位強度の前記参照ボクセルの前記場所をカバーする周辺エリアに関して計算されることを特徴とする請求項１３に記載の方法。
【請求項２１】
前記画像を再構成するステップは、
前記第１の量を前記第２の量の前記反転関数で畳み込んで、前記アップデート値を取得するステップを含むことを特徴とする請求項１３に記載の方法。
【請求項２２】
前記アップデート値は、限定関数によってさらに限定されることを特徴とする請求項１３に記載の方法。
【請求項２３】
前記アップデート値を使って前記物体空間ボクセルセットを更新するステップをさらに含むことを特徴とする請求項１３に記載の方法。
【請求項２４】
前記物体空間ボクセルセットが更新された後、前記解像度グリッドより大きいサイズをもつ異なる解像度グリッドに変更するステップをさらに含むことを特徴とする請求項２３に記載の方法。
【請求項２５】
前記１つまたは複数のプロセッサによって実行されて、請求項１３に記載の方法を実施する前記ソフトウェアを格納するためのものであることを特徴とする請求項１３に記載の１つまたは複数のデータ記憶装置。
【請求項１】
観察物体内に励起エネルギーを送信するための１つまたは複数の送信機と、
前記観察物体への前記送信励起エネルギーに応答して１つまたは複数の検出装置によって受信されるエネルギーをエンコードする投影空間データを生成するための前記１つまたは複数の検出装置と、
前記励起エネルギーを送信するための前記１つまたは複数の送信機および前記投影空間データを生成するための１つまたは複数の受信機を制御するためのコントローラと、
前記投影空間データと予測投影データの間の差を特徴づける第１の量を計算し、
各インパルス応答が単位強度の参照ボクセルに対応する少なくとも１つのインパルス応答に対応する第２の量を計算し、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算し、
前記アップデート値を使って前記観察物体を表す物体空間画像を再構成する
ことによって、
前記投影空間データを受信し、前記投影空間データを処理するための画像再構成装置と
を備えることを特徴とするシステム。
【請求項２】
前記励起エネルギーは、電磁（ＥＭ）エネルギー、Ｘ線エネルギー、粒子線、赤外線エネルギー、光エネルギー、または振動エネルギーの少なくとも１つであることを特徴とする請求項１に記載のシステム。
【請求項３】
前記予測投影データは、物体空間ボクセルセットに基づいて計算され、前記物体空間ボクセルセットは、前記物体空間ボクセルセットが投影空間に投影されると、可変サイズの複数の所定の解像度グリッドをカバーすることを特徴とする請求項１に記載のシステム。
【請求項４】
前記第１の量および前記第２の量の逆関数を使ってアップデート値を計算することは、
前記第１の量を前記第２の量の前記反転関数で畳み込むことを含むことを特徴とする請求項３に記載のシステム。
【請求項５】
前記アップデート値を使って前記観察物体を表す前記物体空間画像を再構成することは、
前記物体空間ボクセルセットを前記アップデート値で更新することを含むことを特徴とする請求項１に記載のシステム。
【請求項６】
前記観察物体を表す前記物体空間画像を受信するための出力装置をさらに備え、前記出力装置は、データ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項１に記載のシステム。
【請求項７】
前記物体空間画像を受信するための調査用コンピュータであって、前記物体空間画像からの診断情報の抽出を実施するように、あるいは前記１つもしくは複数の送信機、前記１つもしくは複数の検出装置、または前記画像再構成装置の少なくとも１つのためのパラメータを微調整するようにプログラムされた、調査用コンピュータをさらに備えることを特徴とする請求項１に記載のシステム。
【請求項８】
１つまたは複数のプロセッサと、
前記１つまたは複数のプロセッサによって実行されるソフトウェアを格納する１つまたは複数のデータ記憶装置と
を備え、前記ソフトウェアは、
物体空間から投影空間に投影される物体を表す入力データを受信するためのソフトウェアコードと、
前記入力データと予測投影データの間の差を特徴づける第１の量を計算するためのソフトウェアコードと、
少なくとも１つのインパルス応答に対応する第２の量を計算するためのソフトウェアコードであって、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応する、ソフトウェアコードと、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算するためのソフトウェアコードと、
前記アップデート値を使って前記観察物体を表す物体空間画像を再構成するためのソフトウェアコードと
を備えることを特徴とするワークステーション。
【請求項９】
前記入力データは、投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから受信されることを特徴とする請求項８に記載のワークステーション。
【請求項１０】
前記入力データは、データバス、ケーブル、有線ネットワーク、または無線ネットワークを介して受信されることを特徴とする請求項８に記載のワークステーション。
【請求項１１】
前記画像を受信するための出力装置をさらに備え、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つであることを特徴とする請求項８に記載のワークステーション。
【請求項１２】
前記ソフトウェアは、前記第１の量を前記第２の量の前記反転関数で畳み込むためのソフトウェアコードをさらに備えることを特徴とする請求項８に記載のワークステーション。
【請求項１３】
１つまたは複数のデータ記憶装置に格納されたソフトウェアコードを実行する１つまたは複数のプロセッサによって実装される方法であって、
投影画像取得システム、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、前記１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーション、または別の１つもしくは複数のプロセッサによって実行されるソフトウェアシミュレーションの少なくとも１つから、物体空間から投影空間に投影される物体を表す入力データを受信するステップと、
解像度グリッド上の前記入力データと前記解像度グリッド上の予測投影データの間の差を特徴づける第１の量を計算するステップと、
少なくとも１つのインパルス応答に対応する第２の量を計算するステップであって、各インパルス応答は、前記物体空間内の場所での単位強度の参照ボクセルに対応する、ステップと、
前記第１の量および前記第２の量の反転関数を使ってアップデート値を計算するステップと、
前記アップデート値を使って前記観察物体を表す画像を再構成するステップと、
前記再構成画像を出力装置に出力するステップであって、前記出力装置は、前記１つもしくは複数のデータ記憶装置、別の１つもしくは複数のデータ記憶装置、表示装置、または印刷装置の少なくとも１つである、ステップと
を含み、
前記入力データを受信するステップ、前記第１の量を計算するステップ、前記第２の量を計算するステップ、前記アップデート値を計算するステップ、前記画像を再構成するステップ、および前記再構成画像を出力するステップは、前記１つまたは複数のデータ記憶装置に格納された前記ソフトウェアコードに従って前記１つまたは複数のプロセッサによって実施されることを特徴とする方法。
【請求項１４】
前記解像度グリッドは、可変サイズの複数の所定の解像度グリッドから取得されることを特徴とする請求項１３に記載の方法。
【請求項１５】
前記予測投影データは、物体空間ボクセルセットに基づいて計算され、前記物体空間ボクセルセットは、前記物体空間ボクセルセットが投影空間に投影されると、可変サイズの前記複数の所定の解像度グリッドをカバーすることを特徴とする請求項１４に記載の方法。
【請求項１６】
前記解像度グリッド上の前記入力データまたは前記物体空間ボクセルセットの少なくとも１つを平滑カーネルで平滑化するステップをさらに含むことを特徴とする請求項１５に記載の方法。
【請求項１７】
前記物体空間ボクセルセットは、値の範囲内に抑制されることを特徴とする請求項１５に記載の方法。
【請求項１８】
前記第１の量は、前記入力データと前記予測投影データとの間の誤差の測度として特徴づけられることを特徴とする請求項１３に記載の方法。
【請求項１９】
前記第２の量は、前記入力データの測定値から取得されるノイズ分散からの寄与率を含むことを特徴とする請求項１３に記載の方法。
【請求項２０】
前記反転関数は、前記物体空間内の単位強度の前記参照ボクセルの前記場所をカバーする周辺エリアに関して計算されることを特徴とする請求項１３に記載の方法。
【請求項２１】
前記画像を再構成するステップは、
前記第１の量を前記第２の量の前記反転関数で畳み込んで、前記アップデート値を取得するステップを含むことを特徴とする請求項１３に記載の方法。
【請求項２２】
前記アップデート値は、限定関数によってさらに限定されることを特徴とする請求項１３に記載の方法。
【請求項２３】
前記アップデート値を使って前記物体空間ボクセルセットを更新するステップをさらに含むことを特徴とする請求項１３に記載の方法。
【請求項２４】
前記物体空間ボクセルセットが更新された後、前記解像度グリッドより大きいサイズをもつ異なる解像度グリッドに変更するステップをさらに含むことを特徴とする請求項２３に記載の方法。
【請求項２５】
前記１つまたは複数のプロセッサによって実行されて、請求項１３に記載の方法を実施する前記ソフトウェアを格納するためのものであることを特徴とする請求項１３に記載の１つまたは複数のデータ記憶装置。
【図１】
【図２】
【図３】
【図４Ａ】
【図４Ｂ】
【図５】
【図６Ａ】
【図６Ｂ】
【図７Ａ】
【図７Ｂ】
【図２】
【図３】
【図４Ａ】
【図４Ｂ】
【図５】
【図６Ａ】
【図６Ｂ】
【図７Ａ】
【図７Ｂ】
【公表番号】特表２０１１−５２６３６７（Ｐ２０１１−５２６３６７Ａ）
【公表日】平成２３年１０月６日（２０１１．１０．６）
【国際特許分類】
【出願番号】特願２０１１−５１６７８４（Ｐ２０１１−５１６７８４）
【出願日】平成２１年６月２９日（２００９．６．２９）
【国際出願番号】ＰＣＴ／ＵＳ２００９／０４９１０２
【国際公開番号】ＷＯ２００９／１５８７１８
【国際公開日】平成２１年１２月３０日（２００９．１２．３０）
【出願人】（５１０３３６２９９）
【Ｆターム（参考）】
【公表日】平成２３年１０月６日（２０１１．１０．６）
【国際特許分類】
【出願日】平成２１年６月２９日（２００９．６．２９）
【国際出願番号】ＰＣＴ／ＵＳ２００９／０４９１０２
【国際公開番号】ＷＯ２００９／１５８７１８
【国際公開日】平成２１年１２月３０日（２００９．１２．３０）
【出願人】（５１０３３６２９９）
【Ｆターム（参考）】
[ Back to top ]