特許第6301081号(P6301081)IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 東芝メディカルシステムズ株式会社の特許一覧

特許6301081X線コンピュータ断層撮影装置、画像再構成方法及び再構成フィルタの構造
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6301081
(24)【登録日】2018年3月9日
(45)【発行日】2018年3月28日
(54)【発明の名称】X線コンピュータ断層撮影装置、画像再構成方法及び再構成フィルタの構造
(51)【国際特許分類】
   A61B 6/03 20060101AFI20180319BHJP
【FI】
   A61B6/03 350S
【請求項の数】28
【全頁数】29
(21)【出願番号】特願2013-158208(P2013-158208)
(22)【出願日】2013年7月30日
(65)【公開番号】特開2014-23936(P2014-23936A)
(43)【公開日】2014年2月6日
【審査請求日】2016年7月11日
(31)【優先権主張番号】13/561,674
(32)【優先日】2012年7月30日
(33)【優先権主張国】US
(73)【特許権者】
【識別番号】594164542
【氏名又は名称】キヤノンメディカルシステムズ株式会社
(74)【代理人】
【識別番号】100108855
【弁理士】
【氏名又は名称】蔵田 昌俊
(74)【代理人】
【識別番号】100109830
【弁理士】
【氏名又は名称】福原 淑弘
(74)【代理人】
【識別番号】100088683
【弁理士】
【氏名又は名称】中村 誠
(74)【代理人】
【識別番号】100103034
【弁理士】
【氏名又は名称】野河 信久
(74)【代理人】
【識別番号】100075672
【弁理士】
【氏名又は名称】峰 隆司
(74)【代理人】
【識別番号】100153051
【弁理士】
【氏名又は名称】河野 直樹
(74)【代理人】
【識別番号】100140176
【弁理士】
【氏名又は名称】砂川 克
(74)【代理人】
【識別番号】100158805
【弁理士】
【氏名又は名称】井関 守三
(74)【代理人】
【識別番号】100172580
【弁理士】
【氏名又は名称】赤穂 隆雄
(74)【代理人】
【識別番号】100179062
【弁理士】
【氏名又は名称】井上 正
(74)【代理人】
【識別番号】100124394
【弁理士】
【氏名又は名称】佐藤 立志
(74)【代理人】
【識別番号】100112807
【弁理士】
【氏名又は名称】岡田 貴志
(74)【代理人】
【識別番号】100111073
【弁理士】
【氏名又は名称】堀内 美保子
(72)【発明者】
【氏名】レリー・ザン
(72)【発明者】
【氏名】アレキサンダー・ザミャチン
【審査官】 伊藤 昭治
(56)【参考文献】
【文献】 特表2008−505676(JP,A)
【文献】 米国特許出願公開第2012/0045108(US,A1)
【文献】 米国特許出願公開第2013/0251286(US,A1)
【文献】 藤田晃史,異なる再構成関数および新しい再構成アルゴリズムを用いた超低線量胸部MDCT画像の検討,日本医学放射線学会雑誌,2003年11月25日,vol. 63, No.9 ,588〜590頁
(58)【調査した分野】(Int.Cl.,DB名)
A61B 6/03
(57)【特許請求の範囲】
【請求項1】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記決定部は、ビュー毎に前記投影データに基づいて前記重み値を決定する、X線コンピュータ断層撮影装置。
【請求項2】
前記再構成フィルタは、下記数式のH(w,view)により規定され、
【数1】
ここで、ωは周波数分散であり、w(view)はビュー単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、kは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
請求項記載のX線コンピュータ断層撮影装置。
【請求項3】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記決定部は、投影線毎に前記投影データに基づいて前記重み値を決定する、線コンピュータ断層撮影装置。
【請求項4】
前記再構成フィルタは、下記数式のH(w,ray)により規定され、
【数2】
ここで、ωは周波数分散であり、w(ray)は投影線単位で決定される重み値であり、αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、kは反復アルゴリズムにおいて反復数をエミュレートするパラメータである、
請求項記載のX線コンピュータ断層撮影装置。
【請求項5】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記発生部は、逆投影が実行された後に、前記再構成フィルタを前記投影データに適用する、線コンピュータ断層撮影装置。
【請求項6】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記ノイズモデルは、ポワソン統計を含む、線コンピュータ断層撮影装置。
【請求項7】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記ノイズモデルは、ガウス分布に従う電子ノイズを含む光子強度計測の複合ポワソン統計によって記述される、線コンピュータ断層撮影装置。
【請求項8】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
前記再構成フィルタを整形するために正規化する正規化部と、
を具備するX線コンピュータ断層撮影装置。
【請求項9】
前記画像が基準画像に基づいてさらに発生される、請求項記載のX線コンピュータ断層撮影装置。
【請求項10】
前記整形部は、前記再構成フィルタとして第1の再構成フィルタと第2の再構成フィルタとを整形し、
前記第1の再構成フィルタは、エッジを保存するための第1の画像を発生するように整形され、第2の再構成フィルタは、ノイズを低減してエッジマップを発生するための第2の画像を発生するように整形され、
前記画像は、前記第1の画像と前記第2の画像と前記エッジマップとに基づいて発生される、
請求項記載のX線コンピュータ断層撮影装置。
【請求項11】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記整形部は、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするように前記再構成フィルタを整形する、線コンピュータ断層撮影装置。
【請求項12】
少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、
前記再構成フィルタに基づいて画像を発生する発生部と、
を具備し、
前記整形部は、前記再構成フィルタをシステムマトリクスに基づいて整形する、線コンピュータ断層撮影装置。
【請求項13】
前記システムマトリクスが前記投影データの不均一サンプリングを指定する、請求項12記載のX線コンピュータ断層撮影装置。
【請求項14】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記決定することは、ビュー毎に前記投影データに基づいて前記重み値を決定する、画像再構成方法。
【請求項15】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記決定することは、投影線毎に前記投影データに基づいて前記重み値を決定する、画像再構成方法。
【請求項16】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記発生することは、逆投影が実行された後に、前記再構成フィルタを前記投影データに適用する、画像再構成方法。
【請求項17】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記ノイズモデルは、ポワソン統計を含む、画像再構成方法。
【請求項18】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記ノイズモデルは、ガウス分布に従う電子ノイズを含む光子強度計測の複合ポワソン統計によって記述される、画像再構成方法。
【請求項19】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生し、
前記再構成フィルタを整形するために正規化する
ことを具備する画像再構成方法。
【請求項20】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記整形することは、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするように前記再構成フィルタを整形する、画像再構成方法。
【請求項21】
少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定し、
前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形し、
前記再構成フィルタに基づいて画像を発生する、
ことを具備し、
前記整形することは、前記再構成フィルタをシステムマトリクスに基づいて整形する、画像再構成方法。
【請求項22】
投影データを重み付けするために少なくとも所定のノイズモデルに従ってビュー毎に決定された重み値を、前記投影データに基づいて決定し、前記決定されたビュー毎の重み値を含むパラメータに基づいて整形された、少なくとも1つの再構成フィルタの構造であって、
X線コンピュータ断層撮影装置の再構成処理に用いられる、再構成フィルタの構造。
【請求項23】
前記再構成フィルタの構造は、前記重み値を含む前記パラメータを各ビューに対応させて配列する請求項22記載の再構成フィルタの構造。
【請求項24】
投影データを重み付けするために少なくとも所定のノイズモデルに従って投影線毎に決定された重み値を、前記投影データに基づいて決定し、前記決定された投影線毎の重み値を含むパラメータに基づいて整形された、少なくとも1つの再構成フィルタの構造であって、
X線コンピュータ断層撮影装置の再構成処理に用いられる、再構成フィルタの構造。
【請求項25】
前記再構成フィルタの構造は、前記重み値を含む前記パラメータを各投影線に対応させて配列する請求項24記載の再構成フィルタの構造。
【請求項26】
投影データを重み付けするために少なくともポワソン統計を含むノイズモデルに従って決定された重み値を、前記投影データに基づいて決定し、前記決定された重み値を含むパラメータに基づいて整形された、少なくとも1つの再構成フィルタの構造であって、
X線コンピュータ断層撮影装置の再構成処理に用いられる、再構成フィルタの構造。
【請求項27】
投影データを重み付けするために少なくともガウス分布に従う電子ノイズを含む光子強度計測の複合ポワソン統計によって記述されるノイズモデルに従って決定された重み値を、前記投影データに基づいて決定し、前記決定された重み値を含むパラメータに基づいて整形された、少なくとも1つの再構成フィルタの構造であって、
X線コンピュータ断層撮影装置の再構成処理に用いられる、再構成フィルタの構造。
【請求項28】
投影データを重み付けするために少なくとも所定のノイズモデルに従って決定された重み値を、前記投影データに基づいて決定し、前記決定された重み値を含むパラメータに基づき、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするように整形された、少なくとも1つの再構成フィルタの構造であって、
X線コンピュータ断層撮影装置の再構成処理に用いられる、再構成フィルタの構造。
【発明の詳細な説明】
【技術分野】
【0001】
本発明の実施形態は、X線コンピュータ断層撮影装置及び画像再構成方法に関する。
【背景技術】
【0002】
フィルタ補正逆投影(FBP)アルゴリズムは、簡単かつ効率的であり、核医学、X線CT、及びMRIにおいて画像を再構成するために使用される。FBPアルゴリズムは、計算効率が高いため、X線CT及び核医学の画像再構成において重宝されている。FBPは、低線量のX線で収集されたデータに適用されると多くのノイズを含む画像を発生してしまう。
【0003】
一般的な傾向として、FBPアルゴリズムは、反復アルゴリズムに徐々に置換されている。FBPアルゴリズムは、ノイズを含む画像を発生してしまう。さらに、FBPアルゴリズムは、ノイズレベルを低減するためにノイズモデルを組み込むことが可能でないことが知られている。この点で、反復アルゴリズムは、投影ノイズモデルを組み込むことにより、FBPアルゴリズムよりもノイズがより少ない画像を発生することができる。しかしながら、反復アルゴリズムは、集中的な計算を必要とする。
【0004】
計算要件にかかわらず、反復アルゴリズムは、最大事後確率(MAP)を使用してノイズと分解能とのバランスがとれた画像を発生することができる。対照的に、FBPアルゴリズムは、MAPまたは事前の画像情報を活用することができないことが知られている。
【発明の概要】
【発明が解決しようとする課題】
【0005】
目的は、FBPアルゴリズムを使用し、計算効率を失わずに画像内のノイズを実質的に低減可能なX線コンピュータ断層撮影装置及び画像再構成方法を提供することにある。
【課題を解決するための手段】
【0006】
本実施形態に係るX線コンピュータ断層撮影装置は、少なくとも所定のノイズモデルに従って、投影データに重み付けするための重み値を決定する決定部と、前記重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する整形部と、前記再構成フィルタに基づいて画像を発生する発生部と、を具備する。
【図面の簡単な説明】
【0007】
図1】本実施形態に係るX線コンピュータ断層撮影装置の構成を示す図。
図2図1の再構成部の構成を示す図。
図3図2の再構成部による処理の典型的な流れを示す図。
図4】本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図。
図5】本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図。
図6A】本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図。
図6B】本実施形態に係る、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図。
図7】本実施形態に係る、FBPアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図。
図8】本実施形態に係る、エッジ保存FBP−MAPアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図。
図9図8のエッジマップの発生処理の典型的な流れを示す図。
図10】本実施形態に係る、事前画像を使用したFBP−MAPアルゴリズムに基づく画像再構成方法の典型的な流れを示す図。
図11】本実施形態に係るノイズ重み付け反復エミュレーションFBPアルゴリズムを用いて画像再構成方法の典型的な流れを示す図。
図12】本実施形態に係るノイズ重み付けFBPアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図。
図13A】従来のFBPフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
図13B】β=0を用いたrFBP−MAP再構成であるrFBP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
図13C】二次平滑化事前分布を使用し、β=0.0005を用いたrFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
図13D】実質的に図13B及び図13Cの重み付け加算である、エッジ保存事前分布を用いたrFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
図13E】エッジ保存FBP−MAP再構成において使用されるエッジマップを示す図。
図13F】vFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図。
【発明を実施するための形態】
【0008】
本実施形態は、ノイズ重み付けフィルタ補正逆投影(filtered back projection)法を使用して投影データを処理するためのX線コンピュータ断層撮影装置及び再構成処理方法に関する。
【0009】
以下、図面を参照しながら、本実施形態に係わるX線コンピュータ断層撮影装置及び画像再構成方法について説明する。
【0010】
図1は、本実施形態に係るX線コンピュータ断層撮影装置の構成を示す図である。図1に示すように、本実施形態に係るX線コンピュータ断層撮影装置は、マルチスライスX線CT装置またはマルチスライスX線CTスキャナである。本実施形態に係るX線コンピュータ断層撮影装置は、ガントリ100を装備している。ガントリ100は、X線管101、円環状フレーム102、及び多列タイプあるいは二次元配列タイプのX線検出器103を含む。X線管101とX線検出器103とは、軸RAの周りを回転するフレーム102上に被検体Sを横断して直径方向に設置される。回転部107は、被検体Sが軸RAに沿って例示されたページ内または例示されたページ外に移動される間に、一回転当たり0.4秒等の高速でフレーム102を回転させる。
【0011】
本実施形態に係るX線コンピュータ断層撮影装置は、高電圧発生装置109を備える。高電圧発生装置109は、X線管101がX線を発生するように、管電圧をX線管101に印加する。例えば、高電圧発生器109はフレーム102上に設置される。X線は、被検体Sに向けて放射される。X線検出器103は、被検体Sを透過したX線を検出するために、被検体Sを横断してX線管101の逆側に位置する。
【0012】
さらに図1に示すように、本実施形態に係るX線コンピュータ断層撮影装置は、X線検出器103から検出された信号を処理するための他のデバイスをさらに含む。データ収集回路(DAS)104は、X線検出器103からの信号出力をチャンネル毎に電圧信号に変換して、当該電圧信号を増幅し、増幅後の電圧信号をデジタル信号に変換する。X線検出器103とDAS104とは、1回転ごとに所定の総投影数(TPRP)を処理するように構成される。
【0013】
上述のデータは、非接触データ伝送部105を介してガントリ100の外部のコンソール内に収納された前処理部106に伝送される。前処理部106は、生データに対して感度補正等の前処理を施す。前処理が施されたデータは、投影データと呼ばれている。記憶部112は、投影データを記憶する。記憶部112は、再構成部114、表示部116、入力部115、及びスキャン計画サポート装置200と共に、データ/制御バスを介してシステム制御部110に接続される。スキャン計画サポート装置200は、スキャン計画を構築するために、撮影技法をサポートするための機能を有している。
【0014】
本実施形態によれば、再構成部114は、ノイズ重み付けを用いたフィルタ補正逆投影(FBP)技法に基づいて、記憶部112に記憶された投影データから画像を再構成する。あるいは、再構成部114は、ノイズ重み付けと基準画像等の事前分布入力(prior in)とも用いたフィルタ補正逆投影に基づいて、投影データから画像を再構成しても良い。再構成部114の上記の2つの実施形態のいずれか1つは、所定の反復再構成アルゴリズムに従って所定数の反復において特定の反復結果をエミュレートする追加の特徴を用いたフィルタ補正逆投影法に基づいて、投影データから画像を再構成する。
【0015】
再構成部114は、ソフトウェアとハードウェアとの組合せにより構成され、特定の実装形態に限定されない。本実施形態に係る再構成部114は、核医学及び磁気共鳴撮像(MRI)等の他のモダリティにも実装可能である。
【0016】
図2は、本実施形態に係る再構成部114の構成を示す図である。図2に示すように再構成114は、ノイズ重み決定部114A、フィルタ整形部114B、及び画像発生部11Cを有する。ノイズ重み決定部114Aは、少なくとも所定のノイズモデルに従って投影データに重み付けを施すための重み値を決定する。フィルタ整形部114Bは、決定された重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する。画像発生部114Cは、整形された再構成フィルタに基づいて画像を発生する。なお、この再構成部114の実装形態は単なる例示であることを言及しておく。以下、ノイズ重み決定部114A、フィルタ整形部114B、及び画像発生部11Cについて詳細に説明する。
【0017】
図3は、再構成部114による処理の典型的な流れを示す図である。ステップS100においてノイズ重み決定部114Aは、少なくとも所定のノイズモデルに従って投影データに重み付けするための重み値を決定する。重み値は、投影データについての所定のデータ単位毎に決定される。データ単位としては、ビューであっても良いし、投影線であっても良い。ステップS120においてフィルタ整形部114Bは、ステップS100で決定された重み値を含むパラメータに基づいて、少なくとも1つの再構成フィルタを整形する。本実施形態に係る整形のために、複数の再構成フィルタがオプションで使用されると同時に、これらのパラメータは重み値に限定されない。ステップS120において画像発生部114Cは、ステップS120において整形された再構成フィルタに基づいて画像を発生する。
【0018】
ステップS120は、所定のパラメータに基づいてフィルタを整形する。上記のように、1つの例示的なパラメータは、投影データに含まれるノイズ成分の分散(以下、ノイズ分散と呼ぶ)である。しかしながら、本実施形態に係る重み値を含むパラメータは、ノイズパラメータに限定されない。例えば、再構成フィルタは、システムマトリクスに基づいて整形されても良い。より詳細には、システムマトリクスは、投影データの不均一サンプリングを指定する。従来のFBPアルゴリズムはオブジェクトが均一にサンプリングされると仮定するが、可変のサンプリング戦略を使用することも可能である。例えば、重要な構造を含むある角度範囲内の信号は、第1の角度間隔でサンプリングされる。対照的に、当該領域外の別の角度範囲内の信号は、第1の角度間隔よりも大きい第2の角度間隔でサンプリングされる。
【0019】
なおステップS100は、投影データの所定のデータ部分の各々に関する重み値を決定する際に、ステップS120の前に実行または完了されなくても良い。すなわち、再構成フィルタが整形される度にオプションで重み値が決定されても良い。
【0020】
さらに図3に示すように、ステップS140において画像発生部114Cは、ステップS120において整形された再構成フィルタを適用しながら、サブステップの所定のセットを実行する。1つの例示的なプロセスにおいては、さらに詳細に説明されるように、所定のFBP技法においてフーリエ変換の後に、整形されたフィルタを周波数領域内で投影データに適用することによって画像が再構成される。別の例示的なプロセスにおいては、所定のFBP技法において畳み込みとしてカスタマイズされたフィルタを空間領域内で投影データに適用することにより画像が再構成される。
【0021】
図4は、本実施形態に係る画像再構成方法におけるビュー毎のノイズ重み値の決定方式を示す図である。ビュー毎の重み付け方式において各ビューは、概して、特定のビューにおける最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、重み値が割り当てられる。図4は、3つの位置、すなわち、位置A、B、及びCにおけるX線源、ならびに位置A’、B’、及びC’における3つの対応する位置におけるX線検出器を例示する。例えば、X線源が位置Aにあるとき、X線検出器は、対応する位置A’にある。位置A’においてX線検出器により収集された投影データはビューからなり、重み値w(θ1)は、当該ビューに関する所定の重み付け方式に従ってビュー角度θ1に応じて決定される。上記のように、ビュー単位の重み値の決定は、1つの例示的なプロセスでは、再構成フィルタの整形に先立って実行および完了される。このビュー単位の重み値の決定は、他の例示的なプロセスでは、再構成フィルタが整形される度に決定される。
【0022】
本実施形態において投影データに含まれるノイズ分散は、ビュー単位の重み付け方式を使用してモデル形成可能である。一般に、ノイズの影響を補償する間に画像を再構成するための典型的な手法は、(1)式により規定される、ノイズ分散に応じて重み付けされた目的関数を最小化することである。
【数1】
【0023】
(1)式中、Aは投影行列であり、Xは列ベクトルとして記述される画像配列であり、Pは列ベクトルとして記述される投影配列であり、Wは重み付けノルムを定義するノイズ重み付け行列である。この目的関数を最小化するために、(2)式に規定される反復ランドウェーバー(Landweber)アルゴリズムが使用される。
【数2】
【0024】
(2)式中、ATは逆投影行列であり、X(k)は第k番目の反復において推定された画像であり、αはステップサイズであり、α>0である。この再帰的表現は、(3)式のように非再帰的表現に書き直すことができる。
【数3】
【0025】
(ATWA)-1=A-1-1(AT-1が存在すると仮定すると、(3)式は、(4)式のように簡素化される。
【数4】
【0026】
(4)式は、ノイズ重み付けFBPアルゴリズムに相当し、周波数領域内における(4)式のフィルタ関数は、(5)式のように表される。
【数5】
【0027】
(5)式中、ωは周波数分散であり、w(view)はビュー単位で決定される重み値である。すなわち、w(view)は、特定のビュー角度におけるノイズ関連の重み関数である。αは反復アルゴリズムにおいてステップサイズをエミュレートするパラメータであり、kは反復アルゴリズムにおいて反復数をエミュレートするパラメータである。vFBPとして示されるビュー毎の重み付けFBPアルゴリズムは、シフト不変の点応答関数(PRF)を有する。
【0028】
上記の実施形態において再構成フィルタは、所定の反復再構成アルゴリズムに従って、所定数k回の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けFBPアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
【0029】
図5は、本実施形態に係る画像再構成方法における投影線毎のノイズ重み値の決定方式を示す図である。投影線毎の重み付け方式では、特定のビュー内の最大ノイズ分散または平均ノイズ分散など、所定の関数に従って、各投影線に重み値が割り当てられる。図5は、所定の位置におけるX線源、ならびに、対応する位置におけるX線検出器を例示する。投影線がX線源から投影されるとき、投影線は対応する位置におけるX線検出器に達する。投影データは、θにおけるビュー角度内の複数の投影線からなる。例えば、重み値w(θ,t1)は、ビュー角度θに応じ、かつ当該投影線に関する所定の重み付け方式に従って、特定の投影線t1について決定される。上記の通り、投影線の重み値の決定は、1実施形態において、再構成フィルタを整形するのに先立って実行および完了される。投影線重み値は、他の実施形態において、再構成フィルタとして決定される。
【0030】
シフト不変PRFに関する要件が除去されると、投影線毎の重み付け方式が構築され、その周波数領域における窓化されたランプフィルタが上記の(5)式を変形することにより、(6)式のように記述される。
【数6】
【0031】
(6)式中、重み値、すなわち、重み係数w(ray)は、投影線で決定される重み計数である。w(ray)は、関連する投影線のノイズ分散に応じて決定される。kは、所定の反復再構成アルゴリズムに従って、特定の反復結果をエミュレートするための所定数の反復を示すパラメータである。H(ω,ray)の逆フーリエ変換がh(t,ray)であるとすると、h(t,ray)はフィルタの空間領域におけるカーネルである。p(r,θ)をビューθ及び投影線rにおける投影であるとし、q(r,θ)をフィルタリング後の投影であるとする。この場合、q(r,θ)は、以下の(7)式のような積分で規定される。
【数7】
【0032】
カーネルh(t,r)はrに依存するため、この積分は畳み込み(コンボルーション)ではない。フィルタリングが(7)式のように空間領域で実施される場合、計算コストは畳み込みと同等である。この新しいrFBPアルゴリズムにおける逆投影ステップは、従来のFBPアルゴリズムの逆投影ステップと同じである。従ってrFBPアルゴリズムは、計算的に効率的である。他方で、画像領域フィルタリングは、投影線毎のノイズ重み付けFBP−MAP(rFBP−MAP)アルゴリズムに関する実装形態で後に説明されるように、伝達関数を用いて2次元フーリエ領域内で実施される。
【0033】
重み係数を割り当てるための典型的な手法は、w(ray)を投影線計測のノイズ分散の逆数にすることである。この手法は、最適化問題に関する目的関数として、尤度関数(すなわち、結合確率密度関数)を使用して補われる。ノイズを含む計測に比してノイズを含まない計測を信頼すべきであるという根本原理を採用する場合、ノイズがより少ない計測により大きなw(ray)が割り当てられ、ノイズがよい多い計測により小さなw(ray)が割り当てられという制限下において任意に重み係数を割り当てることができる。
【0034】
本実施形態に係るrFBPアルゴリズムを実施するための最も容易かつより効率的な様式は、空間領域において投影をフィルタリングするために、(7)式を使用することである。現在、本発明者らは積分カーネルh(t,r)に関する解析式を有していない。各ビューθに関してq(r,θ)を計算するための代替様式は、周波数領域において、以下のように方程式(7)を実施することである。
【0035】
ステップ1:rに関するp(r,θ)の1次元フーリエ変換を見出して、P(ω,θ)を取得する。
【0036】
ステップ2.a:各投影線「ray」に関して、重み係数w(ray)を割り当て、方程式(6)で表されるように周波数領域伝達関数を形成する。実装形態において、ωは周波数指数であり、整数をとる。
【0037】
ステップ2.b:積、Qray(ω,θ)=P(ω,θ)H(ω,ray)を見出す。
【0038】
ステップ2.c:ωに関するQray(ω,θ)の1次元逆フーリエ変換を見出して、q(ray,θ)を取得する。
【0039】
上記の実施形態においてステップ1は、各ビューの列内の全てのrを処理し、一方、ステップ2.aから2.cは一度に1つの投影線のみを処理する。ビュー毎の重み付けFBP(vFBP)アルゴリズムでは、単一の重み値、すなわち、重み係数w(view)がビュー内の全ての投影線に割り当てられる。対照的に、投影線毎の重み付けFBP(rFBP)アルゴリズムでは、別個の重み値、すなわち、重み係数w(ray)がビュー内の投影線の各々に割り当てられる。結果として、ノイズ重み付け方式は、ビュー毎の重み付けFBP(vFBP)アルゴリズムよりも、投影線毎の重み付けFBP(rFBP)アルゴリズムにおいてより正確である。
【0040】
ビュー毎の重み付けFBP(vFBP)アルゴリズム及び投影線毎の重み付けFBP(rFBP)アルゴリズムに関して上で説明されたように、本実施形態は、画像を発生する前に、重み値を含むパラメータに基づいて少なくとも1つの再構成フィルタを整形する。さらに、上記の実施形態では、再構成フィルタは所定の反復再構成アルゴリズムに従って、kで記された所定数の反復において特定の反復結果をエミュレートするように整形される。再構成フィルタにおいて反復をエミュレートするための上記の特徴は、ビュー毎の重み付けFBPアルゴリズムまたは投影線ごとの重み付けFBPアルゴリズムを使用して投影データから画像を再構成するための他の実施形態に含まれるとは限らない。
【0041】
その後、整形された再構成フィルタに基づいて、画像を発生するための工程が実行される。一般に、画像再構成におけるノイズ低減は、画像再構成、特に逆投影が加算処理であり、投影データの適切な重み付けが実質的に逆投影の分散を削減するという概念に基づいている。
【0042】
図6Aは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む画像再構成方法の典型的な流れを示す図である。FBPアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行される前に、再構成フィルタが投影データに適用される。ステップS140Aに先立って再構成フィルタは、ノイズ重み値を含むパラメータに基づいて予め整形されている。ステップS140Aにおいて整形された再構成フィルタが投影データに適用される。ステップS140Aの前に投影データがフーリエ変換に供される場合、フィルタリングされた投影データは、ステップS140Bの前に逆フーリエ変換を受ける。ステップS140Bにおいてフィルタリングされた投影データは、逆投影される。逆投影により、フィルタリングされた投影データに基づいて画像が発生される。
【0043】
図6Bは、予め整形された再構成フィルタに基づいて画像を発生するための追加の工程を含む他の画像再構成方法の典型的な流れを示す図である。FBPアルゴリズムが投影データに関するタスクを実行すると仮定すると、逆投影が実行された後に、再構成フィルタが投影データに適用される。ステップS140Bの前に投影データがフーリエ変換に供される場合、投影データは、ステップS140Bの前に逆フーリエ変換を受ける。ステップS140Bにおいて投影データは、逆投影され、画像が発生される。整形された再構成フィルタは、逆投影されたデータに適用される。ステップS140Aに先立って、再構成フィルタはノイズ重み値を含むパラメータに基づいて予め整形されている。ステップS140Aにおいて、整形された再構成フィルタが投影データに適用される。
【0044】
図7は、FBPアルゴリズムでノイズ重み付けランプフィルタを使用した画像再構成方法の典型的な流れを示す図である。図7の処理においては、逆投影の前段においてノイズ重み付けフィルタリングがフーリエ領域で行わると仮定する。図7に示すように、投影データPDは、複数のビューにおいて収集される。投影データPDは、ビュー又は投影線単位で処理される。投影データPDは、変換されたデータがフーリエ領域または周波数領域に存在するように、フーリエ変換FTに施される。ランプフィルタは、ノイズ重み値を含む所定パラメータのセットに従って整形される。これにより、重み付けランプフィルタWFが形成される。いくつかの実施形態では、複数のランプフィルタは、複数の重み付けランプフィルタWFになるように整形される。これらパラメータは、ノイズ重みパラメータに限定されず、他のパラメータを含んでも良い。この例示的な実施形態では、変換後の投影データは、所定のFBPアルゴリズムにおいてノイズ重み付けランプフィルタNWFに従ってフィルタリングされる。周波数領域内でフィルタリングされたデータは、投影データが逆投影BPにおいて空間領域内で逆投影されるように、逆フーリエ変換IFTを受ける。あるいは、周波数領域内でノイズ重み付けランプフィルタNWFを変換後のデータに適用する代わりに、別のノイズ重み付けランプフィルタNWF2が画像に適用されても良い。この場合、画像に適用されるノイズ重み付けランプフィルタNWF2は、変換後のデータに適用されるノイズ重み付けランプフィルタNWFとは異なる。
【0045】
逆投影BPは、画像発生の間、重み関数を使用して上述の可変角度サンプリングを正規化する。例えば、重要な構造を含む一定の角度範囲内の信号は、第1の角度間隔でサンプリングされる。当該領域外の他の角度範囲内の信号は、第1の角度間隔よりもより大きい第2の角度間隔でサンプリングされる。逆投影BPは、サンプリング角度全体に亘る積分であり、角度サンプリング密度補償は、逆投影積分において本質的にヤコビアン係数である正規化係数を使用することによって達成される。数学的に、逆投影画像は、(8)式のように表現される。
【数8】
【0046】
(8)式中、初めに投影が逆投影される場合、pは角度θにおける投影である。あるいは、FBPアルゴリズムが使用される場合、pは、角度θにおいてフィルタリングされた投影であり、tは、検出器ビンの位置を示す。方程式(8)の離散的な表現は、(9)式のように表される。
【数9】
【0047】
(9)式中、nは検出器の位置指数であり、mは投影角度指数であり、Mは投影が収集されたビューの総数であり、「INT」は最近傍補間を示すために使用される。実際には、「INT」関数は使用せず、2つの近傍検出器ビンの間の線形補間が使用される。
【0048】
均一でない角度サンプリングに関して、(9)式で表されるような簡易的な逆投影が施される場合、当該逆投影は、単位角度当たりのビュー数である密度関数d(0)に依存する。例えば、0<θ<π/4に関して、サンプリング間隔が1°であり、π/4<θ<πに関して3°である場合、密度関数は、(10)式により表現され、逆投影(9)式は、(11)式のように変更される。
【数10】
【0049】
【数11】
【0050】
この場合、サンプリング密度関数は、実際の角度θの代わりに、角度指数inの関数である。
【0051】
図8は、本実施形態に係る、所定のエッジ保存FBP−MAPアルゴリズムでノイズ重み付けランプフィルタの組を使用した画像再構成方法の典型的な流れを示す図である。図8の画像再構成方法は、投影データの組を逆投影して画像を再構成する前段に、フーリエ領域においてノイズ重み付けフィルタリングが実施される。最終的に、2つの画像に基づいて最終画像が発生される。
【0052】
図8に示すように、投影データPDは、複数のビューにおいて収集される。投影データPDは、ビュー単位で処理されても良いし、投影線単位で処理されても良い。投影データPDは、変換後のデータがフーリエ領域または周波数領域にあるようにフーリエ変換FTを受ける。二つのランプフィルタは、ノイズ重み値を含む所定のパラメータのセットに従って個別に整形され、重み付けランプフィルタWF−β0及びWF−βに整形される。ランプフィルタは、正規化によってさらに整形可能である。このようにフィルタは、値βの関数である。値βが大きいほど、画像の平滑化がより強く行われる。例えば、平滑化が行われない場合あるいはボケの無い場合、β値はゼロなるように選択される。変換後の投影データは、所定のFBPアルゴリズムにおいてノイズ重み付けされると良い。ノイズ重みづけされた投影データは、正規化されたランプフィルタWF−β0及びWF−βによってそれぞれフィルタリングされる。その後、周波数領域においてフィルタリングされた投影データは、逆フーリエ変換IFTに施される。逆フーリエ変換後の投影データは、2つの画像X0及びXβが発生されるように空間領域内でそれぞれ逆投影される。平滑化を用いないで逆投影された、あるいは、最小限の平滑化を用いて逆投影された画像X0は、予め用意されたエッジマップEMAP内のエッジマップ情報Bを抽出するために使用される。最小限に平滑化された画像X0、フィルタリングされた画像Xβ、及び抽出されたエッジ情報に基づいて一定の規則に従って最終画像FIMGが発生される。
【0053】
図8に例示される実施形態においてエッジは、下記の(12)式に示す目的関数においてベイジアン正規化項を使用してノイズ除去の間に保存される。(12)式のペナルティ関数g(X)は、注目画素近傍における画素値の段差を計測する。
【数12】
【0054】
(12)式中、Aは投影行列であり、Xは列ベクトルとして書き込まれた画像配列であり、Pは列ベクトルとして書き込まれた投影配列であり、Wは各投影に関する重み関数からなる対角行列であり、bは、(12)式の第1項である忠実度項
【数13】
【0055】
に対するベイジアン項g(X)の重要性を調整する相対的な重み係数である。
【0056】
ペナルティg(X)は、画素値段差の二次関数である。当該画素値段差は、当該画素の画素値とその近傍の画素の平均値との間の差として計算される。大きな画素値段差は大きなペナルティに対応する。当該ペナルティ関数は平滑化画像をもたらす。
【0057】
エッジを過度に平滑化せずに画像ノイズを抑えるために、ペナルティ関数g(X)は、エッジ及び非エッジ領域に関して、異なる特性を有すべきである。そのような能力を有するペナルティ関数として、閾値を有するフーバ関数が用いられる。画素値段差が閾値より小さく、画素値段差がノイズとして分類される場合、フーバ関数は二次式であり、平滑化を施す。画素値段差が閾値より大きく、画素値段差がエッジとして分類される場合、フーバ関数は線形であり、より小さい平滑化を施す。
【0058】
部分毎の一定画像を奨励する別の一般的な事前分布は、TV(全変動)計測である。離散的な一次元関数g(X)の場合、g(X)のTV計測の勾配は3つの値だけをとる。すなわち、g(X)がその左右の近傍より大きい場合、正の値、g(X)がその左右の近傍よりも小さい場合、負の値、g(X)がその左右の近傍の間である場合、0である。勾配タイプの最適化アルゴリズムでは、正の勾配値は、g(X)の値を下方に押さえ込み、負の勾配値はg(X)の値を上方に押し上げ、0の勾配値は何の変更も行われていないことを意味する。従ってTV事前分布は、単調関数を奨励し、振幅を抑え、鮮明なエッジを保存する。
【0059】
図8に例示されるように、実施形態に係る画像再構成方法は、投影データにエッジ保存FBP−MAPアルゴリズムを適用する。βは、1つのノイズ重み付けランプフィルタについてのβ0、すなわち、ゼロに設定される。(12)式を参照すると、β0g(X)はゼロであり、何の平滑化も施されない。その後、何の平滑化も施されずに画像X0が再構成される。他方、他のパラメータセットの場合、βは、他のノイズ重み付けランプフィルタに関するβ、すなわり、非ゼロに設定される。βにより一定量の平滑化が施され、画像Xβが再構成される。これにより、2つの画像Xβ及びX0が再構成される。すなわち、第1の画像Xβは、非ゼロβ値を用いたrFBP−MAP再構成に基づいて再構成され、第2の画像X0はゼロβ値を用いたrFBP再構成に基づいて再構成される。
【0060】
図9は、図8のエッジマップBの発生処理の典型的な流れを示す図である。画像X0はゼロβ値を用いたrFBP再構成法により再構成される。ステップS200において、本質的に、画像勾配のノルムであるソーベルカーネルを使用して、エッジ検出によって画素値段差が評価される。ソーベルカーネルによって生成されたエッジマップは、次いで、ステップS220において、所定の閾値に従って、値0または値1を用いて二値画像に変換される。閾値を超える画像値は、値1に設定され、閾値未満の画像値は値0に設定される。その後、ステップS240において、所定のローパスフィルタによって二値画像をボケさせ、エッジマップBを発生させる。二値画像をボケさせる理由は、二値画像に含まれるエッジ領域と非エッジ領域との境界を滑らかにするためである。最終画像は、下記の(13)式に示すような、XbとX0との一次結合により発生される。
【数14】
【0061】
(14)式中、「.*」は、MATLAB(登録商標)で定義される点毎の乗算を表す。
【0062】
図10は、事前画像を使用したFBP−MAPアルゴリズムに基づく画像再構成方法の典型的な流れを示す図である。一般に、所与の事前画像
【数15】
【0063】
を加えることによって、FBP−MAPアルゴリズムに基づいて画像を再構成するために、目的関数(8)は、下記の(14)式により一般化される。
【数16】
【0064】
反復的な中間解は、(15)式のように表される。
【数17】
【0065】
等価FBP−MAPアルゴリズムは、下記の(16)式と(17)式とに表される2つのフーリエ領域フィルタ関数を有する。
【数18】
【0066】
【数19】
【0067】
図10に示すように、FBP−MAPアルゴリズムは、投影データPDと事前画像PIとの2つの入力を有する。最終画像FIは、投影データPDに由来する第1の画像と、事前画像PIに由来する第2の画像とに基づく加算画像である。FBP−MAP再構成は、(16)式のような検出器フーリエ領域フィルタとその後に行われる逆投影とに基づいて第1の再構成画像を発生する。画像フーリエ領域処理IFDは、(17)式のように画像フーリエ領域フィルタを使用して第2の画像を発生する。画像フーリエ領域フィルタは、2次元であっても良いし3次元であっても良い。
【0068】
図11は、本実施形態に係るノイズ重み付け反復エミュレーションFBPアルゴリズムを用いて画像再構成方法の典型的な流れを示す図である。一般に、再構成フィルタは、所定数の反復において特定の反復結果を模倣(エミュレート)する。再構成フィルタは、所定の反復再構成アルゴリズムに従って、重み値を含むパラメータに基づいて整形される。ノイズ重み付け反復エミュレーションFBPアルゴリズム200Dは、ステップ201において投影データが収集される。投影データは、撮影されているオブジェクトの線積分を表すラドン変換の形式で収集される。投影データは、図1に例示されるようなX線コンピュータ断層撮影装置やPET装置、SPECT装置等の任意のモダリティにより収集される。
【0069】
ステップ202において反復数やノイズ重み値、及びシステムマトリクス等を含むパラメータの組合せに基づいて、少なくとも1つの再構成フィルタが整形される。これらのパラメータは、任意に選択可能である。例えば、パラメータkは、所定の反復再構成アルゴリズムに従って、所定数の反復において特定の反復結果をエミュレートするために選択される。また、投影データに適用される再構成フィルタを整形するための反復パラメータkと共に、予め決定された重み値Wが使用されても良い。その他のフィルタパラメータは、様々なアプリケーションに対処できるように、反復アルゴリズム内のステップサイズと、ベイジアン条件の相対的重み付けとを含む。
【0070】
ステップS202においてパラメータの組合せに従って再構成フィルタが整形された後、ステップ203において、整形された再構成フィルタを投影データに適用する。フィルタが施されたデータは、ステップ204において逆投影される。逆投影により画像が発生する。ステップS205においてエミュレーションを反復するか否かが決定される。エミュレーションを反復する場合、エミュレーションの反復のために再びステップS202に進む。ステップS205においてエミュレーションを反復しない場合、本実施形態に係るノイズ重み付け反復エミュレーションFBPアルゴリズムは終了する。
【0071】
図12は、本実施形態に係るノイズ重み付けFBPアルゴリズムにおける上述のエミュレートされた反復の影響の概念を示す図である。従来のCT撮影においてノイズ制御は、事前フィルタリングまたは事後フィルタリングを使用して実行されている。対照的に本実施形態においてノイズは、ノイズ重み付け及び反復エミュレートされたFBPアルゴリズムを使用して画像再構成において実質的に低減される。画像再構成におけるノイズ低減は、画像再構成(特に、逆投影)が加算処理であり、投影に適切な重み付けを実行することにより逆投影の分散を削減できるという理念に基づいている。
【0072】
図12に示すように、2つの投影、すなわち、計測データが2本の投影線P1及びP2によって表され、解、すなわち画像Iは2本の投影線P1及びP2の交点である。横軸X1及び縦軸Y1は、画像Iにおける2つの変数を表す。1つの計測が画像から行われる場合、2つの画像画素は1本の直線に沿ったいずれかの場所の値をとることができる。2つの独立した計測が行われる場合、それら2本の直線の交点である、それら2つの画素の値が決定される。第3の計測が行われ、すべての計測がノイズを含む場合、それらの3本の直線は、一点において交差しない可能性がある。一般に、画像が多くの画素からなり、多くの矛盾するノイズを伴う計測が行われた場合、それらの計測は単一点で交わらない可能性があるため、撮像問題の解が存在しない可能性がある。
【0073】
数学では、ノイズ計測が矛盾する場合、それらの計測は「投影演算子の範囲内にない」と見なされる。他方で、計測が矛盾しない場合、それらの計測は投影演算子の範囲内にマッピングされる。実際には、投影演算子の範囲内のデータ成分だけが画像領域内に逆投影されるため、このマッピングはFBPアルゴリズムにおいて逆投影によって保証される。すなわち、この投影の矛盾はFBP再構成画像に寄与しない。一般性を失わずに、直線によって示されるように、多くの画素を有する画像の場合ですら、計測は矛盾しないものの、依然として雑音を伴い、投影演算子の範囲内にある状況を仮定することができる。反復アルゴリズムとは異なり、FBPアルゴリズムは、ノイズを伴う投影のセットから、ノイズを伴う単一の解を提供することが可能である。
【0074】
図12に示すように、ノイズがより少ない画像を発生するため、反復アルゴリズムの複数の反復ステップをエミュレートするようにFBPアルゴリズムは修正される。解軌道はk=0の初期値からk=∞の最終的な解に向けて示される。FBPアルゴリズムのその後の反復は、結果として生じる、より少ないノイズを伴って画像を移動させる。反復アルゴリズムは、通常、真の解に達成する前に終了する。これは、計測データが低い光子計数によるノイズを含むため、真の解決策がノイズを伴うものである場合に特に当てはまる。
【0075】
k=0からの1つの解軌道は、真の解決策、すなわち、画像Iに向けた直線であり、この直線は、エミュレートされた反復がノイズ分散によって重み付けされないことを示す。直線上の点NNW1及びNNW2は、それぞれ、第1の所定数の反復および第2の所定数の反復の後の中間結果を示す。この例では、第2の所定数の反復は、第1の所定数の反復よりも大きく、中間結果NNW2は、第2の所定数の反復の後にエミュレートされる。他方で、中間結果NW1、NW2、及びNW3は曲線上に位置し、この曲線は、それらのエミュレートされた反復がノイズ分散によって重み付けされていることを示す。中間結果NW1、NW2、及びNW3は、それぞれ、第3の所定数の反復、第4の所定数の反復、及び第5の所定数の反復の後である。この例では、第5の数は、第3の数の反復及び第4の数の反復より大きく、第3の数は最小である。中間結果NW3は、第5の反復の後にエミュレートされる。この解軌道は、重み係数のセットによって決定される。1つの重み係数がそれぞれの計測に割り当てられる。より大きな係数は、関連する計測がより信頼できるものであることを意味する。一般的な手法は、計測のノイズ分散の逆数に比例する重み係数を割り当てることである。
【0076】
実際には、真のノイズを伴う解の近傍にあるノイズがより少ない疑似解が受け入れられる。疑似解の手順は、ノイズ分散重み値と反復指数kとを含むパラメータの所定のセットによってパラメータ化される。一般に、より小さいkを伴う疑似解は、より平滑な画像に対応し、より大きなkを伴う疑似解は、より鮮明であるが、より多くのノイズを伴う画像に対応する。
【0077】
以下、本実施形態に係る上記の特徴の特定の組合せに関する追加の実施例が記述される。一実施例は、rFBP−MAPによって指定される、最大事後確率(MAP)アルゴリズムを用いた、投影線毎のノイズ重み付けFBPが利用される。rFBP−MAPアルゴリズムに関する導出ステップは、ビューごとの重み付けFBP−MAPアルゴリズムに関する導出ステップと非常に類似している。幾つかの主な導出ステップについて以下に簡単に説明する。目的関数は、(12)式と同じである。ベイジアン条件g(X)が二次式である場合、g(X)の勾配は、幾つかの行列Rに関してRXの形であり、解は、(18)式により表されるランドウェーバーアルゴリズムによって得られる。
【数20】
【0078】
(18)式に示される反復式は、等価の非反復式を有する。初期状態がゼロであるとき、非再帰的数式は、(19)式により表される。
【数21】
【0079】
(ATWA+βR)-1が存在する場合、(19)式中に示される総和は、(20)式により表される閉形式を有する。
【数22】
【0080】
(ATWA)-1が存在する場合、(ATWA+βR)-1=[I+β(ATWA)-1R]-1(ATWA)-1である。十分なサンプリングが逆ラドン行列A-1の存在をもたらすことが仮定される場合、(ATWA)-1=A-1-1(AT-1及び(ATWA)-1TW=A-1となる。従って(20)式は、(21)式のように簡略化可能である。
【数23】
【0081】
不十分なビュー数等の理由でA-1が存在しない場合、FBPアルゴリズムは、投影演算子Aの疑似逆数であり、(ATWA)+TW=A+であることが確認でき、式中、「+」は疑似逆数を示し、A+はランプフィルタリングの後に逆投影ATが続くことである。従って(21)式は、依然として、疑似逆数「+」によって置換することによって逆数「−1」を可とする。以下、この逆数及び疑似逆数を識別しないこととする。「A-1P」の表記は、従来のFBP再構成を表すために使用されることになる。
【0082】
(21)式において、[I−(I−αATWA−αβR)k]は、反復アルゴリズムの第k番目の反復におけるウィンドウ効果を表す。kが無限大になる傾向があるため[I−(I−αATWA−αβR)k]は、識別行列になる。[I+β(ATWA)-1R]-1は、ベイジアン正規化の効果を表す。β=0のとき、ベイジアン修正は効果的ではない。(21)式は、反復アルゴリズムがノイズ重み付け及びベイジアン正規化をどのように処理するかの見識と理解を提供し、(21)式は、同じタスクを実行するためのFBPタイプのアルゴリズムの構築をもたらし得る。
【0083】
最後に、FBPタイプのアルゴリズムが取得されるように、フーリエ領域表記が行列式から導出される。行列A-1は、1次元ランプフィルタリングの後に逆投影が続くか、又は逆投影の後に2次元ランプフィルタリングが続くように処理される。行列式(21)が「逆投影の後にフィルタリングする」アルゴリズムと見なされるとき、(21)式における画像領域フィルタリングは、(22)式に示す伝達関数を用いて2次元フーリエ領域において実行される。
【数24】
【0084】
(22)式中、vx及びvyは、それぞれ、x及びyに関する周波数である。
【数25】
【0085】
は2次元周波数ベクトルである。
【数26】
【0086】
よく知られているように、投影演算子Aが2次元空間内の線積分(すなわち、ラドン変換)であり、ATが演算子(すなわち、逆投影変換)である場合、投影と逆投影との組み合わされた演算子ATAは、元の画像の2次元カーネル1/rとの2次元畳み込みである。この場合、x−yデカルト座標において
【数27】
【0087】
である。1/rの2次元フーリエ変換は、
【数28】
【0088】
である。(22)式においては、Rのフーリエ領域等価(Fourier domain equivalence)はL2Dによって示される。
【0089】
中央断面定理を使用することによって、画像領域の2次元フィルタリングは、投影領域の1次元フィルタリングに変換される。2次元画像のフーリエ領域における伝達関数((22)式)は、下記の(23)式のように1次元投影のフーリエ領域における伝達関数に等しい。
【数29】
【0090】
(23)式中、ωは検出器に沿った変数に関する周波数である。ω=0であるとき、
【数30】
【0091】
を定義する。(23)式において、L1Dは(22)式のL2Dの中央断面である。rFBP−MAPルゴリズムでは、
【数31】
【0092】
は、1次元フーリエ領域で表現される修正ランプフィルタである。反復効果は、
【数32】
【0093】
によって特徴付けられる。ベイジアン効果は、
【数33】
【0094】
によって特徴付けられる。最終的な(k=∞)解は、β=0の場合、ノイズ重み付けwに依存しないが、β≠0の場合、ノイズ重み付けwに依存する。rFBP−MAPアルゴリズムに関する逆投影について修正は必要とされない。
【0095】
(22)式及び(23)式のノイズに依存する重み係数wは、定数ではない。ビュー単にのノイズ重み付けの場合、wはビュー角度に応じ、w=w(view)である。投影線単位のノイズ重み付けの場合、wは投影線に応じ、w=w(ray)である。重み係数を割り当てるための一般的な手法は、w(ray)を投影線計測のノイズ分散の逆数に比例させることである。この手法は、最適化問題に関する目的関数として、尤度関数(すなわち、結合確率密度関数)を使用することによって補われる。根本原理は、ノイズのより大きい計測よりもノイズのより少ない計測を信頼すべきであるということである。kが無限大になる傾向があり、βがゼロになる傾向があるとき、修正ランプフィルタ((23)式)は従来のランプフィルタになり、rFBP−MAPアルゴリズムは、従来のFBPアルゴリズムに低減される。
【0096】
rFBP−MAPアルゴリズムを実施する際、
【数34】
【0097】
の逆フーリエ変換は、1次元修正ランプフィルタの空間領域カーネルである
【数35】
【0098】
になると仮定される。rFBP−MAPアルゴリズムでは、wは投影線の関数であるため、w=w(t,θ)である。p(t,θ)をビューθ及び検出器上の位置tにおける投影にし、q(t,θ)をフィルタリングされた投影にする。この場合、q(t,θ)は、カーネル
【数36】
【0099】
がθ及びtに依存するため、畳み込みではない、以下の(24)式に示す積分により規定される。
【数37】
【0100】
フィルタリング処理が(24)式のように空間領域で実施される場合、計算コストは畳み込みと同じである。従って、rFBP−MAPアルゴリズムは計算的に効率的である。rFBP−MAPアルゴリズムを実施するための最も容易で最も効率的な方法は、空間領域で投影をフィルタリングするために(24)式を使用することである。
【0101】
現在、積分カーネル
【数38】
【0102】
に関する解析式は存在しない。q(t,θ)を計算する代替的な方法は、以下のように、周波数領域内の各ビュー角度θにおいて(24)式を実施することである。
【0103】
ステップ1:tに関するp(t,θ)の1次元フーリエ変換を見出して、P(ω,θ)を取得する。
【0104】
ステップ2:各投影線tについて重み係数w(t)を割り当て、周波数領域伝達関数((23)式)を形成する。ωは離散的な周波数指数であり、0、1、2、...などの整数値をとる。
【0105】
ステップ3:積
【数39】
【0106】
を見出す。
【0107】
ステップ4:ωに関するQ1(ω,θ)の1次元逆フーリエ変換を見出して、q(t,θ)を取得する。
【0108】
ステップ1は、各ビュー内の全てのtを処理し、一方、ステップ2から4は、一度に1つのtだけを処理する。
【0109】
上述のrFBP−MAPアルゴリズムは、平行ビームを用いる撮像ジオメトリのみに限定されない。このrFBP−MAPアルゴリズムは、FBPアルゴリズムが存在する限り、ファンビーム、コーンビーム、平面積分、及び減衰計測等の他の撮像ジオメトリに拡張可能である。FBPアルゴリズムに対する唯一の修正は、投影データフィルタリングである。各計測に重み係数が割り当てられ、周波数領域におけるフィルタ伝達関数は、撮像ジオメトリとは無関係な(23)式として形成される。
【0110】
rFBP−MAPアルゴリズムの効果を例示するために、低線量使用のX線コンピュータ断層撮影装置を使用して、死体の胴体が走査された。画像は、FBPアルゴリズム、ならびにrFBP−MAPアルゴリズムを用いて再構成された。撮像ジオメトリはコーンビームであり、X線源軌道は半径600mmの円である。X線検出器は64列を有し、列の高さは0.5mmであり、各列は896チャネルを有し、ファン角度は49.2°であった。管電圧は120kVであり、電流は60mAであった。360°に亘って1200個のビューが均一にサンプリングされた。
【0111】
FBPアルゴリズムとしてフェルトカンプアルゴリズムを使用し、データは、まず、余弦関数を用いて重み付けされ、次いで、1次元ランプフィルタがコーンビーム投影の各列に適用された。最後に、3次元ボリュームデータを発生するためにコーンビーム逆投影が使用された。rFBP−MAPアルゴリズムを実施するために、1次元ランプフィルタは、(23)式により表されるようにランプフィルタと置換された。反復指数kは、20,000と選択され、ステップサイズαは2.0に設定された。エッジ保存MAP再構成の場合、2つのβ値、すなわち、β=0.0005及びβ=0が使用された。ベイジアンペナルティ関数g(X)は、横断断面内の中央画素値とその4近傍の平均値との間の差の共通二次関数であった。すなわち、Rは、その1次元バージョンが[−1 2 −1]のカーネルを用いた二次導関数である単なるラプラシアンである。ノイズ重み付け関数は、w(t,θ)=exp(−p(t,θ))によって定義された。本実施形態においては、伝送計測が近似的にポワソン分布に従い、線積分p(t,θ)が伝送計測の対数であると仮定された。
【0112】
エッジマップは、2次元ソーベルカーネルを使用して断面毎に取得される。エッジ検出のための閾値は、最大画素値の70%に設定された。次いで、3×3平滑カーネルを二値画像に適用する。平滑カーネルの適用により二値画像に含まれるエッジがボケ、二値画像から(13)式に示すエッジマップBが発生される。
【0113】
画像ボリュームデータは、512×512×512の3次元配列で再構成され、表示のために非中央軸断片が使用される。その他の事前フィルタリング又は事後フィルタリングはこの再構成に適用されなかった。
【0114】
図13Aは、従来のFBPフェルトカンプアルゴリズムに基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。低線量のため、胴体を横切って左から右への著しいストリークアーチファクトが発生している。図13Bは、β=0を用いたrFBP−MAP再構成であるrFBP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。rFBPアルゴリズムにおける投影線単位のノイズ重み付けは、図13Aに出現したストリークアーチファクトを効果的に除去した。図13Cは、二次平滑化事前分布を使用し、β=0.0005を用いたrFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。図13Cに示すように、ノイズ重み付けにより、ストリークアーチファクトは消えている。ボケの事前分布により、画像は特にエッジ領域で過度に平滑に見える。図13Dは、実質的に図13B及び図13Cの重み付け加算である、エッジ保存事前分布を用いたrFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。エッジ検出ソーベルカーネルを用いて、重み付け加算画像BがX0から抽出された。エッジマップは二値画像であった。エッジ領域から非エッジ領域への遷移が平滑であるように、エッジマップをボケさせるためにローパスフィルタが使用された。図13Dの画像は、図13Bよりもノイズがより少ないが、図13Bの主なエッジを平滑でない状態で保持している。図13Eは、エッジ保存FBP−MAP再構成において使用されるエッジマップBを示す図である。
【0115】
計算効率を改善するために、多くのアプリケーションに関して、ビュー単位のノイズ重み付けvFBP−MAPアルゴリズムを使用することが可能である。ノイズ重み付けは、各ビューにおける最大ノイズ分散によって決定される。図13Fは、vFBP−MAP再構成に基づいて再構成された、死体の胴体領域内の横断切断に関する画像を示す図である。vFBP−MAP再構成は、ノイズ重み付け関数wを除いて、全てのパラメータが図13Dにおける再構成のパラメータと同一値に設定された。
【0116】
本実施形態に係る画像再構成においてノイズモデルは、ポワソン統計により規定される。あるいは、ノイズモデルは、ガウス分布により記述される電子ノイズを含む光子強度計測の複合ポワソン統計により規定される。確率論において複合ポワソン分布は、独立した一様に分布したランダム変数の和の確率分布であり、変数の数は、ポワソン分布に従う。ポワソン分布は、入射光子の数の記述に適している。蓄電時間の間に検出器に達する光子数Nは、N〜ポワソン(λ)であり、この場合、λは計数率である。検出された各光子は、ランダムエネルギーXiを有し、その分布はX線管の多色スペクトルと検出器効率とによって説明される。エネルギー積分検出器の場合、積分エネルギーは、Y=ΣNiによって与えられる。積分エネルギーYは、(25)式及び(26)式のように複合ポワソン統計(平均値および分散)に従う。
【数40】
【0117】
【数41】
【0118】
分散利得の概念は、Wによって示され、(27)式のように規定される。
【数42】
【0119】
ポワソン分散に従うランダム変数の場合、W=1である点に留意されたい。複合ポワソン変数の場合、Wは(28)式に従うとする。
【数43】
【0120】
従ってノイズ分散利得は、X線検出器への入射X線のスペクトルの平均有効エネルギーに比例する。一般に、平均スペクトルエネルギーに影響を及ぼす以下の係数が存在する。
【0121】
1)管電圧(kVp)。通常、管電圧は80kVpから140kVpまで様々であり、スペクトルエネルギーの上限を決定する。
【0122】
2)ボウタイフィルタ及びスペクトルフィルタ、ヒール効果。これらは、CTデータ内に常に存在する。ボウタイフィルタ及びヒール効果によりスペクトルは、各画素において異なる。従って各検出器画素に関して、独立して計算を行う必要がある。
【0123】
3)オブジェクトのサイズ及び構造。オブジェクトがビームを硬化させること、すなわち、E[X]を増大することはよく知られている。ビームハードニング効果は、経路長と減衰とに依存する。
【0124】
複合ポワソンモデルにおいてノイズ分散利得に関する平均スペクトルエネルギーの効果が利用される。収集データは、(29)式のように表される。
【数44】
【0125】
(29)式中、gは、検出器利得係数であり、eは分散Veを伴う電子ノイズである。Veを計測するために、X線管のスイッチをオフにしてデータを収集して、分散を計測する。平均及び分散は、各検出器画素に関して独立して時間方向に計測される。収集データの統計は、(30)式により得られる。
【数45】
【0126】
すなわち、(31)式が得られる。
【数46】
【0127】
オブジェクトに関する投影データを収集した後で、平均E[Z]及び分散Var[Z]を計測して、gWを計算することができる。利得gは投影データから容易に取得できない点に留意されたい。未知のスケールgの効果を補償するために、空気のスキャンデータ(WOBJ=W/WAIR)により投影データを較正する。ノイズ分散に関する複合ポワソンモデルは、(32)式により書き直される。
【数47】
【0128】
(32)式中、Wは(28)式により与えられる。
【0129】
本発明の構造および機能の詳細と共に、本発明の多数の特徴および利点が前述の説明に記載されているものの、本開示は、単なる例示であり、特に、部品の形状、サイズ、および構成、ならびにソフトウェア、ハードウェア、またはそれら両方の組合せの形の実装形態に関して詳細な変更を行うことは可能であるが、それらの変更は、添付の請求項が表現される条件の一般的な意味によって示される完全範囲まで本発明の原理内である点を理解されたい。
【0130】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【0131】
100…ガントリ、101…X線管、102…フレーム、103…X線検出器、104…データ収集回路、105…非接触データ伝送部、106…前処理部、107…回転部、108…スリップリング、109…高電圧発生部、110…システム制御部、112…記憶部、114…再構成部、114A…ノイズ重み決定部、114B…フィルタ整形部、114C…画像発生部、115…入力部、116…表示部、117…電流調整部、200…スキャン計画サポート装置
図1
図2
図3
図4
図5
図6A
図6B
図7
図8
図9
図10
図11
図12
図13A
図13B
図13C
図13D
図13E
図13F