(57)【特許請求の範囲】
【請求項１】
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
前記事前知識は、
再構成画像の各画素の領域において定義されるパラメータで、
撮像対象の各組織クラスの存在割合を表現するパラメータと、
各組織クラスのＸ線吸収係数を表現するパラメータと、
各組織クラスの空間的に連続する度合いを表現するパラメータと、
によって特徴付けられる確率分布で表現され、
Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用い、
前記統計推定は、
前記事前知識と尤度から計算される事後分布を用いて、
事後確率最大化による推定(ＭＡＰ推定)、或いは、事後確率の期待値によるベイズ推定により画像再構成及び組織クラス推定を行う、
ことを特徴とするＸ線ＣＴ画像処理方法。
【請求項２】
前記Ｘ線吸収係数を表現するパラメータに関する分布としてガウス分布を用い、ガウス分布における平均の値は組織クラスごとに設定されることを特徴とする請求項１に記載のＸ線ＣＴ画像処理方法。
【請求項３】
前記ＭＡＰ推定は、再構成画像の各画素の領域が属する組織クラスとＸ線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定することを特徴とする請求項１又は２に記載のＸ線ＣＴ画像処理方法。
【請求項４】
前記ベイズ推定は、再構成画像の各画素の領域が属する組織クラスとＸ線吸収係数の両方に関する事後分布の推定を行うことを特徴とする請求項１又は２に記載のＸ線ＣＴ画像処理方法。
【請求項５】
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
３）前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行うステップと、
を備えたことをＸ線ＣＴ画像処理方法。
【請求項６】
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
３）前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定（ベイズ推定）により画像再構成及び組織クラス推定を行うステップと、
を備えたことをＸ線ＣＴ画像処理方法。
【請求項７】
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理プログラムであって、
コンピュータに、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
３）前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのＸ線ＣＴ画像処理プログラム。
【請求項８】
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理プログラムであって、
コンピュータに、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
３）前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定（ベイズ推定）により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのＸ線ＣＴ画像処理プログラム。
【請求項９】
請求項７又は８に記載のＸ線ＣＴ画像処理プログラムを搭載したＸ線ＣＴ画像装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、体内組織構成を考慮したＸ線ＣＴアルゴリズム、そのアルゴリズムを有するプログラムおよび該プログラムが搭載されたＸ線ＣＴ装置に関する。
【背景技術】
【０００２】
従来から、メタルアーティファクトの低減技術や、より低いＸ線被曝量で従来と同程度のＣＴ画像を再構成する技術など、Ｘ線ＣＴの高精度化を可能にするアルゴリズムが研究されている。一般に、Ｘ線ＣＴは、Ｘ線管にＸ線検出器を対向配置するとともに、これらの間にターンテーブルを配置した構成を取る。このターンテーブルの回転軸は、Ｘ線管とＸ線検出器を結ぶＸ線光軸に対して直交する向きとされる。Ｘ線ＣＴ装置においては、測定対象物を３６０°回転させて透過Ｘ線データを収集しなければＣＴ断層画像の再構成を行うことができない。
【０００３】
メタルアーティファクトは、測定対象物にＸ線を強く吸収する金属などの密度の高い物質があると、再構成画像にノイズがのるというものである。すなわち、鉄などの比較的Ｘ線が透過しにくい測定対象物のＣＴ撮像を行った場合、３６０°回転させながら透過Ｘ線データを採取すると、密度の高い物質がＸ線の透過方向に重なった状態での透過Ｘ線データが得られることになって、フィルター補正逆投影法(ＦＢＰ)などの既存手法を用いた場合、再構成画像にはメタルアーティファクトと呼ばれる虚像が生じ、正確な画像の再構成ができないという問題がある。
従って、メタルアーティファクトの低減によって、Ｘ線ＣＴの高精度化が図れることになる。
【０００４】
一方、Ｘ線は強力であれば、再構成画像のＳＮ比が高いことから、Ｘ線被曝量の低減とＳＮ比はトレードオフの関係がある。より低いＸ線被曝量で従来と同程度のＳＮ比の再構成画像を得る技術とは、従来と同じＸ線強度，撮像枚数でも解像度の高い画像が得られるという解像度向上を図る技術や、従来と同じＸ線強度や解像度でも、より撮像枚数を減らして撮像時間を短縮できる撮像時間の短縮を図る技術や、従来と同じ解像度でもＸ線強度を弱めて被曝量を減らせる被曝量の低減を図る技術、従来と同じＸ線強度や解像度でも再構成画像に含まれるアーティファクト(もしくはノイズ)を低減できるノイズ除去を図る技術、従来と同じＸ線強度や解像度でも再構成画像のコントラスト分解能（再構成画像を何階調で表せるか、すなわち、Ｘ線吸収係数の精度）を向上するコントラスト分解能の向上を図る技術がある。
【０００５】
上記の従来技術は、統計推定を行うものと統計推定を行わないものに大きく２つに分類できる。
先ず、統計推定を行わない従来技術としては、画像再構成法の主流であるフィルター補正逆投影法（ＦＢＰ：Filtered
Back Projection）や、ＦＢＰを改良したＰＣＬＩＳ（Projection Completion
Method based on a Linear Interpolation in the Sinogram）が挙げられる。ＦＢＰは、撮像対象に対して色んな方向からＸ線を照射して、その投影像(Ｘ線の吸収像、投影像のセットをシノグラムと呼ぶ)を得た後、得られた投影像を逆方向に投射して画像再構成する際に、シノグラムに周波数領域上でフィルターをかけてボケ除去を行う画像再構成法である。ＰＣＬＩＳは、メタルなどのＸ線を強く吸収する物体によってほとんどＸ線を観測できなかったシノグラムの部分に対してシノグラム上で線形補間を行ってからＦＢＰを行うことでメタルアーティファクトなどの極端な虚像の発生を抑えるものである。
次に、統計推定を行なう従来技術としては、再構成画像に関する事前知識をおかない最尤推定法（ＭＬ)がある。最尤推定法（ＭＬ)は、観測における不確実性を含む物理過程を確率モデルで表現して、その確率モデルを基に最も尤もらしい再構成画像を推定するものである。観測における不確実性を含む物理過程には、観察されるフォトンに関する揺らぎ(ショットノイズ)などが含まれる。
【０００６】
また、撮像対象の物体を異なる位置・角度から撮影した投影画像を複数枚用意できる場合、それら複数の投影画像の情報を処理することで、元の物体の断面像を再構成できることが知られている。投影画像を使用して元の物体の断面像を復元する画像再構成法のうち統計的な解法に関する従来の研究は、最尤推定法（ＭＬ），ＭＡＰ（maximum a posteriori）推定法，ベイズ推定法の３手法が存在する。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】”Suppression of Metal Artifacts in CTUsing a Reconstruction Procedure That Combines MAP and Projection Completion”, IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 28, NO. 2, FEBRUARY2009.
【非特許文献２】” Metalartifact reduction in CT using tissue-class modeling and adaptive prefiltering”, Medical Physics, Vol. 33, No. 8,August 2006.
【非特許文献３】” ベイズアプローチに基づいた断層画像の再構成”，情報処理学会研究報告(2009年)
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上記の従来技術を用いた画像再構成の場合は、従来と同程度の再構成画像を得るためには従来と同程度のＸ線被曝量となり、またメタルアーティファクトの低減が困難であった。
【０００９】
そこで、本発明は、従来技術を用いた画像再構成と比べて、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、メタルアーティファクトを低減可能なＸ線ＣＴ画像処理方法、Ｘ線ＣＴ画像処理プログラムおよび該プログラムが搭載されたＸ線ＣＴ装置を提供することを目的とする。
【課題を解決するための手段】
【００１０】
上記状況に鑑みて、本発明のＸ線ＣＴ画像処理方法は、
Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
（１）観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現され、
（２）上記の再構成画像に関する事前知識は、
再構成画像の各画素の領域において定義されるパラメータで、
撮像対象の各組織クラスの存在割合を表現するパラメータと、
各組織クラスのＸ線吸収係数を表現するパラメータと、
各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現され、
Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用い、
（３）上記の統計推定は、
前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化による推定(ＭＡＰ推定)、或いは、事後確率の期待値によるベイズ推定であり、
画像再構成及び組織クラス推定を行うことを特徴とする。
【００１１】
上記のＸ線ＣＴ画像処理方法によれば、観測されるフォトンに関する確率モデルがポアソン分布などで表現されることにより、従来と比べてより現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。
また、上記のＸ線ＣＴ画像処理方法によれば、人体などの撮像対象の組織分布に関する事前知識について、どのような組織（筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど）がどの程度の割合で分布するか、また、それぞれの組織がどの程度のＸ線吸収係数を有するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいかの分布情報を確率分布の形で表現して、事後確率最大化によるＭＡＰ推定、或いは事後確率の期待値によるベイズ推定統計的推論を用いて、画像再構成や組織推定を行う。組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることが可能である。
【００１２】
ここで、本発明のＸ線ＣＴ画像処理方法において、上記の各組織クラスのＸ線吸収係数を表現するパラメータに関する分布としてガウス分布を用い、かつ、各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布としてボルツマン分布を用いることが好ましい態様である。
【００１３】
各組織クラスのＸ線吸収係数を表現するパラメータに関する分布としてガウス分布を用いることとした理由は、組織クラスごとに決まる特定のＸ線吸収係数（ＣＴ値）をとりやすいことに着目したものである。
各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布としてボルツマン分布を用いることとした理由は、同じ組織クラスが空間的に集まりやすい(組織クラス別に集まりやすさの調整が可能)ことと、標準とされる各組織クラスの割合を表現しやすいことに着目したものである。
【００１４】
また、本発明のＸ線ＣＴ画像処理方法において、上述のＭＡＰ推定は、再構成画像の各画素の領域が属する組織クラスとＸ線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定することが好ましい。
【００１５】
再構成画像の各画素の領域が属する組織クラスとの組合せを推定することは、組織クラス依存性を含めることになる。この組織クラス依存性を含めることで、例えば、空気の組織クラスはピクセル同士の連結がしやすい（空気中に他の組織は入らない）などの組織クラスに依存した知識を含めることができるのである。これにより、より柔軟性をもった分布表現が可能となる。
【００１６】
なお、推定精度を向上すべく、例えば、グラフカットの一種であるα−拡張を用いた最適化アルゴリズムによってＭＡＰ推定することがより好ましい。
【００１７】
また、本発明のＸ線ＣＴ画像処理方法において、上述のベイズ推定は、Ｘ線吸収係数と再構成画像の各画素の領域が属する組織クラスの両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。
なお、ベイズ推定の場合も、ＭＡＰ推定の場合と同様に、組織クラス依存性を含めることができる。これにより、より柔軟性をもった分布表現が可能となる。
【００１８】
また、本発明の他の観点によれば、本発明のＸ線ＣＴ画像処理方法は、Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
３）再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行うステップと、
を備えた構成とされる。
【００１９】
かかる構成によれば、画像の再構成の速度を速めることができ、より高速でＸ線ＣＴ画像処理を行うことができる。
【００２０】
また、本発明の他の観点によれば、本発明のＸ線ＣＴ画像処理方法は、Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理方法であって、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
３）再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定（ベイズ推定）により画像再構成及び組織クラス推定を行うステップと、
を備えた構成とされる。
【００２１】
かかる構成によれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にＸ線ＣＴ画像処理を行うことができる。
【００２２】
また、本発明の他の観点によれば、Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理プログラムであって、
コンピュータに、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
３）再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定（ＭＡＰ推定）により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのプログラムが提供される。
【００２３】
かかるプログラムによれば、画像の再構成の速度を速めることができ、より高速でＸ線ＣＴ画像処理を行うことができる。
【００２４】
また、本発明の他の観点によれば、Ｘ線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うＸ線ＣＴ画像処理プログラムであって、
コンピュータに、
１）投影画像と少なくともＸ線強度を含む計測条件を入力するステップと、
２）投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
３）再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、Ｘ線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
４）前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定（ベイズ推定）により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのプログラムが提供される。
【００２５】
かかるプログラムによれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にＸ線ＣＴ画像処理を行うことができる。
【００２６】
また、本発明は、上記のＸ線ＣＴ画像処理プログラムを搭載したＸ線ＣＴ画像装置を提供することができる。
【発明の効果】
【００２７】
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムによれば、従来技術を用いた画像再構成と比べて、より低いＸ線被曝量で従来と同程度の再構成画像が取得可能で、メタルアーティファクトを低減可能であるといった効果を有する。
【図面の簡単な説明】
【００２８】
【図１】本発明のＸ線ＣＴ画像処理方法（ＭＡＰ推定）の処理フロー
【図２】本発明のＸ線ＣＴ画像処理方法（ベイズ推定）の処理フロー
【図３】Ｘ線ＣＴの説明図
【図４】ＭＡＰ推定を用いたアルゴリズムの擬似コード
【図５】ベイズ推定を用いたアルゴリズムの擬似コード
【図６】事前知識として用いるＣＴ値の分布を示す図
【図７】実験に用いるファントムの説明図
【図８】実施例１のＣａｓｅＡの実験設定での画像再構成を示す図
【図９】実施例１のＣａｓｅＢの実験設定での画像再構成を示す図
【図１０】実施例１のＣａｓｅＡの実験設定での組織クラスの推定結果を示す図
【図１１】実施例１のＣａｓｅＢの実験設定での組織クラスの推定結果を示す図
【図１２】投影画像数とＰＳＮＲの関係を示す図
【図１３】本発明のＸ線ＣＴ画像処理方法の処理時間の説明図
【図１４】実施例２のＣａｓｅＡの実験設定での画像再構成を示す図
【図１５】実施例２のＣａｓｅＢの実験設定での画像再構成を示す図
【図１６】実施例２のＣａｓｅＡの実験設定での組織クラスの推定結果を示す図
【図１７】実施例２のＣａｓｅＢの実験設定での組織クラスの推定結果を示す図
【発明を実施するための最良の形態】
【００２９】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。
【００３０】
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、Ｘ線吸収係数に関する事前知識を用いて組織クラスの統計推定を用いて画像の再構成を行う。特に、放射線による被曝を避けるためにより少ない観測データからの画像再構成を試みる。その一方で、画像再構成は不良設定性を有しているため、適当な事前知識を用いて解に何らかの制約を課すことで解決を試みる。本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、人体は脂肪，筋肉，骨といった限られた物質から構成され、それぞれのＸ線吸収係数の大まかな分布はあらかじめ分かっているものとする。また、それぞれの組織クラスは空間的な連続性を有しており、各組織クラスが人体に占める割合も大まかに分かっているという状況を想定する。
【００３１】
すなわち、再構成画像に関する事前知識について、再構成画像の各画素の領域において定義されるパラメータで、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのＸ線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、これらの知識を階層ベイズモデルにおける事前知識として利用する。
【００３２】
そして、事後分布推定による統計推定を行う。画像再構成において、事後分布を求めることが必要であるが、高次元の隠れ変数に関する積分計算が含まれるため、解析的に実行することは困難であることから、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、事後確率最大化による推定(ＭＡＰ推定)、或いは、事後確率の期待値によるベイズ推定を用いた近似手法を適用することでこの計算困難さを克服することにしたものである（図１，図２のフローチャートを参照）。
【００３３】
（問題の定式化）
先ず、Ｘ線ＣＴ画像のように、様々な方向から投影されたＴ個の投影データを、Ｄ＝｛Ｙ（１），・・・，Ｙ（Ｔ）｝と表すことにする。各々のデータＹ（ｔ）は、ｔ番目の投影によって検出器で検出されるデータの集合であり、Ｙ（ｔ）＝｛ｙ１（１），・・・，ｙＩ（ｔ）｝となる。ｙｉ（ｔ）は、ｉ番目の検出器で検出された光子の数を表す。Ｘ線は物質を透過する際に指数的に減衰することから、下記数式（１）のように表すことができる。
【００３４】
【数１】
0005590548.tif 000002
【００３５】
ここで、ｘｊは観測対象を撮像対象のＸ線吸収係数をラスタスキャンして得られるＪ次元ベクトルｘ＝｛ｘ１，・・・，ｘＪ｝のｊ番目のピクセルのＸ線吸収係数である。ｂｉ（ｔ）は、Ｘ線源から放出される光子数（物体が何も置かれていないときに観測され得る光子数）を表す。また、ｌｉｊ（ｔ）は角度θ（ｔ）から投影された際のｉ番目の検出器で検出される投影線とｊ番目のピクセルが交差する距離に相当するものであり、ｌｉｊ（ｔ）xｊがｔ番目の投影のｉ番目の検出器に入射するＸ線に対するｊ番目のピクセル領域の実効的なＸ線吸収係数を表す（図３を参照）。
【００３６】
（階層ベイズモデル）
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおけるＸ線ＣＴ画像の再構成では、Ｘ線吸収係数ｘは観測データＤの事後分布として推定される。本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、事後分布の平均値としてＸ線吸収係数ｘが得られるとする。事後分布はベイズの定理により下記数式（２）で求めることができる。
【００３７】
【数２】
0005590548.tif 000003
【００３８】
また、人間の体は限られた数の組織から成っており、組織ごとにＸ線吸収係数ｘが定まるという事前知識を用いるため、Ｘ線吸収係数ｘに関して階層的な事前分布を定義することにする。説明の便宜上、以下では、人体の組織は５（＝Ｋ）つの組織（筋肉，脂肪，骨，空気，金属）に分類されると仮定して、事前分布ｐ（ｘ）を隠れ変数ｚ＝｛ｚ１，・・・，ｚＪ｝を用いて下記数式（３）のように表わす。
【００３９】
【数３】
0005590548.tif 000004
【００４０】
ここで、ｚｊはＫ次元の二値変数であり、各要素ｚｊｋはｊ番目のピクセルがｋ番目のクラスに属しているときに１をとり、それ以外の時は０をとるものである。事前分布において、すべてのピクセルはいずれかのクラスに属し、Σｋｚｊｋ＝１を満たすものとする。
【００４１】
次に、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおいて、投影画像に観測されるフォトンに関する物理モデル、事前知識としての組織クラスに関する事前分布、すなわち、撮像対象の各組織クラスの存在割合を表現するパラメータに関する分布、各組織クラスのＸ線吸収係数を表現するパラメータに関する分布、各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布について、以下説明を行う。
【００４２】
（観測されるフォトンに関する物理モデル）
Ｘ線ＣＴにおける観測データの主なノイズの要因は、検出される光子（フォトン）の量子化ノイズであると考えられる。そこで、光子（フォトン）の観測データは、投影毎、検出器毎に独立なポアソン分布に従って生成されるとしてモデル化する。
なお、ポアソン分布に特に限定されるものではなく、物理過程をより良く表現できる他の物理モデルも適用可能である。
上述の数式（１）は、ポアソン分布の平均を表しており、数式（１）の結果を用いて、下記数式（４）で表わす。
【００４３】
【数４】
0005590548.tif 000005
【００４４】
（組織に関する事前分布）
所属する組織クラスの情報が与えられた下で、Ｘ線吸収係数ｘはガウス分布に従うとして、下記数式（５）に示すようにモデル化する。
【００４５】
【数５】
0005590548.tif 000006
【００４６】
ここで、ｖｋとｒ２ｋはガウス分布の平均と分散を表している。平均ｖｋの値は、各々の組織（空気，脂肪，筋肉，骨，金属）ごとにそれぞれ、ｖ１＝0.000，ｖ２＝0.018，ｖ３＝0.022，ｖ４＝0.040，ｖ５＝0.120とする。
また、分散ｒ２ｋは、下記表１に示すように、どの組織に所属するかに関わらず一定とする。なお、各々の平均ｖｋは、対象となる物体の個体差や温度等で変動すると考えられる。これらの不確かさはガウス分布の分散ｒ２ｋによって調整する。
【００４７】
【表１】
0005590548.tif 000007
【００４８】
また、組織のクラスに関する事前分布はボルツマン分布に従うとして、下記数式（６）に示すようにモデル化する。
【００４９】
【数６】
0005590548.tif 000008
【００５０】
上記数式（６）において、Ｚは正規化項であり、エネルギーＥを下記数式（７）のように定義する。
【００５１】
【数７】
0005590548.tif 000009
【００５２】
η（ｊ）は、ｊ番目のピクセルの近傍画素を表している。ＪｋｓｅｌｆとＪｋｉｎｔｅｒは、事前分布をコントロールする非負の定数である。ボルツマン分布は、エネルギーＥ(ｚ) が低ければ低いほどその確率が高くなる。エネルギー項の第一項目は、組織の分布する割合に関わり、第二項目は組織の空間的な連続性に関わる。
【００５３】
（事後分布の推定）
上述した如く、Ｘ線ＣＴ画像の再構成において、事後分布を求めるのであるが、高次元の隠れ変数に関する和計算が含まれるため、解析的に実行することは困難である。このことから、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、事後確率最大化による推定(ＭＡＰ推定)、或いは、事後確率の期待値によるベイズ推定を近似的に行った手法を用いることで計算困難さを回避している。
【００５４】
（ＭＡＰ推定）
ベイズの定理によると、吸収係数ｘと組織クラスｚに関する事後分布（これを前述の事後分布ｐ（ｘ｜Ｄ）と区別すべく、同時事後分布と表現する。）ｐ（ｘ,ｚ｜Ｄ）は、それぞれの事前分布と尤度の積の形で、下記数式（８）のように表せる。このような同時事後分布ｐ（ｘ,ｚ｜Ｄ）の最大値を求めるようにすることで、高次元の隠れ変数に関する和計算が必要なくなるのである。
【００５５】
【数８】
0005590548.tif 000010
【００５６】
上記数式（８）において、変数ｘ、ｚを事後分布最大化基準（ＭＡＰ）により決定する。下記数式（９）に示すように、上記数式（８）の負の対数をとったものを目的関数とし、これを最小化するものとして事後分布最大化を達成するように未知変数を求める。
【００５７】
【数９】
0005590548.tif 000011
【００５８】
しかしながら、連続変数ｘと離散変数ｚに関する同時最適化は困難である。そのため、下記数式（１０）と数式（１１）に示す連続変数ｘと離散変数ｚに関する最適化は交互に最適化を行っている。
【００５９】
【数１０】
0005590548.tif 000012
【００６０】
【数１１】
0005590548.tif 000013
【００６１】
ここで、連続変数ｘ＊に関しては、ＳＣＧ (Scaled conjugate gradient) 法を用いて最適化を行っている。また、離散変数ｚ＊に関して、グラフカットの一種であるα−拡張アルゴリズムにより最適化している。
ＭＡＰ推定のアルゴリズムの擬似コードを、図４に示す。なお、アルゴリズムは事前に決められた所定の収束条件を満たすまで繰り返される。
【００６２】
（グラフカット）
ここでは、グラフカットの簡単に説明する。グラフカットとは、離散変数の最適化を行う際によく用いられるアルゴリズムである。一般の離散変数最適化問題は、ＮＰ困難であることが知られているが、劣モジュラ性と呼ばれる性質を満たす場合にグラフカットにより大域最小解が得られることが示されている。
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、離散変数ｚの最適化としてグラフカットの一種であるα−拡張アルゴリズムを用いることにしている。このアルゴリズムは、大域最小解が得られることが保証されているわけではないが、適用できる範囲が広くコンピュータビジョンなどでよく使用されている。グラフカットは、主として下記数式（１２）の形のエネルギーを最小化するために使われる方法である。
【００６３】
【数１２】
0005590548.tif 000014
【００６４】
ここで、ｚは二値もしくは多値の離散変数であり、ｚのとりうる値の集合をＬとする。ｇｉはｚｉの関数であり、ｈｉｊはｚｉとｚｊの関係を表す項である。Ｅは隣接する変数ｚのセットを表すインデックス集合である。α−拡張アルゴリズムではｈｉｊが下記数式（１３）を満たす必要がある。
【００６５】
【数１３】
0005590548.tif 000015
【００６６】
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおいては、ｈｉｊは上記数式（７）の第二項目に相当し、Ｊｋｉｎｔｅｒ＞０であれば上記数式（１３）が成立するため、α−拡張アルゴリズムを用いることができるのである。
【００６７】
（ベイズ推定）
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおいて用いるベイズ推定では、同時事後確率ｐ（ｘ，ｚ｜ｙ）は試験分布ｑ（ｘ，ｚ）と呼ばれる分布によって評価される。試験分布ｑ（ｘ，ｚ）は、ＫＬ距離と呼ばれる指標を最小化するものとして決定される。ここでのＫＬ距離は下記数式（１４）で定義される。
【００６８】
【数１４】
0005590548.tif 000016
【００６９】
ここで、＜・＞ｑ（ｘ，ｚ）は分布ｑ（ｘ，ｚ）に関して積分計算を行うことを表す演算子である。ＫＬ距離は常に非負であり、ｑ＝ｐのときにのみ０（ゼロ）になる。計算困難性に対処するために、計算可能な分布の集合の中において事後分布に最も近いものを探索する。
試験分布ｑ（ｘ，ｚ）は、上記数式（１４）の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。上記数式（１４）の最小化が可能であるような分布としては、ｑ（ｘ，ｚ）＝Π_iｑ（ｘi|{ｚ}_{Ｎ(i)}）Πiｑ（ｚi）がある。ただし、{ｚ}_{Ｎ(i)}は画素ｉ近傍の画素における組織クラス変数ｚの集合であり、ｑ（ｘi|{ｚ}_{Ｎ(i)}）はガウス分布とする。
ここでは、簡単のため下記数式（１５）に示されるように試験分布ｑ（ｘ，ｚ）に対して因子化仮定と呼ばれる変数ｘとｚの間の独立性の仮定を行う。さらに、ｘｊに関する分布ｑ（ｘｊ）がガウス分布になるものと仮定して、下記数式（１６）に示されるようにモデル化する。これらの仮定の下で、交互最適化を用い最適試験分布を求めている。
【００７０】
【数１５】
0005590548.tif 000017
【００７１】
【数１６】
0005590548.tif 000018
【００７２】
（最適試験分布）
ここで、最適試験分布について説明する。試験分布ｑ（ｘｊ）の平均ｕｊと分散σ２ｊは、ｑ（ｚｊ）を固定した下で、ＤＫＬを最小にするものとして、下記数式（１７）のように得ることができる。
【００７３】
【数１７】
0005590548.tif 000019
【００７４】
これらの最適化は、ＳＣＧ法によって行う。試験分布ｑ（ｘｊ）が推定された後、すなわち平均ｕｊと分散σ２ｊが求まった下で、組織クラスに関する最適試験分布ｑ＊（ｚｊ）は、ＤＫＬを最小にするものとして、下記数式（１８）と数式（１９）に示すように解析的に求まる。なお、下記数式（１９）におけるρｊｋは、下記数式（２０）で示されるものである。
【００７５】
【数１８】
0005590548.tif 000020
【００７６】
【数１９】
0005590548.tif 000021
【００７７】
【数２０】
0005590548.tif 000022
【００７８】
上記数式（２０）において、第一項目と第二項目はクラス事前分布に由来する項であり、第三項目は条件付き事前分布の正規化項から生じたものである。また、最後の項は条件付き事前分布に由来しており、推定された吸収係数ｕｊとクラスごとに定められた吸収係数ｖｊとを比較するものである。このｑ（ｚｊ）に関する最適化は、組織に関して組織のクラスをやわらかいクラスタリングによって決めることに対応している。
また、上記数式（２０）における第二項目は、同じクラスに属する隠れ変数がつながりやすいことを表現しており、再構成画像の滑らかさは隠れ変数を通して実現されている。
【００７９】
ベイズ推定のアルゴリズムの擬似コードを、図５に示す。なお、アルゴリズムは事前に決められた所定の収束条件を満たすまで繰り返される。
【００８０】
以上、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムにおける処理内容について説明した。以下の実施例では、観測対象としてファントムデータを用いた計算機実験を行い、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの画面再構成を、従来技術による画面再構成と比較して評価する。
下記の評価実験から、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムの有用性が理解できる。
【００８１】
実施例１として、Ｘ線被曝を極力抑えた状況を想定し、限られた放射線量の下での実験を行う。実施例１の結果から、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムが、限られたデータからの画像再構成に有効働くことが示される。
そして、実施例２では、投影数は十分得られるが金属の含まれるデータからの断層画像再構成を行う。実施例２の結果から、本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムが、メタルアーティファクト除去に有効に働くことが示される．
【実施例１】
【００８２】
実施例１は、Ｘ線被曝を極力抑えた状況を想定し、限られた放射線量の下での実験を行った。
現在広く用いられている医療用ＣＴでは、検出器の数は７００〜９００程度であり、投影数は８００〜１５００程度である。それらのデータから、２５６×２５６もしくは５１２×５１２のピクセルサイズの断層画像が生成される。Ｘ線吸収係数ｘは、Hounsfield unit (ＨＵ) ＝１０００（ｘ−ｘ０）／ｘ０を用いて表される。ここで、ｘ０は水のＸ線吸収係数である。ＨＵ値は、水のＸ線吸収係数を０、空気のＸ線吸収係数を−１０００として正規化した値である。断層画像の表示を全ての階調で行うと十分なコントラストが得ることができないため、実験ではウインドウ幅を[−500，500] ＨＵ に統一して表示を行った(Ｘ線吸収係数では[0.01，0.03] に相当する。)。
【００８３】
一般に、Ｘ線ＣＴによって再構成される断層画像の精度は、放射線量や放出光子の量によって変化し、再構成画像に生じるノイズと放射線量には下記数式（２１）の関係が成り立つ。
【００８４】
【数２１】
0005590548.tif 000023
【００８５】
ここで、Ｎは画像に生じるノイズを表している。また、Ｂは対象のＸ線透過率，Ｄは入射線量，ｈはスライス厚，ｗはピクセルサイズをそれぞれ表している。放射線の量を増やせば画像に生じるノイズは減少するが、放射線量にしたがって線形に減少するわけではない。これは光子ノイズが、ポワソン分布にしたがって変動するためであると推察する。
【００８６】
本発明のＸ線ＣＴ画像処理方法およびＸ線ＣＴ画像処理プログラムでは、観測対象の物質とＸ線吸収係数を事前知識として用いているが、Ｘ線吸収係数は個体差やＸ線のエネルギーなどによって変化することから、実際の観測対象の吸収係数は予め完全に分かっているわけではない。こうした吸収係数の変動を考慮して、実験では吸収係数の平均値が想定したものから大きくずれた状況を想定する。吸収係数に関する事前分布の分散をｒ２＝１０−５に設定して、観測対象の真のＸ線吸収係数を事前分布の平均ｖｋからずらしたものを二通り(ＣａｓｅＡ, ＣａｓｅＢ)を想定した（図６を参照）。なお、空気と金属のＸ線吸収係数は事前知識とずれはないものと想定し、事前分布の平均値と一致させる。
【００８７】
１つ目の実験設定（ＣａｓｅＡ）では、筋肉，脂肪，骨のＸ線吸収係数を事前分布の平均から標準偏差ｒだけＸ線吸収係数が離れる方向にずらした。これは１５８ＨＵ程度のずれであり、実際のＸ線吸収係数の個体差に比べると十分大きな値であると考えられる。２つ目の実験設定（ＣａｓｅＢ）では、筋肉と脂肪のＸ線吸収係数が近くなる方向に５０ＨＵずつ事前分布の平均から移動させた。こちらも実際の変動よりは十分に大きい値であると考えられる。
【００８８】
設定すべきパラメータとしては、上述したように、Ｘ線吸収係数に関する事前分布の平均ｖｋと分散ｒ２ｋ、組織に関する事前分布のＪｋｓｅｌｆ，Ｊｋｉｎｔｅｒがある。
これらは、実際の観測データセットから決められるべきパラメータである。例えば、肺では太さの異なる血管が張り巡らされており、Ｘ線吸収係数の分布は一様ではない。一方水などの均質な物質では吸収係数は一定になる。これらのことから、吸収係数の分布や組織ごとの繋がりやすさは各々の組織に依存することになる。
本実験では、これらのパラメータに関するロバストさを調べるために、前記表１のように実験設定によらず固定した。
【００８９】
実施例１の実験で用いたファントムを図７（ａ）（ｂ）に示す。それぞれのファントムの大きさは、４７１×３５３ｍｍである。Ｘ線発生源から並行ビームを発生させて、観測対象を通過した光子数を検出器によりカウントした。検出器の数は３６７個であり、１８０°の間に３６回投影を行いデータを獲得した。光子数ｂｉ（ｔ）は１０５とした。再構成する画像の解像度は、２５６×２５６ピクセルとする。比較実験として、フィルター補正逆投影法（ＦＢＰ）と最尤法（ＭＬ）による画像再構成を同時に行った。
ＭＡＰ推定とベイズ推定の画像再構成の結果を、図８ (ＣａｓｅＡの実験設定)，図１０(ＣａｓｅＢの実験設定)に示す。それぞれの図で、(ａ) はＦＢＰによって得られた再構成画像（ＦＢＰ）、 (ｂ) はＭＬによって得られた再構成画像（ＭＬ），(c)は本発明のＸ線ＣＴ画像処理方法（ＭＡＰ推定）によって得られた再構成画像（ＭＡＰ）、(ｄ) は本発明のＸ線ＣＴ画像処理方法（ベイズ推定）によって得られた再構成画像(ＢＡＹＥＳ) である。
【００９０】
ＭＡＰ推定とベイズ推定で得られた再構成画像を、ＰＳＮＲにより評価した。それぞれの結果は、図中に記載の通りである。ＣａｓｅＡ、ＣａｓｅＢいずれも、本発明のＸ線ＣＴ画像処理方法（ベイズ推定）を用いた再構成画像が最も高いＰＳＮＲを示していることがわかる(ＣａｓｅＡ：20.32 ｄＢ, ＣａｓｅＢ：
21.60 ｄＢ)。
図に示されるように、他の手法では投影数が制限されていることにより、ノイズが多く含まれた結果が得られている。また、本発明のＸ線ＣＴ画像処理方法のＭＡＰ推定とベイズ推定によって得られた組織クラス推定の結果を、図９(ＣａｓｅＡの実験設定)と図１１(ＣａｓｅＢの実験設定)に示す。本発明のＸ線ＣＴ画像処理方法により、組織クラスの推定が可能であることが理解できる。図において、上段（ａ）〜（ｅ）はＭＡＰ推定による結果であり、下段（ｆ）〜（ｊ）はベイズ推定による結果である。それぞれの組織に属していると推定されたピクセルは白色で表示されている。
【００９１】
図１２は、投影数と誤差の関係を調べるために、投影枚数ごとの各手法によるＰＳＮＲを計測したものである。対象物体とパラメータは、ＣａｓｅＡの実験設定と同じものを用いるが、計算コストを下げるために再構成画像のサイズは６４×６４ピクセルとしている。各手法での計測を６回ずつ行い、平均値と標準偏差を、エラーバーを用いてプロットした。ほとんどの場合で、本発明のＸ線ＣＴ画像処理方法が従来技術による画像再構成を上回っており、本発明のＸ線ＣＴ画像処理方法が低放射線下でのＸ線ＣＴ画像再構成に有効に作用していることが理解できる。
【００９２】
また、本発明のＸ線ＣＴ画像処理方法のＭＡＰ推定とベイズ推定での結果を比較すると、わずかではあるが、ほとんどの場合でベイズ推定が上回った結果が得られている。パラメータの設定の違いがあるために一概に述べることはできないが、本結果はベイズ推定では再構成画像と所属の組織クラスの不確かさまで考慮した推定を行っているためであると推察する。
なお、本発明のＸ線ＣＴ画像処理方法のＭＡＰ推定とベイズ推定と従来技術（ＭＬ）の処理速度については、図１３に示すように、再構成画像のサイズは６４×６４ピクセルや１２８×１２８ピクセルの場合は大差なく、２５６×２５６ピクセル程度以上の解像度になると、ベイズ推定の処理コストが大きくなることがわかる。
【実施例２】
【００９３】
実施例２では、投影数は十分得られるが金属の含まれるデータからの断層画像再構成を行う。実施例２の実験で用いた金属が含まれるファントムを図７（ｃ）（ｄ）に示す。対象となる画像は人間の頭(それぞれ大きさ４７２×４３６ｍｍ) を模擬したデータであり、歯にあたる部分にそれぞれ(１８，１９，２３ｍｍ) の大きさの金属が埋め込まれているものである。対象物体にＸ線平行ビームを投射して、反対側にある３６７個の検出器で光子の数を検出した。投影は１８０°の間に７２回データをサンプルした。Ｘ線発生器で発生させる光子数ｂｉ（ｔ）は１０６とした（実施例１の実験よりも光子数は多い）。実施例１と同様に、再構成する画像の解像度は、２５６×２５６ピクセルとする。また、比較実験として、フィルター補正逆投影法（ＦＢＰ）と最尤法（ＭＬ）による画像再構成を同時に行った。
ＭＡＰ推定とベイズ推定の画像再構成の結果を、図１４(ＣａｓｅＡの実験設定)，図１６(ＣａｓｅＢの実験設定)に示す。それぞれの図で、(ａ) は真の断層画像であり、（ｂ）はＦＢＰによって得られた再構成画像（ＦＢＰ）、 (ｃ) はＰＣＬＩＳによって得られた再構成画像（ＰＣＬＩＳ），(ｄ) はＭＬによって得られた再構成画像（ＭＬ），(ｅ)は本発明のＸ線ＣＴ画像処理方法（ＭＡＰ推定）によって得られた再構成画像（ＭＡＰ）、(ｆ) は本発明のＸ線ＣＴ画像処理方法（ベイズ推定）によって得られた再構成画像(ＢＡＹＥＳ) である。
【００９４】
ＭＡＰ推定とベイズ推定で得られた再構成画像を、ＰＳＮＲにより評価した。それぞれの結果は、図中に記載の通りである。ＣａｓｅＡ、ＣａｓｅＢいずれも、本発明のＸ線ＣＴ画像処理方法（ベイズ推定）を用いた再構成画像が最も高いＰＳＮＲを示していることがわかる(ＣａｓｅＡ：33.10 ｄＢ, ＣａｓｅＢ:
34.34 ｄＢ)。
また、ＦＢＰやＰＣＬＩＳでは、メタルアーティファクトの影響で、画像に筋が入っていたり、ノイズの多い画像になっている。しかし、本発明のＸ線ＣＴ画像処理方法ではメタルアーティファクトが低減されており、真の画像に近い画像が得られていることがわかる。また、本発明のＸ線ＣＴ画像処理方法であるＭＡＰ推定とベイズ推定によって得られた組織クラス推定の結果を、図１５(ＣａｓｅＡの実験設定)と図１７(ＣａｓｅＢの実験設定)に示す。本発明のＸ線ＣＴ画像処理方法であるＭＡＰ推定とベイズ推定によって得られた図を比較すると、ベイズ推定によって得られた画像(クラス推定を含め) の方がよりよい結果が得られていることがわかる。図において、上段（ａ）〜（ｅ）はＭＡＰ推定による結果であり、下段（ｆ）〜（ｊ）はベイズ推定による結果である。それぞれの組織に属していると推定されたピクセルは白色で表示されている。
【００９５】
上述の説明では、本発明のＸ線ＣＴ画像処理方法では、隠れ変数を導入し組織に関する事前知識を用いることで、少ない投影データにおいて、従来技術より高精細な画像再構成が可能なこと及びメタルアーティファクトを低減できることを示した。また、実施例の実験では、本発明のＸ線ＣＴ画像処理方法による再構成によって精度の高い推定が行われていることを示した。
【００９６】
昨今のＸ線ＣＴの主な利用先である医療診断では、診断補助としてＣＴ画像を組織ごとに領域分割するための研究が多くなされている。これらの方法はいずれもＦＢＰなどを用いて得られたＣＴ画像から領域分割を行っている。領域分割と画像の再構成は互いに深く関わっているため、別々に行うよりも統一的な枠組みの下で推定を行う方が好ましい。本発明のＸ線ＣＴ画像処理方法では、観測データから統一的な枠組みの下で画像再構成と領域分割が実行されているため、副次的にせよ精度の良いセグメント化がなされていると考える。本発明のＸ線ＣＴ画像処理方法は、医療診断補助などへの応用も大いに期待できるものである。
【産業上の利用可能性】
【００９７】
本発明は、医療用のみならず産業用のＸ線ＣＴ装置に有用である。
【図１】
0005590548.tif 000024
【図２】
0005590548.tif 000025
【図３】
0005590548.tif 000026
【図４】
0005590548.tif 000027
【図５】
0005590548.tif 000028
【図６】
0005590548.tif 000029
【図７】
0005590548.tif 000030
【図８】
0005590548.tif 000031
【図９】
0005590548.tif 000032
【図１０】
0005590548.tif 000033
【図１１】
0005590548.tif 000034
【図１２】
0005590548.tif 000035
【図１３】
0005590548.tif 000036
【図１４】
0005590548.tif 000037
【図１５】
0005590548.tif 000038
【図１６】
0005590548.tif 000039
【図１７】
0005590548.tif 000040