(57)【要約】投影空間データを受信することと、１つ又は複数の測定された変換ピクセルを計算することと、１つ又は複数の予測変換ピクセルを計算することと、第１のピクセルイノベーション結果データを計算することと、第２のピクセルイノベーション結果データを計算することと、予備的変換物体更新データを計算することと、変換物体ボクセル密度更新推定値を計算することと、変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算することと、変換密度推定値を使用して観察対象の物体を表す物体空間画像を再構成することと、を実行するための、画像再構成器。【選択図】  図１
2019535405.tif 000002
【特許請求の範囲】
【請求項１】
投影空間データを格納するための非一時的データストレージであって、前記投影空間データが観察対象の物体の密度領域内にあり、１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、非一時的データストレージと、
少なくとも１つのプロセッサを含んでいる画像再構成コンピュータであって、前記少なくとも１つのプロセッサが、
前記密度領域内の前記投影空間データを前記非一時的データストレージから受信することと、
前記密度領域内の前記１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算することと、
前記密度領域内の前記１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して前記変換領域内の１つ又は複数の予測変換ピクセルを計算することと、
前記１つ又は複数の測定された変換ピクセルと前記１つ又は複数の予測変換ピクセルとの間の差異を使用して前記変換領域内の第１のピクセルイノベーション結果データを計算することと、
対応する入力ピクセル値及び参照ピクセル値に対して前記入力変換関数の反転傾斜及び前記参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算することと、
前記第１のピクセルイノベーション結果データと前記ピクセルごとのイノベーションスケーリング行列の対応する要素とのピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算することと、
断層撮影再構成アルゴリズムを使用し、前記第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算することと、
ボクセルごとの更新スケーリング行列の前記対応する要素を使用して前記予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算することであって、前記変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、前記ボクセルごとの更新スケーリング行列の要素に関連付けられる、ことと、
前記変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得することと、
前記変換密度推定値を使用して、前記観察対象の物体を表す物体空間画像を再構成することと、
を実行するよう機能する、画像再構成コンピュータと、
を備えている、システム。
【請求項２】
前記画像再構成コンピュータが、
ピクセルごとの補正ゲイン行列の要素を計算することであって、前記ピクセルごとの補正ゲイン行列の少なくとも１つの要素が、前記第２のピクセルイノベーション結果データ内の対応するピクセルに関連付けられる、ことと、
前記第２のピクセルイノベーション結果データと前記ピクセルごとの補正ゲイン行列とを使用して第３のピクセルイノベーション結果データを計算することであって、
前記予備的変換物体更新データが、前記断層撮影再構成アルゴリズム及び前記第３のピクセルイノベーション結果データを使用して計算される、ことと、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項３】
前記第１のピクセルイノベーション結果データ、前記第２のピクセルイノベーション結果データ、及び前記第３のピクセルイノベーション結果データが、影響関数を使用して処理され、
前記予備的変換物体更新データ及び前記変換物体ボクセル密度更新推定値が、影響関数を使用して処理される、請求項２に記載のシステム。
【請求項４】
前記画像再構成コンピュータが、
正の制約関数を使用して前記変換密度推定値を前記密度領域内に逆変換し、生のボクセル密度推定データを取得することと、
前記生のボクセル密度推定データをフィルタリングして、予備的ボクセル密度推定データを取得することと、
前記予備的ボクセル密度推定データを後処理して、（ａ）等しい又は増加したグリッド分解能、又は（ｂ）等しい又は増加した密度値の分解能のうちの１つ又は複数に基づいて、１つ又は複数のボクセルを作成することと、
投影プロセッサを使用して前記１つ又は複数のボクセルを処理することによって予測投影データを計算することと、
フィルタを前記予測投影データに適用して、平滑化された予測投影データを取得することと、
前記参照変換関数を前記平滑化された予測投影データに適用して、前記１つ又は複数の予測変換ピクセルを取得することと、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項５】
前記投影プロセッサが、少なくとも１つのシステムパラメータ又は物体パラメータを表す、請求項４に記載のシステム。
【請求項６】
前記画像再構成コンピュータが、
１つ又は複数の物体空間のボクセルのセットに基づいて、平滑化された前のボクセル密度推定データを計算することであって、前記１つ又は複数の物体空間のボクセルのセットが、サイズが変化する、単一のグリッド点を含む複数の分解能のグリッドをカバーする、ことと、
前記平滑化された前のボクセル密度推定データを使用して前記変換先行ボクセルデータ推定値を計算することと、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項７】
前記画像再構成コンピュータが、
反復回数を計算して追跡し、事前に決定された制限する反復回数を使用して、ある分解能レベルでそれ以上の計算を終了することと、
反復を計算し、使用しているある分解能レベルで、データのノイズレベル又はその他の基準に関して、物体の投影残差データ又は体積データの誤差の縮小が満たされたときに、反復を終了することと、
ある分解能レベルで、反復を計算し、（ｉ）十分な物体データに、非線形データ変換の小信号処理の基準が適用され、（ｉｉ）その後、終了の前に、反復回数が追加された後に、反復を終了することと、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項８】
投影空間データが、電磁（ＥＭ）エネルギー、Ｘ線エネルギー、粒子線、赤外線エネルギー、光エネルギー、又は振動エネルギーの投影収集プロセスのうちの少なくとも１つから、或いはその他のプロセスの出力から取得される、請求項１に記載のシステム。
【請求項９】
前記画像再構成コンピュータが、
ボクセルのセットについての投影値を計算することと、
逆投影値又はピクセルのセットについての値を表すボクセルを計算することと、
他の情報源を使用して残りのピクセル及びボクセルを計算することと、
ボクセルとピクセルの間の関数関係を使用するか、又は構築することと、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項１０】
前記観察対象の物体を表す前記物体空間画像を受信するための出力デバイスをさらに備えており、前記出力デバイスが、データストレージデバイス、表示デバイス、印刷デバイス、又はコンピュータシステムのうちの少なくとも１つである、請求項１に記載のシステム。
【請求項１１】
前記画像再構成コンピュータが、物体空間ボクセルデータのセット及び予測空間データのセットのうちの少なくとも１つを平滑化するようさらに機能する、請求項１に記載のシステム。
【請求項１２】
前記物体空間画像、前記密度領域からの投影残差、又は前記変換領域からの残差のうちの少なくとも１つを受信するためのインベスティゲータコンピュータをさらに備えており、前記インベスティゲータコンピュータが、少なくとも１つのプロセッサを含んでおり、前記少なくとも１つのプロセッサが、
前記物体空間画像、前記密度領域からの前記投影残差、又は前記変換領域からの残差のうちの少なくとも１つからの診断情報の抽出を実行することと、
前記１つ又は複数の入力投影ピクセル、投影方向、投影の変形、投影生成システムのプロセス、又は前記画像再構成コンピュータのうちの少なくとも１つを処理するためのパラメータを微調整することと、
を実行するよう機能する、請求項１に記載のシステム。
【請求項１３】
再帰的に、収束の後に、外部から制御可能な固定されたパラメータ成分の固定されたセットに関して、第１のイノベーションプロセス、投影残差、又は前記画像再構成コンピュータの物体データの期待値のうちの少なくとも１つから得られる性能の尺度を生成し、
前記画像再構成コンピュータの外部にある前記インベスティゲータコンピュータが、既知の参照物体と前記対応する計算された物体密度との間の前記差異を特徴付ける前記投影残差又は物体の残差のうちの少なくとも１つの特性を使用しながら、パラメータベクトルの制御可能な成分を、初期設定から、外部の目的に対して最低のコストを生成する設定に徐々に変更するようさらに機能する、請求項１２に記載のシステム。
【請求項１４】
前記パラメータベクトルが、少なくとも１つのシステムパラメータ又は物体パラメータを表す、請求項１３に記載のシステム。
【請求項１５】
前記画像再構成コンピュータが、
特定の反復内で、前記第１のピクセルイノベーション結果データの第１のイノベーション残差コスト関数を減らしながら、外部から制御可能な固定されたパラメータのセットに関して、観察対象の物体を再構成すること、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項１６】
前記画像再構成コンピュータが、
投影残差の近傍領域を表すか又は前記近傍領域から得られるパラメータ或いは先験的パラメータを使用してロバストな推定影響関数で前記投影残差を処理することであり、前記パラメータが、局所的な入力測定値、予測値、システムの特性、及び先験的に期待される物体及び測定の特性の関数である、処理すること、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項１７】
前記画像再構成コンピュータが、
選択範囲上の投影残差の平滑化、又は選択範囲上の前記投影残差のスケーリングのうちの少なくとも１つによって、前記投影残差を処理すること
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項１８】
前記入力変換関数が非線形入力変換関数であり、前記参照変換関数が非線形参照変換関数である、請求項１に記載のシステム。
【請求項１９】
前記非線形入力変換関数及び前記非線形参照変換関数が、前記１つ又は複数の入力投影ピクセルの分散安定化関数である、請求項１８に記載のシステム。
【請求項２０】
前記非線形入力変換関数及び前記非線形参照変換関数が、前記１つ又は複数の入力投影ピクセルのポアソン分散安定化平方根関数である、請求項１８に記載のシステム。
【請求項２１】
前記画像再構成コンピュータが、
前記ボクセルごとの更新スケーリング行列を計算することであり、前記ボクセルごとの更新スケーリング行列の要素が、対応する変換密度推定値の正の制約関数のボクセルごとの要素を使用して、ボクセルごとの逆変換の要素の反転傾斜から計算される、計算すること、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項２２】
前記入力変換関数が非線形入力変換関数であり、前記正の制約関数の全部又は一部が前記非線形入力変換関数の逆関数である、請求項２１に記載のシステム。
【請求項２３】
前記正の制約関数が、前記正の制約関数のデータ範囲の少なくとも１つの領域内の外部入力を満たす出力値を計算する、請求項２１に記載のシステム。
【請求項２４】
前記画像再構成コンピュータが、
前記ボクセルごとの更新スケーリング行列を計算することであり、前記ボクセルごとの更新スケーリング行列の要素が、前のボクセル密度推定データを変換先行ボクセルデータ推定値に変換するフィードバック関数の傾斜から計算される、計算すること、
を実行するようにさらに機能する、請求項１に記載のシステム。
【請求項２５】
前記画像再構成コンピュータが、
先行する物体密度のボクセルの予測を使用して、信頼性が低いか、又は不足している投影空間データの一部を再計算又は補正すること、
を実行するようさらに機能する、請求項１に記載のシステム。
【請求項２６】
１つ又は複数のプロセッサを含んでいるコンピュータによって実行可能なソフトウェアコードを格納する非一時的データストレージデバイスであって、前記ソフトウェアコードが、
観察対象の物体の密度領域内の投影空間データを受信することであって、前記投影空間データが、１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、ことと、
前記密度領域内の前記１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算することと、
前記密度領域内の前記１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して前記変換領域内の１つ又は複数の予測変換ピクセルを計算することと、
前記１つ又は複数の測定された変換ピクセルと前記１つ又は複数の予測変換ピクセルとの間の差異を使用して前記変換領域内の第１のピクセルイノベーション結果データを計算することと、
対応する入力ピクセル値及び参照ピクセル値に対して前記入力変換関数の反転傾斜及び前記参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算することと、
前記第１のピクセルイノベーション結果データと前記ピクセルごとのイノベーションスケーリング行列の対応する要素とのピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算することと、
断層撮影再構成アルゴリズムを使用し、前記第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算することと、
ボクセルごとの更新スケーリング行列の前記対応する要素を使用して前記予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算することであって、前記変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、前記ボクセルごとの更新スケーリング行列の要素に関連付けられる、ことと、
前記変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得することと、
前記変換密度推定値を使用して、前記観察対象の物体を表す物体空間画像を再構成することと、
を実行する、非一時的データストレージデバイス。
【請求項２７】
少なくとも１つのプロセッサを含んでいる画像再構成コンピュータによって実行される画像再構成のための方法であって、
観察対象の物体の密度領域内の投影空間データを受信するステップであって、前記投影空間データが、１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、ステップと、
前記密度領域内の前記１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算するステップと、
前記密度領域内の前記１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して前記変換領域内の１つ又は複数の予測変換ピクセルを計算するステップと、
前記１つ又は複数の測定された変換ピクセルと前記１つ又は複数の予測変換ピクセルとの間の差異を使用して前記変換領域内の第１のピクセルイノベーション結果データを計算するステップと、
対応する入力ピクセル値及び参照ピクセル値に対して前記入力変換関数の反転傾斜及び前記参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算するステップと、
前記第１のピクセルイノベーション結果データと前記ピクセルごとのイノベーションスケーリング行列の対応する要素とのピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算するステップと、
断層撮影再構成アルゴリズムを使用し、前記第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算するステップと、
ボクセルごとの更新スケーリング行列の前記対応する要素を使用して前記予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算するステップであって、前記変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、前記ボクセルごとの更新スケーリング行列の要素に関連付けられる、ステップと、
前記変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得するステップと、
前記変換密度推定値を使用して、前記観察対象の物体を表す物体空間画像を再構成するステップと、
を含んでいる、方法。
【発明の詳細な説明】
【関連出願の相互参照】
【０００１】
[0001]本特許出願は、２０１６年１１月１８日に出願された米国特許出願第６２／４２４，１８７号、及び２０１６年１１月２３日に出願された米国特許出願第６２／４２６，０６５号の優先権を請求し、これらの米国特許出願の内容全体は、両方とも参照によって本明細書に組み込まれる。
【発明の分野】
【０００２】
[0002]本発明は、画像再構成の分野に関連している。
【発明の背景】
【０００３】
[0003]断層撮影（ある領域（ピクセルとして表される）のｍ次元の投影に基づくｎ次元空間のその領域内の密度（ボクセルとして表される）の計算（推定）（通常は、０＜ｍ＜ｎ））は、線形のフィルタ補正逆投影（ＦＢＰ：ｆｉｌｔｅｒｅｄ  ｂａｃｋ−ｐｒｏｊｅｃｔｉｏｎ）法及び高速フーリエ変換（ＦＦＴ：ｆａｓｔ  Ｆｏｕｒｉｅｒ  ｔｒａｎｓｆｏｒｍ）法、並びに最新の非線形の反復再構成（ＩＲ：ｉｔｅｒａｔｉｖｅ  ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ）（代数的再構成法（ＡＲＴ：Ａｌｇｅｂｒａｉｃ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  Ｔｅｃｈｎｉｑｕｅ）、同時反復代数的再構成法（ＳＩＲＴ及びＳＡＲＴ：Ｓｉｍｕｌｔａｎｅｏｕｓ  Ｉｔｅｒａｔｉｖｅ  Ａｌｇｅｂｒａｉｃ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  Ｔｅｃｈｎｉｑｕｅ）、及びモデルベース反復再構成（ＭＢＩＲ：Ｍｏｄｅｌ  Ｂａｓｅｄ  Ｉｔｅｒａｔｉｖｅ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ）など）という２つの主要なカテゴリに分類される。しかし、最新の断層撮影の設定からの画像再構成は、通常、非線形の推定の使用に起因して、複雑である。断層撮影の設定は、１つの投影が２，０００ピクセル（１次元）〜４，０００×４，０００ピクセル（２次元）程度を含んでいる場合、１０６〜１０１０個の画素（すなわち、ボクセル）の計算を必要とすることが多い。
【０００４】
[0004]ＦＢＰに関連する重大な欠陥として、限定された定量的精度を実現するために、多数の投影が必要になることが挙げられる。投影の数は、通常、数百又は数千になるが、投影は、可能な限り効率的に使用されているわけではない。
【０００５】
[0005]ＩＲ法の利点として、ＦＢＰに従うときに、再構成誤差を減らすことができることが挙げられる。関連する処理手法（特に、線形手法及びロバスト手法）及びそれらの具体的な利点の概要が、Ｓｕｎｎｅｇａｒｄｈ［１５］に記載されている。さらに、ＩＲは、制約を考慮することができ、特に、密度の推定値が非負になることを保証することができる。最適な物体再構成を見つけるには、すべての再構成誤差の合計である目的関数を最小化する必要がある。最適な再構成を見つけるには、各ボクセルに対する処理が必要である。目的関数は、最低誤差の近くでは、通常、目的関数の二次導関数を表すヘッシアン行列で表され得る。断層撮影では、どの次元においても１，０００ボクセル以上である３次元体積（１，０００×１，０００×１，０００＝１０９ボクセル以上に相当する）の場合、このヘッシアンは、１０９×１０９＝１０１８以上の要素を含むことがある。しかし、このような多数の要素を格納し、操作することは、現在のコンピュータ及び予測可能な将来のシステムの能力の限界を超えている。
【０００６】
[0006]これまで、すべての市販のＩＲシステムの性能は、このヘッシアン行列の特性によって制限されてきた。準ニュートン法を使用する反復は、ヘッシアン行列の評価を回避する。それらの手法の場合、ヘッシアンの最大固有値に対処する最小化ステップを識別する傾向がある。このために、初期ＦＢＰの後に、反復改良が続く。しかし、画像再構成用のヘッシアン行列のサイズが大きいため、改良中に、ヘッシアン行列の構造（及び、ヘッシアン行列の逆行列の構造）が、通常は無視されるか、又は不完全に近似される。さらに、ヘッシアンの固有値の広範な分布のため、現在の最適化手法は、ある反復回数（通常は、数十回〜数千回）を超えると改善を示さなくなる傾向があり、いつかの大きい固有値のみに対処することがある。
【０００７】
[0007]しかし、これらのアルゴリズムのマルチグリッドの変形は、役立つことがあるが、細かいグリッドに関連するヘッシアンのサイズのため、やはり最終的に失敗する。マルチグリッドの分解能とは、本明細書では、反復が実行されるにつれて徐々に細かくなる分解能を使用することを指す。
【０００８】
[0008]さらに、例えば、物体の変形、X線画像化の場合のビームハードニング及び散乱、又は磁気共鳴画像法（ＭＲＩ：ｍａｇｎｅｔｉｃ  ｒｅｓｏｎａｎｃｅ  ｉｍａｇｉｎｇ）の場合の電磁界歪を特徴付ける、数十個又は数百個のシステムパラメータにおけるシステムの不確実性も、測定データの非一貫性を増大させることがある。
【概要】
【０００９】
[0009]本発明の一部の実施形態に従うシステムは、投影空間データを格納するための非一時的データストレージであって、この投影空間データが観察対象の物体の密度領域内にあり、１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、非一時的データストレージと、少なくとも１つのプロセッサを含んでいる画像再構成コンピュータであって、この少なくとも１つのプロセッサが、密度領域内の投影空間データを非一時的データストレージから受信することと、密度領域内の１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算することと、密度領域内の１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して変換領域内の１つ又は複数の予測変換ピクセルを計算することと、１つ又は複数の測定された変換ピクセルと１つ又は複数の予測変換ピクセルの間の差異を使用して変換領域内の第１のピクセルイノベーション結果データを計算することと、対応する入力ピクセル値及び参照ピクセル値に対して入力変換関数の反転傾斜及び参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算することと、第１のピクセルイノベーション結果データ及びピクセルごとのイノベーションスケーリング行列の対応する要素のピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算することと、断層撮影再構成アルゴリズムを使用し、第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算することと、ボクセルごとの更新スケーリング行列の対応する要素を使用して予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算することであって、変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、ボクセルごとの更新スケーリング行列の要素に関連付けられる、ことと、変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得することと、変換密度推定値を使用して、観察対象の物体を表す物体空間画像を再構成することとを実行するよう機能する、画像再構成コンピュータとを備える。
【００１０】
[0010]本発明の一部の実施形態に従って、非一時的データストレージデバイスは、１つ又は複数のプロセッサを含んでいるコンピュータによって実行可能なソフトウェアコードを格納し、このソフトウェアコードは、観察対象の物体の密度領域内の投影空間データを受信することであって、この投影空間データが１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、ことと、密度領域内の１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算することと、密度領域内の１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して変換領域内の１つ又は複数の予測変換ピクセルを計算することと、１つ又は複数の測定された変換ピクセルと１つ又は複数の予測変換ピクセルの間の差異を使用して変換領域内の第１のピクセルイノベーション結果データを計算することと、対応する入力ピクセル値及び参照ピクセル値に対して入力変換関数の反転傾斜及び参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算することと、第１のピクセルイノベーション結果データ及びピクセルごとのイノベーションスケーリング行列の対応する要素のピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算することと、断層撮影再構成アルゴリズムを使用し、第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算することと、ボクセルごとの更新スケーリング行列の対応する要素を使用して予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算することであって、変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、ボクセルごとの更新スケーリング行列の要素に関連付けられる、ことと、変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得することと、変換密度推定値を使用して、観察対象の物体を表す物体空間画像を再構成することとを実行する。
【００１１】
[0011]本発明の一部の実施形態に従って、少なくとも１つのプロセッサを含んでいる画像再構成コンピュータによって実行される画像再構成のための方法は、観察対象の物体の密度領域内の投影空間データを受信することであって、この投影空間データが１つ又は複数の入力投影ピクセル及び１つ又は複数の予測投影空間ピクセルを含んでいる、ことと、密度領域内の１つ又は複数の入力投影ピクセル及び入力変換関数を使用して変換領域内の１つ又は複数の測定された変換ピクセルを計算することと、密度領域内の１つ又は複数の予測投影空間ピクセル及び参照変換関数を使用して変換領域内の１つ又は複数の予測変換ピクセルを計算することと、１つ又は複数の測定された変換ピクセルと１つ又は複数の予測変換ピクセルの間の差異を使用して変換領域内の第１のピクセルイノベーション結果データを計算することと、対応する入力ピクセル値及び参照ピクセル値に対して入力変換関数の反転傾斜及び参照変換関数の反転傾斜を使用してピクセルごとのイノベーションスケーリング行列を計算することと、第１のピクセルイノベーション結果データ及びピクセルごとのイノベーションスケーリング行列の対応する要素のピクセルごとの積を使用して第２のピクセルイノベーション結果データを計算することと、断層撮影再構成アルゴリズムを使用し、第２のピクセルイノベーション結果データに基づいて予備的変換物体更新データを計算することと、ボクセルごとの更新スケーリング行列の対応する要素を使用して予備的変換物体更新データをスケーリングすることによって、変換物体ボクセル密度更新推定値を計算することであって、変換物体ボクセル密度更新推定値の少なくとも１つのボクセルが、ボクセルごとの更新スケーリング行列の要素に関連付けられる、ことと、変換物体ボクセル密度更新推定値を、対応する変換先行ボクセルデータ推定値に加算して、変換密度推定値を取得することと、変換密度推定値を使用して、観察対象の物体を表す物体空間画像を再構成することとを含む。
【図面の簡単な説明】
【００１２】
[0012]ここで、実施形態例が、関連する図面に関して説明される。
【００１３】
【図１】[0013]図１は、本発明の実施形態例に従って、画像化システムを示す図である。
【００１４】
【図２Ａ】[0014]図２Ａは、反復ｉでの、実施形態例に従って再構成された物体の密度を使用する最適化された再帰法によって拡張された高効率のコンピュータ断層撮影（ＣＴ：ｃｏｍｐｕｔｅｒ  ｔｏｍｏｇｒａｐｈｙ）の例示的なデータ入力、参照予測、及び内部フィードバックループ１６５を示す図である（図２と総称される）。
【図２Ｂ】[0014]図２Ｂは、反復ｉでの、実施形態例に従って再構成された物体の密度を使用する最適化された再帰法によって拡張された高効率のコンピュータ断層撮影（ＣＴ：ｃｏｍｐｕｔｅｒ  ｔｏｍｏｇｒａｐｈｙ）の例示的なデータ入力、参照予測、及び内部フィードバックループ１６５を示す図である（図２と総称される）。
【００１５】
【図３】[0015]図３は、造影剤の注入中の例示的なＸ線投影を示す図である。
【００１６】
【図４】[0016]図４は、３Ｄ円錐ビームの再構成に適用されるＬｅｖｅｎｂｅｒｇ-Ｍａｒｑｕａｒｄｔプロセスの実施形態例に従って、６回の投影からの鼓動する水力学的冠動脈枝モデルの再構成を示す図である。
【００１７】
【図５】[0017]図５は、例示的な影響関数を示す図である。
【００１８】
【図６Ａ】[0018]図６Ａは、２００回の投影（図６Ａ）及び５回の投影（図６Ｂ）からの頭部のファントムの例示的な再構成を示す図である（図６と総称される）。
【図６Ｂ】[0018]図６Ｂは、２００回の投影（図６Ａ）及び５回の投影（図６Ｂ）からの頭部のファントムの例示的な再構成を示す図である（図６と総称される）。
【発明を実施するための形態】
【００１９】
[0019]以下では、本発明の実施形態例が詳細に説明される。実施形態例の説明では、明確にするために、特定の用語が採用されている。しかし、本発明は、そのように選択された特定の用語に限定されるよう意図されていない。当業者は、本発明の広範な概念から逸脱することなく、その他の同等の構成要素が採用されてもよく、その他の方法が開発されてもよいということを認識するであろう。本明細書において引用されているすべての参考文献は、それぞれが個別に組み込まれているかのように、参照によって組み込まれる。
【００２０】
[0020]本発明者は、最適化された再帰法によって拡張された高効率のコンピュータ断層撮影の計算を効率的にするための、小信号線形化及び部分問題の簡略化の近似法を発見した。
【００２１】
[0021]本明細書に記載された手法及びシステムは、既知の既存の手法（例えば、米国特許第８，６６０，３２８号及び第８，６６０，３３０号に記載された手法）を拡張することができる。
【００２２】
[0022]図１は、本発明の実施形態例に従って、画像化システムを示している。画像化システムは、コントローラ１３５と、調査器（investigator、インベスティゲータ）１７０（例えば、コンピュータ）との信号通信及びコントローラ１３５との同期のためのインターフェイス１４０（例えば、グラフィックユーザインターフェイス、キーボード、ジョイスティックなど）と、コントローラ１３５からの制御信号に応答した励起エネルギー１２０の生成及び放出のための送信器と、測定されたピクセルデータ１５５を生成するように構成された検出器１４５とを含んでもよい。測定されたピクセルデータ１５５は、デジタル形式でデータストレージデバイス１８０に格納されてもよい。
【００２３】
[0023]物体１２５の内部構造をエンコードする情報を含んでいる測定されたピクセルデータ１５５は、画像再構成器１６５（例えば、コンピュータ）を使用して、再構成済み画像データ２６１に変換され、例えば出力デバイス１７５上で可視化されてもよい。測定されたピクセルデータ１５５は、実験的取得１１０、シミュレーション１８５（例えば、コンピュータ上でのシミュレーション）、及び／又は例えば以前の実験又はシミュレーションから記録された投影ピクセルデータを含んでいるデータストレージデバイス１８０から取得されてもよい。実験的取得１１０、データストレージデバイス１８０、及びシミュレーション１８５は、測定されたピクセルデータ１５５を、画像再構成器１６５に直接、又はネットワーク（例えば、ローカルエリアネットワーク（ＬＡＮ：ｌｏｃａｌ  ａｒｅａ  ｎｅｔｗｏｒｋ）、広域ネットワーク（ＷＡＮ：ｗｉｄｅ  ａｒｅａ  ｎｅｔｗｏｒｋ）、インターネットなど）を介して、リモートから提供してもよい。測定されたピクセルデータ１５５は、任意の形態の励起エネルギーを観察対象の物体１２５に供給した結果を含んでもよく、したがって、Ｘ線コンピュータ断層撮影（ＣＴ）に加えて、例えば、電子顕微鏡、磁気共鳴（ＭＲ：ｍａｇｎｅｔｉｃ  ｒｅｓｏｎａｎｃｅ）、陽電子放射断層撮影（ＰＥＴ：ｐｏｓｉｔｒｏｎ  ｅｍｉｓｓｉｏｎ  ｔｏｍｏｇｒａｐｈｙ）、単一陽電子放射コンピュータ断層撮影（ＳＰＥＣＴ：ｓｉｎｇｌｅ  ｐｏｓｉｔｒｏｎ  ｅｍｉｓｓｉｏｎ  ｃｏｍｐｕｔｅｄ  ｔｏｍｏｇｒａｐｈｙ）、超音波、蛍光発光、多光子顕微鏡（ＭＰＭ：ｍｕｌｔｉ−ｐｈｏｔｏｎ  ｍｉｃｒｏｓｃｏｐｙ）、光干渉断層撮影（ＯＣＴ：ｏｐｔｉｃａｌ  ｃｏｈｅｒｅｎｃｅ  ｔｏｍｏｇｒａｐｈｙ）、コンピュータ断層撮影（ＣＴ）、電磁気（ＥＭ：ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ）エネルギー、Ｘ線エネルギー、粒子線、赤外線エネルギー、光エネルギー、及び／又は振動エネルギーのなどの投影取得プロセスからのデータを含んでもよい。データストレージデバイス１８０は、例えば、非一時的メモリ、ＣＤ−ＲＯＭ、ＤＶＤ、光磁気（ＭＯ：ｍａｇｎｅｔｏ−ｏｐｔｉｃａｌ）ディスク、ハードディスク、フロッピーディスク、ジップディスク、フラッシュドライブ、クラウドデータストレージなどであってもよい。実験的取得１１０及び／又はシミュレーション１８５からの測定されたピクセルデータ１５５は、データストレージデバイス１８０に格納されてもよい。
【００２４】
[0024]画像化システムは、物体１２５を受け取るように構成され、送信器１１５及び検出器１４５との関係において平行移動又は回転するよう機能する、平行移動又は回転するテーブル１３０と、電気的に、又は任意選択的なデータトランスポータ１５０（インターネット、その他のネットワーク、又はデータバスなど）を介して検出器１４５及びコントローラ１３５に結合された、画像再構成器１６５とを含んでもよい。データバスは、例えば、コンピュータ内のコンピュータの構成要素間又はコンピュータ間でデータを転送する電線のサブシステムであってもよい。接続された構成要素間のデータフローは一方向であるが、接続された構成要素間の通信は双方向であってもよい。
【００２５】
[0025]励起エネルギー１２０は、例えばＸ線エネルギー、電磁気（ＥＭ）エネルギー、光エネルギー、赤外線（ＩＲ）エネルギー、粒子エネルギー（例えば、電子ビーム、中性子ビーム、原子ビーム）、超音波などの振動エネルギーなどの、放射エネルギーの形態であってもよい。励起エネルギー１２０は、物体１２５に照射されてもよく、物体１２５は、例えば、ファントム、人間の患者、試料、又はこれらの任意の組み合わせであってもよい。励起エネルギー１２０は、対応するエネルギータイプの送信器１１５によって放出されてもよい。励起エネルギー１２０は、物体１２５を伝搬してもよく、一部が適切な検出器１４５によって受信されてもよい。検出器１４５は、受信されたエネルギーを測定可能な電気信号に変換してもよく、測定された電気信号は、デジタル形式の投影ピクセルデータにさらに変換されてもよい。
【００２６】
[0026]コントローラ１３５は、送信器１１５及び検出器１４５を接続する回路であってもよく、制御信号を送信器１１５及び検出器１４５に送信して、励起エネルギー１２０の送信及び検出器１４５の動作を同期させてもよい。この回路は、アナログ信号、デジタル信号、又は混合信号であってもよい。コントローラ１３５は、制御信号を送信器１１５及び検出器１４５に送信するために、１つ又は複数のプロセッサ及び１つ又は複数のメモリを含んでいるコンピュータであってもよい。
【００２７】
[0027]画像再構成器１６５は、コントローラ１３５に応答し、測定されたピクセルデータ１５５を受信することができ、本発明の実施形態に従う方法によって物体１２５の画像を再構成して、高い計算効率で物体の高忠実度の画像を生成することができる。画像再構成器１６５は、例えば、本発明の実施形態例に従ってコンピュータを動作させるためのソフトウェアコードを格納する１つ又は複数のデータストレージデバイスを含んでいるコンピュータであってもよい。このコンピュータは、１つ又は複数のデータストレージデバイスに格納されたソフトウェアコードを読み込んで実行するために、１つ又は複数のプロセッサ及び非一時的コンピュータ可読媒体を含んでもよい。別の実施形態例では、画像再構成器１６５は、本発明の実施形態例に従って動作可能な１つ又は複数のプログラム格納デバイス及びプログラム実行デバイスを含んでもよい。画像再構成器１６５は、再構成済み画像データ２６１を生成してもよい。画像再構成器１６５は、測定されたピクセルデータ１５５を受信してもよく、例えば非線形変換２００によって、測定されたピクセルデータ１５５を処理し、累積物体密度投影データ２０１を生成し、測定された変換ピクセル２０６、投影データ２０１の非線形変換をさらに計算し、予測変換ピクセルデータ２１１を超える超過を計算して、イノベーションを生成し、再構成済み画像データ２６１を生成してもよく、例えば図２の予測投影データ２７６をさらに生成してもよい。
【００２８】
[0028]出力デバイス１７５は、再構成済み画像データ２６１、再構成済み画像誤差データ１６０、又は再構成画像誤差データ２１６のうちの１つ又は複数を受信してもよい。出力デバイス１７５は、例えば、視覚化デバイス又はデータストレージデバイスであってもよい。視覚化デバイスは、例えば、表示デバイス又は印刷デバイスであってもよい。表示デバイスの例としては、例えば、ブラウン管（ＣＲＴ：ｃａｔｈｏｄｅ  ｒａｙ  ｔｕｂｅ）、発光ダイオード（ＬＥＤ：ｌｉｇｈｔ−ｅｍｉｔｔｉｎｇ  ｄｉｏｄｅ）ディスプレイ、液晶ディスプレイ（ＬＣＤ：ｌｉｑｕｉｄ  ｃｒｙｓｔａｌ  ｄｉｓｐｌａｙ）、デジタル光投影（ＤＬＰ：ｄｉｇｉｔａｌ  ｌｉｇｈｔ  ｐｒｏｊｅｃｔｉｏｎ）モニタ、真空蛍光ディスプレイ（ＶＦＤ：ｖａｃｕｕｍ  ｆｌｏｒｅｓｃｅｎｔ  ｄｉｓｐｌａｙ）、表面伝導型電子放出素子ディスプレイ（ＳＥＤ：ｓｕｒｆａｃｅ−ｃｏｎｄｕｃｔｉｏｎ  ｅｌｅｃｔｒｏｎ−ｅｍｉｔｔｅｒ  ｄｉｓｐｌａｙ）、電界放出ディスプレイ（ＦＥＤ：ｆｉｅｌｄ  ｅｍｉｓｓｉｏｎ  ｄｉｓｐｌａｙ）、ＬＣＯＳ（ｌｉｑｕｉｄ  ｃｒｙｓｔａｌ  ｏｎ  ｓｉｌｉｃｏｎ）ディスプレイなどが挙げられる。印刷デバイスの例としては、例えば、トナーベースプリンタ、液体インクジェットプリンタ、個体インクプリンタ、昇華型プリンタ、並びに感熱式プリンタ及び紫外線（ＵＶ：ｕｌｔｒａｖｉｏｌｅ）プリンタなどのインクレスプリンタなどが挙げられる。印刷デバイスは、３次元（３Ｄ：ｔｈｒｅｅ  ｄｉｍｅｎｓｉｏｎｓ）で印刷してもよい。出力デバイス１７５は、観察対象の物体１２５を表す物体空間画像情報を受信してもよい。出力デバイス１７５の例としては、データストレージデバイス、表示デバイス、印刷デバイス、又は別のコンピュータシステムが挙げられる。
【００２９】
[0029]画像化システムは、インベスティゲータ１７０をさらに含んでよい。インベスティゲータ１７０は、プログラムされたコンピュータであってもよく、再構成済み画像データ２６１、再構成済み画像誤差データ１６０、又は再構成画像誤差データ２１６のうちの１つ又は複数を受信し、その後、アルゴリズム（例えば、事前にプログラムされたルーチン、人工知能、機械学習など）を適用して、物体１２５に関する診断情報を抽出するか、又は送信器１１５、検出器１４５、テーブル１３０、画像再構成器１６５などの制御パラメータを微調整してもよい。一部の実施形態では、インターフェイス１４０又は出力デバイス１７５は、不要であってもよい。一部の実施形態では、インベスティゲータ１７０、コントローラ１３５、及び画像再構成器１６５は、同じコンピュータ上又は別々のコンピュータ上に存在してもよい。インベスティゲータ１７０は、データ２６１、１６０、又は２１６を受信してもよく、データからの診断情報の抽出を実行するようにプログラムされるか、又は例えば１つ又は複数の密度領域の入力ピクセル、投影方向、投影の変形、投影生成システムのプロセス（焦点サイズなど）、又は画像再構成器１６５のうちの少なくとも１つを処理するためのパラメータを微調整するようにプログラムされてもよい。インベスティゲータ１７０は、本明細書では「インベスティゲータコンピュータ」とも呼ばれる。
【００３０】
[0030]本発明の一部の実施形態は、画像再構成器１６５に類似する方法で画像を再構成するように構成された１つ又は複数のプロセッサを備えているワークステーションを提供してもよい。このワークステーションは、画像化システム、データストレージデバイス、又はコンピュータのうちの少なくとも１つから入力データを受信してもよい。入力データは、データバス、ケーブル、有線ネットワーク、無線ネットワークなどを介して受信されてもよい。このワークステーションは、再構成済み画像を受信するための出力デバイス１７５をさらに備えてもよい。この出力デバイスは、データストレージデバイス、表示デバイス、印刷デバイスなどであってもよい。例示的なデータストレージデバイス、表示デバイス、及び印刷デバイスは、前述したとおりである。
【００３１】
[0031]図２は、実施形態に従う系統図を示している。特に、図２Ａは、最適化された再帰法によって拡張された高効率のＣＴに関して、先行推定物体密度２６１の投影２７６と共に、例示的な生の測定データ入力１５５と、変換２００を使用する生の測定値からの投影空間密度データ２０１の計算と、変換投影データ２０６を取得するための変換２０５と、フィルタ２８５を使用する参照変換予測２１１とを示している。イノベーションゲイン調整２２０及び２２５は、図２Ａと図２Ｂの間のブリッジを形成する。図２Ｂは、ボクセル密度推定データ２６１を取得するためのフィードバックセクションと共に、最適化された再帰法によって拡張された高効率のＣＴの例示的な前方処理及びスケーリングされた反転処理を示している。上の処理における一部の、ただし全部ではない変換は、Ｓｕｎｎｅｇａｒｄｈ［１５］における仮定とは対照的に、線形であってもよい。
【００３２】
[0032]一実施形態では、本明細書に記載された発明は、励起エネルギー１２０を観察対象の物体１２５に送信するための１つ又は複数の送信器１１５と、観察対象の物体１２５に送信された励起エネルギー１２０に応答して１つ又は複数の検出器１４５によって受信されたエネルギーをエンコードする投影空間データを生成するための１つ又は複数の検出器１４５と、１つ又は複数の送信器１１５を制御して励起エネルギー１２０を送信し、１つ又は複数の検出器１４５を制御して投影空間データ１５５を生成するためのコントローラ１３５と、投影空間データ１５５を受信し、例えば図２に示され、本明細書において説明されたプロセスによって、投影空間データ１５５を処理するために、少なくとも１つのプロセッサを含んでいる、画像再構成器１６５とを含んでもよい。
【００３３】
[0033]一部の実施形態に従って、画像再構成器１６５は、ボクセルのセットに対して投影値を計算することと、ピクセルのセットに対して投影値を計算することと、逆投影値を計算することと、他の情報源を使用して残りのピクセル及びボクセルを計算することと、特定の物体の先験的に知られた高密度のボクセル及び対応するピクセルの期待される投影範囲、又はボクセル及び対応する不足している投影測定値などの、ボクセルとピクセルの間の関数関係を使用するか、又は構築することとを実行してもよい。例えば、ミッシングウエッジ問題（Ｌ．Ｐａａｖｏｌａｉｎｅｎ他［１６］）又はスパース投影などにおいて、一部の投影値が不足している場合、Ｄｅｍｐｓｔｅｒ  他［１７］に記載されている期待値最大化（ＥＭ：ｅｘｐｅｃｔａｔｉｏｎ  ｍａｘｉｍｉｚａｔｉｏｎ）法を使用して予測を行って、（期待される）不足しているピクセル測定値を置き換えてもよい。
【００３４】
[0034]ＥＭ法は、無限の観測ノイズに相当する不足しているデータのケースを表す。画像の一部の領域では、信頼性が低いと考えられるデータ（ＤＲＬＲ：Ｄａｔａ  ｒｅｇａｒｄｅｄ  ｌｅｓｓ  ｒｅｌｉａｂｌｅ）は、無限の観測ノイズではなく、観測ノイズが増加しているデータとして同等に表されてもよい。ＥＭ法と共に、ＤＲＬＲを含む利用可能なすべてのデータを使用して、それらの領域内のデータを非常に効率的に再計算／推定し、再計算／補正されたＥＭデータ（ＲＥＭＤ：ｒｅ−ｃｏｍｐｕｔｅｄ／ｃｏｒｒｅｃｔｅｄ  ＥＭ  ｄａｔａ）を生成することができる。その後、ＤＲＬＲ及びＲＥＭＤの重み付きの組み合わせを、次のプロセスの反復の入力値として使用して、データの信頼性に依存するデータの重み付けの問題を解決することができる。相対的な重み付けは、例えば、結合されたノイズレベルに基づいてもよい。そのため、画像再構成器１６５は、少なくとも１つの信頼性が低いデータ又は不足している測定データの一部を、重み付き予測投影データ２７６を使用して、再計算又は補正してもよい。
【００３５】
[0035]本発明の実施形態は、例えば、１つ又は複数のプロセッサを含んでいるワークステーションと、１つ又は複数のプロセッサによって実行されるソフトウェアを格納している１つ又は複数の非一時的データストレージデバイス１８０とを含んでもよく、このソフトウェアは、図２に示され、本明細書において説明されているプロセスを実施するためのソフトウェアコードを含んでもよい。
【００３６】
[0036]本発明の実施形態は、１つ又は複数のデータストレージデバイス１８０に格納されたソフトウェアコードを実行する１つ又は複数のプロセッサによって実施される方法を提供してもよく、この方法は、図２に示され、本明細書において説明されているプロセスを実施するためのステップを含んでいる。
【００３７】
[0037]例えば、画像再構成器１６５は、投影空間データ１５５を変換して投影ピクセルデータ２０１を取得することと、投影ピクセルデータ２０１を非線形に変換して測定された変換ピクセルデータ２０６を取得することと、第１のピクセルイノベーション結果データ２１６を計算して、測定された変換ピクセルデータ２０６と予測変換ピクセルデータ２１１の間の差異を特徴付けることと、第１のピクセルイノベーション結果データ２１６をデータストレージデバイス１８０に記録することと、第２のピクセルイノベーション結果データ２２１を計算することであって、第１のピクセルイノベーション結果データ２１６が非線形変換２０５及び２１０の２つのセットの傾斜の反転に基づいて再スケーリングされる、ことと、第２のピクセルイノベーション結果データ２２１に基づいて、先行する反復から第３のピクセルイノベーション結果データ２２６を計算することと、断層撮影再構成アルゴリズムを使用して予備的変換物体密度更新データ２３６を取得して、第２のピクセルイノベーション結果データ２２１を近似的に反転することと、ボクセル密度の非線形変換２９０のセットから得られた傾き及び非線形逆変換２５０のセットの反転傾斜を使用して、予備的変換物体密度更新データ２３６を再スケーリングすることによって、変換物体ボクセル密度更新データ２４１を計算することと、変換物体ボクセル密度更新データ２４１を先行する反復からの先行するボクセルデータ推定値２９１と累算して変換物体ボクセル密度推定値２４６を形成することと、変換物体ボクセル密度推定値２４６を変換することによって、生のボクセル密度データ２５１を計算することと、ローパスフィルタを使用して生のボクセル密度データ２５１を平滑化して、予備的ボクセル密度推定データ２５６を取得することと、物体構造情報を使用して予備的ボクセル密度推定データ２５６を再推定することによって、ボクセル密度推定データ２６１を計算することと、ボクセル密度推定データ２６１をデータストレージデバイス１８０に格納することと、ボクセル密度推定データ２６１を単位遅延時間だけ遅らせて、前のボクセル密度推定データ２６６を取得することと、ローパスフィルタを使用して前のボクセル密度推定データ２６６を平滑化して、平滑化された前のボクセル密度推定データ２７１を取得することと、平滑化された前のボクセル密度推定データ２７１を変換して、先行するボクセルデータ推定値２９１を取得することと、断層撮影投影アルゴリズムを使用して、前のボクセル密度推定データ２６６から予測投影データ２７６を計算することと、ローパスフィルタ２８５を使用して予測投影データ２７６を平滑化して、平滑化された予測投影データ２８６を取得することと、平滑化された予測投影データ２８６を変換して予測変換ピクセルデータ２１１を取得することと、予測投影データ２７６と投影ピクセルデータ２０１を比較することによって、再構成済み画像誤差データ１６０を決定することと、再構成済み画像誤差データ１６０をデータストレージデバイス１８０に格納することとによって、投影空間データを処理してもよく、又はソフトウェアは、これらのことを実行するためのソフトウェアコードを含んでもよく、又は方法は、これらのことを実行するためのステップを含んでもよい。
【００３８】
[0038]本発明の実施形態の説明において、及び図２Ａ及び図２Ｂにおいて、投影空間データ１５５は、「測定されたピクセルデータ１５５」、「物体画像化データ」、「生の測定データ入力１５５」、及び「入力データ１５５」とも呼ばれる。
【００３９】
[0039]再構成済み画像誤差データ１６０は、「残差データ１６０」、「投影残差１６０」、「画像データの残差１６０」、「イノベーションの残差１６０」、「データ１６０」、及び「残差のセット１６０」とも呼ばれる。
【００４０】
[0040]画像再構成器１６５は、本明細書では、「内側フィードバックループ１６５」及び「画像再構成コンピュータ」とも呼ばれる。
【００４１】
[0041]ボックス２００の変換を使用して投影方向に沿って累積物体ボクセル密度によって生成された投影ピクセルデータ２０１は、本明細書では、「投影積分物体密度２０１」、「投影空間密度データ２０１」、「入力信号データ２０１」、「近似的投影データ」、「入力データ２０１」、「投影２０１」、「データ２０１」、「投影データピクセルセットｓ（ｉ）２０１」、及び「入力投影ピクセル」とも呼ばれる。
【００４２】
[0042]変換２０５は、本明細書では「入力変換関数」とも呼ばれる。複数の入力変換関数（又は１つの入力変換関数）２０５は、線形又は非線形であってもよい。
【００４３】
[0043]測定された変換ピクセルデータ２０６は、本明細書では、「測定された変換ピクセル２０６」、「ピクセルデータ２０６」、「変換積分投影物体密度２０６」、「変換投影データ２０６」、「イノベーションデータ２０６」、及び「変換投影物体密度２０６」とも呼ばれる。
【００４４】
[0044]変換２１０は、本明細書では「参照変換関数」とも呼ばれる。複数の参照変換関数（又は１つの参照変換関数）２１０は、線形又は非線形であってもよい。
【００４５】
[0045]予測変換ピクセルデータ２１１は、本明細書では、「変換予測２１１」、「予測データ２１１」、「予測変換ピクセル２１１」、及び「変換予測データ２１１」とも呼ばれる。
【００４６】
[0046]第１のピクセルイノベーション結果データ２１６は、本明細書では、「再構成画像誤差データ２１６」、「画像誤差データ２１６」、「画像データの残差２１６」、「イノベーションの残差２１６」、「データ２１６」、「イノベーション２１６」、「イノベーションデータ２１６」、「残差データ２１６」、及び「残差２１６」とも呼ばれる。
【００４７】
[0047]第２のピクセルイノベーション結果データ２２１は、本明細書では、「誤差データ２２１」、「イノベーションデータ２２１」、「イノベーションの残差２２１」、及び「データ２２１」とも呼ばれる。
【００４８】
[0048]第３のピクセルイノベーション結果データ２２６は、本明細書では、「イノベーションの残差２２６」、「データ２２６」、及び「イノベーションデータ２２６」とも呼ばれる。
【００４９】
[0049]予備的ボクセル更新データ２３６は、本明細書では、「予備的変換物体密度更新データ２３６」、「予備的更新データ２３６」、及び「予備的変換物体更新データ」とも呼ばれる。
【００５０】
[0050]変換物体ボクセル密度更新データ２４１は、本明細書では「更新データ２４１」とも呼ばれる。
【００５１】
[0051]変換物体ボクセル密度更新推定値２４６は、本明細書では、「変換密度推定値２４６」、「変換ボクセル密度推定値２４６」、及び「データ２４６」とも呼ばれる。
【００５２】
[0052]生のボクセル密度データ２５１は、本明細書では「生のボクセル密度推定データ２５１」とも呼ばれる。
【００５３】
[0053]再構成済み画像データ２６１は、本明細書では、「データ２６１」、「推定物体密度２６１」、「ボクセル密度推定データ２６１」、「物体密度データ２６１」、「再構成済み物体２６１」、「物体データ２６１」、「物体空間画像」、及び「画像データ２６１」とも呼ばれる。
【００５４】
[0054]既知の参照モデルと計算された物体密度２６１の間の不一致は、ＭＭＲ−２６１と呼ばれる。
【００５５】
[0055]投影２７６は、本明細書では、「予測投影データ２７６」、「予測空間データ２７６」、「予測投影空間ピクセル２７６」、「データ２７６」、及び「予測フィードバック測定ピクセル２７６」とも呼ばれる。
【００５６】
[0056]平滑化された予測投影データ２８６は、本明細書では「データ２８６」と呼ばれてもよい。
【００５７】
[0057]図２Ａ及び図２Ｂに示されている画像再構成器１６５における高い処理性能は、次の３つのステップのプロセスとして要約することができる。
【００５８】
[0058]第１に、一部の処理ステップでデータ処理の線形化を使用し、線形反転断層撮影手法を効率的に使用して、理想的には物体に関連するすべての空間周波数成分を捕捉できるようにする、フィードバックループ内の投影ピクセルデータ２０１からボクセル密度推定データ２６１へのフィードフォワードデータの反転処理セクション。正のボクセル密度推定データ２６１を実現するために、小信号ボクセルゲイン係数、及び中間の変換物体ボクセル密度推定値２４６の計算を可能にする、非線形変換が含まれる。
【００５９】
[0059]第２に、正の物体ボクセル密度推定データ２６１を計算して正の予測投影データ２７６を生成する、フィードバックループのデータ投影セクション。送信器１１５から開始して、データ２０１までの、対応する物体画像化プロセスを近似することによって、再構成の品質を考慮する投影予測器の成分Ｈは、物体投影プロセスの近似を計算する。反復ごとの最適化は、高効率であり、数値的な山登り法ではなく、フィードバックループの理論的特性に基づいている。
【００６０】
[0060]第３に、未知のシステムの特性をより良く近似するために、画像収集システムにおける不確実性を表す、ベクトルｐに含まれているパラメータのセットが選択され、それらの不確実性の影響を補正するように調整されてもよい。例えば、期待される散乱又はビームハードニングの計算における不確実性は、最初は、投影物体密度を推定するための、測定されたＸ線強度の単純な対数変換によって考慮されず、補正されたベクトルパラメータｐを計算することによって考慮され、調整されてもよい。ベクトルパラメータｐは、本明細書では「パラメータベクトルｐ」とも呼ばれる。ｐの成分の調整は、例えば、測定データの前処理２００において補正を適用して、特定の物体データに関してビームハードニングを考慮するように、対数変換の基数を変更し、２００及び２７５を使用してＸ線散乱を考慮するように、予測信号強度データ値を変更する。性能の尺度が、外部から制御可能な固定されたパラメータ成分の固定されたセットに関して、第１のイノベーションプロセス２１６、投影残差１６０、又は画像再構成コンピュータの物体データの期待値のうちの少なくとも１つから得ることができる。
【００６１】
[0061]同様に、例えば、ラジアル圧縮センシング（ＣＳ：ｃｏｍｐｒｅｓｓｅｄ  ｓｅｎｓｉｎｇ）磁気共鳴画像法（ＭＲＩ）では、定量的に不十分に表された、ケース固有の物理データ収集及び処理特性が含まれることがあり、ベクトルパラメータｐにおいて補正が計算されることがある。この補正プロセスは、物体密度データ２６１を予測投影データ２７６まで進めるループのフィードバックセクションにおいて、発生することがある。このプロセスは、データ１１５の効果を繰り返し、物体画像化データ１５５をボックス２００に転送して、近似的投影データ２０１を生成する。データ収集プロセスの表現及びループのフィードバックセクションにおけるベクトルｐ内の対応する未知のパラメータの調整は、ベクトルｐ内の未知のパラメータの効率的で正確な調整を可能にする。パラメータの変更の効果は、例えば一部分において、ボックス２７５の投影行列Ｈ内の係数で表されてもよい。
【００６２】
[0062]画像再構成器１６５で、測定されたピクセルデータ１５５から予測データ２１１までのプロセスの反復と共に、増加するグリッド分解能及びベクトルパラメータｐ内の固定パラメータを使用して、ＦＢＰの再構成時間に近い画像データ２６１の再構成時間を実現できる。この速度によって、ｐに含まれているシステムパラメータの反復改良が実現可能になる。例えば、画像データの残差１６０又は２１６が、インベスティゲータ１７０に送信される。インベスティゲータ１７０は、例えば高効率のＬｅｖｅｎｂｅｒｇ-Ｍａｒｑｕａｒｄｔプロセスで、例えば画像データの残差１６０、２１６の二乗和、及び利用可能な場合は、既知の物体と計算された物体データの間の誤差の最小化を使用して、ｐ内のパラメータの補正値を計算することができる。Ｌｅｖｅｎｂｅｒｇ-Ｍａｒｑｕａｒｄｔ法は、通常、他の手法によって使用される最急降下法及び関連するアプローチのｋ倍以上の時間がかかる（ｋは、調整するべき未知の数値である）。ｐ内の改良されたパラメータは、画像再構成器１６５に戻される。ｐに適用されるデータ処理のこの第２の層の反復は、物体データ２６１の大きな改善につながる。このプロセスの結果の例が、図４に示されている。
【００６３】
[0063]一実施形態では、図２は画像再構成器１６５の機能を含んでもよい。例えば、図１の測定されたピクセルデータ１５５は、図２に示されているように、入力として供給されてもよい。測定されたピクセルデータ１５５は、実験的取得１１０、シミュレーション１８５、又は以前の実験又はシミュレーション或いはその他のデータソースから記録された投影ピクセルデータを含んでいるデータストレージデバイス１８０から取得されてもよい。測定されたピクセルデータ１５５は、画像再構成器１６５に、直接、又はネットワーク（例えば、ローカルエリアネットワーク（ＬＡＮ）、広域ネットワーク（ＷＡＮ）、インターネットなど）を介してリモートから提供されてもよい。一実施形態では、画像再構成器１６５は、測定されたピクセルデータ１５５を処理して、再構成済み物体２６１を生成してもよい。ボックス２００で、測定されたピクセルデータ１５５が処理されて、投影ピクセルデータ２０１を取得する。測定されたピクセルデータ１５５は、その後のデータ処理の近似の精度に対応できる変換によって処理される。例えば、Ｘ線画像化では、対応する物体密度投影の推定値の予備的近似として取得するために、光子強度が負の対数によって変換されてもよい。この処理は、例えば、ビームハードニング、散乱、又は画像記録システムの欠陥の影響に関して、投影密度の推定値を調整するパラメータベクトルｐの成分で表された変換を改良し、２７５内でのこれらの改良を推定する必要性を減らすことをさらに含んでもよい。ラジアルスキャニング（ＲＳ：ｒａｄｉａｌ  ｓｃａｎｎｉｎｇ）圧縮センシング（ＣＳ）ＭＲＩでは、ボックス２００は、フーリエ変換又は主成分法を使用して、測定されたピクセルデータ１５５を物体密度投影の空間密度データの予備的近似に変換してもよく、ボックス２００は、例えば、パラメータベクトルｐで表された変換を改良し、幾何学的電磁界歪及び設計された励起及び応答の測定条件からの逸脱に関して調整することを、さらに含んでよい。ボックス２００における変換プロセスの詳細は、近似的な物体密度の投影密度を取得するために使用される特定の物理的画像化プロセスを近似するように適応されてもよい。
【００６４】
[0064]図２Ａ及び図２Ｂに示されている画像再構成器１６５における処理を理解するための別の方法は、次のように要約することができる。少なくとも行列Ｚ（ボックス２２０）又は行列Ｌ（ボックス２４０）が決定され、期待される状態データ２７６を使用して状態データ２０１を処理するために使用されてもよく（例えば、変換ピクセルイノベーション２２１を取得するために行列Ｚ（ボックス２２０）を使用するか、又はボクセルの更新２４１を取得するために行列Ｌ（ボックス２４０）を使用する）、この処理は、期待される状態データ２７６、行列Ｚ（ボックス２２０）を決定するためのｆ１１（ボックス２０５）及びｆ１２（ボックス２１０）の関数セット、並びに行列Ｌ（ボックス２４０）を決定するためのｆ２１（ボックス２９０）及びｇ２２（ボックス２５０）の関数セットに依存する。
【００６５】
[0065]ボックス２０５で、物体の内部構造をエンコードする情報を含んでいる投影ピクセルデータ２０１が、ピクセルデータ２０６に非線形（又は線形）に変換されてもよい。ボックス２０５で、要素ｆ１１（例えば、ｆ１１（ｓ）＝ｌｏｇ（ｓ）であり、すべてのピクセル値ｓ＞ｓｏ＞ｅｐｓについて、ピクセルごとにそれぞれｆ１１（ｓｏ）＞ｃｍｉｎである）を使用する通常は非線形な入力変換行列が、ピクセルごとに使用され、変換領域内の１つ又は複数の測定された変換ピクセル値２０６を計算してもよい（例えば、ｆ１１は、例えばビームハードニングを表す、パラメータの関数であってもよい）。ボックス２０５で、データが密度領域から変換領域に変換される。別の実施形態では、ボックス２０５で、ｆ１１が、近似的に、入力データ２０１の分散安定化関数として計算されてもよい。例えば、ボックス２０５で、ｆ１１が、近似的に、ＰＥＴ測定から取得された入力信号データ２０１のポアソン分散安定化平方根関数であってもよい。別の例として、例えば測定ノイズが一定である場合、ｆ１１は、線形関数又は線形関数のセットであってもよい。
【００６６】
[0066]ボックス２１０で、非線形（又は線形）参照変換ｆ１２をピクセルごとに使用して、予測投影データ２７６が１つ又は複数の予測変換ピクセル２１１に変換されてもよい。例えば、非線形参照変換を使用する場合、ｆ１２は対数関数（例えば、ｆ１２＝ｌｏｇ（．））であってもよい。本明細書において使用されるとき、「ピクセル」は、ピクセルの位置又はピクセルの値のことを指してもよい。同様に、本明細書において使用されるとき、「ボクセル」は、ボクセルの位置又はボクセルの値のことを指してもよい。ボックス２１０で、データが密度領域から変換領域に変換される。非線形（又は線形）参照変換ｆ１２は、密度領域からの１つ又は複数の予測投影空間ピクセル２７６（又は、ボックス２８５でフィルタ（ローパスフィルタなど）によって平滑化される場合、平滑化された予測投影データ２８６）を使用して、変換領域内の１つ又は複数の予測変換ピクセル２１１を計算してもよい。例えば、ボックス２８５で、物体の投影２７６に対する２次元の平滑化特性は、例えば３次元においてフィルタ２７０を使用して同じ物体データを平滑化するときの投影の効果の近似である。この概念は、さらに高い次元の再構成にも適用される。言い換えると、投影及び平滑化の線形処理は交換可能である。予測投影空間ピクセル２７６は、前のボクセル密度推定データ２６６を使用して、予測器のＨ（ボックス２７５内）から取得されてよい。非線形（又は線形）参照変換ｆ１２は、ｆ１１のピクセルごとの近似になる傾向がある。ここで、例えば、堅牢（ロバストな推定処理を可能にするために、ｆ１１の要素が、データ２０１のアーチファクトに関する調整を必要とする場合、ｆ１１及びｆ１２は異なってもよい。別の例として、ｆ１２は、線形関数又は線形関数のセットであってもよい。
【００６７】
[0067]ボックス２１５で、変換領域内のイノベーションデータ２１６が計算されてもよい。ボックス２１５の一実施形態では、１つ又は複数の測定された変換ピクセル２０６が、予測変換ピクセル２１１（平滑化された予測投影データ２８６から非線形に変換されたデータ）と比較されて、イノベーションデータ２１６に対応する差異を生成してもよい。言い換えると、ボックス２１５で、対応する測定された変換ピクセル２０６と予測変換ピクセル２１１の間の差異が計算されて、第１のピクセルイノベーション結果データ２１６を生成してもよい。第１のピクセルイノベーション結果データ２１６が、残差分析、及びシステムパラメータの最適化に使用され、最終的に最良の物体密度再構成を生み出すことができる。ボックス２１５で、対応するピクセル位置の、（対応するピクセル位置の近傍領域内の）関連付けられた密度領域のピクセル値のボックス２０５からの非線形入力変換ｆ１１の反転傾斜及びボックス２１０からの非線形参照変換ｆ１２の反転傾斜の、重み付きの組み合わせｚのそれぞれが、個別に計算されてもよい。ノイズが多いか、又は不正確なデータ２７６、及び特にノイズが多いデータ２０１では、特定の密度領域のピクセル値に対応する傾きが、空間的な連続的に隣接する値を含んでもよく、ロバストな関数を使用した再評価を受けてもよい。
【００６８】
[0068]ボックス２２０で、第１のピクセルイノベーション結果データ２１６が、ピクセルごとのマッチング行列Ｚによってスケーリングされて、誤差データ２２１を生成してもよい。ピクセルごとのマッチング行列Ｚは、現在の投影空間データ及びそれらの変換の特性から計算されてもよい。ボックス２２０のスケーリングマッチング行列Ｚは、各ｚが、関連付けられたピクセル位置に対応するように、ｚのセットを編成することによって計算されてもよい。この計算は、対応する入力ピクセル値及び参照ピクセル値に関する非線形入力変換関数２０５の反転傾斜及び非線形参照変換関数２１０の反転傾斜を使用してもよい。行列Ｚの係数は、１つ又は複数の測定された変換ピクセル２０６のピクセル値でのピクセル変換関数ｆ１１及び予測変換ピクセル２１１のピクセル値でのｆ１２の組み合わせられた反転傾斜の表現であってもよい。反転傾斜ｚの組み合わせは、ノイズがない場合には、最適に近い値ｚ−ｏｐｔを生成する可能性が高いプロセスを使用し、ノイズがある場合には、ｚ−ｏｐｔに近い近似値を生成する可能性が高くなるように、反転傾斜を組み合わせる。ボックス２２０で、第２のピクセルイノベーション結果２２１が、第１のピクセルイノベーション結果データ２１６と、ピクセルイノベーションスケーリング行列Ｚの対応する要素の、ピクセルごとの積を使用して計算されてよい。ピクセルごとのイノベーションスケーリング行列Ｚは、対応する非線形入力ピクセル密度変換関数２０５の反転傾斜及び対応する非線形参照ピクセル密度変換関数２１０の反転傾斜をピクセルごとに使用して計算されてもよい。非線形入力ピクセル密度変換関数２０５は、範囲が制限されるか、又はロバスト（例えば、制約された値）であってもよい。非線形参照ピクセル密度変換関数２１０は、範囲が制限されるか、又はロバスト（例えば、制約された値）であってもよい。
【００６９】
[0069]ボックス２２０で、ピクセルの連続的な特性に基づく米国特許第８，６６０，３３０号の手法と同等の行列ゲインの改良及び隣接する第１のピクセルイノベーション結果データ２１６が、任意選択的に統合されてもよい。
【００７０】
[0070]任意選択的なボックス２２５で、ピクセルの連続的な特性に基づく米国特許第８，６６０，３３０号の手法と同等の第２の直列なイノベーションゲイン行列Ｇ及び隣接するピクセルイノベーション結果が、第２のピクセルイノベーション結果データ２２１を第３のピクセルイノベーション結果データ２２６に変換してもよい。ボックス２２５で、ピクセルごとの補正ゲイン行列Ｇの要素が、第２のピクセルイノベーション結果データ２２１及び第３のピクセルイノベーション結果データ２２６内の対応するピクセルに関連付けられてもよい。第２のピクセルイノベーション結果データ２２１は、例えば、米国特許第８，６６０，３３０号のゲイン評価手法に対応する、第１のピクセルイノベーション結果データ２１６又は第２のピクセルイノベーション結果データ２２１の連続的な特性に基づくゲイン評価手法を使用して、第３のピクセルイノベーション結果データ２２６の計算に使用されてもよい。
【００７１】
[0071]ボックス２３０で、イノベーションプロセスのいずれかにおいて連続的な特性を使用して、米国特許第８，６６０，３３０号のゲイン評価手法に対応する生の補正ゲインが決定されてもよい。例えば、ボックス２３０で、重み付け関数を使用して、イノベーションゲインの１つ又は複数の空間的に隣接する係数又は先行する反復からのイノベーションゲインの１つ又は複数の連続的な係数に対する回帰によって、ピクセルごとの補正ゲイン行列２２５が計算される。補正ゲイン行列２２５は、測定及びモデルの欠陥を表す。回帰の前に、イノベーションプロセッサ２３０において、第１のピクセルイノベーション結果データ２１６又は第２のピクセルイノベーション結果データ２２１のうちの少なくとも１つを使用して、生の補正係数が計算され、イノベーションパターンＡ又はＢを識別し、係数の対応する生の補正の増加又は生の補正の減少を伴って変動する。
【００７２】
[0072]ボックス２３５で、第３のピクセルイノベーション結果データ２２６が、断層撮影再構成アルゴリズムに提供され、予備的ボクセル更新データ２３６を取得してもよい。断層撮影再構成アルゴリズムは、例えば、線形再構成アルゴリズム（ＬＲＡ：ｌｉｎｅａｒ  ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  ａｌｇｏｒｉｔｈｍ）であってもよい。例えば、Ｚｅｎｇ他［１］のＤＣ／平均投影値を使用する断層撮影再構成反転アルゴリズムが使用されてもよい。別の例として、断層撮影再構成アルゴリズムは、逆投影法であってもよい。通常、システムパラメータのセットによって特徴付けられる、必要なＤＣ／平均投影値と他の周波数成分の間の関係は、利用できる投影データの量及び方向に依存することがある。例えば、２つの（ノイズがない）直角投影の場合、本発明による再構成は、再構成のための投影のフィルタリングを必要としなくてもよい。
【００７３】
[0073]ボックス２４０で、予備的変換物体密度更新データ２３６が、ボクセルごとの更新スケーリング行列Ｌの要素を使用して再スケーリングされてもよく、各変換物体ボクセル密度更新データ２４１が行列Ｌの要素に関連付けられてもよい。ボクセルごとの更新スケーリング行列Ｌは、変換された物体データ及び正の制約変換物体ボクセルデータ変換の特性から計算されてもよい。行列Ｌの要素の値は、対応する変換ボクセル密度推定値２４６の、ボックス２５０のｇ２２のボクセルごとの要素を使用して、ボクセルごとの逆変換の要素の反転傾斜から計算されてもよい。
【００７４】
[0074]ボックス２４５で、変換密度推定値２４６が、変換物体ボクセル密度更新データ２４１を対応する変換先行ボクセルデータ推定値２９１に加算することによって、計算されてもよい。変換密度推定値２４６は、インベスティゲータ１７０で後で使用するために、データストレージデバイス１８０に格納されてもよい。
【００７５】
[0075]ボックス２５０で、変換密度推定値２４６が、密度領域に逆変換され、生のボクセル密度推定データ２５１を取得してもよい。ボックス２５０で、データが変換領域から密度領域に変換される。変換密度推定値２４６が、行列要素ｇ２２を使用してボクセルごとに非線形に変換され、生のボクセル密度推定データ２５１を生成してもよい。一実施形態では、ボックス２５０で、行列要素ｇ２２は、行列ｇ２２の主要な部分に関して、ｆ１１又はｆ１２行列要素の１つ又は複数の関数タイプの近似逆関数としての関数タイプであってもよく、ｇ２２の行列要素は、変換密度推定値２４６のデータ範囲の少なくとも１つの領域内の外部入力を満たす出力値を計算するために使用されてよい。ボックス２５０の関数は、本明細書では「正の制約関数」と呼ばれてもよい。
【００７６】
[0076]ボックス２５５で、生のボクセル密度推定データ２５１がフィルタリングされて、予備的ボクセル密度推定データ２５６を取得してもよい。ボックス２５５は、ローパスフィルタを使用して、生のボクセル密度推定データ２５１を平滑化してもよい。ボックス２５５は、生のボクセル密度推定データ２５１をフィルタリングして、予備的ボクセル密度推定データ２５６を正として出力してもよい。ボックス２７０及び２８５で、ローパスフィルタが、出力を正に維持する特性を共有してもよく、ボックス２５５の関数を置き換えてもよい。
【００７７】
[0077]ボックス２６０で、物体の予備的ボクセル密度推定データ２５６が、最初の反復では単一のボクセルに対して、その後の反復では多数のボクセルに対して動作するＰによって後処理され、ボクセル密度推定データ２６１の特性に関する適切な先験的情報を使用して、密度値を改良してもよい。この手法は、単一のボクセルから開始することができ、最初に複数のボクセルを必要としない。予備的ボクセル密度推定データ２５６を後処理することによって、（ａ）等しい又は増加したグリッド分解能、又は（ｂ）等しい又は増加した密度値の分解能のうちの１つ又は複数に基づいて、１つ又は複数のボクセルを作成してもよい。ボックス２６５で、ボクセル密度推定データ２６１が、前のボクセル密度推定データ２６６として記録され、表示され、ボックス２７０のフィルタへの入力として送信され、その後、ボックス２９０の変換ｆ２１に送信され、ボックス２７５の予測投影プロセッサＨへの入力としても送信される。
【００７８】
[0078]ボックス２７０で、前のボクセル密度推定データ２６６がフィルタリングされて、平滑化された前のボクセル密度推定データ２７１を取得する。ボックス２７０のフィルタは、ローパスフィルタであってもよい。
【００７９】
[0079]ボックス２９０で、平滑化された前のボクセル密度推定データ２７１が、非線形参照変換ファミリーｆ２１（例えば、ｆ２１＝ｌｏｇ（．））を使用して、先行するボクセルデータの推定データ２９１に変換されてもよい。ボックス２９０で、データが密度領域から変換領域に変換される。非線形参照変換ファミリーｆ２１は、本明細書では、フィードバック関数２９０と呼ばれてもよく、ボックス２９０内のフィードバック関数のセットと呼ばれてもよい。
【００８０】
[0080]ボックス２７５で、前のボクセル密度推定データ２６６が処理され、通常は単一のスパース投影行列Ｈの表現を使用して、予測投影データ２７６を取得してもよい。行列Ｈは、元の物体密度値の画像化を、物体投影密度値として繰り返すように設計されてもよい。しかし、図１の画像化プロセスにおける不確実性が、行列Ｈにおける不確実性をもたらすことがあり、画像化プロセスについてより多くのことが学習されたときに調整するために、ベクトルパラメータｐで表されてもよい。例えば、行列Ｈは、ベクトルｐでパラメータ化された、パターン化されたボクセル外投影成分を使用して、特定のボクセルに関する散乱を表してもよい。しかし、実際の処理では、行列Ｈの明示的な使用が回避されてもよい。
【００８１】
[0081]ボックス２７５の一実施形態では、機能は、例えばＺｈａｎｇ他［２］及びＬｏｎｇ他［３］において説明されている機能に関連していてもよい。
【００８２】
[0082]さらに、行列Ｈで要約されているボックス２７５内の投影プロセスは、例えば、少なくとも１つのシステムパラメータ又は物体パラメータを表してもよい。システムパラメータは、例えば、焦点の形状、焦点のビーム出口の強度及び硬度並びに出口の形状特性、変化する管の供給電圧、ビームに依存するＸ線検出器の特性、又はＸ線散乱のうちの１つ又は複数であってもよい。物体パラメータは、例えば、物体の動きを表してもよい。これらのプロセスの不確実性（例えば、システムパラメータ及び／又は物体パラメータ）は、ベクトルパラメータｐで要約されてもよい。画像再構成では、未知のボクセル値の数が、ベクトルｐ内の未知の係数の数をはるかに上回るため、ｐ自体のほんのわずかの成分が推定の対象になってもよい、ということに注意する。まばらな測定の状況では、ボクセル密度に対して制約を設け、例えばボックス２５５又はボックス２６０を使用することによって、ｐ及び物体のボクセルが推定の対象になってもよい。
【００８３】
[0083]ボックス２７５で、断層撮影投影アルゴリズムが、Ｈ（ｐ）のベクトルパラメータｐと共に使用されてもよい。ｐを変更せずに、図２示されている処理ループが収束した後に、パラメータｐが、性能基準に対して改善された解を見つけるように調整されてもよい。この性能基準は、例えば、イノベーションプロセス２１６の残差の最後のセット、残差のセット１６０、又は物体の計算された値と参照値の間の値の不一致に基づいてもよい。
【００８４】
[0084]ボックス２１０で、非線形変換ｆ１２が予測投影データ２７６に適用されて、変換予測データ２１１を取得してもよい。
【００８５】
[0085]ベクトルパラメータpは、データ２１６、２２１、２２６、又は１６０を、例えば重み付き二乗和の最小化に使用して、システムパラメータを、変化する投影条件に最適に適合させる。投影条件は、例えば、ボックス２７５での投影演算Ｈのパラメータ、ボックス２６０で要約された後処理法Ｐ、並びにボックス２０５での関数ｆ１１、ボックス２１０での関数ｆ１２、ボックス２７０での関数ｆ２１、及びボックス２５０での関数ｇ２２のパラメータで表される。
【００８６】
[0086]データ２４６及び２６１は、外部のシステム変数及び目的を満たすように変更されてもよい。例えば、データ２４６及び２６１に含まれている一部の物体密度変数を先験的に知られた固定値にリセットすると、残りの物体密度変数及びシステムパラメータのより迅速な収束を支援することができる。
【００８７】
[0087]ボックス２７５で、予測投影データ２７６が、１つ又は複数の物体空間のボクセルのセットに基づいて計算されてもよく、物体空間のボクセルのセットは、物体空間のボクセルのセットが投影空間に投影されるときにサイズが変化する、複数の分解能のグリッド（例えば、単一のグリッド点を含む）をカバーしてもよい。
【００８８】
[0088]ボックス２７５で、適切な（１つ又は複数の）ピクセルのセットに関するボクセルの投影値が計算される。残りの投影値は、必要な範囲について、他の情報源を使用して計算されてもよい。
【００８９】
[0089]ボックス２８０で、投影ピクセルデータ２０１及び予測投影データ２７６から残差データ１６０が生成されてもよい。残差データ１６０が測定法のノイズ（例えば、Ｘ線画像化の場合のホワイトノイズ）に近づき、説明できない投影物体密度の痕跡を（例えば、パラメータｐの係数が適切に決定され、図２Ａ及び図２Ｂで表されたループ内のデータの密度が十分の収束したため）無視できる場合、計算された残差データ１６０がインベスティゲータ１７０及び出力デバイス１７５に転送されてもよい。
【００９０】
[0090]ボックス２２５で、各イノベーションでの補正ゲイン行列が、１つ又は複数の生の隣接した先行するイノベーションゲインを重み付けする空間的な連続的重み付け関数を使用して計算されてもよく、それらのイノベーションゲインは、例えば米国特許第８，６６０，３３０号で使用されるループゲイン決定の特徴の少なくとも一部を使用する方法で計算される。図２の現在のフィードバックループへの米国特許第８，６６０，３３０号の特徴の実装は、イノベーション２１６又は２２１が連続して統計的に有意な同一の符号である場合に、ボックス２３０での生のゲイン係数を増やして、パターンＡを作成し、イノベーションが連続して統計的に有意な交互の符号である場合に、ボックス２３０での生のゲイン係数を減らして、パターンＢを作成する。このアプローチの異なるバージョンが可能である。より基本的なプロセスのバージョンは、特定のイノベーションパターンを有する生のゲイン係数の著しく不正な値から開始し、反復ｉごとに、有意な代替のパターンを検出するまで、ゲインを徐々に変更して、関連付けられたパターンの有意性を減らし、関連付けられた生のゲイン係数を前の値にリセットして、このプロセスを再び開始することである。ボックス２３０での先行するイノベーションプロセスを使用するボックス２２５でのゲイン調整又はその他の適切なメカニズムの選択肢は、システムノイズ、例えば密度領域の物体空間内のフィルタ効果、逆投影フィルタ特性の誤差、ロバストな影響関数の使用、単位ＤＣループのフィードバックゲインの逸脱、及びその他のモデルの欠陥などから生じる、反復ループ計算内の準最適なゲイン決定の影響を減らすために提供される。例えば、第１の影響関数のセットを使用して改良された第１のピクセルイノベーション結果データ２１６（Ｓｕｎｎｅｇａｒｄｈ［１５］を比較する）、第２の影響関数のセットを使用して改良された第２のピクセルイノベーション結果データ２２１、及び第３の影響関数のセットによって改良された第３のピクセルイノベーション結果データ２２６に関して、イノベーションが、重み付け関数を使用して補正ゲイン行列の係数を計算してもよく、この重み付け関数は、測定及びモデルの欠陥を考慮するために、１つ又は複数の空間的に隣接するイノベーション行列の係数又は先行する反復からの１つ又は複数の連続的なイノベーション行列の係数に対する回帰によって計算されてもよい。１つ又は複数の連続的なイノベーション行列の係数は、イノベーションゲインの１つ又は複数の連続的な係数と呼ばれてもよい。
【００９１】
[0091]ボックス１７０で、対象の物体が、ｐを制御するインベスティゲータによって、図２のループ内のｐのセットごとの最後の反復によって作成された１６０、２１６、２２１、又は２２６の重み付きの収束したイノベーションの残差を使用して、反復的に改良されてもよい。重み付きのイノベーションの残差の各セットが、図１に示され、図２Ａ及び２Ｂで要約されている内側のフィードバックループ１６５内で、十分適切に収束した反復から取得される場合に、外側のループ内で計算されたパラメータｐの摂動が、最適化に使用されてもよい。
【００９２】
[0092]外部から制御可能なパラメータｐの成分が、初期推定値から、イノベーションの残差１６０、２１６、２２１、２２６、或いは例えば、シミュレーションにおける、又は参照物体の密度に対して物体密度推定値の収束を評価するときなどの、既知のシステムの参照モデルからの逸脱に関して、外部の目的又は最小コストとの最良の一致を形成する推定値に、徐々に変化してもよい。画像再構成コンピュータ［１６５］の外部にあるインベスティゲータコンピュータ［１７０］は、既知の参照物体と対応する計算された物体密度の間の差異を特徴付ける投影残差又は物体の残差のうちの少なくとも１つの特性を使用しながら、パラメータベクトルの制御可能な成分を、初期設定から、外部の目的に対して最低のコストを生成する設定に徐々に変更するよう機能してもよい。
【００９３】
[0093]投影残差は、先験的なパラメータ或いは近傍領域を表すパラメータ及び／又は近傍領域を使用して得られるパラメータを使用した、物体のシミュレーション及び課題の記録によって詳しく調べられた、ロバストな推定影響関数を使用して、処理されてもよい。これらのパラメータは、例えば、局所的な入力測定値、システムの特性、予測値、又は先験的に期待される物体及び測定の特性のうちの１つ又は複数の関数であってもよい。
【００９４】
[0094]投影残差は、選択範囲上の残差の平滑化、選択範囲上の残差のスケーリング、又は選択範囲上の影響関数による重み付けのうちの少なくとも１つを使用して処理されてもよい。
【００９５】
[0095]イノベーションデータ２１６、２２１、２２６、予備的更新データ２３６、及び更新データ２４１が、影響関数を使用して処理されてもよい。測定されたピクセルデータ１５５が、影響関数を使用して処理されてもよい。同様に、既知の物体密度再構成のシミュレーション中に、例えば、不一致の残差ＭＭＲ−２６１が、影響関数を使用して処理されてもよい。
【００９６】
[0096]図２Ａ及び２Ｂのプロセスの実施形態は、例えば、線形再構成アルゴリズム（ＬＲＡ）（例えば、フィルタ補正逆投影（ＦＢＰ）又はその他の資格のある、線形であることが好ましい、再構成アルゴリズム）の能力を使用して、断層撮影の反転問題の固有値に効率的に対処する、最適化された再帰法によって拡張された高効率のコンピュータ断層撮影（ｅＨＥＣＴＯＲ：ｅｘｔｅｎｄｅｄ  ｈｉｇｈ  ｅｆｆｉｃｉｅｎｃｙ  ｃｏｍｐｕｔｅｄ  ｔｏｍｏｇｒａｐｈｙ  ｗｉｔｈ  ｏｐｔｉｍｉｚｅｄ  ｒｅｃｕｒｓｉｏｎｓ）を含んでもよい。このために、ＬＲＡは、２２０、２２５、及び２４０でピクセルごとの小信号データのゲインを使用して線形化されてもよい、非線形構造に埋め込まれてもよい。例えば、ＦＢＰが、単一のステップで、良好だが不完全な物体密度推定値を生成するのと同様に、本明細書では、ＦＢＰが、変換推定問題の線形化モデルに対して再帰的に使用されてもよい。正味の効果は、図２のこの線形化ループが、推定誤差の幾何学的縮小を含む縮小傾向が大きいマッピングを表すということである。この収縮は、測定ノイズ又はモデルの誤差などの残差が、それ以上の反復によって解決できない不一致を引き起こす場合に、終了する。そのような不一致は、投影残差において発生することがあり、ベクトルパラメータｐの再推定のための原動力になる。
【００９７】
[0097]線形化されたｅＨＥＣＴＯＲの場合、理想的には、線形カルマンフィルタが使用されてもよい。しかし、断層撮影再構成の場合、この最適なカルマンフィルタの決定は、最新の計算リソースの能力を超える。ボックス２３５のこのフィルタでは、最適なカルマンゲインを使用する代わりに、ＬＲＡが、現在のメモリレスのカルマンフィルタゲインを近似するために使用される。
【００９８】
[0098]最適なカルマンフィルタを使用して線形フィルタリングをモデル化するときに、一連のイノベーションが、ホワイトノイズプロセスを形成する。イノベーションプロセスの分散は、メッセージモデル及び観測ノイズのランダムな原動力／データから生じる（Ｓａｇｅ他［４］の２６８ページを参照）。しかし、準最適なアプローチである、ｅＨＥＣＴＯＲの非線形ループ設定の場合、中程度の相関関係の連続的なイノベーションが生成される。そのため、ｅＨＥＣＴＯＲのループ構造は、最適な最小二乗アルゴリズム（例えば、最適なカルマンフィルタ）から、いくつかの方法で逸脱している。例えば、ＬＲＡは、最適なフィルタゲインの準最適な近似である。さらに、小信号データの線形化は、最適なフィルタの適用の場合にのみ近似になり、イノベーションプロセスはホワイトプロセスにならない。それでも、中程度の密度データ値の変化の場合、適切に選択されたＬＲＡ（現在のフィルタ補正逆投影断層撮影などで使用されるＬＲＡなど）は、反復ごとに、縮小傾向が大きいマッピングを提供する。例えば、グリッド分解能が徐々に増加する場合、中程度の密度値の変化を期待することができ、数回の反復での非線形フィルタの収束が可能になる。反復の計算コストが、幾何学的合計において蓄積し、最後の反復が、例えば単一の逆投影又は反転のコスト成分に近い、支配的なコスト成分を決定する。
【００９９】
[0099]グリッド分解能が増加するときに、イノベーションプロセスでの空間的な高周波成分が、より顕著になることがある。Ｘ線断層撮影では、例えば、物体密度推定値の収束において、これらの空間的な高周波投影成分が、空間的測定のホワイトノイズに近いノイズを生成する。連続的にグリッド分解能が増加するときに、イノベーションの相関関係は中程度のままである。しかし、これら及び近似のその他の態様の場合、現在のループ設定は、断層撮影再構成の実用的且つ効率的な計算にとって十分であり、図２Ａ及び２Ｂで詳細に示されて
いる画像再構成器１６５において使用される計算の基礎を形成する。
【０１００】
[0100]本発明の一実施形態では、ボックス１６５での反復の連続的プロセス、ボックス２５５でのローパスデータフィルタの無視、小さいイノベーションを伴うグリッド分解能の増加、及び近似的で準最適な、明示的に計算されないＬＲＡのゲインの使用が、状態変数の固定された次元数、明示的に計算されるゲインＫ、及び有限分散を含む適切に指定された前の統計値を使用する最適な線形カルマンフィルタと比較されてもよい。
【０１０１】
[0101]これらのカルマンフィルタの条件のいずれも、図２の計算を正確に表さない。それでも、この計算モデルは、図２の反復プロセスの定性的評価を可能にする。
【０１０２】
[0102]理想的な状況の場合に期待される性能の比較及び推定のために、表１に、Ｓａｇｅ他［４］の２６８ページに示されている最適な離散カルマンフィルタアルゴリズムを編集して示している。
【表１】
2019535405.tif 000003
【０１０３】
[0103]カルマンフィルタの比較の線形化部分は、次の方法でモデル化される。簡単にするために、ノイズ駆動型の動的システムの成分（表１）が存在しない静的システムについて考え、ＶＷ（ｉ＞０）＝０とし、ただし初期分散ＶＸ（０）＝ＶＷ（０）＝ＶＭ≠０が、初期モデル状態共分散行列を表すとする。このモデルの状態ｘｍの推定は、先行する状態の推定と、カルマンゲインＫ（ｉ）によって重み付けされたイノベーションの組み合わせである（Ｓａｇｅ他［４］の２６８ページを参照）。ここで、独立した測定ノイズの寄与の新しいセットではなく、測定成分の同じ固定されたセットを使用しているにもかかわらず、相関関係のない観測ノイズＶＶ（ｉ）が仮定される。連続的な相関関係は、ノイズデータが、その後、相互に一致しないことがあるため、例えば、平滑化処理、グリッド分解能の同様の増加、及び冗長な測定を考慮して、無視される。そのため、連続的な相関関係は、物体密度において（ほとんど）表すことができず（ただし、損傷を与える画像ノイズを除く）、したがって、状態の更新は、ホワイトノイズの存在における計算ステップと同じ計算ステップを使用してもよい。このモデルは、物体１２５が存在しないときのノイズ値の固定されたセットの再構成が、物体１２５の最終的な存在と比較したときに、無視できるボクセル密度を生成する程度まで、正当化することがでる。実際、ボックス２５５、２６０、２７０、及び２８５は、わずかなフィルタ効果を使用して、ボクセルデータ値における高周波成分を抑制し、それに応じて、例えば「物体投影の不一致」に起因する、大きい投影残差値をそのままにすることができる。最初は必ずしも最適でないカルマンフィルタモデルの線形化部分は、物体密度状態ベクトルｘに関して、次のように表すことができる。
ｘ（ｉ）＝ｘ（ｉ−１）＋Ｋ（ｉ）ｘ［ｚ（ｉ）−Ｈ（ｉ）ｘ（ｉ−１）］    ［１．１］
Ｋ（ｉ）＝ＶＸ（ｉ）ＨＴ（ｉ）ＶＶ−１（ｉ）    ［１．２］
ここで、
ＶＷ（０）は、状態の分散ＶＸ（ｉ）に含まれている初期物体密度の分散を表す    ［１．３］
ｚ（ｉ）＝Ｈ（ｉ）ｘ（ｉ）＋ｖ（ｉ）  投影測定ベクトル    ［１．４］
ここで、
Ｈ（ｉ）  物体密度状態ベクトルｘ（ｉ）の観測行列    ［１．５］
ｖ（ｉ）は、観測／測定ノイズである。    ［１．６］
【０１０４】
[0104]システムノイズＶＷ＝０及び一定な相関関係のない観測ノイズＶＶ（ｉ）によって駆動されるメッセージ物体密度モデルの場合、ｉが無限大になるにつれて、カルマンゲインＫ（ｉ）が０に収束する。
Ｋ（ｉ）→０    ［１．７］
観測ノイズの固定パターンを前提として、固定されたグリッド分解能での反復中に、対応するｅＨＥＣＴＯＲにおいてデータ処理の線形化の動作点を決定している間は、次の式をそのままにすることで、十分であることがある。
Ｋ（ｉ＋１）＝Ｋ（ｉ）    ［１．８］
【０１０５】
[0105]物体のボクセル数及び投影ピクセル数が大きいＣＴでは、数値的な理由のため、カルマンゲインフィルタシステム（特に、カルマンゲイン）は、計算不可能になる。しかし、同等の特性を有するＬＲＡの近似を利用できる。最適なカルマンゲインに対するそのような近似は、Ｒａｄｏｎ［５］の設定、Ｗｏｏｄ他［１４］による小規模システムの最小分散の実装、又は最適なカルマンフィルタのような１ステップの反転を提供しないが、大規模システムの場合に、良好な縮小傾向のＬＲＡが得られる。例えば、Ｒａｄｏｎのアプローチ［５］によって動機付けられた近似を使用して、ＦＢＰに類似する縮小傾向のマッピングが得られる（ただし、Ｚｅｎｇ他［１］を参照）。
【０１０６】
[0106]図２は、例えば、対数関数ｆ１１（ボックス２０５）、対数関数ｆ１２（ボックス２１０）、対数関数ｆ２１（ボックス２９０）、指数関数ｇ２２（ボックス２５０）、イノベーションスケーリング行列Ｚ（ボックス２２０）、及び更新スケーリング行列Ｌ（ボックス２４０）を使用する特別な場合も示している。ボックス２２０及び２４０の場合、行列Ｚ及びＬはそれぞれ、対角優位行列によって表されてもよく、各対角要素は、データストリーム内の単一のピクセル又は単一のボクセルなどのデータ要素に対応する。図２Ａのボックス２２０の反復ｉで、行列Ｚ（ｉ）の要素は、第１の成分である、投影データのピクセルセットｓ（ｉ）２０１に関する関数ｆ１１（ｓ（ｉ））（ボックス２０５）の小信号の反転傾斜と、補完的な対応する重み付きの第２の成分である、関数ｆ１２（ｐｆ（ｉ））（ボックス２１０）の小信号の反転傾斜との、ピクセルごとの加重平均を表してもよく、ｐｆ（ｉ）は、予測フィードバック測定ピクセル２７６の値のセットであってもよい。ボックス２４０で、行列Ｌ（ｉ）の要素は、関数ｇ２２（ｄｔｒ（ｉ））のボクセルごとの小信号の反転傾斜、又はその反転傾斜と、関数ｆ２１（２６６）の小信号の傾斜との加重平均を表してもよい。
【０１０７】
[0107]関数ｆ１１（ボックス２０５）及び関数ｆ１２（ボックス２１０）は、同一のセットの対を形成してもよく、関数ｆ２１（ボックス２９０）及び関数ｇ２２（ボックス２５０）は、相互に逆関数であるセットの対を形成してもよい。目的に応じて、関数ｆ１１（ボックス２０５）、ｆ１２（ボックス２１０）、及びｆ２１（ボックス２９０）の機能は、例えば、対数関数（例えば、投影がごく少数である場合）又は平方根関数（例えば、測定された陽電子放射断層撮影（ＰＥＴ）又は単光子断層撮影（ＳＰＥＣＴ）の光子計数信号の分散安定化の場合）であってもよい。関数ｇ２２（ボックス２５０）は、例えば、指数関数又は二次関数の一部であってもよい。関数ｆ１１（ボックス２０５）及び関数ｆ１２（ボックス２１０）は、例えば、１つ又は複数の入力投影ピクセルのポアソン分散安定化平方根関数であってもよい。
【０１０８】
[0108]ＣＴ再構成のための近似的なカルマンフィルタゲイン行列Ｇを得る、米国特許第８，６６０，３３０号の概念の代替及び拡張が使用されてもよい。そのため、図２では、制約を可能にするように設計された、一般化された非線形の近似的カルマンフィルタが要約されている。強化の態様は、ボックス２３５での、現在の非線形フレームワーク内で近似的なカルマンフィルタゲインをエミュレートするための、近似的な逆投影又はその他の反転アルゴリズムの埋め込みである。これらの機能的及び構造的変更は、１ステップの最適な推定器を損なうが、ボクセル密度推定フィルタ内で線形化を使用して、急速に収束する断層撮影再構成を可能にする。一般化では、図２に示されているｆ１１（ボックス２０５）及びｆ１２（ボックス２１０）、並びにｆ２１（ボックス２７０）及びｇ２２（ボックス２５０）の非線形関数の選択に関連する、行列Ｚ（ボックス２２０）及びＬ（ボックス２４０）を追加する。
【０１０９】
[0109]表１から、図２に示されている静的線形モデルに関して、次の式が得られる。
Ｔ（ｉ＋１，ｉ）＝Ｉ    ［２．１］
Γ（ｉ）＝Ｉ    ｉ＝０の場合    ［２．２］
＝０    ｉ＞０の場合
【０１１０】
[0110]前の統計値Ｅ［ｘ（０）］＝μは、インプラントのような、先験的に知られた密度値μｘｉを含む領域などの、ベクトルｙで表されたｘの固定されたサブセクションを含んでもよい。これに対応して、Ｖｙｙ（ｉ）＝０である。それに応じて、反復内でｙがμｅにリセットされ、μｅはｇ２２−１を使用して変換されたμｘｉのサブセットを表す。
【０１１１】
[0111]効果的な線形画像再構成（ＬＲＡ）は、ルートが通常は＜＜１（及び、ｄｅｔ｜｜ＣＸ｜｜＜＜１）である、縮小するマッピングＣＸを生成する。
ＶＸ（ｉ）＝ＣＸＶＸ（ｉ｜ｉ−１）＝ＣＸＶＸ（ｉ−１）    ［２．３］
帰納的分散アルゴリズムから、次式が得られる。
Ｋ（ｉ）Ｈ（ｉ）＝Ｉ−ＣＸ（ｉ）又はＣＸ（ｉ）＝Ｉ−Ｋ（ｉ）Ｈ（ｉ）    ［２．４］
【０１１２】
[0112]小信号表現の場合、新しいカルマンゲインＫ’が指定され、小信号ループゲイン行列Ｚ（ボックス２２０）及びＬ（ボックス２４０）に関連付けられてもよい。小信号ループゲイン行列Ｚ及びＬは、非線形変換の傾斜に基づいてもよい。例えば、ゲイン行列Ｚ（ボックス２２０）は、ボックス２０５（ボックス２０５では、投影ピクセルデータ２０１が入力され、例えば対数的に変換される）の入力変換の動作点及びボックス２１０の予測データ２１１に関連付けられた関数ｆ１１及びｆ１２の傾斜を補償するために使用されてもよい。同様に、例えば、ゲイン行列Ｌ（ボックス２４０）は、出力変換ｇ２２（例えば、指数関数）及びその反転ｆ２１に関連付けられてもよい。線形の場合、カルマンＫ行列（ボックス２３５）は、例えば、フィルタ補正逆投影（ＦＢＰ）又は逆投影フィルタ（ＢＰＦ：ｂａｃｋ−ｐｒｏｊｅｃｔｉｏｎ  ｆｉｌｔｅｒ）近似（Ｚｅｎｇ他［１］を参照）などの、ＬＲＡによって近似されてもよく、Ｈは、散乱などの現象を含む、推定された物体の順投影であってもよい。
【０１１３】
[0113]任意選択的なゲイン行列Ｇ（ボックス２２５）は、Ｋが通常、近似されるため、Ｋ（ボックス２３５）に対する補正を支援してもよく、Ｇ（ボックス２２５）は、エントロピー増加、安定化、及び制約するローパスフィルタ、先験的な制約、及びフィードバックループ内のその他の介入などの影響を補償するのを支援する。ローパスフィルタ２５５及び先験的な物体に関する知識２６０が、特にまばらな投影データの物体密度計算において、ｅＨＥＣＴＯＲの安定性を支援するように、設定されてもよい。
【０１１４】
[0114]すべての行列について、退化した単一のグリッド点のケースに対して、小信号解析が、特定の入力データ２０１から開始し、例えば、スカラーゲイン（行列）Ｇを調整して、ユニティ（ＤＣ）ループゲインを提供してもよい。さらに多くのグリッド点が使用される場合、初期ループゲインは、例えば、米国特許第８，６６０，３３０号又はその他の基準に類似する観測されたイノベーションシーケンスデータに基づいて、行列Ｇを使用して調整されてもよい。フィードバックループの解析の一部の実施形態では（例えば、特定のグリッド分解能での物体再構成問題が、少数の投影を使用する場合などの、まばらな測定データセットを使用するのではない場合）、ボックス２５５での平滑化処理（エントロピー増加処理）の使用が無視されるか、又は影響を最小限に抑えるように調整されてもよい。
【０１１５】
[0115]ノイズが少なく、少量のイノベーションデータ２０６の場合、図２の非線形フィルタリング問題は、例えば表１に示されているカルマンフィルタ方程式を使用して解析されてもよい。
【０１１６】
[0116]簡単にするために、２０１のピクセルデータｓが、中程度のノイズ及び実際の物体密度に近い計算された物体密度の収束に対応する、２７６の予測ピクセル値ｐに近い値を含んでいると仮定する。その場合、例えば、データ２０１の関数ｆ１１のピクセルごとの傾斜値及びデータ２７６の関数ｆ１２のピクセルごとの傾斜値は、おおよそ同じ（ｓ≒ｐ）である。２７６のフィードバックの摂動Δｐ（又は、同様に２０１の測定データの摂動Δｓ）の小信号ループゲイン特性を補償するために、小信号対角ゲイン行列Ｚ（ボックス２２０）が、次の反転傾斜要素を含んでもよい（ｓ≒ｐを仮定する）。
ｚｄｄ＝σｐｄ／σｆ１１（ｐｄ）＝ｚｄ    ［２．５］
ここで、ｚｄは、インデックスｄでのＺ行列の対角非ゼロ要素を表し、σは微分演算子である。傾斜ｚｄは、データ２０１及び２８６のノイズ特性に応じて関数ｆ１２及び関数ｆ１１から得られた重み付きの組み合わせであってもよい。係数ｚｄの使用は、下でさらに説明される影響関数による制約を受けることがある。
【０１１７】
[0117]その後、近似された小信号の摂動がＬＲＡ（ボックス２３５）を通過し、カルマンフィルタゲインを近似し、予備的ボクセル更新を計算してもよい。
【０１１８】
[0118]フィルタの小信号出力（表１のｘ（ｉ）に対応するデータ）の後に、小信号変換行列Ｌ（ボックス２４０）が続き、変換物体密度推定値ｄｔｒ（ｉ）のその後の非線形ゲインを補償してもよい。
【０１１９】
[0119]対角行列Ｌ（ボックス２４０、ボクセル密度Ｖｑに関連する）は、次の要素を含む。
ｌｑｑ＝σｖｑ／σｇ２２（ｖｑ）ｌｑ    ［２．６］
インデックスｑｑをｑに置き換えると、対角要素を示す。係数ｌｑの使用は、下でさらに説明される影響関数による制約を受けることがある。
【０１２０】
[0120]成分２２０、２３５、及び２４０を組み合わせるが、
（ｉ）簡単にするために、Ｇ＝Ｉのままにし、
（ｉｉ）ｄｉｍ（ｚ）≧ｄｉｍ（ｘ）であると仮定し、ｚは冗長であるが、物体密度ｘを計算するために独立して収集された測定値であるとすると、
新しい小信号の、入力データ２０１に依存する、カルマンフィルタに似たゲイン行列Ｋ’が得られる（ｄｉｍ（ｚ）＜ｄｉｍ（ｘ）の場合に使用されるループ内の平滑化処理／エントロピー増加２５５及び２６０の処理を無視する）。
Ｋ’＝ＬＫＺ    ［２．７］
【０１２１】
[0121]関数ｇ２２（ボックス２５０）の成分の項の非線形の寄与の小信号近似（傾斜）から、再び、行列Ｌ−１が得られる。この簡略化されたモデルにおいてボックス２５５及び２６０を無視して、投影行列Ｈ（ボックス２７５）が、ｐｄ（ｉ）≒ｓｄ（ｉ）の範囲内の動作点で、ｆ１１［ｓｄ（ｉ）］（ボックス２０５）とｆ１２［ｐｄ（ｉ）］（ボックス２１０）の重み付きの組み合わせによって形成された線形近似によって変更され、例えば、１／ｐｄ≒１／ｓｄとなる傾斜要素の重み付きの組み合わせを使用するｆ１１及びｆ１２での対数変換の場合に、対角行列Ｚ−１を生成する。したがって、Ｌ−１でのｇ２２の小信号の置換及びｆ１２とｆ１１の組み合わせがＺ−１をもたらし、新しい小信号データ摂動投影行列Ｈ’が得られ、次式で定義される。
Ｈ’＝Ｚ−１ＨＬ−１    ［２．８］
【０１２２】
[0122]新しいカルマンゲインＫ’と観測行列Ｈ’の共有の特性を確認すると、対角（スケーリング）行列Ｌ及びＺを使用して、帰納的アルゴリズムから、次のループゲインが得られる。
Ｋ’Ｈ’＝ＬＫＺＺ−１ＨＬ−１＝ＬＫＨＬ−１＝Ｌ（Ｉ−ＣＸ）Ｌ−１＝（Ｉ−ＬＣＸＬ−１）    ［２．９］
ここで、ｄｅｔ｜｜Ｌ｜｜ｄｅｔ｜｜Ｌ−１｜｜＝１、ＣＸのルート＜＜１であり、適切なＬＲＡの使用を仮定する。小信号近似の場合に、フィードバックループを再帰的に通過することによって、最適なカルマンフィルタゲインに対するＬＲＡの近似の逸脱ＣＸによって制限される幾何学的物体密度縮小係数が得られる。例えば、ＩとＫＨの間の一致に応じて後の分散の削減を示す、表１の線形カルマンフィルタの帰納的分散アルゴリズムを比較する。反対に、一致しない測定成分は、それ以上縮小傾向があるマッピングがない、加重最小二乗推定につながる。収束での残差の構造は、診断システムの性能に役立つ。例えば、ホワイトノイズの投影残差及び投影残差において無視できる物体の特徴（もしあれば）は、収束を示し、再構成において取得されるモデルの精度を定量化することがある。
【０１２３】
[0123]Ｋ（ｉ）（方程式２．９）の有効性を表すＣＸ（ｉ）に対する判断は、次に基づいてもよい。
（ｉ）一連の反復における物体の特徴の残差の縮小、
（ii）測定によって引き起こされる投影の背景ノイズの存在（例えば、Ｘ線断層撮影又は電子線断層撮影の場合の近似的な投影のホワイトノイズ）、物体の投影の特徴の不在（例えば、同じ物体の異なる測定のセットの残差間の差異を調べることによる）、及び
（ｉｉｉ）グリッド分解能の増加でイノベーションを比較するときの、連続的なホワイトイノベーションデータ２１６に近いデータの存在（２１６のより高い分解能でのピクセルに対応するように補間する必要があるより低いグリッド分解能からの遡及的なイノベーションステップ）。
【０１２４】
[0124]ＬＲＡのＫ（ボックス２３５）の不正確さ（例えば、ｄｅｔ｜｜ＣＸ｜｜≠０）及び出力変換（ボックス２５０）後のローパス平滑化／エントロピー増加処理（ボックス２５５）をそれぞれ補償するために、Ｋは、イノベーションゲイン調整行列Ｇ（ボックス２２５）を使用して増大されてもよい。ゲイン調整器Ｇ（ボックス２２５）は、入力データのｆ１１の傾斜（ボックス２０５）及び観測の予測データのｆ１２の傾斜（ボックス２１０）の動作点の範囲内で、ゲイン値を適度に変更する。ゲイン調整行列Ｇは、予測投影密度の一貫した過小推定又は過大推定から生じる、ビームハードニングなどの要因による影響を受けることがある。それらのビームハードニングの影響の補正に固有の方法が、他者によって開発されている。行列Ｌを使用する価値は、単位行列Ｉ（方程式２．９）の２１６〜２１１の小信号ループゲイン全体を支援する関数ｇ２２に関連付けられた小信号ループゲインを補償することであり、それらの差異はボックス２１５において生じ、ノイズ以外の小さい差異を意味し、したがって、反復計算中の物体密度の急速な収束を意味する。
【０１２５】
[0125]主に関数ｆ１１（ボックス２０５）、関数ｆ１２（ボックス２１０）、関数ｆ２１（ボックス２９０）、及び関数ｇ２２（ボックス２５０）の小信号の特性が、ボックス２２０の行列Ｚ及びボックス２４０の行列Ｌに使用されてもよく、小信号行列要素の一部が、０に設定されるか、又は何らかの先験的な情報の存在下で除外されてもよい。例えば、インプラント又は患者の周囲の密度などの物体密度は、ある時点で先験的に知られてもよい。それらの物体密度の推定に関連付けられた要素は、除外されてもよく、対応する物体密度データが、既知のデータ値に置き換えられる。
【０１２６】
[0126]断層撮影再構成の場合、既知ではなく、例えば、インプラントなどの物体の特性、ビームのスペクトル特性、放出されるビームの変化する硬度、物体の動き、散乱、呼吸及び心臓の動き、並びに画像再構成器１６５の外部のその他のシステムの構成要素に起因する、ビームハードニングにおける不確実性を表す他の手段によって観測もされないパラメータは、画像再構成器１６５内で、例えば図２に示されている投影行列Ｈ（ボックス２７５）で表されて、調整されてもよい。そのようなパラメータは、ボックス２００、２０５、２１０、２５０、２５５、２６０、及び２９０の関数を調整してもよい。それらのパラメータは、例えばＬｅｖｅｎｂｅｒｇ-Ｍａｒｑｕａｒｄｔ（ＬＭ）のアプローチを使用して、一緒に推定されるか、又は別々に推定されてもよい。測定値のセットを記録した後に、ＬＭ又は同等のアプローチは残差データから計算し、それらの残差データは、対応する物体の推定が収束したときに、画像再構成器１６５内で変更されたパラメータのセットから計算されており、不確実性を含むパラメータに対する調整が、収束するまで継続する。
【０１２７】
[0127]冠動脈枝の拍動モデルの推定にＬＭ処理を使用する、そのようなＬＭベースの再構成の例が、図３及び図４に示されており、データは、医療用２方向Ｃアームシステムを使用して収集されている。モデルパラメータ計算では、５０個以上のパラメータを使用して、例えばイメージインテンシファイアの静磁場及び変動磁場のゆがみ並びにＣアームの位置の誤差から生じる投影の位置合わせを再推定した。
【０１２８】
[0128]画像再構成器１６５における反復の初期化は、多数のボクセルではなく、例えば投影ごとに単一のボクセル及び単一のピクセルから、低いグリッド分解能で開始してもよい。初期ボクセル密度は、例えば、投影密度の平均から得られる。次に、反復で、グリッド分解能が望ましいレベルに増やされてもよい。分解能の変化は、同じ反復間隔で行われなくてもよく、１回の分解能の増加量が同じでなくてもよい。グリッド分解能が増やされるときにモデルパラメータｐの成分を組み込む数及び順序は、誤差の削減、又は例えばＳｃｈｗａｒｚ［６］及びＡｋａｉｋｅ［７］などにおける情報基準に基づいて、事前に決定されてもよい。パラメータ係数の受容は、例えば、データ２０６又は１６０の測定値、及び再構成でｐの成分を組み込むときの改善の程度に基づいてもよく、パラメータ及びそれらの選択の重要度の決定に使用されてもよい。
【０１２９】
[0129]初期化後に、ボックス２２０及びボックス２４０によって、ループゲインの主要なゲイン調整が決定される。残りのゲイン調整があれば、それらのゲイン調整が、一連の先行するイノベーション２１６のパターンから得ることができる行列Ｇ（ｉ）（ボックス２２５）のゲインを使用して実施されてもよい。固定されたグリッド分解能での有意な振動は、過剰なゲイン、持続する有意な静的イノベーション値の不十分なゲイン、及びランダムなパターンの適切なゲインを示していることがある（例えば、米国特許第８，６６０，３３０号を参照）。
【０１３０】
[0130]発明の手法の比較
【０１３１】
[0131]非線形の反復断層撮影再構成は、多くの人によって、長年にわたって集中的に追及されてきた。例えば、Ｅｌｂａｋｒｉ他［８］及び［９］には、多エネルギーＸ線源スペクトル、エネルギー依存の減衰、及び重複しない材料を考慮するアプローチが記載されている。それらの目的は、既知のエネルギー依存の質量減衰係数を使用して未知の物体密度を計算することである。Ｅｌｂａｋｒｉのアプローチは、この多エネルギーモデルに対してペナルティ付き尤度関数を定式化し、順序付けられたサブセットの反復アルゴリズムを使用して、ボクセルごとに未知の物体密度を推定する。
【０１３２】
[0132]同様の目的に対処する多エネルギーの同時代数的再構成法（ＳＡＲＴ：Ｓｉｍｕｌｔａｎｅｏｕｓ  Ａｌｇｅｂｒａｉｃ  Ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  Ｔｅｃｈｎｉｑｕｅ）（アルゴリズム）のＨｕｍｐｈｒｉｅｓの解析［１０］は、収束条件が減衰係数の特性を含む複数の要因に依存することを示している。さらに、ＳＡＲＴの収束は保証されず、反復のヤコビ行列のスペクトル半径は、必ずしも１未満にならない。
【０１３３】
[0133]米国特許第８，９２３，５８３号は、安定化ペナルティ関数Ｓ（ｘ）を使用してモデルの尤度Ｄ（ｙ；ｘ）を最小化することによって、マルチスペクトル測定からの二重エネルギーＣＴデータを再構成した。米国特許第８，９２３，５８３号には、測定値ｙから物体ｘを再構成するときに、「この方法の欠点は、ｘとｙの間の非線形関係に起因して、関数Ｄ（ｙ；ｘ）をモデル化するのが、通常は非常に複雑であり、計算的に困難であることである」と記載されている。
【０１３４】
[0134]最近、Ｓｔａｕｂ他［１１］は、少数の固有ベクトル表現を物体の再構成に使用して、患者の生体構造の動きをボクセルごとのスケールで使用して追跡する能力を示した。Ｓｔａｕｂ他［１１］は、Ｎｅｌｄｅｒ−Ｍｅａｄシンプレックスアルゴリズムが、物体の動きを表す少数のパラメータを推定するための最もロバストなアプローチであることを発見した。
【０１３５】
[0135]前述した困難を克服する革新的な代替の計算の概要説明を簡単にするために、新しいアプローチのすべてのモデルパラメータを推定の対象にすることができる。ビームハードニング、Ｘ線源スペクトルの変化、機器の動き又はずれなどを特徴付ける外部機器の特性のパラメータ、或いは物体の材料特性、物体の動きなどの物体内部のパラメータが、単一のパラメータベクトルｐにおいて組み合わせられる。システムの変化は、このパラメータの変化によって特徴付けられ、図１に示されているように、１９５を介してインベスティゲータ１７０によって送信されたパラメータデータの計算された推定値によって、修正することができる。再構成パラメータは、例えば、残差二乗和又は残差データの計算されたロバストな推定値のその他の基準を使用してもよく、例えばデータ１６０、データ２１６の影響関数を使用して、例えばシミュレーション、不一致の残差（ＭＭＲ−２６１）の場合に、再構成の性能を決定する。この設定において、パラメータ、機器、及び物体の任意のセットに対して、最適な物体再構成が実行される。パラメータベクトルｐから、再構成誤差を最小化するか、又は再構成の性能を最適化する、大域的に最適なパラメータデータｐのセットが計算される。
【０１３６】
[0136]現在、未知の画像化システムの成分ｐを推定しない場合でも、公開された非線形の反復アルゴリズムの物体密度推定は、前述したように、ゆっくりと収束する。物体の外部又は内部のパラメータを含む、最適なデータセットｐの複数の成分を推定する必要がある適用の場合、それらのシステムは、次によって制限される。
１．単一のモデルパラメータの評価を追加する場合でも、著しく計算量が増加する。
２．収束が達成されたように見える場合でも、強い物体の特徴を示す、反復ごとの投影残差の高い相関関係。この欠点は、例えばイノベーション２１６に関する、誤差密度の縮小傾向が大きいマッピングの欠如に起因しており、最適解から遠いモデルパラメータの効果的な調整を支援しない。
３．明示的な反復的数値最小化（山登り法）を使用して、単一の特定の固定されたパラメータｐに関して性能基準を決定し（一部分において、最小化する解からの遠い距離に起因する）、無限又は単数に近いパラメータ更新行列を作成し、多数の可変モデルパラメータを考慮する能力をさらに制限する、投影残差の高い相関関係。
４．安定化ペナルティ関数ｄ（ｙ：ｘ）を識別するのが困難である（例えば、米国特許第８，９２３，５８３号を参照）。
【０１３７】
[0137]対照的に、本明細書に記載された本発明は、縮小傾向が大きいパラメータ推定手法（方程式２．９）を組み合わせる。本発明の有効性は、例えば次から得ることができる。
１．反復的データ処理の第１のレベルで、線形化されたカルマン投影法によって提供される、数値的ではない、解析的方法を使用し、そのレベル内で、少ない数値計算量（方程式２．９）で、反復から反復への縮小傾向が大きいマッピングを作成する、効率。
２．効果的な線形に近い小信号近似計算を可能にするほど、密度推定値の変化を十分小さく保つ、グリッド分解能の緩やかな増加。
３．ベクトルｐ内のパラメータの数値最適化。各反復内で、上の（１）によって最適解に近づきながら、（イノベーションプロセスを作成する反復の）投影残差の行列が、それらのパラメータの変化を明確に表す（最適化に関連付けられたヘッシアン行列の非特異性）。パラメータは、例えば、図２に概略的に示されている、投影行列Ｈｍｎｏ、変換ｆ１１、ｆ１２、ｆ２１、及びｇ２２、並びにその他の関連するシステム変数の係数に対する影響を表してもよい。
４．ボックス２２５などの安定化ペナルティ関数を設計及び適用することの単純さ。
【０１３８】
[0138]このアプローチ例示的な結果が、図６に示されている。革新的なｅＨＥＣＴＯＲは、ＣＥＲＮのＳＡＲＴの２０１７年８月３０のリリース（Ｂｉｇｕｒｉ他［２１］）と比較した場合、高い画像品質及び１０倍以上速い速度を示す。図６Ａに示されているように、ｅＨＥＣＴＯＲは、ＳＡＲＴより大幅に正確な物体の特徴の再構成及びより小さい残差を示している。さらに、図６Ｂで、ｅＨＥＣＴＯＲの暫定版が、まばらな投影データセット（例えば、５回の投影）から、有用な物体再構成結果を示しており、ＣＥＲＮのＳＡＲＴによる有用な物体再構成結果は、３２回の投影によって取得されている。図６Ａ及び６Ｂは、米国国立衛生研究所（ＮＩＨ：Ｎａｔｉｏｎａｌ  Ｉｎｓｔｉｔｕｔｅｓ  ｏｆ  Ｈｅａｌｔｈ）の国立小児保健発育研究所（ＮＩＣＨＤ：Ｎａｔｉｏｎａｌ  Ｉｎｓｔｉｔｕｔｅ  ｏｆ  Ｃｈｉｌｄ  Ｈｅａｌｔｈ  ａｎｄ  Ｈｕｍａｎ  Ｄｅｖｅｌｏｐｍｅｎｔ）のＥｕｎｉｃｅ  Ｋｅｎｎｅｄｙ  ＳｈｒｉｖｅｒのＢａｓｓｅｒ  Ｌａｂｏｒａｔｏｒｙで実装されたｅＨＥＣＴＯＲ及びＳＡＲＴを示している。これは、ＮＩＨによるｅＨＥＣＴＯＲプロセスの承認ではない。
【０１３９】
[0139]例えば、例えば、図１のインベスティゲータ１７０によって制御されるパラメータベクトルｐのサブセットに対するモデルパラメータの適合最適化は、図２に示されているような、２１６、２２１、２２６、又は１６０の残差データ行列又はＭＭＲ−２６１を使用して、測定と予測の間の誤差を表すそれらの残差を最小化する。２エネルギースキャニングなどの、複数の測定方法を使用する場合、推定問題のコスト関数に対して、重み付けを指定する必要がある。この数値最適化は、反復再構成アルゴリズムの実行時に、特定の解析に基づいて収束される反復で使用されてもよい。その後、ｐを数値的に推定するための好ましい方法は、これらの残差データ行列のいずれかに適用されるＬｅｖｅｎｂｅｒｇ-Ｍａｒｑｕａｒｄｔ手順（ＬＭＰ：Ｌｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔ  ｐｒｏｃｅｄｕｒｅ）である。５０個以上のパラメータを使用する場合でも、収束を簡単に実現することができる。図３の冠動脈の拍動モデルの投影、及び図４の対応する再構成に示されているように、５０個以上のパラメータを使用して、例えば、２つの陰極Ｘ線イメージインテンシファイア及びそれらの個別の地磁場に対する方位の感度を使用することに起因する画像ゆがみ、二重Ｃアームガントリーのアームの独立した動揺、回転位置合わせの中心のわずかな相互のずれ、６つの投影方向に関する冠動脈の拍動モデルの画像フレーム収集のタイミングのオフセット、冠動脈モデルの屈曲、ねじれ、及び脈動、並びにＣアームの方位角の増分の誤差を表す。
【０１４０】
[0140]ｅＨＥＣＴＯＲの効率は、反復ループ内で小信号近似を使用して実現される。特定の分解能での小信号動作の良好な収束は、計算コストがさらに大きい、より高い分解能でのより良い小信号近似を支援する。したがって、小信号近似が「十分に良い」ように見えるとき（例えば、線形モデルからの非線形性のループ効果のずれが１０％未満であるとき）に、より高い分解能への切り替えが発生するべきである。この線形性の目的が、少なくとも、高密度又は重要な領域を表すデータ（例えば、ピーク密度より１％を超えて高い密度値）に適用されるべきである。効果的な収束を実現するための次のような複数のアプローチが、これらの組み合わせを含めて、可能である。
１．分解能の各レベルに必要な期待される反復回数を、理論的及び実験的評価から推定する。期待される収束は、特定の再構成の状況（例えば、頭蓋の再構成）に適用されてもよい。
２．物体データの残差値が反復間で有意に縮小する間、特定の分解能で反復を続行する。特に、残差データの変化がシステムノイズレベルと比較して大きい場合、反復を続行するが、残差データの変化が有意でない場合（例えば、残差データの変動がノイズ又はその他の基準に相当する場合）、反復を終了する。
３．影響関数及びその他の非線形のシステム関数が線形に近い領域（ＮＬＲ：ｎｅａｒ−ｌｉｎｅａｒ  ｒｅｇｉｏｎｓ）内で動作していない間、例えば、影響関数の入力値及び出力値が１０％を超えて異なり、他の非線形関数の入力／出力の関係が、先行する反復と比較したときに１０％を超えて変化する場合、特定の分解能で反復を続行する。
データが、診断上重要な領域の大部分（例えば、領域を含んでいる物体情報の９９％）を十分適切に通って、上で示したすべての関数のＮＬＲを通過した後に、次の分解能に切り替える前に、線形プロセスの誤差の縮小の有効性を使用して、数回（例えば、１回又は２回）の追加の反復を続行する。
【０１４１】
[0141]例えば、図６の計算は、上記の１番及び３番のアプローチの組み合わせを使用している。
【０１４２】
[0142]加えて、最高の分解能レベルでの反復を含む、反復のいずれかの終了は、事前に決定されるか、制限されるか、又は例えば統計的尺度及びその他のコスト基準がしきい値に達することから導かれてよい。
【０１４３】
[0143]残差の空間的又は時間的変化は、再構成中のグリッド分解能の緩やかな増加に関連して、小さくなる傾向があり、効果的な線形化を支援する。さらに、パラメータベクトルタｐの成分は、通常、最大分解能未満で適切に定義される。このようにして、例えば、最も細かいグリッド分解能の半分で、１００個の成分を含むパラメータベクトルｐの１回の更新の計算、及び１００個のパラメータの変動の対応する計算は、（３次元の再構成を仮定して）約１２回のフィルタ補正逆投影に相当する。少数のパラメータベクトルｐの更新は、適切に位置合わせされた再構成を実現するのに十分だった。
【０１４４】
[0144]同様に、吸収に関連する係数の変動を特徴付けるためのパラメータベクトルｐによるシステムのパラメータ化を使用して、物体密度を効率的に推定するｌことができる。例えば、パラメータベクトルｐと共に適切な多項式を使用すると、Ｅｌｂａｋｒｉ他［８］又はＨｕｍｐｈｒｉｅｓの推定目的にとって望ましいように、可変物体密度のビームハードニングの吸収係数が、局所的な物体密度の関数及び投影に対する積分経路長の影響をモデル化することができる。しかし、吸収に関連する係数の効果を計算する手法が、例えば焦点サイズ、投影方向のずれ、及び画像収集中の患者のテーブルの前進などの、画像化システムを特徴付けるパラメータを推定するために使用されてもよい。実際、この方法は、固定された以前の再構成パラメータを使用するより好ましいことがあり、代わりに、システムの自己較正、「生命」システムの性能の検証、及び保守サービスの必要性の削減を可能にする。
【０１４５】
[0145]画像再構成の別の重要な態様は、信号対ノイズ比（ＳＮＲ：ｓｉｇｎａｌ−ｔｏ−ｎｏｉｓｅ  ｒａｔｉｏ）を含めることである。モデルベース反復再構成（ＭＢＩＲ）は、例えば、データを再重み付けする反復手法を使用して、ＳＮＲを考慮するペナルティ付き加重最小二乗法（ＰＷＬＳ：ｐｅｎａｌｉｚｅｄ  ｗｅｉｇｈｔｅｄ  ｌｅａｓｔ  ｓｑｕａｒｅｓ）による密度推定を実現し、それらの重みは、分解された副鼻腔撮影図の逆共分散を近似するように計算される（例えば、米国特許第８，９２３，５８３号を参照）。
【０１４６】
[0146]複数の手法を個別に、又は組み合わせて使用することで、再重み付け及び逆共分散行列の使用を回避して、密度推定の大幅な簡略化が可能になり、例えば米国特許第８，９２３，５８３号に記載されている制限を克服することができる。そのような手法の例としては、以下が挙げられる。
１．分散安定化又は再設計データ変換の使用、
２．表１のゲイン方程式を近似するための、イノベーションの寄与の可変平滑化及び増幅、及び
３．データの信頼性を考慮した、ロバストな推定のための影響関数の使用。
【０１４７】
[0147]分散安定化データ変換
【０１４８】
[0148]独立した測定及びそれらの測定の期待される値に関連する分散の場合、ボックス２０５の分散安定化変換ＳＶが、推定及びシステム性能の評価を簡略化する。一部の測定の分布計算では、既知の特性に基づいて変換を使用することができ、その他の測定の分布計算では、変換が、理論的にではなく、数値的に評価されてもよい。
【０１４９】
[0149]例えばＳＰＥＣＴにおいて、投影２０１を生成するポアソン分布の場合、ボックス２０５の安定化変換ＳＶは、平方根変換の変形である。光子に基づく、ボックス２０１のピクセル強度測定データの近似変換の例は、ボックス２０６で生成されるデータｓ’を使用して、次のようになる。
ｓ’＝ｓｑｒｔ（ｍ）＋ｃ（ｍ，ｑ）    ［４．１］
ここで補正ｃは、基本的な近似からの小さいずれを表す。補正ｃ（ｍ，ｑ）（ｑは定数）は、二次の項に対応するために使用されてもよい。さらに一般的には、分散安定化変換では、定数パラメータｑが、次のように見えてもよい。
ｓ’＝ｆ（ｍ，ｑ）    ［４．２］
【０１５０】
[0150]相関関係のない測定ノイズが存在する状態で、このピクセルごとのスカラー変換のセットに対する解は、短時間で計算することができる。
【０１５１】
[0151]ノイズＶＶに応答したゲインの調整及び平滑化
【０１５２】
[0152]増加した測定ノイズの存在下で、線形再構成アルゴリズム（ＬＲＡ）（例えば、ボックス２３５の逆投影計算）によって近似して、最適なカルマンゲインの特性を複製するために、表１のゲイン方程式は、データ２４６の物体密度を更新するときに、例えばボックス２２０の行列Ｚ又は２２５のＧを使用して、イノベーションデータのゲイン／寄与を削減することを示唆している。さらに、安定性を支援するために、関数ｆ１１（ボックス２０５）及びｆ１２（ボックス２１０）の傾斜が変更されてもよい。
【０１５３】
[0153]同様に、（例えば、非線形入力変換関数ｆ１１（ボックス２０５）に関連して）一部の投影領域内の測定の重要度が他の投影領域内よりも高い／低いと考えられる場合、残差のスケーリングによって、対応するピクセルイノベーション（又は画像データの残差２１６）の重み付けを増やす／減らすことができる。
【０１５４】
[0154]物体密度の変動性に対するノイズＶＶの影響を減らすための代替のアプローチは、例えば推定の次元数の削減に相当する局所的な平滑化によって、投影測定値を局所的に結合することである。
【０１５５】
[0155]影響関数の適用
【０１５６】
[0156]影響関数は、モデルの仮定からの測定の外れ値又は逸脱の影響を減らすことができ、欠陥のないデータと比較した場合に、少ない計算量及びわずかな品質低下で、物体再構成の改善を可能にする。影響関数は、図２のボックス間の接続のいずれかにおいて利用されてもよい。しかし、影響関数は、ボックス２００と２０５の間、又はボックス２１５と２２０の間及びボックス２３５と２４０の間、又はボックス２４０と２４５の間に配置される可能性が最も高い。さらに、影響関数は、可変ＳＮＲを含むデータを考慮するように設計され得る。通常、影響関数は、信号の最大振幅を制限し、重要な入力値を超えて再降下することがある。例えば、Ｂｏｕｍａｎ他［１３］は、エッジを維持する影響関数の特性を研究しており、最大事後確率（ＭＡＰ：ｍａｘｉｍｕｍ  ａ  ｐｏｓｔｅｒｉｏｒｉ）推定及び対数尤度推定においてそれらの影響関数の有用性を示している。しかし、それらの影響関数の目的は、ランダムなフィールド内の特定の特徴を保つことである。これは、物体密度の推定において測定の一貫した構造情報が累算され、ランダム性／エントロピーが残差において表される本発明とは、対照的である。
【０１５７】
[0157]本発明における不完全な測定モデルを前提として、ボックス２１５と２３５の間の第１の影響関数は、例えば、係数ｚを変更するか、又はボックス２２５のゲイン係数を変更することによって、データ２１６に適用されるか、又はボックス２２０に含まれ、ボックス２１５での測定信号と参照信号の比較の後に、ただしボックス２３５で物体密度更新データ信号を生成する本発明の前に、残差データに対して動作してもよい。例えば図２では、影響関数の影響を受ける最有力候補は、データ２１６（変換された測定値と変換予測の間の残差）である。この影響関数を指定するパラメータの選択は、ピクセルデータ２０６及び２１１内の対応する共同の近傍領域の値、データ２１６のピクセルの近傍領域の内の関数値の分布、及び前の反復を基準にしたその分布の変化に依存する。図２内のいずれかの影響関数の位置及び形状は、問題固有である傾向があってよく、適用分野に由来してもよい。
【０１５８】
[0158]図５は、データ２１６に対して、又はボックス２２０内で使用するための影響関数の形状の例を示している。影響関数は、データ２０１、２０６、又は２１１に関連付けられているなどの、データに対する処理の依存を示し、測定値の信頼性を示す。影響関数は、例えばノイズがある測定環境内の、めったにない大きい異常な測定値に対して、その後のデータ処理をロバストにすることができる。参考のために、図５の横座標の間隔［−１，１］は、ガウス分布を有するデータのサンプル値の約６７％を含んでいる。例えば、σ＝１であるイノベーションのガウス分布の場合、この範囲内のサンプル値のみが、イノベーションの最大に近い重み付けをもたらす。中心の傾斜は、１．００であってもよい。漸近傾斜は、０．００であってもよい。
【０１５９】
[0159]影響関数は、ボックス２２０のイノベーションゲイン調整の一部であることが可能である。例えば、Ｘ線断層撮影では、高密度の塊の再構成は、対応する低い光子計数及び関連する投影ピクセルデータ２０１の相対的に大きい変動に起因して、推定値の低いＳＮＲにつながることがある。本発明では、図５に示されているような影響関数を使用することの単純さが、例えば、分散安定化と相まって、米国特許第８，９２３，５８３号において提案されているように、物体ｘを測定値ｙから再構成するときに、物体ｘを測定値ｙから再構成するときに、モデルの尤度Ｄ（ｙ；ｘ）の数値的に要求の厳しい最小化を伴う安定化ペナルティ関数Ｓ（ｘ）を使用する代替的アプローチより、計算的に好ましいということが、理解されてもよい。
【０１６０】
[0160]図１及び図２に示された手法に関して説明された本発明の実施形態例は、例えば、非一時的データストレージデバイス（例えば、ＣＤ−ＲＯＭ、ＤＶＤ、ブルーレイ、光磁気（ＭＯ）ディスク、ハードディスク、フロッピーディスク、ジップディスク、フラッシュドライブなど）に格納されたソフトウェアコードとして提供されてもよい。格納されたソフトウェアコードは、例えば図１及び図２に関して前述された例示的な手法を実行するために、１つ又は複数のプロセッサ、１つ又は複数の非一時的メモリデバイス（例えば、ランダムアクセスメモリ（ＲＡＭ：ｒａｎｄｏｍ−ａｃｃｅｓｓ  ｍｅｍｏｒｙ）デバイス、ダイナミックＲＡＭ（ＤＲＡＭ：ｄｙｎａｍｉｃ  ＲＡＭ）デバイス、フラッシュメモリデバイス、及びスタティックＲＡＭ（ＳＲＡＭ：ｓｔａｔｉｃ  ＲＡＭ）デバイスなど）を含んでいるコンピュータによって読み込み可能であってもよく、実行可能であってもよい。
【０１６１】
[0161]本発明の実施形態例は、図１及び図２に関して前述された例示的な手法を格納して実行するために、１つ又は複数のプログラムストレージ及び実行デバイス（例えば、特定用途向け集積回路（ＡＳＩＣ：ａｐｐｌｉｃａｔｉｏｎ  ｓｐｅｃｉｆｉｃ  ｉｎｔｅｇｒａｔｅｄ  ｃｉｒｃｕｉｔｓ）、フィールドプログラマブルゲートアレイ（ＦＰＧＡ：ｆｉｅｌｄ  ｐｒｏｇｒａｍｍａｂｌｅ  ｇａｔｅ  ａｒｒａｙｓ）、複合プログラム可能論理デバイス（ＣＰＬＤ：ｃｏｍｐｌｅｘ  ｐｒｏｇｒａｍｍａｂｌｅ  ｌｏｇｉｃ  ｄｅｖｉｃｅｓ）、特定用途向け命令セットプロセッサ（ＡＳＩＰ：ａｐｐｌｉｃａｔｉｏｎ  ｓｐｅｃｉｆｉｃ  ｉｎｓｔｒｕｃｔｉｏｎ−ｓｅｔ  ｐｒｏｃｅｓｓｏｒｓ）など）を提供してもよい。
【０１６２】
[0162]本明細書に記載された例及び実施形態は、非限定的な例である。
【０１６３】
[0163]本発明が実施形態例に関して詳細に説明され、当業者にとっては、本発明から逸脱することなく、より広い態様において変更及び修正が行われてもよいということが、前述の説明から明らかになったであろう。したがって、特許請求の範囲において定義されている本発明は、本発明の真の思想に含まれている、そのようなすべての変更及び修正を対象にするよう意図されている。
【０１６４】
[0164]参考文献
【０１６５】
[0165]本明細書では、以下の参考文献が参照されている。
【０１６６】
[00166][1] Zeng, Gengsheng L., and GrantT. Gullberg. "Can the backprojection filtering algorithm be as accurate asthe filtered backprojection algorithm?" NuclearScience Symposium and Medical Imaging Conference, 1994., 1994 IEEE ConferenceRecord. Vol. 3. IEEE, 1994.
【０１６７】
[00167][2] Yingying Zhang-O’Connor, Jeffrey A. Fessler, Fourier-Based Forward andBack-Projectors in Iterative Fan-Beam Tomographic Image Reconstruction, p.582IEEE Trans. Med. Imaging, Vol. 25, No. 5, MAY 2006.
【０１６８】
[00168][3] Long, Yong, Jeffrey A. Fessler,and James M. Balter. "3D forward and back-projection for X-ray CT usingseparable footprints." IEEE transactions on medical imaging 29.11(2010): 1839-1850.
【０１６９】
[00169][4] Sage, A.P., J.L. Melsa,Estimation Theory with Applications to Communications and Control, The KalmanBucy Filter, New York: McGraw-Hill, 1971. ISBN-10: 0070544298, ISBN-13:9780070544291. 
【０１７０】
[00170][5] Radon, J., "Uber die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten," Berichte der Sachsischen Akademie der Wissenschaften, Mathematisch PhysischeKlasse, Leipzig, Germany, vol. 69, pp. 262-277, 1917).
【０１７１】
[00171][6] Schwarz, Gideon,"Estimating the Dimension of a Model." The Annals of Statistics, Vol.6, No. 2, pp. 461-464, 1978.
【０１７２】
[00172][7] Akaike, H., "StochasticTheory of Minimal Realization." IEEE Trans. On Automatic Control, Vol. AC-19,No. 6, pp. 667-674, December, 1974.
【０１７３】
[00173][8] Elbakri, Idris A., and JeffreyA. Fessler. "Statistical image reconstruction for polyenergetic X-raycomputed tomography." IEEE transactions on medical imaging 21.2(2002): 89-99.
【０１７４】
[00174][9] Elbakri, Idris A., and JeffreyA. Fessler. "Segmentation-free statistical image reconstruction forpolyenergetic x-ray computed tomography with experimental validation." Physics in medicine and biology 48.15(2003): 2453-2477.
【０１７５】
[00175][10] Humphries, T. "TechnicalNote: Convergence analysis of a polyenergetic SART algorithm." Medical physics 42.7 (2015): 4007-4014.
【０１７６】
[00176][11] Staub, David, et al. "4DCone-beam CT reconstruction using a motion model based on principal componentanalysis." Medical physics 38.12 (2011): 6697-6709.
【０１７７】
[00177][12] Humphries, T., and A.Faridani. "Reconstruction of CT Images from Sparse-View Polyenergetic DataUsing Total Variation Minimization." 2015IEEE Nuclear Science Symposium, San Diego, California. 2015.
【０１７８】
[00178][13] Bouman, Charles, and KenSauer. "A generalized Gaussian image model for edge-preserving MAPestimation." IEEE Transactions on Image Processing2.3(1993): 296-310.
【０１７９】
[00179][14] Wood, S.L., Morf, M., "Afast implementation of a minimum variance estimator for computerized tomographyimage reconstruction", IEEE, Trans. on Biomed. Eng., Vol. BME-28, No. 2,pp. 56-68, 1981.
【０１８０】
[00180][15] Sunnegardh, Johan, IterativeFiltered Backprojection Methods for Helical Cone-Beam CT, Dissertation No.1264, Linkopings universitet, SE-581 83 Linkoping, Sweden, Aug. 2009.
【０１８１】
[00181][16] L. Paavolainen, Erman Acar,Uygar Tuna, Sari Peltonen, Toshio Moriya, Pan Soonsawad, Varpu Marjomaki, R. Holland Cheng, Ulla Ruotsalainen, Compensation of MissingWedge Effects with Sequential Statistical Reconstruction in ElectronTomography, PLoS One. 2014; 9(10): e108978.
【０１８２】
[00182][17] Dempster, A.P., Laird, N.M.,Rubin, D.B., Maximum Likelihood from incomplete data via the EM algorithm, J.R.Stat. Soc. Series B, vol. 39, pp. 1-38, 1977.
【０１８３】
[00183][18] U.S. Patent No. 8,660,328. 
【０１８４】
[00184][19] U.S. Patent No. 8,660,330.
【０１８５】
[00185][20] U.S. Patent No. 8,923,583.
【０１８６】
[00186][21] Ander Biguri, Manjit Dosanjh,Steven Hancock, Manuchehr Soleimani, TIGRE: a MATLAB-GPU toolbox for CBCT imagereconstruction, Biomedical Physics & Engineering Express, Volume 2, Number5, 8 September 2016.
【図１】
2019535405.tif 000004
【図２Ａ】
2019535405.tif 000005
【図２Ｂ】
2019535405.tif 000006
【図３】
2019535405.tif 000007
【図４】
2019535405.tif 000008
【図５】
2019535405.tif 000009
【図６Ａ】
2019535405.tif 000010
【図６Ｂ】
2019535405.tif 000011
【国際調査報告】
2019535405.tif 000012