(57)【要約】  Ｘ線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減可能なＸ線ＣＴ画像処理方法を提供する。投影Ｘ線の観測に関する確率分布を設定し統計推定を行う。確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現される。多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有すること、検出器におけるＸ線の検出時間が一瞬ではなく所定の幅を有すること、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有することによって、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射するとして表現される。そして、多重Ｘ線和を表現する観測過程と物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、統計推定が事後分布の期待値によるベイズ推定を行う。
2014185078.tif 000002
【特許請求の範囲】
【請求項１】
Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理方法であって、
Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記統計推定が前記Ｘ線吸収係数及び前記物質に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定するステップと、
を備え、
前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、
前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、
ことを特徴とするＸ線ＣＴ画像処理方法。
【請求項２】
前記多重Ｘ線和に、Ｘ線波長に関する和を含める場合、
各画素の波長依存性を有するＸ線吸収係数に関する事前分布は、
前記物質毎およびＸ線波長毎のＸ線吸収係数の分布を表現するパラメータを用いて表現される、
ことを特徴とする請求項１に記載のＸ線ＣＴ画像処理方法。
【請求項３】
前記多重Ｘ線和に、Ｘ線波長に関する和を含める場合、
各画素の波長依存性を有するＸ線吸収係数は、
波長に依存しない各画素のＸ線吸収係数密度と、画素に依存せず前記物質毎に定められるＸ線波長間のＸ線吸収係数比の積として表現され、
前記物質毎に定められる、Ｘ線波長間のＸ線吸収係数比を表現するパラメータを用いて表現され、
各画素の波長依存性を有するＸ線吸収係数の事前分布が前記物質を条件とした各画素の前記Ｘ線吸収係数密度に対する条件付き事前分布と前記物質に対する事前分布によって表現される、
ことを特徴とする請求項１に記載のＸ線ＣＴ画像処理方法。
【請求項４】
前記物質に対する事前分布が、それぞれの物質の出現する度合いを示すパラメータとそれぞれの物質が空間的に連続する度合いを示すパラメータとで表現される事前分布で設定されることを特徴とする請求項３に記載のＸ線ＣＴ画像処理方法。
【請求項５】
前記物質毎に定められている前記Ｘ線吸収係数密度の条件付き事前分布が各物質のサブクラス毎に定められている前記Ｘ線吸収係数密度の条件付き事前分布と各物質のサブクラスの事前分布の積を各物質のとりえる全てのサブクラスに関して和をとった混合確率分布で表現され、前記各物質のサブクラスの事前分布が、各物質のサブクラスの出現する度合いを示すパラメータと前記各物質のサブクラスが空間的に連続する度合いを示すパラメータとで設定されることを特徴とする請求項３に記載のＸ線ＣＴ画像処理方法。
【請求項６】
Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理装置であって、
Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力する入力手段と、
前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記Ｘ線吸収係数及び前記物質に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定する推定手段と、
を備え、
前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として設定され、
前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、ことを特徴とするＸ線ＣＴ画像処理装置。
【請求項７】
Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理プログラムであって、
コンピュータに、
Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
前記Ｘ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、前記統計推定が前記Ｘ線吸収係数及び前記物質(及び物質のサブクラス)に関する事後分布の期待値によって前記Ｘ線吸収係数及び前記物質を推定するステップと、
を実行させ、
前記Ｘ線投影画像の観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、
前記多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである、ことを特徴とするＸ線ＣＴ画像処理プログラム。
【請求項８】
請求項７に記載のＸ線ＣＴ画像処理プログラムが搭載されたＸ線ＣＴ画像装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、体内組織構成を考慮したＸ線ＣＴ（Computed Tomography）画像処理方法、Ｘ線ＣＴ画像処理プログラム及びＸ線ＣＴ画像装置に関する。
【背景技術】
【０００２】
Ｘ線ＣＴ（Computed Tomography）は、観測対象にＸ線を照射することで得られる投影データから観測対象内部をＸ線吸収係数の分布として画像化する技術である。
一般に、Ｘ線ＣＴは、Ｘ線管にＸ線検出器を対向配置するとともに、これらの間にターンテーブルを配置した構成を取る。このターンテーブルの回転軸は、Ｘ線管とＸ線検出器を結ぶＸ線光軸に対して直交する向きとされる。Ｘ線ＣＴ装置においては、測定対象物に対して複数の角度からの透過Ｘ線データを収集してＣＴ断層画像の再構成を行う。一般的には、１８０°以上、回転させて透過Ｘ線データを収集する。【０００３】
従来から、より低いＸ線被曝量で従来と同程度のＣＴ画像を再構成する技術、メタルアーティファクトの低減技術など、Ｘ線ＣＴ再構成の高精度化を可能にするアルゴリズムが研究されている。上記のメタルアーティファクトは、測定対象物にＸ線を強く吸収する金属などの密度の高い高吸収物質があったときに、正確な再構成ができず再構成画像に人為的なノイズがのるというものである。Ｘ線ＣＴ再構成は、様々な投影角度において得られる投影像から、測定対象物を透過したＸ線の減衰の度合いを求め、内部の物質のＸ線吸収係数の分布を推定する手法である。しかし、高吸収物質はＸ線を強く減衰させるため、ほとんど透過Ｘ線が検出されず観測ノイズに紛れてしまう。特に既存手法であるフィルター補正逆投影法（ＦＢＰ：Filtered Back Projection）を用いた場合、ＦＢＰ法ではノイズの影響を受けやすいことから、再構成画像にはメタルアーティファクトと呼ばれる虚像が生じ、正確な画像の再構成ができないという問題がある。このようにＦＢＰ法はノイズの影響を受けやすいのであるが、主なノイズにはポアソン分布に従うショットノイズが含まれるため照射Ｘ線の強度を高くすることで信号雑音比を向上することができる。これによって相対的なノイズの影響を小さくし、画像再構成の精度を向上させることができるが、その一方で放射線被曝量が増大させてしまい、発癌などの確率的健康リスクを生じさせてしまう。【０００４】
このように、Ｘ線被曝量の低減とＳＮ比はトレードオフの関係があり、より低いＸ線被曝量で従来と同程度のＳＮ比の再構成画像を得る技術として、統計推定法が提案されている。統計推定法では、観測過程が確率的なノイズによる不確実性を含むことや、撮像対象の物体に関する事前知識を利用することで、低被曝環境下における推定の不良設定性を軽減することができる。 結果として、従来と同じＸ線強度や解像度でも再構成画像に含まれるアーティファクトを低減できる。【０００５】
現在、画像再構成法の主流であるフィルター補正逆投影法（ＦＢＰ：Filtered Back Projection）や、ＦＢＰを改良したＰＣＬＩＳ（Projection Completion Method based on a Linear Interpolation in the Sinogram）は、統計推定を行わないものである。
統計推定法は、最尤推定法（Maximum Likelihood Estimation: ＭＬＥ），ＭＡＰ（Maximum A Posteriori）推定法，ベイズ推定法の３手法に大別される。
最尤推定法（ＭＬＥ)は、推定対象である再構成画像に関する事前知識をおかない統計推定法で、観測における不確実性を含む物理過程を確率モデルで表現して、その確率モデルを基に最も尤もらしい再構成画像を推定するものである。観測における不確実性を含む物理過程には、観察される透過Ｘ線のフォトン数に関する揺らぎ(ショットノイズ)などが含まれる。
ＭＡＰ推定法，ベイズ推定法は、推定対象に関する事前知識を確率モデルの形式で表現し、推定に取り入れるものである。すなわち、ＭＡＰ推定法やベイズ推定法は、観測に関する確率モデルと推定対象に関する確率モデルの両方を勘案して、推定が行われる。 Ｘ線CTにおいては、推定対象はＸ線吸収係数である。具体的には、ＭＡＰ推定法では、観測に関する確率モデルと推定対象となる物体のＸ線吸収係数に関する確率モデルから計算されるＸ線吸収係数の事後分布の最大値を推定値とし、ベイズ推定法では同じＸ線吸収係数の事後分布の期待値を推定値とする。ＭＡＰ推定法では、事後分布が最大となる一点のみを取り上げて推定を行うが、ベイズ推定法では事後分布の期待値をとるため、事後分布が大きくなるＸ線吸収係数の集団に対する平均によって推定を行う。 事後分布が最大となる点は、観測ノイズによって敏感に変化しうるが、事後分布の期待値は観測ノイズに対してそれほど敏感に影響されないため、ベイズ推定法のほうがより安定し、精度の高い推定となる。【０００６】
ここで、通常のＸ線ＣＴアルゴリズムの多くは画像再構成にあたって、照射するＸ線が線スペクトルであると仮定し単一の吸収係数を推定している。しかしながら、実際のＸ線は単一エネルギーではなく、Ｘ線管に加えられた最高管電圧で決まるエネルギーを最大として低い方へ連続的に分布するエネルギースペクトルで構成されている。そのため、物質にＸ線を照射するとエネルギーの高いＸ線はほとんど吸収されずに透過するが、エネルギーの低いＸ線は減弱が大きいため、物質を透過した後にはＸ線のスペクトル分布は高エネルギー側にシフトするビームハードニングと呼ばれる現象が生じる。
物質のＸ線吸収係数のエネルギー依存性（波長依存性）の影響によって生じる現象は、Ｘ線源が線スペクトルであると仮定し、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しない観測過程のモデルでは表現できず、正しくない推定が行われてしまう。このように物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルによる推定でビームハードニングアーティファクトが生じる。【０００７】
特に、金属などの極端にＸ線吸収係数が高い高吸収物体があるとストリークアーティファクトと呼ばれる筋状のアーティファクトや高吸収物体周辺のＸ線吸収係数が実際より小さく見積もられるダークバンドアーティファクトが生じることが知られているが、これはビームハードニングの影響を考慮していないことが大きく影響していると考える。
これについて、以下に説明する。
金属など非常に高いＸ線吸収係数をもつ高吸収物体を透過すると、Ｘ線は大きく減弱する。その性質が観測モデルに表現されていれば、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルであっても、高吸収物体近辺に大きなＸ線吸収係数をもつ物体が存在するであろうことを推定できる。しかしながら、高吸収物体でも高エネルギーのＸ線は透過しやすい。そのため、実際にはＸ線源の主要なエネルギーをもつフォトンの多くが高吸収物体を透過中に吸収されていたとしても多少のＸ線が観測される。高吸収物体が離れた位置に複数ある場合、その複数の高吸収物体のいずれも透過するＸ線は、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルでは高いＸ線吸収係数をもつ物体が重なると透過Ｘ線がさらに大きく減衰することを予想する。 しかし、実際にはある一定量の高エネルギーのＸ線は高吸収物体を複数透過したとしても観測されてしまう。このため、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮しないモデルで、単純に高吸収物体が存在すると仮定すると実際の観測と矛盾が生じてしまう。このようにして、物質のＸ線吸収係数のエネルギー依存性を考慮しない観測モデルと実際の観測とができるだけ矛盾しないようにＸ線吸収係数の推定を行うと、高吸収物体周辺のＸ線吸収係数が実際より小さく見積もられてしまうダークバンドアーティファクトが生じる。【０００８】
メタルアーティファクトを軽減する方法として、投影データを金属領域、非金属領域に分離し、それぞれの領域において画像再構成し、得た画像を重ね合わせるという方法が既に知られている。このようにすることで、金属領域における投影データの矛盾を正し、メタルアーティファクトの発生を軽減することができる。また、本願発明者らは、再構成画像となるＸ線吸収係数に関する事前知識をおいて統計推定を行うＸ線ＣＴ画像処理方法であって、事前知識を、再構成画像の各画素の領域において定義されるパラメータで、撮像対象となる人体の各組織の存在割合を表現するパラメータと、各組織のＸ線吸収係数を表現するパラメータと、各組織の空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現するものを用いることで、メタルアーティファクトを低減可能なＸ線ＣＴ画像処理方法を提案している（特許文献１を参照）。
しかし、これらの方法は、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮するわけではないため、実際の観測過程と物質のＸ線吸収係数のエネルギー依存性を考慮しない観測過程とのズレのために生じるビームハードニングアーティファクトを除去するには至っていない。【０００９】
そして、ビームハードニングアーティファクトを軽減する手法として、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮したモデルが知られている（非特許文献１〜４）。Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮するモデルでは、Ｘ線のエネルギー毎のＸ線吸収係数を推定する必要がある。これは、より正確な物理過程に基づくが、推定問題としては不良設定性の問題を生じさせる。たとえば、単純のためＸ線のエネルギーを離散化し、二種類のエネルギーにおけるＸ線吸収係数を推定する問題を考えたとする。このとき、観測される投影像は従来と同じであっても、推定しなければならないＸ線吸収係数は従来に比べて２倍に増えてしまう。このように解となるＸ線吸収係数の自由度が高まるため、実際の観測とそれほど矛盾を生じない解がいくつも生じてしまう。この場合、ちょっとした観測ノイズに推定結果が敏感に影響され安定しない。 そこで、解となるＸ線吸収係数の自由度を抑える制約が必要となる。 先ず、非特許文献１〜３では、各画素におけるエネルギー依存の吸収係数を既知の物質の吸収係数関数の重み付き和で表現することで、物質のエネルギー依存性（波長依存性）を表現している。Ｘ線吸収係数のエネルギー依存性（波長依存性）は、物質毎に異なるが高エネルギーになるに従って減衰する傾向など各物質間で共通する性質が存在する。Ｘ線のエネルギーを入力変数とし、そのエネルギーでのＸ線吸収係数を出力する関数を吸収係数関数と呼ぶことにすると、異なる物質であっても類似した吸収係数関数をもつため、どの物質の吸収係数関数も少数の、典型的には二つの物質の吸収係数関数の重み付き和（基底関数表現）でよく表現できる。例えば、非特許文献２では、物質クラスの決定をＦＢＰで得られた画像をもとに行い、その知識をもとに画像再構成を行っている。非特許文献４では、画素毎のＸ線吸収係数は、既知の物質種毎のＸ線吸収係数の組み合わせで表現できるという仮定で推定を行っている。
このように、Ｘ線吸収係数のエネルギー依存性（波長依存性）を考慮する方法は、その不良設定性を抑制するため、既知の物質のＸ線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らすことが行われてきた。しかしながら、その一方で、強すぎる制約が課せられる問題が指摘されている。【００１０】
一方、ハードウェア的には、ＤＥＣＴ（Dual Energy Computed Tomography）の登場により、異なるＸ線管電圧で撮影し、Ｘ線吸収係数の違いから物質のエネルギー依存性（波長依存性）を推定できるようになった。このＤＥＣＴを用いた様々な画像再構成法が知られている（非特許文献５〜８）。しかし、ＤＥＣＴを用いた画像再構成法は、一般的に単色Ｘ線ＣＴより被曝量が大きく、同程度の被曝量に抑えようとすると推定精度が悪くなってしまうなどの問題がある。また、非特許文献５は、事後的に画像処理によって補正する方法であり、その補正によって更なる虚像を生じさせる可能性もあり限界がある。【００１１】
物質のＸ線吸収係数のエネルギー依存性を考慮しない観測モデルは、検出器に入射するＸ線が複数の波長からなるのではなく、ある代表的な波長のみからなることを想定することになる。 これは、検出器に複数のＸ線スペクトルが多重に入射する過程を、単純化したものと考えることができる。このようにＸ線が多重に検出器に入射する現象は、スペクトル空間だけではなく、時間的、空間的にも生じる。すなわち、検出器はある一定の検出時間内に多重のＸ線が入射し、ある一定の面積の検出面内に多重のＸ線が入射する。これら投影したＸ線が多重に検出器に入射する観測過程において、検出器に入射するＸ線を一本の光線で代表させた場合、Ｘ線吸収係数のエネルギー依存性を考慮しない時にビームハードニングアーティファクトが生じるのと同様にＸ線吸収係数の推定にバイアスが生じることが知られている（非特許文献９）。つまり、投影したＸ線が多重に検出器に入射する観測過程においては、一本の光線ではなく（Ｘ線の多重性を無視しないで）、複数の光線を想定した多重のＸ線の和と捉えることが望ましい。
【先行技術文献】
【特許文献】
【００１２】
【特許文献１】特開２０１１−１５６３０２号公報
【非特許文献】
【００１３】
【非特許文献１】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
【非特許文献２】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
【非特許文献３】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
【非特許文献４】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)
【非特許文献５】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
【非特許文献６】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
【非特許文献７】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
【非特許文献８】山崎陽一, 戸田尚宏, 低被爆環境下におけるDual-Energy Ｘ線CT雑音低減アルゴリズム, IEICE Technical Report MBE2008-113, (2009), 101-106
【非特許文献９】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
【発明の概要】
【発明が解決しようとする課題】
【００１４】
上述の如く、波長が短くエネルギーの強い硬Ｘ線は、波長が長くエネルギーの弱い軟Ｘ線が物体を透過する力に比べて、物体を透過しやすいため、物体に入射したＸ線のスペクトル分布は物体を通過後、軟Ｘ線が減衰し硬Ｘ線側にパワースペクトルのピークが偏ったものとなり、ビームハードニングと呼ばれる現象を生じている。
理想的には、投影によって得られるＸ線のスペクトル分布を分光法やフィルターを用いて推定し、各スペクトル強度のＸ線を測定し、その測定をもとに各スペクトルのＸ線吸収係数の推定を行えば問題は解決するのであるが、実際に分光すると各スペクトルで得られる信号のＳＮ比が下がるため推定精度は悪化するという問題がある。 そのため従来から、情報処理アルゴリズムでビームハードニングアーティファクトの除去を行うことが要望されている。【００１５】
様々なビームハードニングアーティファクトの除去を行う情報処理アルゴリズムが提案されている中で、物質のＸ線吸収係数の波長依存性を考慮する方法があり、特に、組織のＸ線吸収係数の波長依存性を利用した基底関数表現を用いることで、推定するパラメータの自由度を減らす技術がある（非特許文献１，２を参照）。【００１６】
しかし、非特許文献１，２では、組織の種類を離散変数で表現した組織クラスを決定論的に推定するために、想定する組織クラス数が適切でない場合や、推定された組織クラスが誤っていた場合、解が大きく異なってしまい、既知の組織クラスで表現されない組織に対して適切な表現ができないといった問題がある。
また、推定する組織クラスに関して制約がおかれなかったため、不良設定性の解消が十分でなく、組織クラスの推定が不安定になりがちであるといった問題がある。
特に、非特許文献１では、波長間でＸ線吸収係数の比が可変であり、異なる組織クラスを重み付けする。組織クラス数が少ないと表現できないＸ線吸収係数が現れる可能性があり、組織クラス数が多いと拘束の意味をなさずに不良設定性が生じてしまう。【００１７】
そこで、本発明は、従来技術を用いた画像再構成と比べて、物質の種類を表す離散変数である物質クラスを介してＸ線吸収係数を「確率的」に表現し、物質クラスに確率的な拘束を設けることで不良設定性を抑制する一方で、各画素がどの物質クラスに所属するかを決定論的ではなく確率的に推定することでより柔軟なＸ線吸収係数の表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、アーティファクトを低減可能なＸ線ＣＴ画像処理方法、Ｘ線ＣＴ画像処理プログラムおよび該プログラムが搭載されたＸ線ＣＴ画像処理装置を提供することを目的とする。  ここで、本発明において物質クラスは、必ずしも生体の組織とは限らないので組織クラスを包含する意味で用いている。【００１８】
また、検出器に入射するＸ線は、一本の光線ではなく、複数の光線を想定した多重のＸ線の和と捉えることが必要である。上述の如く、このＸ線の多重性を無視して、一本の光線で代表させることによりバイアスが生じる。
そこで、本発明は、投影したＸ線が多重に検出器に入射する観測に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減可能なＸ線ＣＴ画像処理方法、Ｘ線ＣＴ画像処理プログラムおよび該プログラムが搭載されたＸ線ＣＴ画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【００１９】
上記状況に鑑みて、本発明のＸ線ＣＴ画像処理方法は、Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行う。
上記の確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現する。ここで、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。
そして、多重Ｘ線和を表現する観測過程と物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、統計推定がＸ線吸収係数及び物質に関する事後分布の期待値によって線吸収係数及び物質を推定する。【００２０】
ここで、多重Ｘ線和で表現される観測過程は、以下の（ａ），（ｂ）の特徴を備える。
(ａ) 各画素の波長依存性を有するＸ線吸収係数は、各画素のＸ線吸収係数密度と、物質毎に定められているＸ線波長間のＸ線吸収係数比の積で表現される。Ｘ線吸収係数の波長依存性を考えない場合のＸ線吸収係数の推定は、波長依存性を有するＸ線吸収係数において想定するＸ線波長域が１つしかないと考えてＸ線吸収係数密度のみを推定することと等価となる。そのため、波長依存性を有するＸ線吸収係数の各画素のＸ線吸収係数密度と、物質毎に定められているＸ線波長間のＸ線吸収係数比の積による表現は、Ｘ線吸収係数の波長依存性を考えない場合のＸ線吸収係数の表現の拡張となっている。
(ｂ)投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現するものであり、
多重Ｘ線和は、
検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有すること、
検出器におけるＸ線の検出時間が一瞬ではなく所定の幅を有すること、
Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有すること、
によって、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射することを想定する。そして、想定した複数の光線のそれぞれで投影を計算し、それらの投影によって励起された（電気）信号の和が検出器で観測されるとする。【００２１】
また、撮像対象がどのようなＸ線吸収係数をとる傾向があるかを表すＸ線吸収係数に関する事前分布は、物質クラスという変数を介して決められる。 物質クラスは、撮像対象の内部がどの物質から構成されるかを表す撮像対象内部の各画素で定義される変数である。物質クラスは、どの物質クラスがどの程度生じやすいか、またそれぞれの物質クラスの空間的滑らかさの度合いについて確率的な制約を加えることにより波長依存のＸ線吸収係数の自由度を確率的に制約する。 
また、物質クラスが定まった時に、どのようなＸ線吸収係数密度をとりやすいかという確率分布を２つのガンマ分布の和で表現される場合には、各ガンマ分布のいずれに所属するかを表現する物質クラスのサブクラス変数を導入し、物質クラスとサブクラスの組み合わせに対してどの物質クラスとサブクラスの組み合わせがどの程度生じやすいか、またそれぞれの物質クラスとサブクラスの組み合わせがどの程度、空間的に滑らかさであるか、の確率的な制約を加えることにより波長依存のＸ線吸収係数の自由度を確率的に制約する。【００２２】
本発明のＸ線ＣＴ画像処理方法によれば、ビームハードニングアーティファクトをはじめとする多重Ｘ線和の観測過程を無視することで生じるアーティファクトを除去するために、物質のＸ線吸収係数の波長依存性など多重Ｘ線和の観測過程を考慮した確率分布を設定しＸ線ＣＴ画像処理を行う。ここで、物質のＸ線吸収係数の波長依存性を考慮する際には、波長ごとのＸ線吸収係数を全て未知変数として推定すると、問題の不良設定性のために信頼のおける解が得られないことから、いくつかの既知の物質が特定の波長依存性をもったＸ線吸収係数を有することを表現する事前分布を利用することにし、再構成画像の各画素が既知の物質クラスのうちのいずれの物質クラスから構成されるかを推定するという、より易しい推定問題に変更することで、より信頼のおける解を得る。すなわち、波長ごとに正の実数変数であるＸ線吸収係数を推定する代わりに、波長に依存しない正の実数変数であるＸ線吸収係数密度と離散変数である物質クラスを推定する。この各画素の物質クラスは直接的に観測できない隠れ状態変数、すなわち物質クラスの変数を隠れた確率変数として扱う。【００２３】
この隠れ状態変数の物質クラスに対して、空間的滑らかさの確率的な制約を課すことで、波長依存のＸ線吸収係数の自由度を確率的に制約することとした。その上で隠れた変数をもつ確率分布の推論として、ベイズ推論を近似的に行うことで安定した精度の高い推定（収束が保証され、観測ノイズに対してあまり敏感に反応しない安定した推定）を行うのである。また、本発明におけるベイズ推定によれば、Ｘ線吸収係数と再構成画像の各画素の領域が属する物質クラスの両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。ベイズ推定により、各画素における物質クラスの所属確率と各画素のＸ線吸収係数密度を確率的に表現できるので、決定論的に物質クラスを選択する場合に比べてより柔軟性をもったきめ細やかなＸ線吸収係数の表現ができる。【００２４】
本発明のＸ線ＣＴ画像処理方法では、物質クラスの変数を確率変数として扱うことで、既知の物質クラスで表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、本発明のＸ線ＣＴ画像処理方法では、物質クラスに対して事前知識として空間的滑らかさなどの制約を課すことで不良設定性を抑制することができる。 また、本発明のＸ線ＣＴ画像処理方法では、近似ベイズ推論を行うことで、収束性が保証され、かつ、推定の不安定性(オーバーフィッティング)を抑えることができる。
すなわち、本発明のＸ線ＣＴ画像処理方法を用いることにより、物理過程により忠実でかつ精度の高い安定した推定が行えるのである。【００２５】
本発明において、物質クラスは、組織クラスを包含する意味で用いている。物質クラスは、離散変数で表現される。物質クラスの存在割合（密度）とは、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス（筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど）がどの程度の割合で分布するのを示すもので、非負のスカラーのパラメータである。
物質クラス毎に波長間で固定されたＸ線吸収係数比を表現するパラメータとは、確率的に表現されたＸ線吸収係数において、波長間で吸収係数の比は固定と仮定し、それを事前知識の確率分布に用いて、波長依存のＸ線吸収係数の自由度を制約するものである。
また、物質クラスが空間的に連続する度合いを表現するパラメータは、例えば、空気の物質クラスは画素同士の連結がしやすい（空気中に他の組織は入らない）など空間的に連続する度合いが高いという制約を課すことができるので、物質クラスが空間的に連続する度合いを事前知識の確率分布に用いて、波長依存のＸ線吸収係数の自由度を制約するものである。【００２６】
上記のＸ線ＣＴ画像処理方法によれば、人体などの撮像対象の物質クラス分布に関する事前知識について、どのような物質クラス（筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど）がどの程度の割合で分布するか、また、それぞれの物質クラスがどの程度のＸ線吸収係数を有するか、また、それぞれの物質クラスがどの程度、空間的に連続して分布しやすいかの分布情報を確率分布の形で表現して、事後確率の期待値によるベイズ推定を用いて、画像再構成や物質推定を行う。
なお、物質クラスに関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と物質クラスの推定精度のさらなる向上を図ることが可能である。【００２７】
本発明のＸ線ＣＴ画像処理方法における多重Ｘ線和にＸ線波長に関する和を含める場合を想定する。
この場合、各画素の波長依存性を有するＸ線吸収係数に対する事前分布を、物質クラス毎およびＸ線波長毎のＸ線吸収係数の分布を表現するパラメータを介して特徴づけることによって、どのようなＸ線吸収係数が生じやすいかを確率的に限定するものである。【００２８】
他の観点から、本発明のＸ線ＣＴ画像処理方法における多重Ｘ線和にＸ線波長に関する和を含める場合、各画素の波長依存性を有するＸ線吸収係数は、物質クラスに依存しない各画素のＸ線吸収係数密度と、物質クラス毎に定められているＸ線波長間のＸ線吸収係数比の積で表現される。
そして、各画素の波長依存性を有するＸ線吸収係数は、物質クラス毎に定められている、Ｘ線波長間のＸ線吸収係数比を表現するパラメータで特徴付けられる。
そして、各画素の波長依存性を有するＸ線吸収係数の事前分布が物質クラスを条件とした各画素のＸ線吸収係数密度に対する条件付き事前分布と物質クラスに対する事前分布によって表現される。
すなわち、Ｘ線吸収係数を各画素のＸ線吸収係数密度と物質クラス依存波長依存の定数の積で表現する場合、Ｘ線吸収係数の事前分布は、各画素のＸ線吸収係数密度に対する事前分布と物質クラスに対する事前分布から一意に決定されることになる。【００２９】
ここで、本発明のＸ線ＣＴ画像処理方法において、上記の物質クラス毎に定められているＸ線吸収係数密度の分布が、混合ガンマ分布で表現される。また、物質クラスが、出現する度合いを表現するパラメータと物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現される。
これは、Ｘ線吸収係数の事前分布の表現についての一例である。
物質クラスを条件としたＸ線吸収係数密度に対して混合ガンマ分布を仮定し、物質クラス変数に対してボルツマン分布を仮定することによって、Ｘ線吸収係数の事前分布が決定されるようにする。【００３０】
物質クラス毎に定められているＸ線吸収係数密度の分布として（混合）ガンマ分布で表現されるとした理由は、その分布形状が単峰な分布を表現でき、かつその確率変数のとる値が必ず正の値となることを保証できることに着目したものである。また、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程において検出器において観測される信号の分布としてガウス分布で表現されるとした理由は、光子散乱や量子化雑音などの雑音により、理想値からずれた値を示すことと、投影したＸ線が多重である場合にも計算が用意であることを考慮して、これらの雑音をガウス分布でよく表現されると仮定したものである。
また、物質クラスが、出現する度合いを表現するパラメータと前記物質クラスが空間的に連続する度合いを表現するパラメータとで特徴付けられるボルツマン分布で表現されるとした理由は、同じ組織が空間的に集まりやすい(組織別に集まりやすさの調整が可能)ことと、標準とされる各組織の割合を表現しやすいことに着目したものである。【００３１】
また、本発明のＸ線ＣＴ画像処理方法において、ベイズ推定は、物質クラス、Ｘ線吸収係数密度、物質クラスのサブクラス変数、それぞれに関する事後分布の推定が行われる。
ここで、サブクラス変数は、物質クラスが定まった時に、どのようなＸ線吸収係数密度をとりやすいかという確率分布を２つのガンマ分布の和で表現される場合に、いずれのガンマ分布に所属するかを指定する離散変数である。
物質クラスが与えられた条件下でのＸ線吸収係数の事前分布は、特定の値(たとえば平均値)を取りやすいと考えられるものの、その平均値に対して対称な確率分布となるとは限らなかったり、確率分布の尖度が大小様々であったりすることが予想できることから、このような複雑な事前分布を表現できるように２つの確率分布の重み付き和である混合確率分布で表現したものである。さらにＸ線吸収係数は正の値しか取らないため、この混合確率分布をガンマ分布の混合確率分布で表現されると仮定したものである。物質クラスのサブクラス変数は、２つ分布の混合確率分布においてどちらの確率分布に所属するかを指定するために便宜的に導入される０又は１の二値変数である。【００３２】
本発明のＸ線ＣＴ画像処理装置は、Ｘ線投影画像の観測に関する確率分布を設定し統計推定を行うＸ線ＣＴ画像処理装置であって、Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力する入力手段と、観測に関する確率分布と、物質クラスという再構成画像の各画素にどの物質が存在するかを表す変数を介して表現されるＸ線吸収係数に関する事前分布と、物質クラスに関する事前分布のもとで、Ｘ線吸収係数及び物質クラスに関する事後分布の期待値によるベイズ推定でＸ線吸収係数及び物質クラスを推定する推定手段を備える。ここで、観測に関する確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を通じてどのような投影像が得られやすいかを表す確率分布として設定される。また、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。【００３３】
上記のＸ線ＣＴ画像処理装置によれば、物質クラスの存在割合（密度）を表現するパラメータを用いて、物質クラスの変数を確率変数として扱うことで、予め想定していた物質クラスから外れる物質に対しても、予め想定した物質クラスの重み付き和によって表現することで、近似表現が可能になる。また、物質クラスが空間的に連続する度合いを表現するパラメータを用いることにより、物質クラスに対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定精度が向上することが期待できる。【００３４】
また、本発明の他の観点によれば、Ｘ線吸収係数に関する確率分布をおいて統計推定を行うＸ線ＣＴ画像処理プログラムであって、コンピュータに、以下の１）〜２）のステップを実行させるためのプログラムである。
１）Ｘ線投影画像と少なくともＸ線強度を含む計測条件を入力するステップ
２）上記のＸ線投影画像の観測に関する確率分布と、物質に関するパラメータをもつＸ線吸収係数に関する事前分布のもとで、Ｘ線吸収係数及び物質に関する事後分布の期待値でＸ線吸収係数及び物質を推定するステップ
ここで、上記の確率分布は、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程として表現され、多重Ｘ線和は、検出器におけるＸ線の検出面が空間的に１点ではなく所定の広がりを有し、Ｘ線の検出時間が一瞬ではなく所定の幅を有し、かつ、Ｘ線のスペクトル分布が線スペクトルではなく所定の広がりを有し、空間的な位置と方向および波長が異なる複数のＸ線が検出器に多重に入射し、それらの多重のＸ線の重み付き和が検出器で観測されることである。【００３５】
上記のＸ線ＣＴ画像処理プログラムによれば、物質の存在割合（密度）を表現するパラメータを用いて、物質の変数を確率変数として扱うことで、既知の物質で表現されない物質に対しても、その重み付き和による柔軟な表現が可能になる。また、物質が空間的に連続する度合いを表現するパラメータを用いることにより、物質に対して確率分布として空間的滑らかさの制約を課すことで不良設定性を抑制することができる。 ベイズ推定を近似的に行うことで、収束性が保証され、かつ、推定の不安定性を抑えることができる。【００３６】
また、本発明は、上記のＸ線ＣＴ画像処理プログラムを搭載したＸ線ＣＴ画像処理装置を提供することができる。
【発明の効果】
【００３７】
本発明によれば、従来技術を用いた画像再構成と比べて、Ｘ線吸収係数を確率的に表現することで、より柔軟な表現を可能とし、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、ビームハードニングアーティファクトを低減できるといった効果を有する。
また、本発明によれば、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程に関する確率分布を設定し統計推定を行うことで、一本の光線で代表させることによるバイアスを低減できるといった効果を有する。
【図面の簡単な説明】
【００３８】
【図１】本発明のＸ線ＣＴ画像処理方法の処理フロー
【図２】Ｘ線ＣＴの説明図
【図３】更新の説明図
【図４】金属と骨のＸ線吸収係数の確率分布図
【図５】実施例１の計算機実験での画像再構成を示す図
【図６】実施例２の計算機実験での画像再構成を示す図
【図７】平行ビームにおける投影の幾何の説明図
【図８】平行ビームにおける投影と画素化の説明図
【図９】回転時の投影の幾何の説明図
【発明を実施するための最良の形態】
【００３９】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。【００４０】
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、Ｘ線吸収係数に関する事前分布を設定し統計推定を用いて画像の再構成を行う。特に、Ｘ線スペクトルのエネルギー分布（波長分布）を考慮し、エネルギー依存（波長依存）のＸ線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングアーティファクトを解消する。
推定する変数増加による不良設定性に対処するため、適当な事前知識を用いて解に何らかの制約を課すことで解決を試みる必要がある。本発明では基底関数表現を用いて、その基底を既知の物質のエネルギー依存（波長依存）のＸ線吸収係数とし、その重みはどの物質クラスに所属するかを表現する離散変数と、Ｘ線吸収係数密度を表現する非負のスカラーの連続変数の積で与えられるとし、物質クラスを条件としたＸ線吸収係数密度の事前分布と物質クラスの事前分布を課して解に制約を課す。【００４１】
また、本発明では、人体組織を対象とする再構成画像の場合、人体組織は脂肪，筋肉，骨といった限られた物質から構成され、それぞれの物質のＸ線吸収係数の大まかな分布はあらかじめ分かっているので、その知識を事前分布として表現して推定に利用する。 この推定のためには、各画素がどの物質クラスに所属するかの推定を行う必要があるが、それぞれの物質は空間的な連続性を有しており、各物質が人体に占める割合も大まかに分かっているという状況を想定できることから、物質クラスを表す離散変数に空間的な滑らかさ制約や占有割合などの事前知識を導入することにより不良設定性の抑制を行う。【００４２】
投影Ｘ線の観測に関する観測過程については、投影したＸ線が多重に検出器に入射する多重Ｘ線和の観測過程を表現するものであり、Ｘ線が有限の広がりを持つ波長域（スペクトル分布）を有することによって、波長が異なる複数のＸ線が検出器に入射することを表現する。
再構成画像に関する確率分布については、再構成画像の各画素の領域において、Ｘ線吸収係数がどのような値をとりやすいかを表現するかを特徴付けるものであって、各物質クラスと各サブクラス変数の組み合わせがＸ線波長毎にどのようなＸ線吸収係数密度を取りやすいかを表現する混合ガンマ分布と、各物質クラスと各サブクラス変数の組み合わせがどの程度、出現しやすいかを表現するパラメータと、各物質クラスと各サブクラス変数の組み合わせがどの程度、空間的に滑らかに連続しやすいかを表現するパラメータからなるボルツマン分布で表現する。【００４３】
そして、Ｘ線吸収係数と物質クラスに関する事後分布の期待値推定による統計推定を行う。画像再構成において、この事後分布を求めることが必要であるが、各画素で定義される物質クラスという高次元の隠れ離散状態変数に関する和計算が含まれるため、解析的に計算を実行することは困難である。よって、本発明では、近似手法を適用することでこの計算困難さを克服する（図１のフローチャートを参照）。【００４４】
（問題の定式化）
Ｘ線ＣＴ画像のように、様々な方向から投影されたＴ個の投影データを、Ｄ＝｛Ｙ（１），・・・，Ｙ（Ｔ）｝と表すことにする。各々のデータＹ（ｔ）は、ｔ番目の投影によって検出器で検出されるデータの集合であり、Ｙ（ｔ）＝｛ｙ１（１），・・・，ｙＩ（ｔ）｝となる。但し、Ｔは投影数、Ｉは検出器の数、ｙｉ（ｔ）はｉ番目の検出器で検出された光子の数を表す。
そして、Ｘ線は物質を透過する際に指数的に減衰することから、下記数式（１）のように表すことができる。【００４５】
【数１】
2014185078.tif 000003
【００４６】
ここで、ｘｊは観測対象を撮像対象のＸ線吸収係数をラスタスキャンして得られるＪ次元ベクトルｘ＝｛ｘ１，・・・，ｘＪ｝のｊ番目の画素のＸ線吸収係数である。ｂｉ（ｔ）は、Ｘ線源から放出される光子数（物体が何も置かれていないときに観測され得る光子数）を表す。また、ｌｉｊ（ｔ）は角度θ（ｔ）から投影された際のｉ番目の検出器で検出される投影線とｊ番目の画素が交差する距離に相当するものであり、ｌｉｊ（ｔ）xｊがｔ番目の投影のｉ番目の検出器に入射するＸ線に対するｊ番目の画素領域の実効的なＸ線吸収係数を表す（図２を参照）。【００４７】
ここで、Ｘ線のエネルギースペクトルをＥ通りに離散化して表現することを考える。入射Ｘ線の内、エネルギーｅ（ｅ＝１，・・・，Ｅ）に属するＸ線の光子数をＮｅ（ｔ）とすると、物質に透過した後に観測される光子数の平均値は、下記数式（２）のように表すことができる。
ただし、ｘｊｅはエネルギーｅのＸ線に対するｊ番目の画素のＸ線吸収係数である。【００４８】
【数２】
2014185078.tif 000004
【００４９】
（Ｘ線吸収係数のエネルギー依存性について）
ｊ番目の画素におけるエネルギーｅでのＸ線吸収係数ｘｊｅを、エネルギー毎、画素毎に自由に決定すると不良設定性が生じる。実際には、いずれの物質も高エネルギーになるにつれてＸ線吸収係数が小さくなるという物質間で共通した性質を有する。つまり、異なるエネルギー間のＸ線吸収係数の比に特定の傾向が存在する。そこで、物質種を区別する物質クラスを導入して、Ｘ線吸収係数を制約した形で表現する。
観測対象の物質がＣ種類の物質クラスに分類されると仮定する。ここで、ｊ番目の画素の物質クラスを変数ｚｊを用いて表現する。ただし、ｚｊは、ｚｊ＝｛ｚｊ１，・・・, ｚｊＣ｝と表される。各要素ｚｊｃは０もしくは１をとる二値変数であり、ｊ番目の画素が物質クラスｃに属している場合に１をとり、他の場合は０をとる。この物質クラスに依存して一意に決まるエネルギー間のＸ線吸収係数の比を表現する変数として、変化率ｒｃｅ（＞０）を定義する。ただし、変化率ｒｃｅは、ｅ＝１の時ｒｃｅ１＝１と正規化されているものとする。ｅ＝１におけるｊ番目の画素のＸ線吸収係数をｘｊ（＞０）で表現する。以上の表現を用いて、Ｘ線吸収係数ｘｊｅは下記数式（３）のように表すことができる。【００５０】
【数３】
2014185078.tif 000005
【００５１】
上記数式（３）の表現を用いた場合、各画素で推定する必要がある変数は、スカラーの連続変数となるｘｊとＣ種類のいずれかの値をとる離散変数ｚｊの２つだけとなる。【００５２】
（観測モデルについて）
各検出器で観測される光子数の理想値は、上述の数式（１）に示した通りである。しかし、実際には光子散乱や量子化雑音などの雑音によって理想値からずれた値を示すことがわかっている。これらの雑音が、ガウス分布でよく表現されると仮定すると、各検出器における観測データは下記数式（４）のように表すことができる。但し、Ｘ＝｛ｘ１，・・・，ｘＪ｝、Ｚ＝｛ｚ１，・・・，ｚＪ｝である。
また、観測が、検出器毎、投影毎に独立であると仮定すると、投影データの集合Ｄ（観測データ）は下記数式（５）のように表すことができる。但し、Ｔは投影数、Ｉは検出器の数である。ここで、ガウス分布の標準偏差を表すσは、ショット雑音の性質を考慮して、下記数式のように観測ｙｉ（ｔ）に依存する関数として良い。【００５３】
【数４】
2014185078.tif 000006
【００５４】
【数５】
2014185078.tif 000007
【００５５】
（事前分布について）
事前分布については、観測対象の物質種ごとのＸ線吸収係数に関する事前知識を表現する。画素ｊでのＸ線吸収係数ｘｊは画素ｊにおける物質クラスｚｊｃのみに依存して決定されると仮定すると、事前分布ｐ（Ｘ｜Ｚ）は下記数式（６）のように表すことができる。【００５６】
【数６】
2014185078.tif 000008
【００５７】
物質クラスｚｊｃが与えられたもとでは、Ｘ線吸収係数ｘｊは特定の値を取りやすいことが予想できる。ここでは事前分布ｐ（ｘｊ｜ｚｊ）は、分布の非対称性や様々な異なる尖度を表現できるよう２つの分布の混合分布として下記数式（７）のように表すことができる。但し、変数ｂｊは二値変数である。また、ベクトルＢをＢ＝｛ｂ１，・・・，ｂＪ｝で定義する。【００５８】
【数７】
2014185078.tif 000009
【００５９】
画素毎の物質クラスｚｊｃと変数ｂｊｄが与えられたとき、Ｘ線吸収係数ｘｊはガンマ分布に従うと仮定すると下記数式（８）のように表すことができる。但し、ｕｃｄ，ｖｃｄはガンマ分布のパラメータであり、平均ｕｃｄ／ｖｃｄ、分散ｕｃｄ／ｖ２ｃｄである。これらのパラメータは物質クラスに関する経験的な知識により定める。【００６０】
【数８】
2014185078.tif 000010
【００６１】
また、事前分布において変数ｚｊｃは下記数式（９）を満たす、０もしくは１をとる二値変数とする。ガンマ分布を用いることでＸ線吸収係数が負にならない制約を表現できる。【００６２】
【数９】
2014185078.tif 000011
【００６３】
ここで、変数Ｚ，Ｂの事前分布に関しては変数間の独立性を仮定せず、下記数式（１０）で表されるように、ボルツマン分布に従うと仮定する。但し、Ａは正規化項であり、エネルギーＨ（Ｚ，Ｂ）は下記数式（１１）のように定義する。下記数式（１１）では、２次元平面の近傍として縦横に隣接する４画素を考えており、Ｇｃｄｓｅｌｆ，Ｇｃｄｉｎｔｅｒはそれぞれ非負の定数である。【００６４】
【数１０】
2014185078.tif 000012
【００６５】
【数１１】
2014185078.tif 000013
【００６６】
（事後分布について）
本発明のＸ線ＣＴ画像の再構成では、Ｘ線吸収係数Ｘ、物質クラスＺ，変数Ｂは観測データＤの事後分布として推定される。Ｘ線吸収係数Ｘ、物質クラスＺ、変数Ｂの事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）はベイズの定理により下記数式（１２）のように表すことができる。【００６７】
【数１２】
2014185078.tif 000014
【００６８】
（ベイズ推定）
ベイズ推定では、上記数式（１２）で表される事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）の期待値が必要になるが、この事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）を解析的に計算するのは困難であることに鑑みて、試験分布ｑ（Ｘ，Ｚ，Ｂ）により事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）を近似する。試験分布ｑ（Ｘ，Ｚ，Ｂ）は、上記数式（１２）の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。そして、試験分布ｑ（Ｘ，Ｚ，Ｂ）と事後分布ｐ（Ｘ，Ｚ，Ｂ｜Ｄ）の誤差をＫＬ（Kullback−Leibler）距離で評価し、ＫＬ距離を最小化する試験分布ｑ（Ｘ，Ｚ，Ｂ）を算定する。ここで、ＫＬ距離は下記数式（１３）のように表すことができる。
但し、＜・＞ｑ（Ｘ，Ｚ，Ｂ）は分布ｑ（Ｘ，Ｚ，Ｂ）に関して積分計算を行うことを表す演算子である。ＫＬ距離は常に非負であり、ｑ＝ｐのときにのみ０（ゼロ）になる。＜・＞ｑ（Ｘ，Ｚ，Ｂ）はｑ（Ｘ，Ｚ，Ｂ）に関して期待値をとることを表している。【００６９】
【数１３】
2014185078.tif 000015
【００７０】
この最適化を容易にするため、試験分布ｑ（Ｘ，Ｚ，Ｂ）がその周辺分布ｑ（Ｘ）とｑ（Ｚ，Ｂ）の積で与えられる、すなわち変数Ｘと｛Ｚ，Ｂ｝の間に独立性があることを仮定する。また、ｑ（Ｘ），ｑ（Ｚ，Ｂ）の各要素について画素間で独立であると仮定し、さらにｑ（Ｘ）の各要素ｑ（ｘｊ）はガンマ分布に従うとする。
このときｑ（ｚｊ、ｂｊ）は、確率変数ｚｊ, ｂｊが離散変数であるため必然的に多項分布となる。このことから、試験分布ｑ（Ｘ，Ｚ，Ｂ）は下記数式（１４），（１５）のように表すことができる。但し、試験分布ｑ（Ｘ），ｑ（Ｚ，Ｂ）について同時に最適化するのは困難であるため、下記数式（１６），（１７）のように交互に最適化を行うものとする。【００７１】
【数１４】
2014185078.tif 000016
【００７２】
【数１５】
2014185078.tif 000017
【００７３】
【数１６】
2014185078.tif 000018
【００７４】
【数１７】
2014185078.tif 000019
【００７５】
ここで、ｑ（ｘｊ）のパラメータαｊ，βｊの最適化にはＳＣＧ（Scaled Conjugate Gradient）法を用いた。試験分布ｑ（Ｘ）とｊ番目の画素以外の       ｑ（ｚｊ´，ｂｊ´) (但し、ｊ´はｊと異なる) が固定のもとで、ＫＬ距離を最小化するｑ（ｚｊ，ｂｊ)は解析的に求めることができる。この解析解はｊの近傍画素のｑ（ｚｊ´，ｂｊ´）（ｊ´∈η（ｊ）)に依存する。そこで、近傍画素の分布を変更せずにｑ（Ｚ，Ｂ）を更新するために、図３に示すように市松模様となる画素の部分集合を交互に最適化している。すなわち、パラメータｑ（Ｚ，Ｂ）は全画素を同時に更新するのではなく、図３に示す市松模様のパターンで画素の部分集合を交互に最適化している。図３（ａ）のグレー部分の画素に該当するｑ（ｚｊ，ｂｊ）を更新したのち、図３（ｂ）のグレー部分の画素に該当するｑ（ｚｊ，ｂｊ）を更新する。【００７６】
以上、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおける処理内容について説明した。以下の実施例では、本発明の有効性を調べるために計算機実験を行い、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの画面再構成を、従来技術による画面再構成と比較して評価する。
下記の評価実験から、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの有用性が理解できるであろう。【００７７】
実施例では、再構成画像と真の画像と誤差を測れるように、適当な真の画像を用意し、それに対してＸ線投影をシミュレートして人工的に投影像（サイノグラム）を生成した。Ｘ線投影は、物質のＸ線吸収係数のエネルギー依存性（波長依存性）を考慮するために入射するＸ線のエネルギーを６０ｋｅＶ，７０ｋｅＶ，８０ｋｅＶ，９０ｋｅＶ，１００ｋｅＶの５通り（それぞれのフォトン数は、６×１０４，２．５×１０４，１．５×１０４，７×１０３，５×１０３）に分割し、それぞれで物質固有のＸ線吸収係数に従って減衰させた。検出器においてはその減衰したＸ線の光子にポアソンノイズに従って生成したノイズ（雑音）が加えられたものが観測されるとした。【００７８】
以下の実施例１，２において、検出器の個数は９５（個）、画像の解像度は６４×６４（画素）、投影は１°から１８０°の間で１°間隔としている。再構成画像の推定の良さを定量的に評価するために、下記数式（１８）で表されるＰＳＮＲ（Peak Signal to Noise Ratio) を用いた。【００７９】
【数１８】
2014185078.tif 000020
【００８０】
上記数式（１８）において、ＭＡＸは画素値の最大値であり、ＭＳＥは真の画像と再構成画像との平均２乗誤差である。ＰＳＮＲが高いほどＭＳＥが小さい。すなわち、より真の画像に近いといえる。以下、ＦＢＰ法、Ｘ線エネルギー分布を考慮しない事前分布をガウス分布としてベイズ推定により行う方法（以下、「従来方法」と略する）、Ｘ線エネルギーを３つに離散化した本発明のＸ線ＣＴ画像処理方法（以下、「本発明方法」と略する）の３通りで画像再構成を行い、それぞれ比較した結果を示す。
なお、金属と骨の事前分布は図４に示す通りである。図４は、金属と骨の事前分布と真の画像の６０ｋｅＶにおけるＸ線吸収係数を示している。金属と骨のＸ線吸収係数は、それぞれ０．７，０．４を平均とするガンマ分布で表現される確率分布である。真の画像のＸ線吸収係数は、それぞれのＸ線エネルギーにおいて、事前分布の平均値より５％小さい値であるとしている。
【実施例１】
【００８１】
実施例１では、金属を４つ含むファントムに対する画像再構成についての計算機実験について図５を参照して説明する。
図５（ａ）は真の画像、（ｂ）はＦＢＰ法による再構成画像、（ｃ）は単純ベイズ法による再構成画像、（ｄ）は本発明方法による再構成画像である。
実施例１では、図５（ａ）に示す真の画像に対して、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射Ｘ線の光子数の合計を１．１２×１０５とし、本発明のＸ線ＣＴ画像処理方法において、Ｘ線のエネルギーを６０ｋｅＶ，８０ｋｅＶ，１００ｋｅＶの３つに離散化し、それぞれに属する光子数を７．５×１０４，３×１０４，１０．７×１０４であるとした。繰り返し回数は比較のため共に２５回と設定した。推定結果を、図５（ｂ）〜（ｄ）に示す。【００８２】
また、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれにおけるＰＳＮＲを下記表１に示す。ＰＳＮＲでは、ＦＢＰ法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。【００８３】
【表１】
2014185078.tif 000021
【実施例２】
【００８４】
実施例２では、金属インプラントを含む歯のファントムに対する画像再構成についての計算機実験について図６を参照して説明する。
図６（ａ）は真の画像、（ｂ）はＦＢＰ法による再構成画像、（ｃ）は単純ベイズ法による再構成画像、（ｄ）は本発明方法による再構成画像である。
実施例２では、図６（ａ）に示す真の画像に対して、ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれの方法を用いて画像再構成実験を行った。入射Ｘ線の光子数の合計を１．１２×１０５とし、本発明のＸ線ＣＴ画像処理方法において、Ｘ線のエネルギーを６０ｋｅＶ，８０ｋｅＶ，１００ｋｅＶの３つに離散化し、それぞれに属する光子数を７．５×１０４，３×１０４，１０．７×１０４であるとした。繰り返し回数は比較のため共に２５回と設定した。推定結果を、図６（ｂ）〜（ｄ）に示す。【００８５】
ＦＢＰ法、単純ベイズ法、本発明方法のそれぞれにおけるＰＳＮＲを下記表２に示す。実施例１の実験結果と同様に、ＰＳＮＲでは、ＦＢＰ法、単純ベイズ法と比べて本発明方法が大きく、本発明方法がより精度の高い推定を行えることが理解できる。【００８６】
【表２】
2014185078.tif 000022
【００８７】
以上のように、Ｘ線スペクトルのエネルギー分布を考慮し、エネルギー依存性（波長依存性）のあるＸ線吸収係数を統計モデルに基づいて推定することにより、ビームハードニングの影響による様々なアーティファクトを軽減できることを示した。実施例１，２の計算機実験としてはＸ線エネルギーを３つに分割した場合のみを行ったが、このように実際のエネルギーの分割数より少ない数で離散表現しても、アーティファクトが大きく軽減できることがわかる。さらにエネルギーの分割数を増やすことでより良い推定を行うことが可能であろう。
【実施例３】
【００８８】
次に、投影Ｘ線の観測に関する観測過程について、多重Ｘ線和が、検出器におけるＸ線の検出面が空間的に有限の広がりを有すること、検出器におけるＸ線の検出時間が有限の幅を有すること、Ｘ線が有限の広がりを持つ波長域（スペクトル分布）を有すること、によって、空間的な位置と方向、波長が異なる複数のＸ線が検出器に入射することを表現できることを説明する。【００８９】
まず、平行ビームを例にとり、図７を参照して説明する。図７のように２次元平面で平行ビームによる投影を想定する。ここで、Ｘ線の放出される方向をｖ軸、それに直交する方向をｕ軸ととる。この座標系で以下のように変数を定義する。
位置（ｕ，ｖ）における物質の波長λのＸ線の吸収係数をｘ（ｕ，ｖ，λ）とする。但し、物質は、ｖ∈Ｌに含まれるとする。
ｉ番目の検出器の検出面Ｓｉに、単位時間あたり入射する波長λのフォトン数をｙｉ（λ）とする。平行ビームのＸ線源が、単位時間、単位面積あたり放出する波長λのフォトン数をｄ（λ）、波長λのフォトン１個が励起する電子の平均的な数をρ（λ）とする。
このとき単位時間にi番目の検出器で観測される電流Ｙｉの期待値は、下記数式（１９）のように表される。さらに、ある一定時間の電流を貯めてそれを観測する場合は、時間に関する積分を行うことになる。【００９０】
【数１９】
2014185078.tif 000023
【００９１】
上記数式（１９）は逆投影法の仮定する線積分（上記のｖに関する積分）に帰着できないことから、検出面の空間的な広がりや、Ｘ線のスペクトル分布の拡がりが存在する場合は線積分を仮定する推定手法を用いるとアーティファクトが生じることがわかる。同様に、検出器が受光している時間における回転を考慮する場合も線積分に帰着できない。上記の多重積分の形では計算が難しいので積分を和で近似することを考える。これが多重Ｘ線和を考えることになることを示す。まず、上記数式（１９）において特定の波長λにのみ注目し、下記数式（２０）におけるｕ, ｖに関する積分を和に置き換える。【００９２】
【数２０】
2014185078.tif 000024
【００９３】
計算機上で上記のｕ, ｖに関する積分を近似的に表現するために離散化を行う（図８を参照）。  空間をｕ軸、ｖ軸にそって離散化すると、上記数式（２０）におけるｖに関する積分は下記数式（２１）のように近似できる。【００９４】
【数２１】
2014185078.tif 000025
【００９５】
ここで、ｘｊはｊ番目の画素のＸ線吸収係数であり、Ｉ（ｕ）はｕ軸上のｕの位置にある画素のインデックスの集合、ΔＬは画素の１辺の長さを表す。画素の１辺の長さをΔＬ＝１とすると、上記数式（２１）は下記数式（２２）として表現できる。【００９６】
【数２２】
2014185078.tif 000026
【００９７】
したがって、画素の１辺の長さが１のときは、ｙｉは下記数式（２３）となり、さらに同じ検出器に入る光線に関する積分（ｕに関する積分）を和で近似すると、下記数式（２４）のように表現できる。【００９８】
【数２３】
2014185078.tif 000027
【００９９】
【数２４】
2014185078.tif 000028
【０１００】
上記数式２４において、１／Ｎは積分を短冊で切って近似するときの短冊の幅を示している。
これはインデックスｋごとに線積分で表現される投影の和を考えていることに相当する。
以上から、検出器の空間的拡がりを考慮する場合においては多重のＸ線を考慮する必要があることがわかる。【０１０１】
上記では、Ｘ線が画素に対してまっすぐ入射する場合を想定している。しかしながら、一般のＸ線ＣＴ装置などのように、Ｘ線源と検出器が対向しながら回転するような装置（図９を参照）では、画素に対して斜めに入射する状況を想定する必要がある。
この場合、投影についての下記数式（２４）を、下記数式（２５）のように一般化することで表現できる。【０１０２】
【数２５】
2014185078.tif 000029
【０１０３】
【数２６】
2014185078.tif 000030
【０１０４】
ここで、ΔＮ＝１／Ｎであり、ｒ（ｋ）は光線ｋの検出平面上の位置を表し、｛ｋ｜ｒ（ｋ）∈Ｓｉが検出器ｉに到達する光線の集合を表す。Ｉｋは光線ｋが交差する画素のインデックスの集合を表し、ｌｋｊは光線ｋが画素ｊと交差する長さを表している。上記数式（２５）は上記数式（２６）のｌｋｊを１とした場合と考えることができる。
上記数式（２６）は、平行ビームのほか、ファンビーム、コーンビームでも光線の向きを考えて、ｌｋｊを計算することで同様の式が成立する。
上記数式（２６）からわかることは、検出面の拡がりを考えた投影の観測モデルを考えた時には、検出面内を通る光線ごとの投影を考える必要がある点である。【０１０５】
上記数式（２６）の中のｌｋｊｘｊの和を、検出器ｉや光線ｋに関してまとめて表記するならば、Ｌｘのような行列演算で表現できる。検出面の拡がりを考えなかった場合には、行列Ｌは（検出器の数）×（再構成画像の画素数）であったものが、検出面の拡がりを考えた場合には、行列Ｌは（検出器を通る光線の数）×（再構成画像の画素数）となる。光線の数は、多ければ積分を精度よく近似できる一方で、計算量が増す問題がある。【０１０６】
次に、投影中の回転（時間に関する積分）を想定する。まず時刻ｔにおいて観測されたフォトン数をｙｉ，ｔとすると、下記数式（２７）で表現できる。【０１０７】
【数２７】
2014185078.tif 000031
【０１０８】
検出器の検出時間（露光時間）の間に計測されたフォトン数は、ｙｉ，ｔを時刻ｔに関して積分したものとなるので、下記数式（２８）で表現できる。【０１０９】
【数２８】
2014185078.tif 000032
【０１１０】
上記数式（２８）について、時刻に関する積分を和で置き換えると、下記数式（２９）のように表現できる。【０１１１】
【数２９】
2014185078.tif 000033
【０１１２】
前述のように、検出面の拡がりを考慮すると上記数式（２８）のように各時刻ｔで検出器に入射するフォトンｙｉ，ｔが複数の光線束の和で表現される。時間に関する積分は、さらなる和を含めるため新たな光線束に関する和を考えることに等価である。そこで、光線インデックスの集合Ｓｉを時刻に依存するようＳｉ，ｔとし、それぞれの光線ｋが通る画素インデックスの集合Ｉｋを改めて定義しなおすと、検出時間（露光時間）Ｔで計測されたフォトン数は、下記数式（３０）で表現できる。【０１１３】
【数３０】
2014185078.tif 000034
【０１１４】
上記数式（３０）を波長に関する積分を和で置き換えると、下記数式（３１）で表現できる。【０１１５】
【数３１】
2014185078.tif 000035
【０１１６】
以上をまとめると、検出面の拡がりと投影中の回転,Ｘ線の波長依存性を考慮した場合、下記数式（３２）で表現できる。ここで、ｂｍ＝ｄｍρｍである。【０１１７】
【数３２】
2014185078.tif 000036
【０１１８】
上記数式（３２）のように、ｍやｋに関する和が指数関数の外側でとられることを多重Ｘ線和と呼ぶ。これは、指数関数の外側に和（もしくは積分）があることを想定しない線積分では表現できない。このような多重Ｘ線和のある観測過程一般に対して、本明細書で開示したベイズ推定に基づくアルゴリズムを適用することができる。
【産業上の利用可能性】
【０１１９】
本発明は、ビームハードニングアーティファクトを低減するＸ線ＣＴ装置として、医療用のみならず産業用のＸ線ＣＴ装置に有用である。また、ベルトコンベアによる移動や回転が無視できないベルトコンベア型Ｘ線ＣＴ装置や歯科用Ｘ線ＣＴ装置としても有用である。
【図１】
2014185078.tif 000037
【図２】
2014185078.tif 000038
【図３】
2014185078.tif 000039
【図４】
2014185078.tif 000040
【図５】
2014185078.tif 000041
【図６】
2014185078.tif 000042
【図７】
2014185078.tif 000043
【図８】
2014185078.tif 000044
【図９】
2014185078.tif 000045
【手続補正書】
【提出日】平成27年8月6日(2015.8.6)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
観測対象にＸ線を投影し、その投影されたＸ線をＸ線検出器において検出することで得られる投影Ｘ線データをもとに観測対象内部をＸ線吸収係数の空間分布として画像化するＸ線ＣＴ画像処理方法において、
１）Ｘ線検出器によって検出される投影Ｘ線が、空間的に異なる位置に入射する複数の投影Ｘ線、時間的に異なる時間に入射する複数の投影Ｘ線、および、異なる波長の複数の投影Ｘ線からなる、複数の投影Ｘ線の和であるとした観測過程の場合に、前記観測過程における投影Ｘ線の確率分布として設定するステップと、
２）観測対象内部を構成する物質のＸ線吸収係数および物質に関するパラメータに関する事前確率分布を設定するステップと、
３）前記投影Ｘ線データと前記観測過程における投影Ｘ線の確率分布と前記事前確率分布から、前記Ｘ線吸収係数及び前記物質に関する事後確率分布を求め、前記Ｘ線吸収係数の空間分布を推定するＸ線吸収係数推定ステップ、
を備えたことを特徴とするＸ線ＣＴ画像処理法。
【請求項２】
前記Ｘ線吸収係数推定ステップは、
前記観測過程における投影Ｘ線の確率分布と前記事前確率分布とから、前記事後確率分布とその間の分布間距離が小さくなるように、前記Ｘ線吸収係数及び前記物質に関する試験分布を算出するステップと、
前記Ｘ線吸収係数に関する試験分布の期待値を求めることである、
ことを特徴とする請求項１のＸ線ＣＴ画像処理法。
【請求項３】
前記複数の投影Ｘ線の和に、異なるＸ線波長の投影Ｘ線に関する和を含める場合では、
各画素の波長依存性を有するＸ線吸収係数に関する事前確率分布は、
前記物質毎およびＸ線波長毎のＸ線吸収係数の分布を表現するパラメータを用いて設定される、
ことを特徴とする請求項１又は２に記載のＸ線ＣＴ画像処理法。
【請求項４】
前記複数の投影Ｘ線の和に、異なるＸ線波長の投影Ｘ線に関する和を含める場合では、
各画素の波長依存性を有するＸ線吸収係数は、
波長に依存しない各画素のＸ線吸収係数密度と、画素に依存せず前記物質毎に定められるＸ線波長間のＸ線吸収係数比の積として設定され、
前記物質毎に定められる、Ｘ線波長間のＸ線吸収係数比を表現するパラメータを用いて設定され、
各画素の波長依存性を有するＸ線吸収係数の事前確率分布が前記物質を条件とした各画素の前記Ｘ線吸収係数密度に対する条件付き事前確率分布と前記物質に対する事前確率分布によって設定される、
ことを特徴とする請求項１〜３の何れかに記載のＸ線ＣＴ画像処理法。
【請求項５】
前記物質に対する事前確率分布が、それぞれの物質の出現する度合いを示すパラメータとそれぞれの物質が空間的に連続する度合いを示すパラメータとで表現される事前確率分布で設定されることを特徴とする請求項４に記載のＸ線ＣＴ画像処理法。
【請求項６】
前記物質毎に定められている前記Ｘ線吸収係数密度の条件付き事前確率分布が各物質のサブクラス毎に定められている前記Ｘ線吸収係数密度の条件付き事前確率分布と各物質のサブクラスの事前確率分布の積を各物質のとりえる全てのサブクラスに関して和をとった混合確率分布で設定され、前記各物質のサブクラスの事前確率分布が、各物質のサブクラスの出現する度合いを示すパラメータと前記各物質のサブクラスが空間的に連続する度合いを示すパラメータとで設定されることを特徴とする請求項４に記載のＸ線ＣＴ画像処理法。
【請求項７】
観測対象にＸ線を投影し、その投影されたＸ線をＸ線検出器において検出することで得られる投影Ｘ線データをもとに観測対象内部をＸ線吸収係数の空間分布として画像化するＸ線ＣＴ装置において、
Ｘ線検出器によって検出される投影Ｘ線が、空間的に異なる位置に入射する投影Ｘ線、時間的に異なる時間に入射する複数の投影Ｘ線、および、異なる波長の複数の投影Ｘ線から選ばれる少なくとも一つからなる、複数のＸ線の和であるとした観測過程の場合に、前記観測過程における投影Ｘ線の確率分布を設定する手段、
観測対象内部を構成する物質のＸ線吸収係数および物質に関するパラメータに関する事前確率分布を設定する手段、
前記投影Ｘ線データと前記投影Ｘ線データの尤もらしさを表す確率分布と前記事前確率分布から、前記Ｘ線吸収係数及び前記物質に関する事後確率分布を求め、前記Ｘ線吸収係数の空間分布を推定する手段、
を有することを特徴とするＸ線ＣＴ装置。
【請求項８】
観測対象にＸ線を投影し、その投影されたＸ線をＸ線検出器において検出することで得られる投影Ｘ線データをもとに観測対象内部をＸ線吸収係数の空間分布として画像化するプログラムであって、
コンピュータに、
１）Ｘ線検出器によって検出される投影Ｘ線が、Ｘ線検出器の検出面に対し、空間的に異なる位置に入射する複数の投影Ｘ線、時間的に異なる時間に入射する複数の投影Ｘ線、および、異なる波長の複数の投影Ｘ線から選ばれる少なくとも一つからなる、複数のＸ線の和であるとした観測過程の場合に、前記観測過程における投影Ｘ線の確率分布を設定するステップと、
２）観測対象内部を構成する物質のＸ線吸収係数および物質に関するパラメータに関する事前確率分布を設定するステップと、
３）前記投影Ｘ線データと前記観測過程における投影Ｘ線の確率分布と前記事前確率分布から、前記Ｘ線吸収係数及び前記物質に関する事後確率分布を求め、前記Ｘ線吸収係数の空間分布を推定するステップ、
を実行させるＸ線ＣＴ画像処理プログラム。
【請求項９】
前記Ｘ線吸収係数推定ステップは、
前記観測過程における投影Ｘ線の確率分布と前記事前確率分布とから、前記事後確率分布とその間の分布間距離が小さくなるように、前記Ｘ線吸収係数及び前記物質に関する試験分布を算出するステップと、
前記Ｘ線吸収係数に関する試験分布の期待値を求めることである、
ことを特徴とする請求項８のＸ線ＣＴ画像処理プログラム。
【請求項１０】
請求項８又は９のＸ線ＣＴ画像処理プログラムが搭載されたＸ線ＣＴ装置。
【国際調査報告】
2014185078.tif 000046
2014185078.tif 000047