TOP > 国内特許検索 > X線CT画像処理方法,X線CT画像処理プログラム及びX線CT画像装置 > 明細書

明細書 :X線CT画像処理方法,X線CT画像処理プログラム及びX線CT画像装置

発行国 日本国特許庁(JP)
公報種別 再公表特許(A1)
発行日 平成29年2月23日(2017.2.23)
発明の名称または考案の名称 X線CT画像処理方法,X線CT画像処理プログラム及びX線CT画像装置
国際特許分類 A61B   6/03        (2006.01)
FI A61B 6/03 350J
国際予備審査の請求 未請求
全頁数 31
出願番号 特願2015-516922 (P2015-516922)
国際出願番号 PCT/JP2014/002571
国際公開番号 WO2014/185078
国際出願日 平成26年5月15日(2014.5.15)
国際公開日 平成26年11月20日(2014.11.20)
優先権出願番号 2013102646
優先日 平成25年5月15日(2013.5.15)
優先権主張国 日本国(JP)
指定国 AP(BW , GH , GM , KE , LR , LS , MW , MZ , NA , RW , SD , SL , SZ , TZ , UG , ZM , ZW) , EA(AM , AZ , BY , KG , KZ , RU , TJ , TM) , EP(AL , AT , BE , BG , CH , CY , CZ , DE , DK , EE , ES , FI , FR , GB , GR , HR , HU , IE , IS , IT , LT , LU , LV , MC , MK , MT , NL , NO , PL , PT , RO , RS , SE , SI , SK , SM , TR) , OA(BF , BJ , CF , CG , CI , CM , GA , GN , GQ , GW , KM , ML , MR , NE , SN , TD , TG) , AE , AG , AL , AM , AO , AT , AU , AZ , BA , BB , BG , BH , BN , BR , BW , BY , BZ , CA , CH , CL , CN , CO , CR , CU , CZ , DE , DK , DM , DO , DZ , EC , EE , EG , ES , FI , GB , GD , GE , GH , GM , GT , HN , HR , HU , ID , IL , IN , IR , IS , JP , KE , KG , KN , KP , KR , KZ , LA , LC , LK , LR , LS , LT , LU , LY , MA , MD , ME , MG , MK , MN , MW , MX , MY , MZ , NA , NG , NI , NO , NZ , OM , PA , PE , PG , PH , PL , PT , QA , RO , RS , RU , RW , SA , SC , SD , SE , SG , SK , SL , SM , ST , SV , SY , TH , TJ , TM , TN , TR , TT , TZ , UA , UG , US
発明者または考案者 【氏名】前田 新一
【氏名】吉川 大悟
【氏名】田中 匠
【氏名】石井 信
出願人 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
個別代理人の代理人 【識別番号】110000822、【氏名又は名称】特許業務法人グローバル知財
審査請求 未請求
テーマコード 4C093
Fターム 4C093AA22
4C093CA13
4C093EA02
4C093FC24
4C093FD03
要約 X線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いX線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減可能なX線CT画像処理方法を提供する。投影X線の観測に関する確率分布を設定し統計推定を行う。確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程として表現される。多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有すること、検出器におけるX線の検出時間が一瞬ではなく所定の幅を有すること、X線のスペクトル分布が線スペクトルではなく所定の広がりを有することによって、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射するとして表現される。そして、多重X線和を表現する観測過程と物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、統計推定が事後分布の期待値によるベイズ推定を行う。
特許請求の範囲 【請求項1】
X線投影画像の観測に関する確率分布を設定し統計推定を行うX線CT画像処理方法であって、
X線投影画像と少なくともX線強度を含む計測条件を入力するステップと、
前記X線投影画像の観測に関する確率分布と、物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、前記統計推定が前記X線吸収係数及び前記物質に関する事後分布の期待値によって前記X線吸収係数及び前記物質を推定するステップと、
を備え、
前記X線投影画像の観測に関する確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程として表現され、
前記多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである、
ことを特徴とするX線CT画像処理方法。
【請求項2】
前記多重X線和に、X線波長に関する和を含める場合、
各画素の波長依存性を有するX線吸収係数に関する事前分布は、
前記物質毎およびX線波長毎のX線吸収係数の分布を表現するパラメータを用いて表現される、
ことを特徴とする請求項1に記載のX線CT画像処理方法。
【請求項3】
前記多重X線和に、X線波長に関する和を含める場合、
各画素の波長依存性を有するX線吸収係数は、
波長に依存しない各画素のX線吸収係数密度と、画素に依存せず前記物質毎に定められるX線波長間のX線吸収係数比の積として表現され、
前記物質毎に定められる、X線波長間のX線吸収係数比を表現するパラメータを用いて表現され、
各画素の波長依存性を有するX線吸収係数の事前分布が前記物質を条件とした各画素の前記X線吸収係数密度に対する条件付き事前分布と前記物質に対する事前分布によって表現される、
ことを特徴とする請求項1に記載のX線CT画像処理方法。
【請求項4】
前記物質に対する事前分布が、それぞれの物質の出現する度合いを示すパラメータとそれぞれの物質が空間的に連続する度合いを示すパラメータとで表現される事前分布で設定されることを特徴とする請求項3に記載のX線CT画像処理方法。
【請求項5】
前記物質毎に定められている前記X線吸収係数密度の条件付き事前分布が各物質のサブクラス毎に定められている前記X線吸収係数密度の条件付き事前分布と各物質のサブクラスの事前分布の積を各物質のとりえる全てのサブクラスに関して和をとった混合確率分布で表現され、前記各物質のサブクラスの事前分布が、各物質のサブクラスの出現する度合いを示すパラメータと前記各物質のサブクラスが空間的に連続する度合いを示すパラメータとで設定されることを特徴とする請求項3に記載のX線CT画像処理方法。
【請求項6】
X線投影画像の観測に関する確率分布を設定し統計推定を行うX線CT画像処理装置であって、
X線投影画像と少なくともX線強度を含む計測条件を入力する入力手段と、
前記X線投影画像の観測に関する確率分布と、物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、前記X線吸収係数及び前記物質に関する事後分布の期待値によって前記X線吸収係数及び前記物質を推定する推定手段と、
を備え、
前記X線投影画像の観測に関する確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程として設定され、
前記多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである、ことを特徴とするX線CT画像処理装置。
【請求項7】
X線投影画像の観測に関する確率分布を設定し統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
X線投影画像と少なくともX線強度を含む計測条件を入力するステップと、
前記X線投影画像の観測に関する確率分布と、物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、前記統計推定が前記X線吸収係数及び前記物質(及び物質のサブクラス)に関する事後分布の期待値によって前記X線吸収係数及び前記物質を推定するステップと、
を実行させ、
前記X線投影画像の観測に関する確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程として表現され、
前記多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである、ことを特徴とするX線CT画像処理プログラム。
【請求項8】
請求項7に記載のX線CT画像処理プログラムが搭載されたX線CT画像装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、体内組織構成を考慮したX線CT(Computed Tomography)画像処理方法、X線CT画像処理プログラム及びX線CT画像装置に関する。
【背景技術】
【0002】
X線CT(Computed Tomography)は、観測対象にX線を照射することで得られる投影データから観測対象内部をX線吸収係数の分布として画像化する技術である。
一般に、X線CTは、X線管にX線検出器を対向配置するとともに、これらの間にターンテーブルを配置した構成を取る。このターンテーブルの回転軸は、X線管とX線検出器を結ぶX線光軸に対して直交する向きとされる。X線CT装置においては、測定対象物に対して複数の角度からの透過X線データを収集してCT断層画像の再構成を行う。一般的には、180°以上、回転させて透過X線データを収集する。
【0003】
従来から、より低いX線被曝量で従来と同程度のCT画像を再構成する技術、メタルアーティファクトの低減技術など、X線CT再構成の高精度化を可能にするアルゴリズムが研究されている。上記のメタルアーティファクトは、測定対象物にX線を強く吸収する金属などの密度の高い高吸収物質があったときに、正確な再構成ができず再構成画像に人為的なノイズがのるというものである。X線CT再構成は、様々な投影角度において得られる投影像から、測定対象物を透過したX線の減衰の度合いを求め、内部の物質のX線吸収係数の分布を推定する手法である。しかし、高吸収物質はX線を強く減衰させるため、ほとんど透過X線が検出されず観測ノイズに紛れてしまう。特に既存手法であるフィルター補正逆投影法(FBP:Filtered Back Projection)を用いた場合、FBP法ではノイズの影響を受けやすいことから、再構成画像にはメタルアーティファクトと呼ばれる虚像が生じ、正確な画像の再構成ができないという問題がある。このようにFBP法はノイズの影響を受けやすいのであるが、主なノイズにはポアソン分布に従うショットノイズが含まれるため照射X線の強度を高くすることで信号雑音比を向上することができる。これによって相対的なノイズの影響を小さくし、画像再構成の精度を向上させることができるが、その一方で放射線被曝量が増大させてしまい、発癌などの確率的健康リスクを生じさせてしまう。
【0004】
このように、X線被曝量の低減とSN比はトレードオフの関係があり、より低いX線被曝量で従来と同程度のSN比の再構成画像を得る技術として、統計推定法が提案されている。統計推定法では、観測過程が確率的なノイズによる不確実性を含むことや、撮像対象の物体に関する事前知識を利用することで、低被曝環境下における推定の不良設定性を軽減することができる。 結果として、従来と同じX線強度や解像度でも再構成画像に含まれるアーティファクトを低減できる。
【0005】
現在、画像再構成法の主流であるフィルター補正逆投影法(FBP:Filtered Back Projection)や、FBPを改良したPCLIS(Projection Completion Method based on a Linear Interpolation in the Sinogram)は、統計推定を行わないものである。
統計推定法は、最尤推定法(Maximum Likelihood Estimation: MLE),MAP(Maximum A Posteriori)推定法,ベイズ推定法の3手法に大別される。
最尤推定法(MLE)は、推定対象である再構成画像に関する事前知識をおかない統計推定法で、観測における不確実性を含む物理過程を確率モデルで表現して、その確率モデルを基に最も尤もらしい再構成画像を推定するものである。観測における不確実性を含む物理過程には、観察される透過X線のフォトン数に関する揺らぎ(ショットノイズ)などが含まれる。
MAP推定法,ベイズ推定法は、推定対象に関する事前知識を確率モデルの形式で表現し、推定に取り入れるものである。すなわち、MAP推定法やベイズ推定法は、観測に関する確率モデルと推定対象に関する確率モデルの両方を勘案して、推定が行われる。 X線CTにおいては、推定対象はX線吸収係数である。具体的には、MAP推定法では、観測に関する確率モデルと推定対象となる物体のX線吸収係数に関する確率モデルから計算されるX線吸収係数の事後分布の最大値を推定値とし、ベイズ推定法では同じX線吸収係数の事後分布の期待値を推定値とする。MAP推定法では、事後分布が最大となる一点のみを取り上げて推定を行うが、ベイズ推定法では事後分布の期待値をとるため、事後分布が大きくなるX線吸収係数の集団に対する平均によって推定を行う。 事後分布が最大となる点は、観測ノイズによって敏感に変化しうるが、事後分布の期待値は観測ノイズに対してそれほど敏感に影響されないため、ベイズ推定法のほうがより安定し、精度の高い推定となる。
【0006】
ここで、通常のX線CTアルゴリズムの多くは画像再構成にあたって、照射するX線が線スペクトルであると仮定し単一の吸収係数を推定している。しかしながら、実際のX線は単一エネルギーではなく、X線管に加えられた最高管電圧で決まるエネルギーを最大として低い方へ連続的に分布するエネルギースペクトルで構成されている。そのため、物質にX線を照射するとエネルギーの高いX線はほとんど吸収されずに透過するが、エネルギーの低いX線は減弱が大きいため、物質を透過した後にはX線のスペクトル分布は高エネルギー側にシフトするビームハードニングと呼ばれる現象が生じる。
物質のX線吸収係数のエネルギー依存性(波長依存性)の影響によって生じる現象は、X線源が線スペクトルであると仮定し、物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮しない観測過程のモデルでは表現できず、正しくない推定が行われてしまう。このように物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮しないモデルによる推定でビームハードニングアーティファクトが生じる。
【0007】
特に、金属などの極端にX線吸収係数が高い高吸収物体があるとストリークアーティファクトと呼ばれる筋状のアーティファクトや高吸収物体周辺のX線吸収係数が実際より小さく見積もられるダークバンドアーティファクトが生じることが知られているが、これはビームハードニングの影響を考慮していないことが大きく影響していると考える。
これについて、以下に説明する。
金属など非常に高いX線吸収係数をもつ高吸収物体を透過すると、X線は大きく減弱する。その性質が観測モデルに表現されていれば、物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮しないモデルであっても、高吸収物体近辺に大きなX線吸収係数をもつ物体が存在するであろうことを推定できる。しかしながら、高吸収物体でも高エネルギーのX線は透過しやすい。そのため、実際にはX線源の主要なエネルギーをもつフォトンの多くが高吸収物体を透過中に吸収されていたとしても多少のX線が観測される。高吸収物体が離れた位置に複数ある場合、その複数の高吸収物体のいずれも透過するX線は、物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮しないモデルでは高いX線吸収係数をもつ物体が重なると透過X線がさらに大きく減衰することを予想する。 しかし、実際にはある一定量の高エネルギーのX線は高吸収物体を複数透過したとしても観測されてしまう。このため、物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮しないモデルで、単純に高吸収物体が存在すると仮定すると実際の観測と矛盾が生じてしまう。このようにして、物質のX線吸収係数のエネルギー依存性を考慮しない観測モデルと実際の観測とができるだけ矛盾しないようにX線吸収係数の推定を行うと、高吸収物体周辺のX線吸収係数が実際より小さく見積もられてしまうダークバンドアーティファクトが生じる。
【0008】
メタルアーティファクトを軽減する方法として、投影データを金属領域、非金属領域に分離し、それぞれの領域において画像再構成し、得た画像を重ね合わせるという方法が既に知られている。このようにすることで、金属領域における投影データの矛盾を正し、メタルアーティファクトの発生を軽減することができる。また、本願発明者らは、再構成画像となるX線吸収係数に関する事前知識をおいて統計推定を行うX線CT画像処理方法であって、事前知識を、再構成画像の各画素の領域において定義されるパラメータで、撮像対象となる人体の各組織の存在割合を表現するパラメータと、各組織のX線吸収係数を表現するパラメータと、各組織の空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現するものを用いることで、メタルアーティファクトを低減可能なX線CT画像処理方法を提案している(特許文献1を参照)。
しかし、これらの方法は、X線吸収係数のエネルギー依存性(波長依存性)を考慮するわけではないため、実際の観測過程と物質のX線吸収係数のエネルギー依存性を考慮しない観測過程とのズレのために生じるビームハードニングアーティファクトを除去するには至っていない。
【0009】
そして、ビームハードニングアーティファクトを軽減する手法として、X線吸収係数のエネルギー依存性(波長依存性)を考慮したモデルが知られている(非特許文献1~4)。X線吸収係数のエネルギー依存性(波長依存性)を考慮するモデルでは、X線のエネルギー毎のX線吸収係数を推定する必要がある。これは、より正確な物理過程に基づくが、推定問題としては不良設定性の問題を生じさせる。たとえば、単純のためX線のエネルギーを離散化し、二種類のエネルギーにおけるX線吸収係数を推定する問題を考えたとする。このとき、観測される投影像は従来と同じであっても、推定しなければならないX線吸収係数は従来に比べて2倍に増えてしまう。このように解となるX線吸収係数の自由度が高まるため、実際の観測とそれほど矛盾を生じない解がいくつも生じてしまう。この場合、ちょっとした観測ノイズに推定結果が敏感に影響され安定しない。 そこで、解となるX線吸収係数の自由度を抑える制約が必要となる。 先ず、非特許文献1~3では、各画素におけるエネルギー依存の吸収係数を既知の物質の吸収係数関数の重み付き和で表現することで、物質のエネルギー依存性(波長依存性)を表現している。X線吸収係数のエネルギー依存性(波長依存性)は、物質毎に異なるが高エネルギーになるに従って減衰する傾向など各物質間で共通する性質が存在する。X線のエネルギーを入力変数とし、そのエネルギーでのX線吸収係数を出力する関数を吸収係数関数と呼ぶことにすると、異なる物質であっても類似した吸収係数関数をもつため、どの物質の吸収係数関数も少数の、典型的には二つの物質の吸収係数関数の重み付き和(基底関数表現)でよく表現できる。例えば、非特許文献2では、物質クラスの決定をFBPで得られた画像をもとに行い、その知識をもとに画像再構成を行っている。非特許文献4では、画素毎のX線吸収係数は、既知の物質種毎のX線吸収係数の組み合わせで表現できるという仮定で推定を行っている。
このように、X線吸収係数のエネルギー依存性(波長依存性)を考慮する方法は、その不良設定性を抑制するため、既知の物質のX線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らすことが行われてきた。しかしながら、その一方で、強すぎる制約が課せられる問題が指摘されている。
【0010】
一方、ハードウェア的には、DECT(Dual Energy Computed Tomography)の登場により、異なるX線管電圧で撮影し、X線吸収係数の違いから物質のエネルギー依存性(波長依存性)を推定できるようになった。このDECTを用いた様々な画像再構成法が知られている(非特許文献5~8)。しかし、DECTを用いた画像再構成法は、一般的に単色X線CTより被曝量が大きく、同程度の被曝量に抑えようとすると推定精度が悪くなってしまうなどの問題がある。また、非特許文献5は、事後的に画像処理によって補正する方法であり、その補正によって更なる虚像を生じさせる可能性もあり限界がある。
【0011】
物質のX線吸収係数のエネルギー依存性を考慮しない観測モデルは、検出器に入射するX線が複数の波長からなるのではなく、ある代表的な波長のみからなることを想定することになる。 これは、検出器に複数のX線スペクトルが多重に入射する過程を、単純化したものと考えることができる。このようにX線が多重に検出器に入射する現象は、スペクトル空間だけではなく、時間的、空間的にも生じる。すなわち、検出器はある一定の検出時間内に多重のX線が入射し、ある一定の面積の検出面内に多重のX線が入射する。これら投影したX線が多重に検出器に入射する観測過程において、検出器に入射するX線を一本の光線で代表させた場合、X線吸収係数のエネルギー依存性を考慮しない時にビームハードニングアーティファクトが生じるのと同様にX線吸収係数の推定にバイアスが生じることが知られている(非特許文献9)。つまり、投影したX線が多重に検出器に入射する観測過程においては、一本の光線ではなく(X線の多重性を無視しないで)、複数の光線を想定した多重のX線の和と捉えることが望ましい。
【先行技術文献】
【0012】

【特許文献1】特開2011-156302号公報
【0013】

【非特許文献1】Joseph A. O ’Sullivan, and Jasenka Benac, Alternating Minimization Algorithms for Transmission Tomography, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 26, NO. 3 (2007), 283-297
【非特許文献2】Idris A. Elbakri and Jeffrey A. Fessler, Senior Member, IEEE, Statistical Image Reconstruction for Polyenergetic X-Ray Computed Tomography, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 21, NO. 2 (2002), 89-99
【非特許文献3】Joseph A. O’Sullivan, Bruce R. Whiting, Donald L. Snyder, and Orville A. Earl, Image Reconstruction from Data Acquired With an X-Ray Computerized Tomographic System Having Energy-Integrating Detectors
【非特許文献4】Chye Hwang Yan, Robert T. Whalen, Gary S. Beaupre, Shin Y. Yen, and Sandy Napel, Reconstruction Algorithm for Polychromatic CT Imaging: Application to Beam Hardening Correction, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 19, NO. 1, (2000)
【非特許文献5】A J Coleman and M Sinclair, A beam-hardening correction using dual-energy computed tomography, Phys. Med. Biol., 1985, Vol 30, No. 11, (1985), 1251-1256
【非特許文献6】Jeffrey A. Fessler, Idris Elbakri, Predrag Sukovic, Neal H. Clinthorne, Maximum-likelihood dualenergy tomographic image reconstruction, Medical Imaging 2002: Image Processing, Milan Sonka, J. Michael Fitzpatrick,38 Editors, oceedings of SPIE Vol. 4684 (2002), 38-49
【非特許文献7】Joonki Noh, Jeffrey A. Fessler Paul E. Kinahan, Statistical Sinogram Restoration in Dual-Energy CT for PET Attenuation Correction, IEEE Trans Med Imaging, (2009) 1688-1702
【非特許文献8】山崎陽一, 戸田尚宏, 低被爆環境下におけるDual-Energy X線CT雑音低減アルゴリズム, IEICE Technical Report MBE2008-113, (2009), 101-106
【非特許文献9】Ruhrnschopf, E. P., G. Schwierz, et al. (1981). Nonlinearity and Inhomogenity Effects in Computerized Tomography Due to the Exponential Attenuation of Radiation. Mathematical Aspects of Computerized Tomography. G. T. Herman and F. Natterer, Springer Berlin Heidelberg. 8:252-269
【発明の概要】
【発明が解決しようとする課題】
【0014】
上述の如く、波長が短くエネルギーの強い硬X線は、波長が長くエネルギーの弱い軟X線が物体を透過する力に比べて、物体を透過しやすいため、物体に入射したX線のスペクトル分布は物体を通過後、軟X線が減衰し硬X線側にパワースペクトルのピークが偏ったものとなり、ビームハードニングと呼ばれる現象を生じている。
理想的には、投影によって得られるX線のスペクトル分布を分光法やフィルターを用いて推定し、各スペクトル強度のX線を測定し、その測定をもとに各スペクトルのX線吸収係数の推定を行えば問題は解決するのであるが、実際に分光すると各スペクトルで得られる信号のSN比が下がるため推定精度は悪化するという問題がある。 そのため従来から、情報処理アルゴリズムでビームハードニングアーティファクトの除去を行うことが要望されている。
【0015】
様々なビームハードニングアーティファクトの除去を行う情報処理アルゴリズムが提案されている中で、物質のX線吸収係数の波長依存性を考慮する方法があり、特に、組織のX線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らす技術がある(非特許文献1,2を参照)。
【0016】
しかし、非特許文献1,2では、組織の種類を離散変数で表現した組織クラスを決定論的に推定するために、想定する組織クラス数が適切でない場合や、推定された組織クラスが誤っていた場合、解が大きく異なってしまい、既知の組織クラスで表現されない組織に対して適切な表現ができないといった問題がある。
また、推定する組織クラスに関して制約がおかれなかったため、不良設定性の解消が十分でなく、組織クラスの推定が不安定になりがちであるといった問題がある。
特に、非特許文献1では、波長間でX線吸収係数の比が可変であり、異なる組織クラスを重み付けする。組織クラス数が少ないと表現できないX線吸収係数が現れる可能性があり、組織クラス数が多いと拘束の意味をなさずに不良設定性が生じてしまう。
【0017】
そこで、本発明は、従来技術を用いた画像再構成と比べて、物質の種類を表す離散変数である物質クラスを介してX線吸収係数を「確率的」に表現し、物質クラスに確率的な拘束を設けることで不良設定性を抑制する一方で、各画素がどの物質クラスに所属するかを決定論的ではなく確率的に推定することでより柔軟なX線吸収係数の表現を可能とし、より低いX線被曝量で従来と同程度の再構成画像が取得可能で、アーティファクトを低減可能なX線CT画像処理方法、X線CT画像処理プログラムおよび該プログラムが搭載されたX線CT画像処理装置を提供することを目的とする。 ここで、本発明において物質クラスは、必ずしも生体の組織とは限らないので組織クラスを包含する意味で用いている。
【0018】
また、検出器に入射するX線は、一本の光線ではなく、複数の光線を想定した多重のX線の和と捉えることが必要である。上述の如く、このX線の多重性を無視して、一本の光線で代表させることによりバイアスが生じる。
そこで、本発明は、投影したX線が多重に検出器に入射する観測に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減可能なX線CT画像処理方法、X線CT画像処理プログラムおよび該プログラムが搭載されたX線CT画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【0019】
上記状況に鑑みて、本発明のX線CT画像処理方法は、X線投影画像の観測に関する確率分布を設定し統計推定を行う。
上記の確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程を表現する。ここで、多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである。
そして、多重X線和を表現する観測過程と物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、統計推定がX線吸収係数及び物質に関する事後分布の期待値によって線吸収係数及び物質を推定する。
【0020】
ここで、多重X線和で表現される観測過程は、以下の(a),(b)の特徴を備える。
(a) 各画素の波長依存性を有するX線吸収係数は、各画素のX線吸収係数密度と、物質毎に定められているX線波長間のX線吸収係数比の積で表現される。X線吸収係数の波長依存性を考えない場合のX線吸収係数の推定は、波長依存性を有するX線吸収係数において想定するX線波長域が1つしかないと考えてX線吸収係数密度のみを推定することと等価となる。そのため、波長依存性を有するX線吸収係数の各画素のX線吸収係数密度と、物質毎に定められているX線波長間のX線吸収係数比の積による表現は、X線吸収係数の波長依存性を考えない場合のX線吸収係数の表現の拡張となっている。
(b)投影したX線が多重に検出器に入射する多重X線和の観測過程を表現するものであり、
多重X線和は、
検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有すること、
検出器におけるX線の検出時間が一瞬ではなく所定の幅を有すること、
X線のスペクトル分布が線スペクトルではなく所定の広がりを有すること、
によって、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射することを想定する。そして、想定した複数の光線のそれぞれで投影を計算し、それらの投影によって励起された(電気)信号の和が検出器で観測されるとする。
【0021】
また、撮像対象がどのようなX線吸収係数をとる傾向があるかを表すX線吸収係数に関する事前分布は、物質クラスという変数を介して決められる。 物質クラスは、撮像対象の内部がどの物質から構成されるかを表す撮像対象内部の各画素で定義される変数である。物質クラスは、どの物質クラスがどの程度生じやすいか、またそれぞれの物質クラスの空間的滑らかさの度合いについて確率的な制約を加えることにより波長依存のX線吸収係数の自由度を確率的に制約する。
また、物質クラスが定まった時に、どのようなX線吸収係数密度をとりやすいかという確率分布を2つのガンマ分布の和で表現される場合には、各ガンマ分布のいずれに所属するかを表現する物質クラスのサブクラス変数を導入し、物質クラスとサブクラスの組み合わせに対してどの物質クラスとサブクラスの組み合わせがどの程度生じやすいか、またそれぞれの物質クラスとサブクラスの組み合わせがどの程度、空間的に滑らかさであるか、の確率的な制約を加えることにより波長依存のX線吸収係数の自由度を確率的に制約する。
【0022】
本発明のX線CT画像処理方法によれば、ビームハードニングアーティファクトをはじめとする多重X線和の観測過程を無視することで生じるアーティファクトを除去するために、物質のX線吸収係数の波長依存性など多重X線和の観測過程を考慮した確率分布を設定しX線CT画像処理を行う。ここで、物質のX線吸収係数の波長依存性を考慮する際には、波長ごとのX線吸収係数を全て未知変数として推定すると、問題の不良設定性のために信頼のおける解が得られないことから、いくつかの既知の物質が特定の波長依存性をもったX線吸収係数を有することを表現する事前分布を利用することにし、再構成画像の各画素が既知の物質クラスのうちのいずれの物質クラスから構成されるかを推定するという、より易しい推定問題に変更することで、より信頼のおける解を得る。すなわち、波長ごとに正の実数変数であるX線吸収係数を推定する代わりに、波長に依存しない正の実数変数であるX線吸収係数密度と離散変数である物質クラスを推定する。この各画素の物質クラスは直接的に観測できない隠れ状態変数、すなわち物質クラスの変数を隠れた確率変数として扱う。
【0023】
この隠れ状態変数の物質クラスに対して、空間的滑らかさの確率的な制約を課すことで、波長依存のX線吸収係数の自由度を確率的に制約することとした。その上で隠れた変数をもつ確率分布の推論として、ベイズ推論を近似的に行うことで安定した精度の高い推定(収束が保証され、観測ノイズに対してあまり敏感に反応しない安定した推定)を行うのである。また、本発明におけるベイズ推定によれば、X線吸収係数と再構成画像の各画素の領域が属する物質クラスの両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。ベイズ推定により、各画素における物質クラスの所属確率と各画素のX線吸収係数密度を確率的に表現できるので、決定論的に物質クラスを選択する場合に比べてより柔軟性をもったきめ細やかなX線吸収係数の表現ができる。
【0024】
本発明のX線CT画像処理方法では、物質クラスの変数を確率変数として扱うことで、既知の物質クラスで表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、本発明のX線CT画像処理方法では、物質クラスに対して事前知識として空間的滑らかさなどの制約を課すことで不良設定性を抑制することができる。 また、本発明のX線CT画像処理方法では、近似ベイズ推論を行うことで、収束性が保証され、かつ、推定の不安定性(オーバーフィッティング)を抑えることができる。
すなわち、本発明のX線CT画像処理方法を用いることにより、物理過程により忠実でかつ精度の高い安定した推定が行えるのである。
【0025】
本発明において、物質クラスは、組織クラスを包含する意味で用いている。物質クラスは、離散変数で表現される。物質クラスの存在割合(密度)とは、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス(筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど)がどの程度の割合で分布するのを示すもので、非負のスカラーのパラメータである。
物質クラス毎に波長間で固定されたX線吸収係数比を表現するパラメータとは、確率的に表現されたX線吸収係数において、波長間で吸収係数の比は固定と仮定し、それを事前知識の確率分布に用いて、波長依存のX線吸収係数の自由度を制約するものである。
また、物質クラスが空間的に連続する度合いを表現するパラメータは、例えば、空気の物質クラスは画素同士の連結がしやすい(空気中に他の組織は入らない)など空間的に連続する度合いが高いという制約を課すことができるので、物質クラスが空間的に連続する度合いを事前知識の確率分布に用いて、波長依存のX線吸収係数の自由度を制約するものである。
【0026】
上記のX線CT画像処理方法によれば、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス(筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど)がどの程度の割合で分布するか、また、それぞれの物質クラスがどの程度のX線吸収係数を有するか、また、それぞれの物質クラスがどの程度、空間的に連続して分布しやすいかの分布情報を確率分布の形で表現して、事後確率の期待値によるベイズ推定を用いて、画像再構成や物質推定を行う。
なお、物質クラスに関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と物質クラスの推定精度のさらなる向上を図ることが可能である。
【0027】
本発明のX線CT画像処理方法における多重X線和にX線波長に関する和を含める場合を想定する。
この場合、各画素の波長依存性を有するX線吸収係数に対する事前分布を、物質クラス毎およびX線波長毎のX線吸収係数の分布を表現するパラメータを介して特徴づけることによって、どのようなX線吸収係数が生じやすいかを確率的に限定するものである。
【0028】
他の観点から、本発明のX線CT画像処理方法における多重X線和にX線波長に関する和を含める場合、各画素の波長依存性を有するX線吸収係数は、物質クラスに依存しない各画素のX線吸収係数密度と、物質クラス毎に定められているX線波長間のX線吸収係数比の積で表現される。
そして、各画素の波長依存性を有するX線吸収係数は、物質クラス毎に定められている、X線波長間のX線吸収係数比を表現するパラメータで特徴付けられる。
そして、各画素の波長依存性を有するX線吸収係数の事前分布が物質クラスを条件とした各画素のX線吸収係数密度に対する条件付き事前分布と物質クラスに対する事前分布によって表現される。
すなわち、X線吸収係数を各画素のX線吸収係数密度と物質クラス依存波長依存の定数の積で表現する場合、X線吸収係数の事前分布は、各画素のX線吸収係数密度に対する事前分布と物質クラスに対する事前分布から一意に決定されることになる。
【0029】
ここで、本発明のX線CT画像処理方法において、上記の物質クラス毎に定められているX線吸収係数密度の分布が、混合ガンマ分布で表現される。また、物質クラスが、出現する度合いを表現するパラメータと物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現される。
これは、X線吸収係数の事前分布の表現についての一例である。
物質クラスを条件としたX線吸収係数密度に対して混合ガンマ分布を仮定し、物質クラス変数に対してボルツマン分布を仮定することによって、X線吸収係数の事前分布が決定されるようにする。
【0030】
物質クラス毎に定められているX線吸収係数密度の分布として(混合)ガンマ分布で表現されるとした理由は、その分布形状が単峰な分布を表現でき、かつその確率変数のとる値が必ず正の値となることを保証できることに着目したものである。また、投影したX線が多重に検出器に入射する多重X線和の観測過程において検出器において観測される信号の分布としてガウス分布で表現されるとした理由は、光子散乱や量子化雑音などの雑音により、理想値からずれた値を示すことと、投影したX線が多重である場合にも計算が用意であることを考慮して、これらの雑音をガウス分布でよく表現されると仮定したものである。
また、物質クラスが、出現する度合いを表現するパラメータと前記物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現されるとした理由は、同じ組織が空間的に集まりやすい(組織別に集まりやすさの調整が可能)ことと、標準とされる各組織の割合を表現しやすいことに着目したものである。
【0031】
また、本発明のX線CT画像処理方法において、ベイズ推定は、物質クラス、X線吸収係数密度、物質クラスのサブクラス変数、それぞれに関する事後分布の推定が行われる。
ここで、サブクラス変数は、物質クラスが定まった時に、どのようなX線吸収係数密度をとりやすいかという確率分布を2つのガンマ分布の和で表現される場合に、いずれのガンマ分布に所属するかを指定する離散変数である。
物質クラスが与えられた条件下でのX線吸収係数の事前分布は、特定の値(たとえば平均値)を取りやすいと考えられるものの、その平均値に対して対称な確率分布となるとは限らなかったり、確率分布の尖度が大小様々であったりすることが予想できることから、このような複雑な事前分布を表現できるように2つの確率分布の重み付き和である混合確率分布で表現したものである。さらにX線吸収係数は正の値しか取らないため、この混合確率分布をガンマ分布の混合確率分布で表現されると仮定したものである。物質クラスのサブクラス変数は、2つ分布の混合確率分布においてどちらの確率分布に所属するかを指定するために便宜的に導入される0又は1の二値変数である。
【0032】
本発明のX線CT画像処理装置は、X線投影画像の観測に関する確率分布を設定し統計推定を行うX線CT画像処理装置であって、X線投影画像と少なくともX線強度を含む計測条件を入力する入力手段と、観測に関する確率分布と、物質クラスという再構成画像の各画素にどの物質が存在するかを表す変数を介して表現されるX線吸収係数に関する事前分布と、物質クラスに関する事前分布のもとで、X線吸収係数及び物質クラスに関する事後分布の期待値によるベイズ推定でX線吸収係数及び物質クラスを推定する推定手段を備える。ここで、観測に関する確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程を通じてどのような投影像が得られやすいかを表す確率分布として設定される。また、多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである。
【0033】
上記のX線CT画像処理装置によれば、物質クラスの存在割合(密度)を表現するパラメータを用いて、物質クラスの変数を確率変数として扱うことで、予め想定していた物質クラスから外れる物質に対しても、予め想定した物質クラスの重み付き和によって表現することで、近似表現が可能になる。また、物質クラスが空間的に連続する度合いを表現するパラメータを用いることにより、物質クラスに対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定精度が向上することが期待できる。
【0034】
また、本発明の他の観点によれば、X線吸収係数に関する確率分布をおいて統計推定を行うX線CT画像処理プログラムであって、コンピュータに、以下の1)~2)のステップを実行させるためのプログラムである。
1)X線投影画像と少なくともX線強度を含む計測条件を入力するステップ
2)上記のX線投影画像の観測に関する確率分布と、物質に関するパラメータをもつX線吸収係数に関する事前分布のもとで、X線吸収係数及び物質に関する事後分布の期待値でX線吸収係数及び物質を推定するステップ
ここで、上記の確率分布は、投影したX線が多重に検出器に入射する多重X線和の観測過程として表現され、多重X線和は、検出器におけるX線の検出面が空間的に1点ではなく所定の広がりを有し、X線の検出時間が一瞬ではなく所定の幅を有し、かつ、X線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のX線が検出器に多重に入射し、それらの多重のX線の重み付き和が検出器で観測されることである。
【0035】
上記のX線CT画像処理プログラムによれば、物質の存在割合(密度)を表現するパラメータを用いて、物質の変数を確率変数として扱うことで、既知の物質で表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、物質が空間的に連続する度合いを表現するパラメータを用いることにより、物質に対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定の不安定性を抑えることができる。
【0036】
また、本発明は、上記のX線CT画像処理プログラムを搭載したX線CT画像処理装置を提供することができる。
【発明の効果】
【0037】
本発明によれば、従来技術を用いた画像再構成と比べて、X線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いX線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減できるといった効果を有する。
また、本発明によれば、投影したX線が多重に検出器に入射する多重X線和の観測過程に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減できるといった効果を有する。
【図面の簡単な説明】
【0038】
【図1】本発明のX線CT画像処理方法の処理フロー
【図2】X線CTの説明図
【図3】更新の説明図
【図4】金属と骨のX線吸収係数の確率分布図
【図5】実施例1の計算機実験での画像再構成を示す図
【図6】実施例2の計算機実験での画像再構成を示す図
【図7】平行ビームにおける投影の幾何の説明図
【図8】平行ビームにおける投影と画素化の説明図
【図9】回転時の投影の幾何の説明図

【発明を実施するための最良の形態】
【0039】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。
【0040】
本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、X線吸収係数に関する事前分布を設定し統計推定を用いて画像の再構成を行う。特に、X線スペクトルのエネルギー分布(波長分布)を考慮し、エネルギー依存(波長依存)のX線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングアーティファクトを解消する。
推定する変数増加による不良設定性に対処するため、適当な事前知識を用いて解に何らかの制約を課すことで解決を試みる必要がある。本発明では基底関数表現を用いて、その基底を既知の物質のエネルギー依存(波長依存)のX線吸収係数とし、その重みはどの物質クラスに所属するかを表現する離散変数と、X線吸収係数密度を表現する非負のスカラーの連続変数の積で与えられるとし、物質クラスを条件としたX線吸収係数密度の事前分布と物質クラスの事前分布を課して解に制約を課す。
【0041】
また、本発明では、人体組織を対象とする再構成画像の場合、人体組織は脂肪,筋肉,骨といった限られた物質から構成され、それぞれの物質のX線吸収係数の大まかな分布はあらかじめ分かっているので、その知識を事前分布として表現して推定に利用する。 この推定のためには、各画素がどの物質クラスに所属するかの推定を行う必要があるが、それぞれの物質は空間的な連続性を有しており、各物質が人体に占める割合も大まかに分かっているという状況を想定できることから、物質クラスを表す離散変数に空間的な滑らかさ制約や占有割合などの事前知識を導入することにより不良設定性の抑制を行う。
【0042】
投影X線の観測に関する観測過程については、投影したX線が多重に検出器に入射する多重X線和の観測過程を表現するものであり、X線が有限の広がりを持つ波長域(スペクトル分布)を有することによって、波長が異なる複数のX線が検出器に入射することを表現する。
再構成画像に関する確率分布については、再構成画像の各画素の領域において、X線吸収係数がどのような値をとりやすいかを表現するかを特徴付けるものであって、各物質クラスと各サブクラス変数の組み合わせがX線波長毎にどのようなX線吸収係数密度を取りやすいかを表現する混合ガンマ分布と、各物質クラスと各サブクラス変数の組み合わせがどの程度、出現しやすいかを表現するパラメータと、各物質クラスと各サブクラス変数の組み合わせがどの程度、空間的に滑らかに連続しやすいかを表現するパラメータからなるボルツマン分布で表現する。
【0043】
そして、X線吸収係数と物質クラスに関する事後分布の期待値推定による統計推定を行う。画像再構成において、この事後分布を求めることが必要であるが、各画素で定義される物質クラスという高次元の隠れ離散状態変数に関する和計算が含まれるため、解析的に計算を実行することは困難である。よって、本発明では、近似手法を適用することでこの計算困難さを克服する(図1のフローチャートを参照)。
【0044】
(問題の定式化)
X線CT画像のように、様々な方向から投影されたT個の投影データを、D={Y(1),・・・,Y(T)}と表すことにする。各々のデータY(t)は、t番目の投影によって検出器で検出されるデータの集合であり、Y(t)={y(1),・・・,y(t)}となる。但し、Tは投影数、Iは検出器の数、y(t)はi番目の検出器で検出された光子の数を表す。
そして、X線は物質を透過する際に指数的に減衰することから、下記数式(1)のように表すことができる。
【0045】
【数1】
JP2014185078A1_000003t.gif

【0046】
ここで、xは観測対象を撮像対象のX線吸収係数をラスタスキャンして得られるJ次元ベクトルx={x,・・・,x}のj番目の画素のX線吸収係数である。b(t)は、X線源から放出される光子数(物体が何も置かれていないときに観測され得る光子数)を表す。また、lij(t)は角度θ(t)から投影された際のi番目の検出器で検出される投影線とj番目の画素が交差する距離に相当するものであり、lij(t)xがt番目の投影のi番目の検出器に入射するX線に対するj番目の画素領域の実効的なX線吸収係数を表す(図2を参照)。
【0047】
ここで、X線のエネルギースペクトルをE通りに離散化して表現することを考える。入射X線の内、エネルギーe(e=1,・・・,E)に属するX線の光子数をN(t)とすると、物質に透過した後に観測される光子数の平均値は、下記数式(2)のように表すことができる。
ただし、xjeはエネルギーeのX線に対するj番目の画素のX線吸収係数である。
【0048】
【数2】
JP2014185078A1_000004t.gif

【0049】
(X線吸収係数のエネルギー依存性について)
j番目の画素におけるエネルギーeでのX線吸収係数xjeを、エネルギー毎、画素毎に自由に決定すると不良設定性が生じる。実際には、いずれの物質も高エネルギーになるにつれてX線吸収係数が小さくなるという物質間で共通した性質を有する。つまり、異なるエネルギー間のX線吸収係数の比に特定の傾向が存在する。そこで、物質種を区別する物質クラスを導入して、X線吸収係数を制約した形で表現する。
観測対象の物質がC種類の物質クラスに分類されると仮定する。ここで、j番目の画素の物質クラスを変数zを用いて表現する。ただし、zは、z={zj1,・・・, zjC}と表される。各要素zjcは0もしくは1をとる二値変数であり、j番目の画素が物質クラスcに属している場合に1をとり、他の場合は0をとる。この物質クラスに依存して一意に決まるエネルギー間のX線吸収係数の比を表現する変数として、変化率rce(>0)を定義する。ただし、変化率rceは、e=1の時rce1=1と正規化されているものとする。e=1におけるj番目の画素のX線吸収係数をx(>0)で表現する。以上の表現を用いて、X線吸収係数xjeは下記数式(3)のように表すことができる。
【0050】
【数3】
JP2014185078A1_000005t.gif

【0051】
上記数式(3)の表現を用いた場合、各画素で推定する必要がある変数は、スカラーの連続変数となるxとC種類のいずれかの値をとる離散変数zの2つだけとなる。
【0052】
(観測モデルについて)
各検出器で観測される光子数の理想値は、上述の数式(1)に示した通りである。しかし、実際には光子散乱や量子化雑音などの雑音によって理想値からずれた値を示すことがわかっている。これらの雑音が、ガウス分布でよく表現されると仮定すると、各検出器における観測データは下記数式(4)のように表すことができる。但し、X={x,・・・,x}、Z={z,・・・,z}である。
また、観測が、検出器毎、投影毎に独立であると仮定すると、投影データの集合D(観測データ)は下記数式(5)のように表すことができる。但し、Tは投影数、Iは検出器の数である。ここで、ガウス分布の標準偏差を表すσは、ショット雑音の性質を考慮して、下記数式のように観測y(t)に依存する関数として良い。
【0053】
【数4】
JP2014185078A1_000006t.gif

【0054】
【数5】
JP2014185078A1_000007t.gif

【0055】
(事前分布について)
事前分布については、観測対象の物質種ごとのX線吸収係数に関する事前知識を表現する。画素jでのX線吸収係数xは画素jにおける物質クラスzjcのみに依存して決定されると仮定すると、事前分布p(X|Z)は下記数式(6)のように表すことができる。
【0056】
【数6】
JP2014185078A1_000008t.gif

【0057】
物質クラスzjcが与えられたもとでは、X線吸収係数xは特定の値を取りやすいことが予想できる。ここでは事前分布p(x|z)は、分布の非対称性や様々な異なる尖度を表現できるよう2つの分布の混合分布として下記数式(7)のように表すことができる。但し、変数bは二値変数である。また、ベクトルBをB={b,・・・,b}で定義する。
【0058】
【数7】
JP2014185078A1_000009t.gif

【0059】
画素毎の物質クラスzjcと変数bjdが与えられたとき、X線吸収係数xはガンマ分布に従うと仮定すると下記数式(8)のように表すことができる。但し、ucd,vcdはガンマ分布のパラメータであり、平均ucd/vcd、分散ucd/vcdである。これらのパラメータは物質クラスに関する経験的な知識により定める。
【0060】
【数8】
JP2014185078A1_000010t.gif

【0061】
また、事前分布において変数zjcは下記数式(9)を満たす、0もしくは1をとる二値変数とする。ガンマ分布を用いることでX線吸収係数が負にならない制約を表現できる。
【0062】
【数9】
JP2014185078A1_000011t.gif

【0063】
ここで、変数Z,Bの事前分布に関しては変数間の独立性を仮定せず、下記数式(10)で表されるように、ボルツマン分布に従うと仮定する。但し、Aは正規化項であり、エネルギーH(Z,B)は下記数式(11)のように定義する。下記数式(11)では、2次元平面の近傍として縦横に隣接する4画素を考えており、Gcdself,Gcdinterはそれぞれ非負の定数である。
【0064】
【数10】
JP2014185078A1_000012t.gif

【0065】
【数11】
JP2014185078A1_000013t.gif

【0066】
(事後分布について)
本発明のX線CT画像の再構成では、X線吸収係数X、物質クラスZ,変数Bは観測データDの事後分布として推定される。X線吸収係数X、物質クラスZ、変数Bの事後分布p(X,Z,B|D)はベイズの定理により下記数式(12)のように表すことができる。
【0067】
【数12】
JP2014185078A1_000014t.gif

【0068】
(ベイズ推定)
ベイズ推定では、上記数式(12)で表される事後分布p(X,Z,B|D)の期待値が必要になるが、この事後分布p(X,Z,B|D)を解析的に計算するのは困難であることに鑑みて、試験分布q(X,Z,B)により事後分布p(X,Z,B|D)を近似する。試験分布q(X,Z,B)は、上記数式(12)の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そして、試験分布q(X,Z,B)と事後分布p(X,Z,B|D)の誤差をKL(Kullback-Leibler)距離で評価し、KL距離を最小化する試験分布q(X,Z,B)を算定する。ここで、KL距離は下記数式(13)のように表すことができる。
但し、<・>q(X,Z,B)は分布q(X,Z,B)に関して積分計算を行うことを表す演算子である。KL距離は常に非負であり、q=pのときにのみ0(ゼロ)になる。<・>q(X,Z,B)はq(X,Z,B)に関して期待値をとることを表している。
【0069】
【数13】
JP2014185078A1_000015t.gif

【0070】
この最適化を容易にするため、試験分布q(X,Z,B)がその周辺分布q(X)とq(Z,B)の積で与えられる、すなわち変数Xと{Z,B}の間に独立性があることを仮定する。また、q(X),q(Z,B)の各要素について画素間で独立であると仮定し、さらにq(X)の各要素q(x)はガンマ分布に従うとする。
このときq(z、b)は、確率変数z, bが離散変数であるため必然的に多項分布となる。このことから、試験分布q(X,Z,B)は下記数式(14),(15)のように表すことができる。但し、試験分布q(X),q(Z,B)について同時に最適化するのは困難であるため、下記数式(16),(17)のように交互に最適化を行うものとする。
【0071】
【数14】
JP2014185078A1_000016t.gif

【0072】
【数15】
JP2014185078A1_000017t.gif

【0073】
【数16】
JP2014185078A1_000018t.gif

【0074】
【数17】
JP2014185078A1_000019t.gif

【0075】
ここで、q(x)のパラメータα,βの最適化にはSCG(Scaled Conjugate Gradient)法を用いた。試験分布q(X)とj番目の画素以外の q(z´,b´) (但し、j´はjと異なる) が固定のもとで、KL距離を最小化するq(z,b)は解析的に求めることができる。この解析解はjの近傍画素のq(z´,b´)(j´∈η(j))に依存する。そこで、近傍画素の分布を変更せずにq(Z,B)を更新するために、図3に示すように市松模様となる画素の部分集合を交互に最適化している。すなわち、パラメータq(Z,B)は全画素を同時に更新するのではなく、図3に示す市松模様のパターンで画素の部分集合を交互に最適化している。図3(a)のグレー部分の画素に該当するq(z,b)を更新したのち、図3(b)のグレー部分の画素に該当するq(z,b)を更新する。
【0076】
以上、本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおける処理内容について説明した。以下の実施例では、本発明の有効性を調べるために計算機実験を行い、本発明のX線CT画像処理方法およびX線CT画像処理プログラムの画面再構成を、従来技術による画面再構成と比較して評価する。
下記の評価実験から、本発明のX線CT画像処理方法およびX線CT画像処理プログラムの有用性が理解できるであろう。
【0077】
実施例では、再構成画像と真の画像と誤差を測れるように、適当な真の画像を用意し、それに対してX線投影をシミュレートして人工的に投影像(サイノグラム)を生成した。X線投影は、物質のX線吸収係数のエネルギー依存性(波長依存性)を考慮するために入射するX線のエネルギーを60keV,70keV,80keV,90keV,100keVの5通り(それぞれのフォトン数は、6×10,2.5×10,1.5×10,7×10,5×10)に分割し、それぞれで物質固有のX線吸収係数に従って減衰させた。検出器においてはその減衰したX線の光子にポアソンノイズに従って生成したノイズ(雑音)が加えられたものが観測されるとした。
【0078】
以下の実施例1,2において、検出器の個数は95(個)、画像の解像度は64×64(画素)、投影は1°から180°の間で1°間隔としている。再構成画像の推定の良さを定量的に評価するために、下記数式(18)で表されるPSNR(Peak Signal to Noise Ratio) を用いた。
【0079】
【数18】
JP2014185078A1_000020t.gif

【0080】
上記数式(18)において、MAXは画素値の最大値であり、MSEは真の画像と再構成画像との平均2乗誤差である。PSNRが高いほどMSEが小さい。すなわち、より真の画像に近いといえる。以下、FBP法、X線エネルギー分布を考慮しない事前分布をガウス分布としてベイズ推定により行う方法(以下、「従来方法」と略する)、X線エネルギーを3つに離散化した本発明のX線CT画像処理方法(以下、「本発明方法」と略する)の3通りで画像再構成を行い、それぞれ比較した結果を示す。
なお、金属と骨の事前分布は図4に示す通りである。図4は、金属と骨の事前分布と真の画像の60keVにおけるX線吸収係数を示している。金属と骨のX線吸収係数は、それぞれ0.7,0.4を平均とするガンマ分布で表現される確率分布である。真の画像のX線吸収係数は、それぞれのX線エネルギーにおいて、事前分布の平均値より5%小さい値であるとしている。
【実施例1】
【0081】
実施例1では、金属を4つ含むファントムに対する画像再構成についての計算機実験について図5を参照して説明する。
図5(a)は真の画像、(b)はFBP法による再構成画像、(c)は単純ベイズ法による再構成画像、(d)は本発明方法による再構成画像である。
実施例1では、図5(a)に示す真の画像に対して、FBP法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射X線の光子数の合計を1.12×10とし、本発明のX線CT画像処理方法において、X線のエネルギーを60keV,80keV,100keVの3つに離散化し、それぞれに属する光子数を7.5×10,3×10,10.7×10であるとした。繰り返し回数は比較のため共に25回と設定した。推定結果を、図5(b)~(d)に示す。
【0082】
また、FBP法、単純ベイズ法、本発明方法のそれぞれにおけるPSNRを下記表1に示す。PSNRでは、FBP法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。
【0083】
【表1】
JP2014185078A1_000021t.gif

【実施例2】
【0084】
実施例2では、金属インプラントを含む歯のファントムに対する画像再構成についての計算機実験について図6を参照して説明する。
図6(a)は真の画像、(b)はFBP法による再構成画像、(c)は単純ベイズ法による再構成画像、(d)は本発明方法による再構成画像である。
実施例2では、図6(a)に示す真の画像に対して、FBP法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射X線の光子数の合計を1.12×10とし、本発明のX線CT画像処理方法において、X線のエネルギーを60keV,80keV,100keVの3つに離散化し、それぞれに属する光子数を7.5×10,3×10,10.7×10であるとした。繰り返し回数は比較のため共に25回と設定した。推定結果を、図6(b)~(d)に示す。
【0085】
FBP法、単純ベイズ法、本発明方法のそれぞれにおけるPSNRを下記表2に示す。実施例1の実験結果と同様に、PSNRでは、FBP法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。
【0086】
【表2】
JP2014185078A1_000022t.gif

【0087】
以上のように、X線スペクトルのエネルギー分布を考慮し、エネルギー依存性(波長依存性)のあるX線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングの影響による様々なアーティファクトを軽減できることを示した。実施例1,2の計算機実験としてはX線エネルギーを3つに分割した場合のみを行ったが、このように実際のエネルギーの分割数より少ない数で離散表現しても、アーティファクトが大きく軽減できることがわかる。さらにエネルギーの分割数を増やすことでより良い推定を行うことが可能であろう。
【実施例3】
【0088】
次に、投影X線の観測に関する観測過程について、多重X線和が、検出器におけるX線の検出面が空間的に有限の広がりを有すること、検出器におけるX線の検出時間が有限の幅を有すること、X線が有限の広がりを持つ波長域(スペクトル分布)を有すること、によって、空間的な位置と方向、波長が異なる複数のX線が検出器に入射することを表現できることを説明する。
【0089】
まず、平行ビームを例にとり、図7を参照して説明する。図7のように2次元平面で平行ビームによる投影を想定する。ここで、X線の放出される方向をv軸、それに直交する方向をu軸ととる。この座標系で以下のように変数を定義する。
位置(u,v)における物質の波長λのX線の吸収係数をx(u,v,λ)とする。但し、物質は、v∈Lに含まれるとする。
i番目の検出器の検出面Sに、単位時間あたり入射する波長λのフォトン数をy(λ)とする。平行ビームのX線源が、単位時間、単位面積あたり放出する波長λのフォトン数をd(λ)、波長λのフォトン1個が励起する電子の平均的な数をρ(λ)とする。
このとき単位時間にi番目の検出器で観測される電流Yの期待値は、下記数式(19)のように表される。さらに、ある一定時間の電流を貯めてそれを観測する場合は、時間に関する積分を行うことになる。
【0090】
【数19】
JP2014185078A1_000023t.gif

【0091】
上記数式(19)は逆投影法の仮定する線積分(上記のvに関する積分)に帰着できないことから、検出面の空間的な広がりや、X線のスペクトル分布の拡がりが存在する場合は線積分を仮定する推定手法を用いるとアーティファクトが生じることがわかる。同様に、検出器が受光している時間における回転を考慮する場合も線積分に帰着できない。上記の多重積分の形では計算が難しいので積分を和で近似することを考える。これが多重X線和を考えることになることを示す。まず、上記数式(19)において特定の波長λにのみ注目し、下記数式(20)におけるu, vに関する積分を和に置き換える。
【0092】
【数20】
JP2014185078A1_000024t.gif

【0093】
計算機上で上記のu, vに関する積分を近似的に表現するために離散化を行う(図8を参照)。 空間をu軸、v軸にそって離散化すると、上記数式(20)におけるvに関する積分は下記数式(21)のように近似できる。
【0094】
【数21】
JP2014185078A1_000025t.gif

【0095】
ここで、xはj番目の画素のX線吸収係数であり、I(u)はu軸上のuの位置にある画素のインデックスの集合、ΔLは画素の1辺の長さを表す。画素の1辺の長さをΔL=1とすると、上記数式(21)は下記数式(22)として表現できる。
【0096】
【数22】
JP2014185078A1_000026t.gif

【0097】
したがって、画素の1辺の長さが1のときは、yは下記数式(23)となり、さらに同じ検出器に入る光線に関する積分(uに関する積分)を和で近似すると、下記数式(24)のように表現できる。
【0098】
【数23】
JP2014185078A1_000027t.gif

【0099】
【数24】
JP2014185078A1_000028t.gif

【0100】
上記数式24において、1/Nは積分を短冊で切って近似するときの短冊の幅を示している。
これはインデックスkごとに線積分で表現される投影の和を考えていることに相当する。
以上から、検出器の空間的拡がりを考慮する場合においては多重のX線を考慮する必要があることがわかる。
【0101】
上記では、X線が画素に対してまっすぐ入射する場合を想定している。しかしながら、一般のX線CT装置などのように、X線源と検出器が対向しながら回転するような装置(図9を参照)では、画素に対して斜めに入射する状況を想定する必要がある。
この場合、投影についての下記数式(24)を、下記数式(25)のように一般化することで表現できる。
【0102】
【数25】
JP2014185078A1_000029t.gif

【0103】
【数26】
JP2014185078A1_000030t.gif

【0104】
ここで、ΔN=1/Nであり、r(k)は光線kの検出平面上の位置を表し、{k|r(k)∈Sが検出器iに到達する光線の集合を表す。Iは光線kが交差する画素のインデックスの集合を表し、lkjは光線kが画素jと交差する長さを表している。上記数式(25)は上記数式(26)のlkjを1とした場合と考えることができる。
上記数式(26)は、平行ビームのほか、ファンビーム、コーンビームでも光線の向きを考えて、lkjを計算することで同様の式が成立する。
上記数式(26)からわかることは、検出面の拡がりを考えた投影の観測モデルを考えた時には、検出面内を通る光線ごとの投影を考える必要がある点である。
【0105】
上記数式(26)の中のlkjの和を、検出器iや光線kに関してまとめて表記するならば、Lxのような行列演算で表現できる。検出面の拡がりを考えなかった場合には、行列Lは(検出器の数)×(再構成画像の画素数)であったものが、検出面の拡がりを考えた場合には、行列Lは(検出器を通る光線の数)×(再構成画像の画素数)となる。光線の数は、多ければ積分を精度よく近似できる一方で、計算量が増す問題がある。
【0106】
次に、投影中の回転(時間に関する積分)を想定する。まず時刻tにおいて観測されたフォトン数をyi,tとすると、下記数式(27)で表現できる。
【0107】
【数27】
JP2014185078A1_000031t.gif

【0108】
検出器の検出時間(露光時間)の間に計測されたフォトン数は、yi,tを時刻tに関して積分したものとなるので、下記数式(28)で表現できる。
【0109】
【数28】
JP2014185078A1_000032t.gif

【0110】
上記数式(28)について、時刻に関する積分を和で置き換えると、下記数式(29)のように表現できる。
【0111】
【数29】
JP2014185078A1_000033t.gif

【0112】
前述のように、検出面の拡がりを考慮すると上記数式(28)のように各時刻tで検出器に入射するフォトンyi,tが複数の光線束の和で表現される。時間に関する積分は、さらなる和を含めるため新たな光線束に関する和を考えることに等価である。そこで、光線インデックスの集合Sを時刻に依存するようSi,tとし、それぞれの光線kが通る画素インデックスの集合Iを改めて定義しなおすと、検出時間(露光時間)Tで計測されたフォトン数は、下記数式(30)で表現できる。
【0113】
【数30】
JP2014185078A1_000034t.gif

【0114】
上記数式(30)を波長に関する積分を和で置き換えると、下記数式(31)で表現できる。
【0115】
【数31】
JP2014185078A1_000035t.gif

【0116】
以上をまとめると、検出面の拡がりと投影中の回転,X線の波長依存性を考慮した場合、下記数式(32)で表現できる。ここで、b=dρである。
【0117】
【数32】
JP2014185078A1_000036t.gif

【0118】
上記数式(32)のように、mやkに関する和が指数関数の外側でとられることを多重X線和と呼ぶ。これは、指数関数の外側に和(もしくは積分)があることを想定しない線積分では表現できない。このような多重X線和のある観測過程一般に対して、本明細書で開示したベイズ推定に基づくアルゴリズムを適用することができる。
【産業上の利用可能性】
【0119】
本発明は、ビームハードニングアーティファクトを低減するX線CT装置として、医療用のみならず産業用のX線CT装置に有用である。また、ベルトコンベアによる移動や回転が無視できないベルトコンベア型X線CT装置や歯科用X線CT装置としても有用である。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8