TOP > 国内特許検索 > X線CT画像処理方法,X線CTプログラムおよび該プログラムが搭載されたX線CT装置 > 明細書

明細書 :X線CT画像処理方法,X線CTプログラムおよび該プログラムが搭載されたX線CT装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第5590548号 (P5590548)
公開番号 特開2011-156302 (P2011-156302A)
登録日 平成26年8月8日(2014.8.8)
発行日 平成26年9月17日(2014.9.17)
公開日 平成23年8月18日(2011.8.18)
発明の名称または考案の名称 X線CT画像処理方法,X線CTプログラムおよび該プログラムが搭載されたX線CT装置
国際特許分類 A61B   6/03        (2006.01)
A61B   6/00        (2006.01)
FI A61B 6/03 360D
A61B 6/03 360J
A61B 6/00 350D
請求項の数または発明の数 9
全頁数 32
出願番号 特願2010-022623 (P2010-022623)
出願日 平成22年2月3日(2010.2.3)
審査請求日 平成25年1月30日(2013.1.30)
特許権者または実用新案権者 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
発明者または考案者 【氏名】前田 新一
【氏名】石井 信
【氏名】福田 航
【氏名】兼村 厚範
個別代理人の代理人 【識別番号】110000822、【氏名又は名称】特許業務法人グローバル知財
審査官 【審査官】亀澤 智博
参考文献・文献 特開2008-062035(JP,A)
特表2009-518143(JP,A)
特表2005-536308(JP,A)
国際公開第2009/008102(WO,A1)
特開2007-111526(JP,A)
特開2007-044520(JP,A)
特開2004-174253(JP,A)
特開平7-057533(JP,A)
調査した分野 A61B 6/00 - 6/14
特許請求の範囲 【請求項1】
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
前記事前知識は、
再構成画像の各画素の領域において定義されるパラメータで、
撮像対象の各組織クラスの存在割合を表現するパラメータと、
各組織クラスのX線吸収係数を表現するパラメータと、
各組織クラスの空間的に連続する度合いを表現するパラメータと、
によって特徴付けられる確率分布で表現され、
X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用い、
前記統計推定は、
前記事前知識と尤度から計算される事後分布を用いて、
事後確率最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定により画像再構成及び組織クラス推定を行う、
ことを特徴とするX線CT画像処理方法。
【請求項2】
前記X線吸収係数を表現するパラメータに関する分布としてガウス分布を用い、ガウス分布における平均の値は組織クラスごとに設定されることを特徴とする請求項1に記載のX線CT画像処理方法。
【請求項3】
前記MAP推定は、再構成画像の各画素の領域が属する組織クラスとX線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定することを特徴とする請求項1又は2に記載のX線CT画像処理方法。
【請求項4】
前記ベイズ推定は、再構成画像の各画素の領域が属する組織クラスとX線吸収係数の両方に関する事後分布の推定を行うことを特徴とする請求項1又は2に記載のX線CT画像処理方法。
【請求項5】
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
3)前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップと、
を備えたことをX線CT画像処理方法。
【請求項6】
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
3)前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定(ベイズ推定)により画像再構成及び組織クラス推定を行うステップと、
を備えたことをX線CT画像処理方法。
【請求項7】
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
3)前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのX線CT画像処理プログラム。
【請求項8】
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)前記投影画像に観測されるフォトンに関して物理過程を確率分布で表現するステップと、
3)前記事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定(ベイズ推定)により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのX線CT画像処理プログラム。
【請求項9】
請求項7又は8に記載のX線CT画像処理プログラムを搭載したX線CT画像装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、体内組織構成を考慮したX線CTアルゴリズム、そのアルゴリズムを有するプログラムおよび該プログラムが搭載されたX線CT装置に関する。
【背景技術】
【0002】
従来から、メタルアーティファクトの低減技術や、より低いX線被曝量で従来と同程度のCT画像を再構成する技術など、X線CTの高精度化を可能にするアルゴリズムが研究されている。一般に、X線CTは、X線管にX線検出器を対向配置するとともに、これらの間にターンテーブルを配置した構成を取る。このターンテーブルの回転軸は、X線管とX線検出器を結ぶX線光軸に対して直交する向きとされる。X線CT装置においては、測定対象物を360°回転させて透過X線データを収集しなければCT断層画像の再構成を行うことができない。
【0003】
メタルアーティファクトは、測定対象物にX線を強く吸収する金属などの密度の高い物質があると、再構成画像にノイズがのるというものである。すなわち、鉄などの比較的X線が透過しにくい測定対象物のCT撮像を行った場合、360°回転させながら透過X線データを採取すると、密度の高い物質がX線の透過方向に重なった状態での透過X線データが得られることになって、フィルター補正逆投影法(FBP)などの既存手法を用いた場合、再構成画像にはメタルアーティファクトと呼ばれる虚像が生じ、正確な画像の再構成ができないという問題がある。
従って、メタルアーティファクトの低減によって、X線CTの高精度化が図れることになる。
【0004】
一方、X線は強力であれば、再構成画像のSN比が高いことから、X線被曝量の低減とSN比はトレードオフの関係がある。より低いX線被曝量で従来と同程度のSN比の再構成画像を得る技術とは、従来と同じX線強度,撮像枚数でも解像度の高い画像が得られるという解像度向上を図る技術や、従来と同じX線強度や解像度でも、より撮像枚数を減らして撮像時間を短縮できる撮像時間の短縮を図る技術や、従来と同じ解像度でもX線強度を弱めて被曝量を減らせる被曝量の低減を図る技術、従来と同じX線強度や解像度でも再構成画像に含まれるアーティファクト(もしくはノイズ)を低減できるノイズ除去を図る技術、従来と同じX線強度や解像度でも再構成画像のコントラスト分解能(再構成画像を何階調で表せるか、すなわち、X線吸収係数の精度)を向上するコントラスト分解能の向上を図る技術がある。
【0005】
上記の従来技術は、統計推定を行うものと統計推定を行わないものに大きく2つに分類できる。
先ず、統計推定を行わない従来技術としては、画像再構成法の主流であるフィルター補正逆投影法(FBP:Filtered
Back Projection)や、FBPを改良したPCLIS(Projection Completion
Method based on a Linear Interpolation in the Sinogram)が挙げられる。FBPは、撮像対象に対して色んな方向からX線を照射して、その投影像(X線の吸収像、投影像のセットをシノグラムと呼ぶ)を得た後、得られた投影像を逆方向に投射して画像再構成する際に、シノグラムに周波数領域上でフィルターをかけてボケ除去を行う画像再構成法である。PCLISは、メタルなどのX線を強く吸収する物体によってほとんどX線を観測できなかったシノグラムの部分に対してシノグラム上で線形補間を行ってからFBPを行うことでメタルアーティファクトなどの極端な虚像の発生を抑えるものである。
次に、統計推定を行なう従来技術としては、再構成画像に関する事前知識をおかない最尤推定法(ML)がある。最尤推定法(ML)は、観測における不確実性を含む物理過程を確率モデルで表現して、その確率モデルを基に最も尤もらしい再構成画像を推定するものである。観測における不確実性を含む物理過程には、観察されるフォトンに関する揺らぎ(ショットノイズ)などが含まれる。
【0006】
また、撮像対象の物体を異なる位置・角度から撮影した投影画像を複数枚用意できる場合、それら複数の投影画像の情報を処理することで、元の物体の断面像を再構成できることが知られている。投影画像を使用して元の物体の断面像を復元する画像再構成法のうち統計的な解法に関する従来の研究は、最尤推定法(ML),MAP(maximum a posteriori)推定法,ベイズ推定法の3手法が存在する。
【先行技術文献】
【0007】

【非特許文献1】”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.
【非特許文献2】” Metalartifact reduction in CT using tissue-class modeling and adaptive prefiltering”, Medical Physics, Vol. 33, No. 8,August 2006.
【非特許文献3】” ベイズアプローチに基づいた断層画像の再構成”,情報処理学会研究報告(2009年)
【発明の概要】
【発明が解決しようとする課題】
【0008】
しかしながら、上記の従来技術を用いた画像再構成の場合は、従来と同程度の再構成画像を得るためには従来と同程度のX線被曝量となり、またメタルアーティファクトの低減が困難であった。
【0009】
そこで、本発明は、従来技術を用いた画像再構成と比べて、より低いX線被曝量で従来と同程度の再構成画像が取得可能で、メタルアーティファクトを低減可能なX線CT画像処理方法、X線CT画像処理プログラムおよび該プログラムが搭載されたX線CT装置を提供することを目的とする。
【課題を解決するための手段】
【0010】
上記状況に鑑みて、本発明のX線CT画像処理方法は、
X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
(1)観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現され、
(2)上記の再構成画像に関する事前知識は、
再構成画像の各画素の領域において定義されるパラメータで、
撮像対象の各組織クラスの存在割合を表現するパラメータと、
各組織クラスのX線吸収係数を表現するパラメータと、
各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現され、
X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用い、
(3)上記の統計推定は、
前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定であり、
画像再構成及び組織クラス推定を行うことを特徴とする。
【0011】
上記のX線CT画像処理方法によれば、観測されるフォトンに関する確率モデルがポアソン分布などで表現されることにより、従来と比べてより現実の観測過程に近い物理過程を表現することができ、フォトンノイズなどによる観測の不確定性を表現できる。
また、上記のX線CT画像処理方法によれば、人体などの撮像対象の組織分布に関する事前知識について、どのような組織(筋肉などの通常細胞、脂肪などの軟細胞、骨、メタルなど)がどの程度の割合で分布するか、また、それぞれの組織がどの程度のX線吸収係数を有するか、また、それぞれの組織がどの程度、空間的に連続して分布しやすいかの分布情報を確率分布の形で表現して、事後確率最大化によるMAP推定、或いは事後確率の期待値によるベイズ推定統計的推論を用いて、画像再構成や組織推定を行う。組織分布に関する事前知識は、ある固定した平均的なパラメータを用いても良いが、体格や既往歴、性別、年齢などの個人差や撮像部位に応じて適切に変化させることで、画像再構成と組織推定の精度のさらなる向上を図ることが可能である。
【0012】
ここで、本発明のX線CT画像処理方法において、上記の各組織クラスのX線吸収係数を表現するパラメータに関する分布としてガウス分布を用い、かつ、各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布としてボルツマン分布を用いることが好ましい態様である。
【0013】
各組織クラスのX線吸収係数を表現するパラメータに関する分布としてガウス分布を用いることとした理由は、組織クラスごとに決まる特定のX線吸収係数(CT値)をとりやすいことに着目したものである。
各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布としてボルツマン分布を用いることとした理由は、同じ組織クラスが空間的に集まりやすい(組織クラス別に集まりやすさの調整が可能)ことと、標準とされる各組織クラスの割合を表現しやすいことに着目したものである。
【0014】
また、本発明のX線CT画像処理方法において、上述のMAP推定は、再構成画像の各画素の領域が属する組織クラスとX線吸収係数の両方に関する事後分布において最も尤もらしい組合せを推定することが好ましい。
【0015】
再構成画像の各画素の領域が属する組織クラスとの組合せを推定することは、組織クラス依存性を含めることになる。この組織クラス依存性を含めることで、例えば、空気の組織クラスはピクセル同士の連結がしやすい(空気中に他の組織は入らない)などの組織クラスに依存した知識を含めることができるのである。これにより、より柔軟性をもった分布表現が可能となる。
【0016】
なお、推定精度を向上すべく、例えば、グラフカットの一種であるα-拡張を用いた最適化アルゴリズムによってMAP推定することがより好ましい。
【0017】
また、本発明のX線CT画像処理方法において、上述のベイズ推定は、X線吸収係数と再構成画像の各画素の領域が属する組織クラスの両方に関する事後分布を推定するにあたり、事後分布をよく近似可能な試験分布を用いて近似する。
なお、ベイズ推定の場合も、MAP推定の場合と同様に、組織クラス依存性を含めることができる。これにより、より柔軟性をもった分布表現が可能となる。
【0018】
また、本発明の他の観点によれば、本発明のX線CT画像処理方法は、X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
3)再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップと、
を備えた構成とされる。
【0019】
かかる構成によれば、画像の再構成の速度を速めることができ、より高速でX線CT画像処理を行うことができる。
【0020】
また、本発明の他の観点によれば、本発明のX線CT画像処理方法は、X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理方法であって、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
3)再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定(ベイズ推定)により画像再構成及び組織クラス推定を行うステップと、
を備えた構成とされる。
【0021】
かかる構成によれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にX線CT画像処理を行うことができる。
【0022】
また、本発明の他の観点によれば、X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
3)再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率最大化推定(MAP推定)により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのプログラムが提供される。
【0023】
かかるプログラムによれば、画像の再構成の速度を速めることができ、より高速でX線CT画像処理を行うことができる。
【0024】
また、本発明の他の観点によれば、X線吸収係数に関する事前知識を用いて、画像再構成及び組織クラスの統計推定を行うX線CT画像処理プログラムであって、
コンピュータに、
1)投影画像と少なくともX線強度を含む計測条件を入力するステップと、
2)投影画像に観測されるフォトンに関してポアソン分布などの物理過程を確率分布で表現するステップと、
3)再構成画像に関する事前知識として、再構成画像の各画素の領域において定義されるパラメータであって、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、X線吸収係数の事前分布は、組織クラスを隠れ変数とする事前分布を用い、組織クラスの事前分布は、各組織クラスの存在割合を表現するパラメータに関する項と、各組織クラスの空間的に連続する度合いを表現するパラメータに関する項との和で表されるエネルギー関数をもつボルツマン分布を用いるステップと、
4)前記事前知識と尤度から計算される事後分布を用いて、事後確率の期待値推定(ベイズ推定)により画像再構成及び組織クラス推定を行うステップと、
を実行させるためのプログラムが提供される。
【0025】
かかるプログラムによれば、推定の不確実性を考慮に入れた推論が可能で、またパラメータの学習が容易であり、より正確にX線CT画像処理を行うことができる。
【0026】
また、本発明は、上記のX線CT画像処理プログラムを搭載したX線CT画像装置を提供することができる。
【発明の効果】
【0027】
本発明のX線CT画像処理方法およびX線CT画像処理プログラムによれば、従来技術を用いた画像再構成と比べて、より低いX線被曝量で従来と同程度の再構成画像が取得可能で、メタルアーティファクトを低減可能であるといった効果を有する。
【図面の簡単な説明】
【0028】
【図1】本発明のX線CT画像処理方法(MAP推定)の処理フロー
【図2】本発明のX線CT画像処理方法(ベイズ推定)の処理フロー
【図3】X線CTの説明図
【図4】MAP推定を用いたアルゴリズムの擬似コード
【図5】ベイズ推定を用いたアルゴリズムの擬似コード
【図6】事前知識として用いるCT値の分布を示す図
【図7】実験に用いるファントムの説明図
【図8】実施例1のCaseAの実験設定での画像再構成を示す図
【図9】実施例1のCaseBの実験設定での画像再構成を示す図
【図10】実施例1のCaseAの実験設定での組織クラスの推定結果を示す図
【図11】実施例1のCaseBの実験設定での組織クラスの推定結果を示す図
【図12】投影画像数とPSNRの関係を示す図
【図13】本発明のX線CT画像処理方法の処理時間の説明図
【図14】実施例2のCaseAの実験設定での画像再構成を示す図
【図15】実施例2のCaseBの実験設定での画像再構成を示す図
【図16】実施例2のCaseAの実験設定での組織クラスの推定結果を示す図
【図17】実施例2のCaseBの実験設定での組織クラスの推定結果を示す図

【発明を実施するための最良の形態】
【0029】
以下、本発明の実施形態について、図面を参照しながら詳細に説明していく。ただし、本発明の範囲は、以下の実施例や図示例に限定されるものではなく、幾多の変更または変形が可能である。
【0030】
本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、X線吸収係数に関する事前知識を用いて組織クラスの統計推定を用いて画像の再構成を行う。特に、放射線による被曝を避けるためにより少ない観測データからの画像再構成を試みる。その一方で、画像再構成は不良設定性を有しているため、適当な事前知識を用いて解に何らかの制約を課すことで解決を試みる。本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、人体は脂肪,筋肉,骨といった限られた物質から構成され、それぞれのX線吸収係数の大まかな分布はあらかじめ分かっているものとする。また、それぞれの組織クラスは空間的な連続性を有しており、各組織クラスが人体に占める割合も大まかに分かっているという状況を想定する。
【0031】
すなわち、再構成画像に関する事前知識について、再構成画像の各画素の領域において定義されるパラメータで、撮像対象の各組織クラスの存在割合を表現するパラメータと、各組織クラスのX線吸収係数を表現するパラメータと、各組織クラスの空間的に連続する度合いを表現するパラメータと、によって特徴付けられる確率分布で表現し、これらの知識を階層ベイズモデルにおける事前知識として利用する。
【0032】
そして、事後分布推定による統計推定を行う。画像再構成において、事後分布を求めることが必要であるが、高次元の隠れ変数に関する積分計算が含まれるため、解析的に実行することは困難であることから、本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、事後確率最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定を用いた近似手法を適用することでこの計算困難さを克服することにしたものである(図1,図2のフローチャートを参照)。
【0033】
(問題の定式化)
先ず、X線CT画像のように、様々な方向から投影されたT個の投影データを、D={Y(1),・・・,Y(T)}と表すことにする。各々のデータY(t)は、t番目の投影によって検出器で検出されるデータの集合であり、Y(t)={y(1),・・・,y(t)}となる。y(t)は、i番目の検出器で検出された光子の数を表す。X線は物質を透過する際に指数的に減衰することから、下記数式(1)のように表すことができる。
【0034】
【数1】
JP0005590548B2_000002t.gif

【0035】
ここで、xは観測対象を撮像対象のX線吸収係数をラスタスキャンして得られるJ次元ベクトルx={x,・・・,x}のj番目のピクセルのX線吸収係数である。b(t)は、X線源から放出される光子数(物体が何も置かれていないときに観測され得る光子数)を表す。また、lij(t)は角度θ(t)から投影された際のi番目の検出器で検出される投影線とj番目のピクセルが交差する距離に相当するものであり、lij(t)xがt番目の投影のi番目の検出器に入射するX線に対するj番目のピクセル領域の実効的なX線吸収係数を表す(図3を参照)。
【0036】
(階層ベイズモデル)
本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおけるX線CT画像の再構成では、X線吸収係数xは観測データDの事後分布として推定される。本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、事後分布の平均値としてX線吸収係数xが得られるとする。事後分布はベイズの定理により下記数式(2)で求めることができる。
【0037】
【数2】
JP0005590548B2_000003t.gif

【0038】
また、人間の体は限られた数の組織から成っており、組織ごとにX線吸収係数xが定まるという事前知識を用いるため、X線吸収係数xに関して階層的な事前分布を定義することにする。説明の便宜上、以下では、人体の組織は5(=K)つの組織(筋肉,脂肪,骨,空気,金属)に分類されると仮定して、事前分布p(x)を隠れ変数z={z,・・・,z}を用いて下記数式(3)のように表わす。
【0039】
【数3】
JP0005590548B2_000004t.gif

【0040】
ここで、zはK次元の二値変数であり、各要素zjkはj番目のピクセルがk番目のクラスに属しているときに1をとり、それ以外の時は0をとるものである。事前分布において、すべてのピクセルはいずれかのクラスに属し、Σjk=1を満たすものとする。
【0041】
次に、本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおいて、投影画像に観測されるフォトンに関する物理モデル、事前知識としての組織クラスに関する事前分布、すなわち、撮像対象の各組織クラスの存在割合を表現するパラメータに関する分布、各組織クラスのX線吸収係数を表現するパラメータに関する分布、各組織クラスの空間的に連続する度合いを表現するパラメータに関する分布について、以下説明を行う。

【0042】
(観測されるフォトンに関する物理モデル)
X線CTにおける観測データの主なノイズの要因は、検出される光子(フォトン)の量子化ノイズであると考えられる。そこで、光子(フォトン)の観測データは、投影毎、検出器毎に独立なポアソン分布に従って生成されるとしてモデル化する。
なお、ポアソン分布に特に限定されるものではなく、物理過程をより良く表現できる他の物理モデルも適用可能である。
上述の数式(1)は、ポアソン分布の平均を表しており、数式(1)の結果を用いて、下記数式(4)で表わす。
【0043】
【数4】
JP0005590548B2_000005t.gif

【0044】
(組織に関する事前分布)
所属する組織クラスの情報が与えられた下で、X線吸収係数xはガウス分布に従うとして、下記数式(5)に示すようにモデル化する。
【0045】
【数5】
JP0005590548B2_000006t.gif

【0046】
ここで、vとrはガウス分布の平均と分散を表している。平均vの値は、各々の組織(空気,脂肪,筋肉,骨,金属)ごとにそれぞれ、v=0.000,v=0.018,v=0.022,v=0.040,v=0.120とする。
また、分散rは、下記表1に示すように、どの組織に所属するかに関わらず一定とする。なお、各々の平均vは、対象となる物体の個体差や温度等で変動すると考えられる。これらの不確かさはガウス分布の分散rによって調整する。
【0047】
【表1】
JP0005590548B2_000007t.gif

【0048】
また、組織のクラスに関する事前分布はボルツマン分布に従うとして、下記数式(6)に示すようにモデル化する。
【0049】
【数6】
JP0005590548B2_000008t.gif

【0050】
上記数式(6)において、Zは正規化項であり、エネルギーEを下記数式(7)のように定義する。
【0051】
【数7】
JP0005590548B2_000009t.gif

【0052】
η(j)は、j番目のピクセルの近傍画素を表している。JselfとJinterは、事前分布をコントロールする非負の定数である。ボルツマン分布は、エネルギーE(z) が低ければ低いほどその確率が高くなる。エネルギー項の第一項目は、組織の分布する割合に関わり、第二項目は組織の空間的な連続性に関わる。
【0053】
(事後分布の推定)
上述した如く、X線CT画像の再構成において、事後分布を求めるのであるが、高次元の隠れ変数に関する和計算が含まれるため、解析的に実行することは困難である。このことから、本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、事後確率最大化による推定(MAP推定)、或いは、事後確率の期待値によるベイズ推定を近似的に行った手法を用いることで計算困難さを回避している。
【0054】
(MAP推定)
ベイズの定理によると、吸収係数xと組織クラスzに関する事後分布(これを前述の事後分布p(x|D)と区別すべく、同時事後分布と表現する。)p(x,z|D)は、それぞれの事前分布と尤度の積の形で、下記数式(8)のように表せる。このような同時事後分布p(x,z|D)の最大値を求めるようにすることで、高次元の隠れ変数に関する和計算が必要なくなるのである。
【0055】
【数8】
JP0005590548B2_000010t.gif

【0056】
上記数式(8)において、変数x、zを事後分布最大化基準(MAP)により決定する。下記数式(9)に示すように、上記数式(8)の負の対数をとったものを目的関数とし、これを最小化するものとして事後分布最大化を達成するように未知変数を求める。
【0057】
【数9】
JP0005590548B2_000011t.gif

【0058】
しかしながら、連続変数xと離散変数zに関する同時最適化は困難である。そのため、下記数式(10)と数式(11)に示す連続変数xと離散変数zに関する最適化は交互に最適化を行っている。
【0059】
【数10】
JP0005590548B2_000012t.gif

【0060】
【数11】
JP0005590548B2_000013t.gif

【0061】
ここで、連続変数xに関しては、SCG (Scaled conjugate gradient) 法を用いて最適化を行っている。また、離散変数zに関して、グラフカットの一種であるα-拡張アルゴリズムにより最適化している。
MAP推定のアルゴリズムの擬似コードを、図4に示す。なお、アルゴリズムは事前に決められた所定の収束条件を満たすまで繰り返される。
【0062】
(グラフカット)
ここでは、グラフカットの簡単に説明する。グラフカットとは、離散変数の最適化を行う際によく用いられるアルゴリズムである。一般の離散変数最適化問題は、NP困難であることが知られているが、劣モジュラ性と呼ばれる性質を満たす場合にグラフカットにより大域最小解が得られることが示されている。
本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、離散変数zの最適化としてグラフカットの一種であるα-拡張アルゴリズムを用いることにしている。このアルゴリズムは、大域最小解が得られることが保証されているわけではないが、適用できる範囲が広くコンピュータビジョンなどでよく使用されている。グラフカットは、主として下記数式(12)の形のエネルギーを最小化するために使われる方法である。
【0063】
【数12】
JP0005590548B2_000014t.gif

【0064】
ここで、zは二値もしくは多値の離散変数であり、zのとりうる値の集合をLとする。gはzの関数であり、hijはzとzの関係を表す項である。Eは隣接する変数zのセットを表すインデックス集合である。α-拡張アルゴリズムではhijが下記数式(13)を満たす必要がある。
【0065】
【数13】
JP0005590548B2_000015t.gif

【0066】
本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおいては、hijは上記数式(7)の第二項目に相当し、Jinter>0であれば上記数式(13)が成立するため、α-拡張アルゴリズムを用いることができるのである。
【0067】
(ベイズ推定)
本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおいて用いるベイズ推定では、同時事後確率p(x,z|y)は試験分布q(x,z)と呼ばれる分布によって評価される。試験分布q(x,z)は、KL距離と呼ばれる指標を最小化するものとして決定される。ここでのKL距離は下記数式(14)で定義される。
【0068】
【数14】
JP0005590548B2_000016t.gif

【0069】
ここで、<・>q(x,z)は分布q(x,z)に関して積分計算を行うことを表す演算子である。KL距離は常に非負であり、q=pのときにのみ0(ゼロ)になる。計算困難性に対処するために、計算可能な分布の集合の中において事後分布に最も近いものを探索する。
試験分布q(x,z)は、上記数式(14)の最小化が可能であるようなもの、あるいは近似的な最小化が可能であるようなものであれば、任意に選ぶことができる。上記数式(14)の最小化が可能であるような分布としては、q(x,z)=Π_iq(xi|{z}_{N(i)})Πiq(zi)がある。ただし、{z}_{N(i)}は画素i近傍の画素における組織クラス変数zの集合であり、q(xi|{z}_{N(i)})はガウス分布とする。
ここでは、簡単のため下記数式(15)に示されるように試験分布q(x,z)に対して因子化仮定と呼ばれる変数xとzの間の独立性の仮定を行う。さらに、xに関する分布q(x)がガウス分布になるものと仮定して、下記数式(16)に示されるようにモデル化する。これらの仮定の下で、交互最適化を用い最適試験分布を求めている。
【0070】
【数15】
JP0005590548B2_000017t.gif

【0071】
【数16】
JP0005590548B2_000018t.gif

【0072】
(最適試験分布)
ここで、最適試験分布について説明する。試験分布q(x)の平均uと分散σは、q(z)を固定した下で、DKLを最小にするものとして、下記数式(17)のように得ることができる。
【0073】
【数17】
JP0005590548B2_000019t.gif

【0074】
これらの最適化は、SCG法によって行う。試験分布q(x)が推定された後、すなわち平均uと分散σが求まった下で、組織クラスに関する最適試験分布q(z)は、DKLを最小にするものとして、下記数式(18)と数式(19)に示すように解析的に求まる。なお、下記数式(19)におけるρjkは、下記数式(20)で示されるものである。
【0075】
【数18】
JP0005590548B2_000020t.gif

【0076】
【数19】
JP0005590548B2_000021t.gif

【0077】
【数20】
JP0005590548B2_000022t.gif

【0078】
上記数式(20)において、第一項目と第二項目はクラス事前分布に由来する項であり、第三項目は条件付き事前分布の正規化項から生じたものである。また、最後の項は条件付き事前分布に由来しており、推定された吸収係数uとクラスごとに定められた吸収係数vとを比較するものである。このq(z)に関する最適化は、組織に関して組織のクラスをやわらかいクラスタリングによって決めることに対応している。
また、上記数式(20)における第二項目は、同じクラスに属する隠れ変数がつながりやすいことを表現しており、再構成画像の滑らかさは隠れ変数を通して実現されている。
【0079】
ベイズ推定のアルゴリズムの擬似コードを、図5に示す。なお、アルゴリズムは事前に決められた所定の収束条件を満たすまで繰り返される。
【0080】
以上、本発明のX線CT画像処理方法およびX線CT画像処理プログラムにおける処理内容について説明した。以下の実施例では、観測対象としてファントムデータを用いた計算機実験を行い、本発明のX線CT画像処理方法およびX線CT画像処理プログラムの画面再構成を、従来技術による画面再構成と比較して評価する。
下記の評価実験から、本発明のX線CT画像処理方法およびX線CT画像処理プログラムの有用性が理解できる。
【0081】
実施例1として、X線被曝を極力抑えた状況を想定し、限られた放射線量の下での実験を行う。実施例1の結果から、本発明のX線CT画像処理方法およびX線CT画像処理プログラムが、限られたデータからの画像再構成に有効働くことが示される。
そして、実施例2では、投影数は十分得られるが金属の含まれるデータからの断層画像再構成を行う。実施例2の結果から、本発明のX線CT画像処理方法およびX線CT画像処理プログラムが、メタルアーティファクト除去に有効に働くことが示される.
【実施例1】
【0082】
実施例1は、X線被曝を極力抑えた状況を想定し、限られた放射線量の下での実験を行った。
現在広く用いられている医療用CTでは、検出器の数は700~900程度であり、投影数は800~1500程度である。それらのデータから、256×256もしくは512×512のピクセルサイズの断層画像が生成される。X線吸収係数xは、Hounsfield unit (HU) =1000(x-x)/xを用いて表される。ここで、xは水のX線吸収係数である。HU値は、水のX線吸収係数を0、空気のX線吸収係数を-1000として正規化した値である。断層画像の表示を全ての階調で行うと十分なコントラストが得ることができないため、実験ではウインドウ幅を[-500,500] HU に統一して表示を行った(X線吸収係数では[0.01,0.03] に相当する。)。
【0083】
一般に、X線CTによって再構成される断層画像の精度は、放射線量や放出光子の量によって変化し、再構成画像に生じるノイズと放射線量には下記数式(21)の関係が成り立つ。
【0084】
【数21】
JP0005590548B2_000023t.gif

【0085】
ここで、Nは画像に生じるノイズを表している。また、Bは対象のX線透過率,Dは入射線量,hはスライス厚,wはピクセルサイズをそれぞれ表している。放射線の量を増やせば画像に生じるノイズは減少するが、放射線量にしたがって線形に減少するわけではない。これは光子ノイズが、ポワソン分布にしたがって変動するためであると推察する。
【0086】
本発明のX線CT画像処理方法およびX線CT画像処理プログラムでは、観測対象の物質とX線吸収係数を事前知識として用いているが、X線吸収係数は個体差やX線のエネルギーなどによって変化することから、実際の観測対象の吸収係数は予め完全に分かっているわけではない。こうした吸収係数の変動を考慮して、実験では吸収係数の平均値が想定したものから大きくずれた状況を想定する。吸収係数に関する事前分布の分散をr=10-5に設定して、観測対象の真のX線吸収係数を事前分布の平均vからずらしたものを二通り(CaseA, CaseB)を想定した(図6を参照)。なお、空気と金属のX線吸収係数は事前知識とずれはないものと想定し、事前分布の平均値と一致させる。
【0087】
1つ目の実験設定(CaseA)では、筋肉,脂肪,骨のX線吸収係数を事前分布の平均から標準偏差rだけX線吸収係数が離れる方向にずらした。これは158HU程度のずれであり、実際のX線吸収係数の個体差に比べると十分大きな値であると考えられる。2つ目の実験設定(CaseB)では、筋肉と脂肪のX線吸収係数が近くなる方向に50HUずつ事前分布の平均から移動させた。こちらも実際の変動よりは十分に大きい値であると考えられる。
【0088】
設定すべきパラメータとしては、上述したように、X線吸収係数に関する事前分布の平均vと分散r、組織に関する事前分布のJself,Jinterがある。
これらは、実際の観測データセットから決められるべきパラメータである。例えば、肺では太さの異なる血管が張り巡らされており、X線吸収係数の分布は一様ではない。一方水などの均質な物質では吸収係数は一定になる。これらのことから、吸収係数の分布や組織ごとの繋がりやすさは各々の組織に依存することになる。
本実験では、これらのパラメータに関するロバストさを調べるために、前記表1のように実験設定によらず固定した。
【0089】
実施例1の実験で用いたファントムを図7(a)(b)に示す。それぞれのファントムの大きさは、471×353mmである。X線発生源から並行ビームを発生させて、観測対象を通過した光子数を検出器によりカウントした。検出器の数は367個であり、180°の間に36回投影を行いデータを獲得した。光子数bi(t)は105とした。再構成する画像の解像度は、256×256ピクセルとする。比較実験として、フィルター補正逆投影法(FBP)と最尤法(ML)による画像再構成を同時に行った。
MAP推定とベイズ推定の画像再構成の結果を、図8 (CaseAの実験設定),図10(CaseBの実験設定)に示す。それぞれの図で、(a) はFBPによって得られた再構成画像(FBP)、 (b) はMLによって得られた再構成画像(ML),(c)は本発明のX線CT画像処理方法(MAP推定)によって得られた再構成画像(MAP)、(d) は本発明のX線CT画像処理方法(ベイズ推定)によって得られた再構成画像(BAYES) である。
【0090】
MAP推定とベイズ推定で得られた再構成画像を、PSNRにより評価した。それぞれの結果は、図中に記載の通りである。CaseA、CaseBいずれも、本発明のX線CT画像処理方法(ベイズ推定)を用いた再構成画像が最も高いPSNRを示していることがわかる(CaseA:20.32 dB, CaseB:
21.60 dB)。
図に示されるように、他の手法では投影数が制限されていることにより、ノイズが多く含まれた結果が得られている。また、本発明のX線CT画像処理方法のMAP推定とベイズ推定によって得られた組織クラス推定の結果を、図9(CaseAの実験設定)と図11(CaseBの実験設定)に示す。本発明のX線CT画像処理方法により、組織クラスの推定が可能であることが理解できる。図において、上段(a)~(e)はMAP推定による結果であり、下段(f)~(j)はベイズ推定による結果である。それぞれの組織に属していると推定されたピクセルは白色で表示されている。
【0091】
図12は、投影数と誤差の関係を調べるために、投影枚数ごとの各手法によるPSNRを計測したものである。対象物体とパラメータは、CaseAの実験設定と同じものを用いるが、計算コストを下げるために再構成画像のサイズは64×64ピクセルとしている。各手法での計測を6回ずつ行い、平均値と標準偏差を、エラーバーを用いてプロットした。ほとんどの場合で、本発明のX線CT画像処理方法が従来技術による画像再構成を上回っており、本発明のX線CT画像処理方法が低放射線下でのX線CT画像再構成に有効に作用していることが理解できる。
【0092】
また、本発明のX線CT画像処理方法のMAP推定とベイズ推定での結果を比較すると、わずかではあるが、ほとんどの場合でベイズ推定が上回った結果が得られている。パラメータの設定の違いがあるために一概に述べることはできないが、本結果はベイズ推定では再構成画像と所属の組織クラスの不確かさまで考慮した推定を行っているためであると推察する。
なお、本発明のX線CT画像処理方法のMAP推定とベイズ推定と従来技術(ML)の処理速度については、図13に示すように、再構成画像のサイズは64×64ピクセルや128×128ピクセルの場合は大差なく、256×256ピクセル程度以上の解像度になると、ベイズ推定の処理コストが大きくなることがわかる。
【実施例2】
【0093】
実施例2では、投影数は十分得られるが金属の含まれるデータからの断層画像再構成を行う。実施例2の実験で用いた金属が含まれるファントムを図7(c)(d)に示す。対象となる画像は人間の頭(それぞれ大きさ472×436mm) を模擬したデータであり、歯にあたる部分にそれぞれ(18,19,23mm) の大きさの金属が埋め込まれているものである。対象物体にX線平行ビームを投射して、反対側にある367個の検出器で光子の数を検出した。投影は180°の間に72回データをサンプルした。X線発生器で発生させる光子数b(t)は106とした(実施例1の実験よりも光子数は多い)。実施例1と同様に、再構成する画像の解像度は、256×256ピクセルとする。また、比較実験として、フィルター補正逆投影法(FBP)と最尤法(ML)による画像再構成を同時に行った。
MAP推定とベイズ推定の画像再構成の結果を、図14(CaseAの実験設定),図16(CaseBの実験設定)に示す。それぞれの図で、(a) は真の断層画像であり、(b)はFBPによって得られた再構成画像(FBP)、 (c) はPCLISによって得られた再構成画像(PCLIS),(d) はMLによって得られた再構成画像(ML),(e)は本発明のX線CT画像処理方法(MAP推定)によって得られた再構成画像(MAP)、(f) は本発明のX線CT画像処理方法(ベイズ推定)によって得られた再構成画像(BAYES) である。
【0094】
MAP推定とベイズ推定で得られた再構成画像を、PSNRにより評価した。それぞれの結果は、図中に記載の通りである。CaseA、CaseBいずれも、本発明のX線CT画像処理方法(ベイズ推定)を用いた再構成画像が最も高いPSNRを示していることがわかる(CaseA:33.10 dB, CaseB:
34.34 dB)。
また、FBPやPCLISでは、メタルアーティファクトの影響で、画像に筋が入っていたり、ノイズの多い画像になっている。しかし、本発明のX線CT画像処理方法ではメタルアーティファクトが低減されており、真の画像に近い画像が得られていることがわかる。また、本発明のX線CT画像処理方法であるMAP推定とベイズ推定によって得られた組織クラス推定の結果を、図15(CaseAの実験設定)と図17(CaseBの実験設定)に示す。本発明のX線CT画像処理方法であるMAP推定とベイズ推定によって得られた図を比較すると、ベイズ推定によって得られた画像(クラス推定を含め) の方がよりよい結果が得られていることがわかる。図において、上段(a)~(e)はMAP推定による結果であり、下段(f)~(j)はベイズ推定による結果である。それぞれの組織に属していると推定されたピクセルは白色で表示されている。
【0095】
上述の説明では、本発明のX線CT画像処理方法では、隠れ変数を導入し組織に関する事前知識を用いることで、少ない投影データにおいて、従来技術より高精細な画像再構成が可能なこと及びメタルアーティファクトを低減できることを示した。また、実施例の実験では、本発明のX線CT画像処理方法による再構成によって精度の高い推定が行われていることを示した。
【0096】
昨今のX線CTの主な利用先である医療診断では、診断補助としてCT画像を組織ごとに領域分割するための研究が多くなされている。これらの方法はいずれもFBPなどを用いて得られたCT画像から領域分割を行っている。領域分割と画像の再構成は互いに深く関わっているため、別々に行うよりも統一的な枠組みの下で推定を行う方が好ましい。本発明のX線CT画像処理方法では、観測データから統一的な枠組みの下で画像再構成と領域分割が実行されているため、副次的にせよ精度の良いセグメント化がなされていると考える。本発明のX線CT画像処理方法は、医療診断補助などへの応用も大いに期待できるものである。
【産業上の利用可能性】
【0097】
本発明は、医療用のみならず産業用のX線CT装置に有用である。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16