(57)【要約】【課題】ＦＢＰアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減すること。【解決手段】ノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する。フィルタ整形部１１４Ｂは、重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。画像発生部１１４Ｃは、再構成フィルタに基づいて画像を発生する。【選択図】図２
2014023936.tif 000002
【特許請求の範囲】
【請求項１】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備するＸ線コンピュータ断層撮影装置。
【請求項２】
前記決定部は、ビュー毎に前記投影データに基づいて前記重み値を決定する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項３】
前記再構成フィルタは、下記数式のＨ（ｗ，ｖｉｅｗ）により規定され、
【数１】
2014023936.tif 000050
ここで、ωは周波数分散であり、ｗ（ｖｉｅｗ）はビュー単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
請求項２記載のＸ線コンピュータ断層撮影装置。
【請求項４】
前記決定部は、投影線毎に前記投影データに基づいて前記重み値を決定する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項５】
前記再構成フィルタは、下記数式のＨ（ｗ，ｒａｙ）により規定され、
【数２】
2014023936.tif 000051
ここで、ωは周波数分散であり、ｗ（ｒａｙ）は投影線単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
請求項４記載のＸ線コンピュータ断層撮影装置。
【請求項６】
前記発生部は、逆投影が実行される前に、前記再構成フィルタを前記投影データに適用する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項７】
前記発生部は、逆投影が実行された後に、前記再構成フィルタを前記投影データに適用する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項８】
前記ノイズモデルは、ポワソン統計を含む、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項９】
前記ノイズモデルは、ガウス分布に従う電子ノイズを含む光子強度計測の複合ポワソン統計によって記述される、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項１０】
前記再構成フィルタを整形するために正規化する正規化部をさらに備える、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項１１】
前記画像が基準画像に基づいてさらに発生される、請求項１０記載のＸ線コンピュータ断層撮影装置。
【請求項１２】
前記整形部は、前記再構成フィルタとして第１の再構成フィルタと第２の再構成フィルタとを成形し、
前記第１の再構成フィルタは、エッジを保存するための第１の画像を発生するように整形され、第２の再構成フィルタは、ノイズを低減してエッジマップを発生するための第２の画像を発生するように整形され、
前記画像は、前記第１の画像と前記第２の画像と前記エッジマップとに基づいて発生される、
請求項１０記載のＸ線コンピュータ断層撮影装置。
【請求項１３】
前記整形部は、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするように前記再構成フィルタを整形する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項１４】
前記成形部は、前記再構成フィルタをシステムマトリクスに基づいて整形する、請求項１記載のＸ線コンピュータ断層撮影装置。
【請求項１５】
前記システムマトリクスが前記投影データの不均一サンプリングを指定する、請求項１４記載のＸ線コンピュータ断層撮影装置。
【請求項１６】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備する画像再構成方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、Ｘ線コンピュータ断層撮影装置及び画像再構成方法に関する。
【背景技術】
【０００２】
フィルタ補正逆投影（ＦＢＰ）アルゴリズムは、簡単かつ効率的であり、核医学、Ｘ線ＣＴ、及びＭＲＩにおいて画像を再構成するために使用される。ＦＢＰアルゴリズムは、計算効率が高いため、Ｘ線ＣＴ及び核医学の画像再構成において重宝されている。ＦＢＰは、低線量のＸ線で収集されたデータに適用されると多くのノイズを含む画像を発生してしまう。
【０００３】
一般的な傾向として、ＦＢＰアルゴリズムは、反復アルゴリズムに徐々に置換されている。ＦＢＰアルゴリズムは、ノイズを含む画像を発生してしまう。さらに、ＦＢＰアルゴリズムは、ノイズレベルを低減するためにノイズモデルを組み込むことが可能でないことが知られている。この点で、反復アルゴリズムは、投影ノイズモデルを組み込むことにより、ＦＢＰアルゴリズムよりもノイズがより少ない画像を発生することができる。しかしながら、反復アルゴリズムは、集中的な計算を必要とする。
【０００４】
計算要件にかかわらず、反復アルゴリズムは、最大事後確率（ＭＡＰ）を使用してノイズと分解能とのバランスがとれた画像を発生することができる。対照的に、ＦＢＰアルゴリズムは、ＭＡＰまたは事前の画像情報を活用することができないことが知られている。
【発明の概要】
【発明が解決しようとする課題】
【０００５】
目的は、ＦＢＰアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減可能なＸ線コンピュータ断層撮影装置及び画像再構成方法を提供することにある。
【課題を解決するための手段】
【０００６】
本実施形態に係るＸ線コンピュータ断層撮影装置は、少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、前記重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する整形部と、前記再構成フィルタに基づいて画像を発生する発生部と、を具備する。
【図面の簡単な説明】
【０００７】
【図１】本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図。
【図２】図１の再構成部の構成を示す図。
【図３】図２の再構成部による処理の典型的な流れを示す図。
【図４】本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図。
【図５】本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図。
【図６Ａ】本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図。
【図６Ｂ】本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図。
【図７】本実施形態に係る、ＦＢＰアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図。
【図８】本実施形態に係る、エッジ保存ＦＢＰ−ＭＡＰアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図。
【図９】図８のエッジマップの発生処理の典型的な流れを示す図。
【図１０】本実施形態に係る、事前画像を使用したＦＢＰ−ＭＡＰアルゴリズムに基づく画像再構成方法の典型的な流れを示す図。
【図１１】本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムを用いて画像再構成方法の典型的な流れを示す図。
【図１２】本実施形態に係るノイズ重み付けＦＢＰアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図。
【図１３Ａ】従来のＦＢＰフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【図１３Ｂ】β＝０を用いたｒＦＢＰ−ＭＡＰ再構成であるｒＦＢＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【図１３Ｃ】二次平滑化事前分布を使用し、β＝０．０００５を用いたｒＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【図１３Ｄ】実質的に図１３Ｂ及び図１３Ｃの重み付け加算である、エッジ保存事前分布を用いたｒＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【図１３Ｅ】エッジ保存ＦＢＰ−ＭＡＰ再構成において使用されるエッジマップを示す図。
【図１３Ｆ】ｖＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【発明を実施するための形態】
【０００８】
本実施形態は、ノイズ重み付けフィルタ補正逆投影（filtered back projection）法を使用して投影データを処理するためのＸ線コンピュータ断層撮影装置及び再構成処理方法に関する。
【０００９】
以下、図面を参照しながら、本実施形態に係わるＸ線コンピュータ断層撮影装置及び画像再構成方法について説明する。
【００１０】
図１は、本実施形態に係るＸ線コンピュータ断層撮影装置の構成を示す図である。図１に示すように、本実施形態に係るＸ線コンピュータ断層撮影装置は、マルチスライスＸ線ＣＴ装置またはマルチスライスＸ線ＣＴスキャナである。本実施形態に係るＸ線コンピュータ断層撮影装置は、ガントリ１００を装備している。ガントリ１００は、Ｘ線管１０１、円環状フレーム１０２、及び多列タイプあるいは二次元配列タイプのＸ線検出器１０３を含む。Ｘ線管１０１とＸ線検出器１０３とは、軸ＲＡの周りを回転するフレーム１０２上に被検体Ｓを横断して直径方向に設置される。回転部１０７は、被検体Ｓが軸ＲＡに沿って例示されたページ内または例示されたページ外に移動される間に、一回転当たり０．４秒等の高速でフレーム１０２を回転させる。
【００１１】
本実施形態に係るＸ線コンピュータ断層撮影装置は、高電圧発生装置１０９を備える。高電圧発生装置１０９は、Ｘ線管１０１がＸ線を発生するように、管電圧をＸ線管１０１に印加する。例えば、高電圧発生器１０９はフレーム１０２上に設置される。Ｘ線は、被検体Ｓに向けて放射される。Ｘ線検出器１０３は、被検体Ｓを透過したＸ線を検出するために、被検体Ｓを横断してＸ線管１０１の逆側に位置する。
【００１２】
さらに図１に示すように、本実施形態に係るＸ線コンピュータ断層撮影装置は、Ｘ線検出器１０３から検出された信号を処理するための他のデバイスをさらに含む。データ収集回路（ＤＡＳ）１０４は、Ｘ線検出器１０３からの信号出力をチャンネル毎に電圧信号に変換して、当該電圧信号を増幅し、増幅後の電圧信号をデジタル信号に変換する。Ｘ線検出器１０３とＤＡＳ１０４とは、１回転ごとに所定の総投影数（ＴＰＲＰ）を処理するように構成される。
【００１３】
上述のデータは、非接触データ伝送部１０５を介してガントリ１００の外部のコンソール内に収納された前処理部１０６に伝送される。前処理部１０６は、生データに対して感度補正等の前処理を施す。前処理が施されたデータは、投影データと呼ばれている。記憶部１１２は、投影データを記憶する。記憶部１１２は、再構成部１１４、表示部１１６、入力部１１５、及びスキャン計画サポート装置２００と共に、データ／制御バスを介してシステム制御部１１０に接続される。スキャン計画サポート装置２００は、スキャン計画を構築するために、撮影技法をサポートするための機能を有している。
【００１４】
本実施形態によれば、再構成部１１４は、ノイズ重み付けを用いたフィルタ補正逆投影（ＦＢＰ）技法に基づいて、記憶部１１２に記憶された投影データから画像を再構成する。あるいは、再構成部１１４は、ノイズ重み付けと基準画像等の事前分布入力（prior in）とも用いたフィルタ補正逆投影に基づいて、投影データから画像を再構成しても良い。再構成部１１４の上記の２つの実施形態のいずれか１つは、所定の反復再構成アルゴリズムに従って所定数の反復において特定の反復結果をエミュレートする追加の特徴を用いたフィルタ補正逆投影法に基づいて、投影データから画像を再構成する。
【００１５】
再構成部１１４は、ソフトウェアとハードウェアとの組合せにより構成され、特定の実装形態に限定されない。本実施形態に係る再構成部１１４は、核医学及び磁気共鳴撮像（ＭＲＩ）等の他のモダリティにも実装可能である。
【００１６】
図２は、本実施形態に係る再構成部１１４の構成を示す図である。図２に示すように再構成１１４は、ノイズ重み決定部１１４Ａ、フィルタ整形部１１４Ｂ、及び画像発生部１１Ｃを有する。ノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って投影データに重み付けを施すための重み値を決定する。フィルタ整形部１１４Ｂは、決定された重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。画像発生部１１４Ｃは、整形された再構成フィルタに基づいて画像を発生する。なお、この再構成部１１４の実装形態は単なる例示であることを言及しておく。以下、ノイズ重み決定部１１４Ａ、フィルタ整形部１１４Ｂ、及び画像発生部１１Ｃについて詳細に説明する。
【００１７】
図３は、再構成部１１４による処理の典型的な流れを示す図である。ステップＳ１００においてノイズ重み決定部１１４Ａは、少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定する。重み値は、投影データについての所定のデータ単位毎に決定される。データ単位としては、ビューであっても良いし、投影線であっても良い。ステップＳ１２０においてフィルタ整形部１１４Ｂは、ステップＳ１００で決定された重み値を含むパラメータに基づいて、少なくとも１つの再構成フィルタを整形する。本実施形態に係る整形のために、複数の再構成フィルタがオプションで使用されると同時に、これらのパラメータは重み値に限定されない。ステップＳ１２０において画像発生部１１４Ｃは、ステップＳ１２０において整形された再構成フィルタに基づいて画像を発生する。
【００１８】
ステップＳ１２０は、所定のパラメータに基づいてフィルタを整形する。上記のように、１つの例示的なパラメータは、投影データに含まれるノイズ成分の分散（以下、ノイズ分散と呼ぶ）である。しかしながら、本実施形態に係る重み値を含むパラメータは、ノイズパラメータに限定されない。例えば、再構成フィルタは、システムマトリクスに基づいて整形されても良い。より詳細には、システムマトリクスは、投影データの不均一サンプリングを指定する。従来のＦＢＰアルゴリズムはオブジェクトが均一にサンプリングされると仮定するが、可変のサンプリング戦略を使用することも可能である。例えば、重要な構造を含むある角度範囲内の信号は、第１の角度間隔でサンプリングされる。対照的に、当該領域外の別の角度範囲内の信号は、第１の角度間隔よりも大きい第２の角度間隔でサンプリングされる。
【００１９】
なおステップＳ１００は、投影データの所定のデータ部分の各々に関する重み値を決定する際に、ステップＳ１２０の前に実行または完了されなくても良い。すなわち、再構成フィルタが整形される度にオプションで重み値が決定されても良い。
【００２０】
さらに図３に示すように、ステップＳ１４０において画像発生部１１４Ｃは、ステップＳ１２０において整形された再構成フィルタを適用しながら、サブステップの所定のセットを実行する。１つの例示的なプロセスにおいては、さらに詳細に説明されるように、所定のＦＢＰ技法においてフーリエ変換の後に、整形されたフィルタを周波数領域内で投影データに適用することによって画像が再構成される。別の例示的なプロセスにおいては、所定のＦＢＰ技法において畳み込みとしてカスタマイズされたフィルタを空間領域内で投影データに適用することにより画像が再構成される。
【００２１】
図４は、本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図である。ビュー毎の重み付け方式において各ビューは、概して、特定のビューにおける最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、重み値が割り当てられる。図４は、３つの位置、すなわち、位置Ａ、Ｂ、及びＣにおけるＸ線源、ならびに位置Ａ’、Ｂ’、及びＣ’における３つの対応する位置におけるＸ線検出器を例示する。例えば、Ｘ線源が位置Ａにあるとき、Ｘ線検出器は、対応する位置Ａ’にある。位置Ａ’においてＸ線検出器により収集された投影データはビューからなり、重み値ｗ（θ1）は、当該ビューに関する所定の重み付け方式に従ってビュー角度θ1に応じて決定される。上記のように、ビュー単位の重み値の決定は、１つの例示的なプロセスでは、再構成フィルタの整形に先立って実行および完了される。このビュー単位の重み値の決定は、他の例示的なプロセスでは、再構成フィルタが整形される度に決定される。
【００２２】
本実施形態において投影データに含まれるノイズ分散は、ビュー単位の重み付け方式を使用してモデル形成可能である。一般に、ノイズの影響を補償する間に画像を再構成するための典型的な手法は、（１）式により規定される、ノイズ分散に応じて重み付けされた目的関数を最小化することである。
【数１】
2014023936.tif 000003
【００２３】
（１）式中、Ａは投影行列であり、Ｘは列ベクトルとして記述される画像配列であり、Ｐは列ベクトルとして記述される投影配列であり、Ｗは重み付けノルムを定義するノイズ重み付け行列である。この目的関数を最小化するために、（２）式に規定される反復ランドウェーバー（Landweber）アルゴリズムが使用される。
【数２】
2014023936.tif 000004
【００２４】
（２）式中、ＡTは逆投影行列であり、Ｘ(k)は第ｋ番目の反復において推定された画像であり、αはステップサイズであり、α＞０である。この再帰的表現は、（３）式のように非再帰的表現に書き直すことができる。
【数３】
2014023936.tif 000005
【００２５】
（ＡTＷＡ）-1＝Ａ-1Ｗ-1（ＡT）-1が存在すると仮定すると、（３）式は、（４）式のように簡素化される。
【数４】
2014023936.tif 000006
【００２６】
（４）式は、ノイズ重み付けＦＢＰアルゴリズムに相当し、周波数領域内における（４）式のフィルタ関数は、（５）式のように表される。
【数５】
2014023936.tif 000007
【００２７】
（５）式中、ωは周波数分散であり、ｗ（ｖｉｅｗ）はビュー単位で決定される重み値である。すなわち、ｗ（ｖｉｅｗ）は、特定のビュー角度におけるノイズ関連の重み関数である。αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、ｋは反復アルゴリズムにおいて反復数をエミュレートするパラメータである。ｖＦＢＰとして示されるビュー毎の重み付けＦＢＰアルゴリズムは、シフト不変の点応答関数（ＰＲＦ）を有する。
【００２８】
上記の実施形態において再構成フィルタは、所定の反復再構成アルゴリズムに従って、所定数ｋ回の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けＦＢＰアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
【００２９】
図５は、本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図である。投影線毎の重み付け方式では、特定のビュー内の最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、各投影線に重み値が割り当てられる。図５は、所定の位置におけるＸ線源、ならびに、対応する位置におけるＸ線検出器を例示する。投影線がＸ線源から投影されるとき、投影線は対応する位置におけるＸ線検出器に達する。投影データは、θにおけるビュー角度内の複数の投影線からなる。例えば、重み値ｗ（θ，ｔ1）は、ビュー角度θに応じ、かつ当該投影線に関する所定の重み付け方式に従って、特定の投影線ｔ1について決定される。上記の通り、投影線の重み値の決定は、１実施形態において、再構成フィルタを整形するのに先立って実行および完了される。投影線重み値は、他の実施形態において、再構成フィルタとして決定される。
【００３０】
シフト不変ＰＲＦに関する要件が除去されると、投影線毎の重み付け方式が構築され、その周波数領域における窓化されたランプフィルタが上記の（５）式を変形することにより、（６）式のように記述される。
【数６】
2014023936.tif 000008
【００３１】
（６）式中、重み値、すなわち、重み係数ｗ（ｒａｙ）は、投影線で決定される重み計数である。ｗ（ｒａｙ）は、関連する投影線のノイズ分散に応じて決定される。ｋは、所定の反復再構成アルゴリズムに従って、特定の反復結果をエミュレートするための所定数の反復を示すパラメータである。Ｈ（ω，ｒａｙ）の逆フーリエ変換がｈ（ｔ，ｒａｙ）であるとすると、ｈ（ｔ，ｒａｙ）はフィルタの空間領域におけるカーネルである。ｐ（ｒ，θ）をビューθ及び投影線ｒにおける投影であるとし、ｑ（ｒ，θ）をフィルタリング後の投影であるとする。この場合、ｑ（ｒ，θ）は、以下の（７）式のような積分で規定される。
【数７】
2014023936.tif 000009
【００３２】
カーネルｈ（ｔ，ｒ）はｒに依存するため、この積分は畳み込み（コンボルーション）ではない。フィルタリングが（７）式のように空間領域で実施される場合、計算コストは畳み込みと同等である。この新しいｒＦＢＰアルゴリズムにおける逆投影ステップは、従来のＦＢＰアルゴリズムの逆投影ステップと同じである。従ってｒＦＢＰアルゴリズムは、計算的に効率的である。他方で、画像領域フィルタリングは、投影線毎のノイズ重み付けＦＢＰ−ＭＡＰ（ｒＦＢＰ−ＭＡＰ）アルゴリズムに関する実装形態で後に説明されるように、伝達関数を用いて２次元フーリエ領域内で実施される。
【００３３】
重み係数を割り当てるための典型的な手法は、ｗ（ｒａｙ）を投影線計測のノイズ分散の逆数にすることである。この手法は、最適化問題に関する目的関数として、尤度関数（すなわち、結合確率密度関数）を使用して補われる。ノイズを含む計測に比してノイズを含まない計測を信頼すべきであるという根本原理を採用する場合、ノイズがより少ない計測により大きなｗ（ｒａｙ）が割り当てられ、ノイズがよい多い計測により小さなｗ（ｒａｙ）が割り当てられという制限下において任意に重み係数を割り当てることができる。
【００３４】
本実施形態に係るｒＦＢＰアルゴリズムを実施するための最も容易かつより効率的な様式は、空間領域において投影をフィルタリングするために、（７）式を使用することである。現在、本発明者らは積分カーネルｈ（ｔ，ｒ）に関する解析式を有していない。各ビューθに関してｑ（ｒ，θ）を計算するための代替様式は、周波数領域において、以下のように方程式（７）を実施することである。
【００３５】
ステップ１：ｒに関するｐ（ｒ，θ）の１次元フーリエ変換を見出して、Ｐ（ω，θ）を取得する。
【００３６】
ステップ２．ａ：各投影線「ｒａｙ」に関して、重み係数ｗ（ｒａｙ）を割り当て、方程式（６）で表されるように周波数領域伝達関数を形成する。実装形態において、ωは周波数指数であり、整数をとる。
【００３７】
ステップ２．ｂ：積、Ｑray（ω，θ）＝Ｐ（ω，θ）Ｈ（ω，ｒａｙ）を見出す。
【００３８】
ステップ２．ｃ：ωに関するＱray（ω，θ）の１次元逆フーリエ変換を見出して、ｑ（ｒａｙ，θ）を取得する。
【００３９】
上記の実施形態においてステップ１は、各ビューの列内の全てのｒを処理し、一方、ステップ２．ａから２．ｃは一度に１つの投影線のみを処理する。ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズムでは、単一の重み値、すなわち、重み係数ｗ（ｖｉｅｗ）がビュー内の全ての投影線に割り当てられる。対照的に、投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムでは、別個の重み値、すなわち、重み係数ｗ（ｒａｙ）がビュー内の投影線の各々に割り当てられる。結果として、ノイズ重み付け方式は、ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズムよりも、投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムにおいてより正確である。
【００４０】
ビュー毎の重み付けＦＢＰ（ｖＦＢＰ）アルゴリズム及び投影線毎の重み付けＦＢＰ（ｒＦＢＰ）アルゴリズムに関して上で説明されたように、本実施形態は、画像を発生する前に、重み値を含むパラメータに基づいて少なくとも１つの再構成フィルタを整形する。さらに、上記の実施形態では、再構成フィルタは所定の反復再構成アルゴリズムに従って、ｋで記された所定数の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けＦＢＰアルゴリズムまたは投影線ごとの重み付けＦＢＰアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
【００４１】
その後、整形された再構成フィルタに基づいて、画像を発生するための工程が実行される。一般に、画像再構成におけるノイズ低減は、画像再構成、特に逆投影が加算処理であり、投影データの適切な重み付けが実質的に逆投影の分散を削減するという概念に基づいている。
【００４２】
図６Ａは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図である。ＦＢＰアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行される前に、再構成フィルタが投影データに適用される。ステップＳ１４０Ａに先立って再構成フィルタは、ノイズ重み値を含むパラメータに基づいて予め整形されている。ステップＳ１４０Ａにおいて整形された再構成フィルタが投影データに適用される。ステップＳ１４０Ａの前に投影データがフーリエ変換に供される場合、フィルタリングされた投影データは、ステップＳ１４０Ｂの前に逆フーリエ変換を受ける。ステップＳ１４０Ｂにおいてフィルタリングされた投影データは、逆投影される。逆投影により、フィルタリングされた投影データに基づいて画像が発生される。
【００４３】
図６Ｂは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図である。ＦＢＰアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行された後に、再構成フィルタが投影データに適用される。ステップＳ１４０Ｂの前に投影データがフーリエ変換に供される場合、投影データは、ステップＳ１４０Ｂの前に逆フーリエ変換を受ける。ステップＳ１４０Ｂにおいて投影データは、逆投影され、画像が発生される。整形された再構成フィルタは、逆投影されたデータに適用される。ステップＳ１４０Ａに先立って、再構成フィルタはノイズ重み値を含むパラメータに基づいて予め整形されている。ステップＳ１４０Ａにおいて、整形された再構成フィルタが投影データに適用される。
【００４４】
図７は、ＦＢＰアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図である。図７の処理においては、逆投影の前段においてノイズ重み付けフィルタリングがフーリエ領域で行わると仮定する。図７に示すように、投影データＰＤは、複数のビューにおいて収集される。投影データＰＤは、ビュー又は投影線単位で処理される。投影データＰＤは、変換されたデータがフーリエ領域または周波数領域に存在するように、フーリエ変換ＦＴに施される。ランプフィルタは、ノイズ重み値を含む所定パラメータのセットに従って整形される。これにより、重み付けランプフィルタＷＦが形成される。いくつかの実施形態では、複数のランプフィルタは、複数の重み付けランプフィルタＷＦになるように整形される。これらパラメータは、ノイズ重みパラメータに限定されず、他のパラメータを含んでも良い。この例示的な実施形態では、変換後の投影データは、所定のＦＢＰアルゴリズムにおいてノイズ重み付けランプフィルタＮＷＦに従ってフィルタリングされる。周波数領域内でフィルタリングされたデータは、投影データが逆投影ＢＰにおいて空間領域内で逆投影されるように、逆フーリエ変換ＩＦＴを受ける。あるいは、周波数領域内でノイズ重み付けランプフィルタＮＷＦを変換後のデータに適用する代わりに、別のノイズ重み付けランプフィルタＮＷＦ２が画像に適用されても良い。この場合、画像に適用されるノイズ重み付けランプフィルタＮＷＦ２は、変換後のデータに適用されるノイズ重み付けランプフィルタＮＷＦとは異なる。
【００４５】
逆投影ＢＰは、画像発生の間、重み関数を使用して上述の可変角度サンプリングを正規化する。例えば、重要な構造を含む一定の角度範囲内の信号は、第１の角度間隔でサンプリングされる。当該領域外の他の角度範囲内の信号は、第１の角度間隔よりもより大きい第２の角度間隔でサンプリングされる。逆投影ＢＰは、サンプリング角度全体に亘る積分であり、角度サンプリング密度補償は、逆投影積分において本質的にヤコビアン係数である正規化係数を使用することによって達成される。数学的に、逆投影画像は、（８）式のように表現される。
【数８】
2014023936.tif 000010
【００４６】
（８）式中、初めに投影が逆投影される場合、ｐは角度θにおける投影である。あるいは、ＦＢＰアルゴリズムが使用される場合、ｐは、角度θにおいてフィルタリングされた投影であり、ｔは、検出器ビンの位置を示す。方程式（８）の離散的な表現は、（９）式のように表される。
【数９】
2014023936.tif 000011
【００４７】
（９）式中、ｎは検出器の位置指数であり、ｍは投影角度指数であり、Ｍは投影が収集されたビューの総数であり、「ＩＮＴ」は最近傍補間を示すために使用される。実際には、「ＩＮＴ」関数は使用せず、２つの近傍検出器ビンの間の線形補間が使用される。
【００４８】
均一でない角度サンプリングに関して、（９）式で表されるような簡易的な逆投影が施される場合、当該逆投影は、単位角度当たりのビュー数である密度関数ｄ（０）に依存する。例えば、０＜θ＜π／４に関して、サンプリング間隔が１°であり、π／４＜θ＜πに関して３°である場合、密度関数は、（１０）式により表現され、逆投影（９）式は、（１１）式のように変更される。
【数１０】
2014023936.tif 000012
【００４９】
【数１１】
2014023936.tif 000013
【００５０】
この場合、サンプリング密度関数は、実際の角度θの代わりに、角度指数ｉｎの関数である。
【００５１】
図８は、本実施形態に係る、所定のエッジ保存ＦＢＰ−ＭＡＰアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図である。図８の画像再構成方法は、投影データの組を逆投影して画像を再構成する前段に、フーリエ領域においてノイズ重み付けフィルタリングが実施される。最終的に、２つの画像に基づいて最終画像が発生される。
【００５２】
図８に示すように、投影データＰＤは、複数のビューにおいて収集される。投影データＰＤは、ビュー単位で処理されても良いし、投影線単位で処理されても良い。投影データＰＤは、変換後のデータがフーリエ領域または周波数領域にあるようにフーリエ変換ＦＴを受ける。二つのランプフィルタは、ノイズ重み値を含む所定のパラメータのセットに従って個別に整形され、重み付けランプフィルタＷＦ−β0及びＷＦ−βに整形される。ランプフィルタは、正規化によってさらに整形可能である。このようにフィルタは、値βの関数である。値βが大きいほど、画像の平滑化がより強く行われる。例えば、平滑化が行われない場合あるいはボケの無い場合、β値はゼロなるように選択される。変換後の投影データは、所定のＦＢＰアルゴリズムにおいてノイズ重み付けされると良い。ノイズ重みづけされた投影データは、正規化されたランプフィルタＷＦ−β0及びＷＦ−βによってそれぞれフィルタリングされる。その後、周波数領域においてフィルタリングされた投影データは、逆フーリエ変換ＩＦＴに施される。逆フーリエ変換後の投影データは、２つの画像Ｘ0及びＸβが発生されるように空間領域内でそれぞれ逆投影される。平滑化を用いないで逆投影された、あるいは、最小限の平滑化を用いて逆投影された画像Ｘ0は、予め用意されたエッジマップＥＭＡＰ内のエッジマップ情報Ｂを抽出するために使用される。最小限に平滑化された画像Ｘ0、フィルタリングされた画像Ｘβ、及び抽出されたエッジ情報に基づいて一定の規則に従って最終画像ＦＩＭＧが発生される。
【００５３】
図８に例示される実施形態においてエッジは、下記の（１２）式に示す目的関数においてベイジアン正規化項を使用してノイズ除去の間に保存される。（１２）式のペナルティ関数ｇ（Ｘ）は、注目画素近傍における画素値の段差を計測する。
【数１２】
2014023936.tif 000014
【００５４】
（１２）式中、Ａは投影行列であり、Ｘは列ベクトルとして書き込まれた画像配列であり、Ｐは列ベクトルとして書き込まれた投影配列であり、Ｗは各投影に関する重み関数からなる対角行列であり、ｂは、（１２）式の第１項である忠実度項
【数１３】
2014023936.tif 000015
【００５５】
に対するベイジアン項ｇ（Ｘ）の重要性を調整する相対的な重み係数である。
【００５６】
ペナルティｇ（Ｘ）は、画素値段差の二次関数である。当該画素値段差は、当該画素の画素値とその近傍の画素の平均値との間の差として計算される。大きな画素値段差は大きなペナルティに対応する。当該ペナルティ関数は平滑化画像をもたらす。
【００５７】
エッジを過度に平滑化せずに画像ノイズを抑えるために、ペナルティ関数ｇ（Ｘ）は、エッジ及び非エッジ領域に関して、異なる特性を有すべきである。そのような能力を有するペナルティ関数として、閾値を有するフーバ関数が用いられる。画素値段差が閾値より小さく、画素値段差がノイズとして分類される場合、フーバ関数は二次式であり、平滑化を施す。画素値段差が閾値より大きく、画素値段差がエッジとして分類される場合、フーバ関数は線形であり、より小さい平滑化を施す。
【００５８】
部分毎の一定画像を奨励する別の一般的な事前分布は、ＴＶ（全変動）計測である。離散的な一次元関数ｇ（Ｘ）の場合、ｇ（Ｘ）のＴＶ計測の勾配は３つの値だけをとる。すなわち、ｇ（Ｘ）がその左右の近傍より大きい場合、正の値、ｇ（Ｘ）がその左右の近傍よりも小さい場合、負の値、ｇ（Ｘ）がその左右の近傍の間である場合、０である。勾配タイプの最適化アルゴリズムでは、正の勾配値は、ｇ（Ｘ）の値を下方に押さえ込み、負の勾配値はｇ（Ｘ）の値を上方に押し上げ、０の勾配値は何の変更も行われていないことを意味する。従ってＴＶ事前分布は、単調関数を奨励し、振幅を抑え、鮮明なエッジを保存する。
【００５９】
図８に例示されるように、実施形態に係る画像再構成方法は、投影データにエッジ保存ＦＢＰ−ＭＡＰアルゴリズムを適用する。βは、１つのノイズ重み付けランプフィルタについてのβ0、すなわち、ゼロに設定される。（１２）式を参照すると、β0ｇ（Ｘ）はゼロであり、何の平滑化も施されない。その後、何の平滑化も施されずに画像Ｘ0が再構成される。他方、他のパラメータセットの場合、βは、他のノイズ重み付けランプフィルタに関するβ、すなわり、非ゼロに設定される。βにより一定量の平滑化が施され、画像Ｘβが再構成される。これにより、２つの画像Ｘβ及びＸ0が再構成される。すなわち、第１の画像Ｘβは、非ゼロβ値を用いたｒＦＢＰ−ＭＡＰ再構成に基づいて再構成され、第２の画像Ｘ0はゼロβ値を用いたｒＦＢＰ再構成に基づいて再構成される。
【００６０】
図９は、図８のエッジマップＢの発生処理の典型的な流れを示す図である。画像Ｘ0はゼロβ値を用いたｒＦＢＰ再構成法により再構成される。ステップＳ２００において、本質的に、画像勾配のノルムであるソーベルカーネルを使用して、エッジ検出によって画素値段差が評価される。ソーベルカーネルによって生成されたエッジマップは、次いで、ステップＳ２２０において、所定の閾値に従って、値０または値１を用いて二値画像に変換される。閾値を超える画像値は、値１に設定され、閾値未満の画像値は値０に設定される。その後、ステップＳ２４０において、所定のローパスフィルタによって二値画像をボケさせ、エッジマップＢを発生させる。二値画像をボケさせる理由は、二値画像に含まれるエッジ領域と非エッジ領域との境界を滑らかにするためである。最終画像は、下記の（１３）式に示すような、ＸbとＸ0との一次結合により発生される。
【数１４】
2014023936.tif 000016
【００６１】
（１４）式中、「．*」は、ＭＡＴＬＡＢ（登録商標）で定義される点毎の乗算を表す。
【００６２】
図１０は、事前画像を使用したＦＢＰ−ＭＡＰアルゴリズムに基づく画像再構成方法の典型的な流れを示す図である。一般に、所与の事前画像
【数１５】
2014023936.tif 000017
【００６３】
を加えることによって、ＦＢＰ−ＭＡＰアルゴリズムに基づいて画像を再構成するために、目的関数（８）は、下記の（１４）式により一般化される。
【数１６】
2014023936.tif 000018
【００６４】
反復的な中間解は、（１５）式のように表される。
【数１７】
2014023936.tif 000019
【００６５】
等価ＦＢＰ−ＭＡＰアルゴリズムは、下記の（１６）式と（１７）式とに表される２つのフーリエ領域フィルタ関数を有する。
【数１８】
2014023936.tif 000020
【００６６】
【数１９】
2014023936.tif 000021
【００６７】
図１０に示すように、ＦＢＰ−ＭＡＰアルゴリズムは、投影データＰＤと事前画像ＰＩとの２つの入力を有する。最終画像ＦＩは、投影データＰＤに由来する第１の画像と、事前画像ＰＩに由来する第２の画像とに基づく加算画像である。ＦＢＰ−ＭＡＰ再構成は、（１６）式のような検出器フーリエ領域フィルタとその後に行われる逆投影とに基づいて第１の再構成画像を発生する。画像フーリエ領域処理ＩＦＤは、（１７）式のように画像フーリエ領域フィルタを使用して第２の画像を発生する。画像フーリエ領域フィルタは、２次元であっても良いし３次元であっても良い。
【００６８】
図１１は、本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムを用いて画像再構成方法の典型的な流れを示す図である。一般に、再構成フィルタは、所定数の反復において特定の反復結果を模倣（エミュレート）する。再構成フィルタは、所定の反復再構成アルゴリズムに従って、重み値を含むパラメータに基づいて整形される。ノイズ重み付け反復エミュレーションＦＢＰアルゴリズム２００Ｄは、ステップ２０１において投影データが収集される。投影データは、撮影されているオブジェクトの線積分を表すラドン変換の形式で収集される。投影データは、図１に例示されるようなＸ線コンピュータ断層撮影装置やＰＥＴ装置、ＳＰＥＣＴ装置等の任意のモダリティにより収集される。
【００６９】
ステップ２０２において反復数やノイズ重み値、及びシステムマトリクス等を含むパラメータの組合せに基づいて、少なくとも１つの再構成フィルタが整形される。これらのパラメータは、任意に選択可能である。例えば、パラメータｋは、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするために選択される。また、投影データに適用される再構成フィルタを整形するための反復パラメータｋと共に、予め決定された重み値Ｗが使用されても良い。その他のフィルタパラメータは、様々なアプリケーションに対処できるように、反復アルゴリズム内のステップサイズと、ベイジアン条件の相対的重み付けとを含む。
【００７０】
ステップＳ２０２においてパラメータの組合せに従って再構成フィルタが整形された後、ステップ２０３において、整形された再構成フィルタを投影データに適用する。フィルタが施されたデータは、ステップ２０４において逆投影される。逆投影により画像が発生する。ステップＳ２０５においてエミュレーションを反復するか否かが決定される。エミュレーションを反復する場合、エミュレーションの反復のために再びステップＳ２０２に進む。ステップＳ２０５においてエミュレーションを反復しない場合、本実施形態に係るノイズ重み付け反復エミュレーションＦＢＰアルゴリズムは終了する。
【００７１】
図１２は、本実施形態に係るノイズ重み付けＦＢＰアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図である。従来のＣＴ撮影においてノイズ制御は、事前フィルタリングまたは事後フィルタリングを使用して実行されている。対照的に本実施形態においてノイズは、ノイズ重み付け及び反復エミュレートされたＦＢＰアルゴリズムを使用して画像再構成において実質的に低減される。画像再構成におけるノイズ低減は、画像再構成（特に、逆投影）が加算処理であり、投影に適切な重み付けを実行することにより逆投影の分散を削減できるという理念に基づいている。
【００７２】
図１２に示すように、２つの投影、すなわち、計測データが２本の投影線Ｐ１及びＰ２によって表され、解、すなわち画像Ｉは２本の投影線Ｐ１及びＰ２の交点である。横軸Ｘ1及び縦軸Ｙ1は、画像Ｉにおける２つの変数を表す。１つの計測が画像から行われる場合、２つの画像画素は１本の直線に沿ったいずれかの場所の値をとることができる。２つの独立した計測が行われる場合、それら２本の直線の交点である、それら２つの画素の値が決定される。第３の計測が行われ、すべての計測がノイズを含む場合、それらの３本の直線は、一点において交差しない可能性がある。一般に、画像が多くの画素からなり、多くの矛盾するノイズを伴う計測が行われた場合、それらの計測は単一点で交わらない可能性があるため、撮像問題の解が存在しない可能性がある。
【００７３】
数学では、ノイズ計測が矛盾する場合、それらの計測は「投影演算子の範囲内にない」と見なされる。他方で、計測が矛盾しない場合、それらの計測は投影演算子の範囲内にマッピングされる。実際には、投影演算子の範囲内のデータ成分だけが画像領域内に逆投影されるため、このマッピングはＦＢＰアルゴリズムにおいて逆投影によって保証される。すなわち、この投影の矛盾はＦＢＰ再構成画像に寄与しない。一般性を失わずに、直線によって示されるように、多くの画素を有する画像の場合ですら、計測は矛盾しないものの、依然として雑音を伴い、投影演算子の範囲内にある状況を仮定することができる。反復アルゴリズムとは異なり、ＦＢＰアルゴリズムは、ノイズを伴う投影のセットから、ノイズを伴う単一の解を提供することが可能である。
【００７４】
図１２に示すように、ノイズがより少ない画像を発生するため、反復アルゴリズムの複数の反復ステップをエミュレートするようにＦＢＰアルゴリズムは修正される。解軌道はｋ＝０の初期値からｋ＝∞の最終的な解に向けて示される。ＦＢＰアルゴリズムのその後の反復は、結果として生じる、より少ないノイズを伴って画像を移動させる。反復アルゴリズムは、通常、真の解に達成する前に終了する。これは、計測データが低い光子計数によるノイズを含むため、真の解決策がノイズを伴うものである場合に特に当てはまる。
【００７５】
ｋ＝０からの１つの解軌道は、真の解決策、すなわち、画像Ｉに向けた直線であり、この直線は、エミュレートされた反復がノイズ分散によって重み付けされないことを示す。直線上の点ＮＮＷ１及びＮＮＷ２は、それぞれ、第１の所定数の反復および第２の所定数の反復の後の中間結果を示す。この例では、第２の所定数の反復は、第１の所定数の反復よりも大きく、中間結果ＮＮＷ２は、第２の所定数の反復の後にエミュレートされる。他方で、中間結果ＮＷ１、ＮＷ２、及びＮＷ３は曲線上に位置し、この曲線は、それらのエミュレートされた反復がノイズ分散によって重み付けされていることを示す。中間結果ＮＷ１、ＮＷ２、及びＮＷ３は、それぞれ、第３の所定数の反復、第４の所定数の反復、及び第５の所定数の反復の後である。この例では、第５の数は、第３の数の反復及び第４の数の反復より大きく、第３の数は最小である。中間結果ＮＷ３は、第５の反復の後にエミュレートされる。この解軌道は、重み係数のセットによって決定される。１つの重み係数がそれぞれの計測に割り当てられる。より大きな係数は、関連する計測がより信頼できるものであることを意味する。一般的な手法は、計測のノイズ分散の逆数に比例する重み係数を割り当てることである。
【００７６】
実際には、真のノイズを伴う解の近傍にあるノイズがより少ない疑似解が受け入れられる。疑似解の手順は、ノイズ分散重み値と反復指数ｋとを含むパラメータの所定のセットによってパラメータ化される。一般に、より小さいｋを伴う疑似解は、より平滑な画像に対応し、より大きなｋを伴う疑似解は、より鮮明であるが、より多くのノイズを伴う画像に対応する。
【００７７】
以下、本実施形態に係る上記の特徴の特定の組合せに関する追加の実施例が記述される。一実施例は、ｒＦＢＰ−ＭＡＰによって指定される、最大事後確率（ＭＡＰ）アルゴリズムを用いた、投影線毎のノイズ重み付けＦＢＰが利用される。ｒＦＢＰ−ＭＡＰアルゴリズムに関する導出ステップは、ビューごとの重み付けＦＢＰ−ＭＡＰアルゴリズムに関する導出ステップと非常に類似している。幾つかの主な導出ステップについて以下に簡単に説明する。目的関数は、（１２）式と同じである。ベイジアン条件ｇ（Ｘ）が二次式である場合、ｇ（Ｘ）の勾配は、幾つかの行列Ｒに関してＲＸの形であり、解は、（１８）式により表されるランドウェーバーアルゴリズムによって得られる。
【数２０】
2014023936.tif 000022
【００７８】
（１８）式に示される反復式は、等価の非反復式を有する。初期状態がゼロであるとき、非再帰的数式は、（１９）式により表される。
【数２１】
2014023936.tif 000023
【００７９】
（ＡTＷＡ＋βＲ）-1が存在する場合、（１９）式中に示される総和は、（２０）式により表される閉形式を有する。
【数２２】
2014023936.tif 000024
【００８０】
（ＡTＷＡ）-1が存在する場合、（ＡTＷＡ＋βＲ）-1＝［Ｉ＋β（ＡTＷＡ）-1Ｒ］-1（ＡTＷＡ）-1である。十分なサンプリングが逆ラドン行列Ａ-1の存在をもたらすことが仮定される場合、（ＡTＷＡ）-1＝Ａ-1Ｗ-1（ＡT）-1及び（ＡTＷＡ）-1ＡTＷ＝Ａ-1となる。従って（２０）式は、（２１）式のように簡略化可能である。
【数２３】
2014023936.tif 000025
【００８１】
不十分なビュー数等の理由でＡ-1が存在しない場合、ＦＢＰアルゴリズムは、投影演算子Ａの疑似逆数であり、（ＡTＷＡ）+ＡTＷ＝Ａ+であることが確認でき、式中、「＋」は疑似逆数を示し、Ａ+はランプフィルタリングの後に逆投影ＡTが続くことである。従って（２１）式は、依然として、疑似逆数「＋」によって置換することによって逆数「−１」を可とする。以下、この逆数及び疑似逆数を識別しないこととする。「Ａ-1Ｐ」の表記は、従来のＦＢＰ再構成を表すために使用されることになる。
【００８２】
（２１）式において、［Ｉ−（Ｉ−αＡTＷＡ−αβＲ）k］は、反復アルゴリズムの第ｋ番目の反復におけるウィンドウ効果を表す。ｋが無限大になる傾向があるため［Ｉ−（Ｉ−αＡTＷＡ−αβＲ）k］は、識別行列になる。［Ｉ＋β（ＡTＷＡ）-1Ｒ］-1は、ベイジアン正規化の効果を表す。β＝０のとき、ベイジアン修正は効果的ではない。（２１）式は、反復アルゴリズムがノイズ重み付け及びベイジアン正規化をどのように処理するかの見識と理解を提供し、（２１）式は、同じタスクを実行するためのＦＢＰタイプのアルゴリズムの構築をもたらし得る。
【００８３】
最後に、ＦＢＰタイプのアルゴリズムが取得されるように、フーリエ領域表記が行列式から導出される。行列Ａ-1は、１次元ランプフィルタリングの後に逆投影が続くか、又は逆投影の後に２次元ランプフィルタリングが続くように処理される。行列式（２１）が「逆投影の後にフィルタリングする」アルゴリズムと見なされるとき、（２１）式における画像領域フィルタリングは、（２２）式に示す伝達関数を用いて２次元フーリエ領域において実行される。
【数２４】
2014023936.tif 000026
【００８４】
（２２）式中、ｖx及びｖyは、それぞれ、ｘ及びｙに関する周波数である。
【数２５】
2014023936.tif 000027
【００８５】
は２次元周波数ベクトルである。
【数２６】
2014023936.tif 000028
【００８６】
よく知られているように、投影演算子Ａが２次元空間内の線積分（すなわち、ラドン変換）であり、ＡTが演算子（すなわち、逆投影変換）である場合、投影と逆投影との組み合わされた演算子ＡTＡは、元の画像の２次元カーネル１／ｒとの２次元畳み込みである。この場合、ｘ−ｙデカルト座標において
【数２７】
2014023936.tif 000029
【００８７】
である。１／ｒの２次元フーリエ変換は、
【数２８】
2014023936.tif 000030
【００８８】
である。（２２）式においては、Ｒのフーリエ領域等価（Fourier domain equivalence）はＬ2Dによって示される。
【００８９】
中央断面定理を使用することによって、画像領域の２次元フィルタリングは、投影領域の１次元フィルタリングに変換される。２次元画像のフーリエ領域における伝達関数（（２２）式）は、下記の（２３）式のように１次元投影のフーリエ領域における伝達関数に等しい。
【数２９】
2014023936.tif 000031
【００９０】
（２３）式中、ωは検出器に沿った変数に関する周波数である。ω＝０であるとき、
【数３０】
2014023936.tif 000032
【００９１】
を定義する。（２３）式において、Ｌ1Dは（２２）式のＬ2Dの中央断面である。ｒＦＢＰ−ＭＡＰルゴリズムでは、
【数３１】
2014023936.tif 000033
【００９２】
は、１次元フーリエ領域で表現される修正ランプフィルタである。反復効果は、
【数３２】
2014023936.tif 000034
【００９３】
によって特徴付けられる。ベイジアン効果は、
【数３３】
2014023936.tif 000035
【００９４】
によって特徴付けられる。最終的な（ｋ＝∞）解は、β＝０の場合、ノイズ重み付けｗに依存しないが、β≠０の場合、ノイズ重み付けｗに依存する。ｒＦＢＰ−ＭＡＰアルゴリズムに関する逆投影について修正は必要とされない。
【００９５】
（２２）式及び（２３）式のノイズに依存する重み係数ｗは、定数ではない。ビュー単にのノイズ重み付けの場合、ｗはビュー角度に応じ、ｗ＝ｗ（ｖｉｅｗ）である。投影線単位のノイズ重み付けの場合、ｗは投影線に応じ、ｗ＝ｗ（ｒａｙ）である。重み係数を割り当てるための一般的な手法は、ｗ（ｒａｙ）を投影線計測のノイズ分散の逆数に比例させることである。この手法は、最適化問題に関する目的関数として、尤度関数（すなわち、結合確率密度関数）を使用することによって補われる。根本原理は、ノイズのより大きい計測よりもノイズのより少ない計測を信頼すべきであるということである。ｋが無限大になる傾向があり、βがゼロになる傾向があるとき、修正ランプフィルタ（（２３）式）は従来のランプフィルタになり、ｒＦＢＰ−ＭＡＰアルゴリズムは、従来のＦＢＰアルゴリズムに低減される。
【００９６】
ｒＦＢＰ−ＭＡＰアルゴリズムを実施する際、
【数３４】
2014023936.tif 000036
【００９７】
の逆フーリエ変換は、１次元修正ランプフィルタの空間領域カーネルである
【数３５】
2014023936.tif 000037
【００９８】
になると仮定される。ｒＦＢＰ−ＭＡＰアルゴリズムでは、ｗは投影線の関数であるため、ｗ＝ｗ（ｔ，θ）である。ｐ（ｔ，θ）をビューθ及び検出器上の位置ｔにおける投影にし、ｑ（ｔ，θ）をフィルタリングされた投影にする。この場合、ｑ（ｔ，θ）は、カーネル
【数３６】
2014023936.tif 000038
【００９９】
がθ及びｔに依存するため、畳み込みではない、以下の（２４）式に示す積分により規定される。
【数３７】
2014023936.tif 000039
【０１００】
フィルタリング処理が（２４）式のように空間領域で実施される場合、計算コストは畳み込みと同じである。従って、ｒＦＢＰ−ＭＡＰアルゴリズムは計算的に効率的である。ｒＦＢＰ−ＭＡＰアルゴリズムを実施するための最も容易で最も効率的な方法は、空間領域で投影をフィルタリングするために（２４）式を使用することである。
【０１０１】
現在、積分カーネル
【数３８】
2014023936.tif 000040
【０１０２】
に関する解析式は存在しない。ｑ（ｔ，θ）を計算する代替的な方法は、以下のように、周波数領域内の各ビュー角度θにおいて（２４）式を実施することである。
【０１０３】
ステップ１：ｔに関するｐ（ｔ，θ）の１次元フーリエ変換を見出して、Ｐ（ω，θ）を取得する。
【０１０４】
ステップ２：各投影線ｔについて重み係数ｗ（ｔ）を割り当て、周波数領域伝達関数（（２３）式）を形成する。ωは離散的な周波数指数であり、０、１、２、．．．などの整数値をとる。
【０１０５】
ステップ３：積
【数３９】
2014023936.tif 000041
【０１０６】
を見出す。
【０１０７】
ステップ４：ωに関するＱ1（ω，θ）の１次元逆フーリエ変換を見出して、ｑ（ｔ，θ）を取得する。
【０１０８】
ステップ１は、各ビュー内の全てのｔを処理し、一方、ステップ２から４は、一度に１つのｔだけを処理する。
【０１０９】
上述のｒＦＢＰ−ＭＡＰアルゴリズムは、平行ビームを用いる撮像ジオメトリのみに限定されない。このｒＦＢＰ−ＭＡＰアルゴリズムは、ＦＢＰアルゴリズムが存在する限り、ファンビーム、コーンビーム、平面積分、及び減衰計測等の他の撮像ジオメトリに拡張可能である。ＦＢＰアルゴリズムに対する唯一の修正は、投影データフィルタリングである。各計測に重み係数が割り当てられ、周波数領域におけるフィルタ伝達関数は、撮像ジオメトリとは無関係な（２３）式として形成される。
【０１１０】
ｒＦＢＰ−ＭＡＰアルゴリズムの効果を例示するために、低線量使用のＸ線コンピュータ断層撮影装置を使用して、死体の胴体が走査された。画像は、ＦＢＰアルゴリズム、ならびにｒＦＢＰ−ＭＡＰアルゴリズムを用いて再構成された。撮像ジオメトリはコーンビームであり、Ｘ線源軌道は半径６００ｍｍの円である。Ｘ線検出器は６４列を有し、列の高さは０．５ｍｍであり、各列は８９６チャネルを有し、ファン角度は４９．２°であった。管電圧は１２０ｋＶであり、電流は６０ｍＡであった。３６０°に亘って１２００個のビューが均一にサンプリングされた。
【０１１１】
ＦＢＰアルゴリズムとしてフェルトカンプアルゴリズムを使用し、データは、まず、余弦関数を用いて重み付けされ、次いで、１次元ランプフィルタがコーンビーム投影の各列に適用された。最後に、３次元ボリュームデータを発生するためにコーンビーム逆投影が使用された。ｒＦＢＰ−ＭＡＰアルゴリズムを実施するために、１次元ランプフィルタは、（２３）式により表されるようにランプフィルタと置換された。反復指数ｋは、２０，０００と選択され、ステップサイズαは２．０に設定された。エッジ保存ＭＡＰ再構成の場合、２つのβ値、すなわち、β＝０．０００５及びβ＝０が使用された。ベイジアンペナルティ関数ｇ（Ｘ）は、横断断面内の中央画素値とその４近傍の平均値との間の差の共通二次関数であった。すなわち、Ｒは、その１次元バージョンが［−１  ２  −１］のカーネルを用いた二次導関数である単なるラプラシアンである。ノイズ重み付け関数は、ｗ（ｔ，θ）＝ｅｘｐ（−ｐ（ｔ，θ））によって定義された。本実施形態においては、伝送計測が近似的にポワソン分布に従い、線積分ｐ（ｔ，θ）が伝送計測の対数であると仮定された。
【０１１２】
エッジマップは、２次元ソーベルカーネルを使用して断面毎に取得される。エッジ検出のための閾値は、最大画素値の７０％に設定された。次いで、３×３平滑カーネルを二値画像に適用する。平滑カーネルの適用により二値画像に含まれるエッジがボケ、二値画像から（１３）式に示すエッジマップＢが発生される。
【０１１３】
画像ボリュームデータは、５１２×５１２×５１２の３次元配列で再構成され、表示のために非中央軸断片が使用される。その他の事前フィルタリング又は事後フィルタリングはこの再構成に適用されなかった。
【０１１４】
図１３Ａは、従来のＦＢＰフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。低線量のため、胴体を横切って左から右への著しいストリークアーチファクトが発生している。図１３Ｂは、β＝０を用いたｒＦＢＰ−ＭＡＰ再構成であるｒＦＢＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。ｒＦＢＰアルゴリズムにおける投影線単位のノイズ重み付けは、図１３Ａに出現したストリークアーチファクトを効果的に除去した。図１３Ｃは、二次平滑化事前分布を使用し、β＝０．０００５を用いたｒＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。図１３Ｃに示すように、ノイズ重み付けにより、ストリークアーチファクトは消えている。ボケの事前分布により、画像は特にエッジ領域で過度に平滑に見える。図１３Ｄは、実質的に図１３Ｂ及び図１３Ｃの重み付け加算である、エッジ保存事前分布を用いたｒＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。エッジ検出ソーベルカーネルを用いて、重み付け加算画像ＢがＸ0から抽出された。エッジマップは二値画像であった。エッジ領域から非エッジ領域への遷移が平滑であるように、エッジマップをボケさせるためにローパスフィルタが使用された。図１３Ｄの画像は、図１３Ｂよりもノイズがより少ないが、図１３Ｂの主なエッジを平滑でない状態で保持している。図１３Ｅは、エッジ保存ＦＢＰ−ＭＡＰ再構成において使用されるエッジマップＢを示す図である。
【０１１５】
計算効率を改善するために、多くのアプリケーションに関して、ビュー単位のノイズ重み付けｖＦＢＰ−ＭＡＰアルゴリズムを使用することが可能である。ノイズ重み付けは、各ビューにおける最大ノイズ分散によって決定される。図１３Ｆは、ｖＦＢＰ−ＭＡＰ再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。ｖＦＢＰ−ＭＡＰ再構成は、ノイズ重み付け関数ｗを除いて、全てのパラメータが図１３Ｄにおける再構成のパラメータと同一値に設定された。
【０１１６】
本実施形態に係る画像再構成においてノイズモデルは、ポワソン統計により規定される。あるいは、ノイズモデルは、ガウス分布により記述される電子ノイズを含む光子強度計測の複合ポワソン統計により規定される。確率論において複合ポワソン分布は、独立した一様に分布したランダム変数の和の確率分布であり、変数の数は、ポワソン分布に従う。ポワソン分布は、入射光子の数の記述に適している。蓄電時間の間に検出器に達する光子数Ｎは、Ｎ〜ポワソン（λ）であり、この場合、λは計数率である。検出された各光子は、ランダムエネルギーＸiを有し、その分布はＸ線管の多色スペクトルと検出器効率とによって説明される。エネルギー積分検出器の場合、積分エネルギーは、Ｙ＝ΣNＸiによって与えられる。積分エネルギーＹは、（２５）式及び（２６）式のように複合ポワソン統計（平均値および分散）に従う。
【数４０】
2014023936.tif 000042
【０１１７】
【数４１】
2014023936.tif 000043
【０１１８】
分散利得の概念は、Ｗによって示され、（２７）式のように規定される。
【数４２】
2014023936.tif 000044
【０１１９】
ポワソン分散に従うランダム変数の場合、Ｗ＝１である点に留意されたい。複合ポワソン変数の場合、Ｗは（２８）式に従うとする。
【数４３】
2014023936.tif 000045
【０１２０】
従ってノイズ分散利得は、Ｘ線検出器への入射Ｘ線のスペクトルの平均有効エネルギーに比例する。一般に、平均スペクトルエネルギーに影響を及ぼす以下の係数が存在する。
【０１２１】
１）管電圧（ｋＶｐ）。通常、管電圧は８０ｋＶｐから１４０ｋＶｐまで様々であり、スペクトルエネルギーの上限を決定する。
【０１２２】
２）ボウタイフィルタ及びスペクトルフィルタ、ヒール効果。これらは、ＣＴデータ内に常に存在する。ボウタイフィルタ及びヒール効果によりスペクトルは、各画素において異なる。従って各検出器画素に関して、独立して計算を行う必要がある。
【０１２３】
３）オブジェクトのサイズ及び構造。オブジェクトがビームを硬化させること、すなわち、Ｅ［Ｘ］を増大することはよく知られている。ビームハードニング効果は、経路長と減衰とに依存する。
【０１２４】
複合ポワソンモデルにおいてノイズ分散利得に関する平均スペクトルエネルギーの効果が利用される。収集データは、（２９）式のように表される。
【数４４】
2014023936.tif 000046
【０１２５】
（２９）式中、ｇは、検出器利得係数であり、ｅは分散Ｖeを伴う電子ノイズである。Ｖeを計測するために、Ｘ線管のスイッチをオフにしてデータを収集して、分散を計測する。平均及び分散は、各検出器画素に関して独立して時間方向に計測される。収集データの統計は、（３０）式により得られる。
【数４５】
2014023936.tif 000047
【０１２６】
すなわち、（３１）式が得られる。
【数４６】
2014023936.tif 000048
【０１２７】
オブジェクトに関する投影データを収集した後で、平均Ｅ［Ｚ］及び分散Ｖａｒ［Ｚ］を計測して、ｇＷを計算することができる。利得ｇは投影データから容易に取得できない点に留意されたい。未知のスケールｇの効果を補償するために、空気のスキャンデータ（ＷOBJ＝Ｗ／ＷAIR）により投影データを較正する。ノイズ分散に関する複合ポワソンモデルは、（３２）式により書き直される。
【数４７】
2014023936.tif 000049
【０１２８】
（３２）式中、Ｗは（２８）式により与えられる。
【０１２９】
本発明の構造および機能の詳細と共に、本発明の多数の特徴および利点が前述の説明に記載されているものの、本開示は、単なる例示であり、特に、部品の形状、サイズ、および構成、ならびにソフトウェア、ハードウェア、またはそれら両方の組合せの形の実装形態に関して詳細な変更を行うことは可能であるが、それらの変更は、添付の請求項が表現される条件の一般的な意味によって示される完全範囲まで本発明の原理内である点を理解されたい。
【０１３０】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１３１】
１００…ガントリ、１０１…Ｘ線管、１０２…フレーム、１０３…Ｘ線検出器、１０４…データ収集回路、１０５…非接触データ伝送部、１０６…前処理部、１０７…回転部、１０８…スリップリング、１０９…高電圧発生部、１１０…システム制御部、１１２…記憶部、１１４…再構成部、１１４Ａ…ノイズ重み決定部、１１４Ｂ…フィルタ整形部、１１４Ｃ…画像発生部、１１５…入力部、１１６…表示部、１１７…電流調整部、２００…スキャン計画サポート装置
【図１】
2014023936.tif 000052
【図２】
2014023936.tif 000053
【図３】
2014023936.tif 000054
【図４】
2014023936.tif 000055
【図５】
2014023936.tif 000056
【図６Ａ】
2014023936.tif 000057
【図６Ｂ】
2014023936.tif 000058
【図７】
2014023936.tif 000059
【図８】
2014023936.tif 000060
【図９】
2014023936.tif 000061
【図１０】
2014023936.tif 000062
【図１１】
2014023936.tif 000063
【図１２】
2014023936.tif 000064
【図１３Ａ】
2014023936.tif 000065
【図１３Ｂ】
2014023936.tif 000066
【図１３Ｃ】
2014023936.tif 000067
【図１３Ｄ】
2014023936.tif 000068
【図１３Ｅ】
2014023936.tif 000069
【図１３Ｆ】
2014023936.tif 000070