画像処理装置、画像処理システム、画像処理方法およびプログラム
【課題】良好な再構成画像を生成する。
【解決手段】画像処理装置３は、光検出器２への投影線上でボクセルから放出され検出された光子の計数結果に基づいて投影線毎に投影値を算出する投影データ生成部３１と、投影値の集合を複数のサブセットに分割し、サブセット毎に、投影値を使用した反復演算を実行することによりボクセル毎の濃度を推定する濃度推定部３２とを備える。各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【解決手段】画像処理装置３は、光検出器２への投影線上でボクセルから放出され検出された光子の計数結果に基づいて投影線毎に投影値を算出する投影データ生成部３１と、投影値の集合を複数のサブセットに分割し、サブセット毎に、投影値を使用した反復演算を実行することによりボクセル毎の濃度を推定する濃度推定部３２とを備える。各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被測定物から放出された光を検出し、その検出結果に基づいて被測定物の内部構造を表す２次元画像または３次元画像を再構成するための画像処理技術に関する。
【背景技術】
【０００２】
人体などの被測定物を透過したＸ線、あるいは被測定物から放出されたγ線といった放射線を検出し、その検出結果に基づいて被測定物の内部構造の２次元画像または３次元画像を再構成するための画像処理技術の研究開発が進められている。このような画像処理技術は、Ｘ線ＣＴ装置（Ｘ線コンピュータ断層撮影装置）や陽電子放射断層撮像装置（ＰＥＴ：Positron Emission Tomography）装置で使用されている。
【０００３】
画像再構成のアルゴリズムとしては、ＭＬ−ＥＭ（最尤推定−期待値最大化：Maximum Likelihood - Expectation Maximization）法やＯＳ−ＥＭ（Ordered Subsets - Expectation Maximization）法などの逐次近似的な画像再構成法が広く知られている。ＭＬ−ＥＭ法では、後に詳述するように、ｋ回目（ｋは２以上の整数）の濃度値を使用してｋ＋１回目の濃度値を算出するという逐次演算が実行される。ＭＬ−ＥＭ法は、各回で全ての投影データを使用した逐次演算を実行するので膨大な計算量を必要とし、濃度値が収束するまでに非常に時間がかかるという問題がある。ＯＳ−ＥＭ法は、投影データをいくつかのサブセット（部分集合）に分割し、サブセットごとに逐次演算を実行することで計算時間の短縮を図るアルゴリズムである。よって、ＭＬ−ＥＭ法は、サブセットの数を１つとした場合のＯＳ−ＥＭ法と同じである。また、サブセットの数が多いほどに、濃度値の収束に要する時間も短くなる。
【０００４】
ＭＬ−ＥＭ法は、たとえば、非特許文献１や非特許文献２に開示されており、ＯＳ−ＥＭ法は、たとえば、非特許文献３に開示されている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, Vol. MI-1, No. 2, pp. 113-122.
【非特許文献２】K. Lange and R. Carson, "EM reconstruction algorithms for emission and transmission tomography," J. Comput. Assist. Tomogra., 1984, Vol. 8, No. 2, pp. 306-316.
【非特許文献３】H. M. Hudson and R. S. Larkin, "Accelerated image reconstruction using ordered subsets of projection data," IEEE Trans. Med. Imaging, 1994, Vol. 13, pp. 601-609.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
ＯＳ−ＥＭ法に基づいた逐次演算により濃度値が一様にかつ精度良く収束するためには、異なるサブセット間で投影方向の統計誤差の偏りが低いことが要求される。特に、互いに対向する一対の検出プレートを有する対向型ＰＥＴ装置では、後述するように、一対の検出プレートの検出器要素の対で投影方向が決まるので、逐次演算で考慮しなければならない投影方向の範囲（バラツキ）が広い。それ故、サブセットの構成法が悪ければ、或るサブセットで投影方向の偏りが生じてしまい、濃度値の精度や収束速度が低くなることがある。
【０００７】
上記に鑑みて本発明の目的は、サブセット間で投影方向の統計誤差の偏りを極力小さくし得、これにより良好な再構成画像を得ることを可能にする画像処理装置、画像処理システム、画像処理方法およびプログラムを提供する点にある。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理装置が提供される。この画像処理装置は、光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【０００９】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理システムが提供される。この画像処理システムは、互いに対向する一対の検出プレートを有する光検出器と、前記一対の検出プレートを結ぶ投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【００１０】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理方法が提供される。この画像処理方法は、（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出するステップと、（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定するステップと、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【００１１】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成するために、コンピュータ読み取り可能な記録媒体から読み出されて画像再構成処理をコンピュータに実行させるプログラムが提供される。この画像再構成処理は、（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する処理と、（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する処理と、を含み、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【発明の効果】
【００１２】
上述の通り、本発明による画像処理装置、画像処理システム、画像処理方法およびプログラムでは、投影値の集合を複数のサブセットに分割し、サブセット毎に、投影値を使用した反復演算を実行することによりボクセル毎の濃度が推定される。各サブセットは、投影線の方位が均等に分布するように構成されているので、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができる。したがって、各ボクセルの濃度を正確に算出することができ、良好な再構成画像を得ることが可能である。
【図面の簡単な説明】
【００１３】
【図１】本発明に係る一実施形態の画像再構成システム（画像処理システム）の概略構成を示す機能ブロック図である。
【図２】投影線を説明するための模式図である。
【図３】（Ａ）〜（Ｃ）は、Ｘ軸方向への傾きを有する投影方向の３つのパターンを概略的に例示する図である。
【図４】サブセットの構造を模式的に示す図である。
【図５】画像処理装置の動作手順の一例を概略的に示すフローチャートである。
【図６】画像再構成処理の手順の一例を概略的に示すフローチャートである。
【図７】画像処理装置の動作手順の他の例を概略的に示すフローチャートである。
【図８】画像再構成処理の手順の他の例を概略的に示すフローチャートである。
【発明を実施するための形態】
【００１４】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
【００１５】
図１は、本発明に係る一実施形態の画像再構成システム（画像処理システム）１の概略構成を示す機能ブロック図である。この画像再構成システム１は、光検出器である対向型ＰＥＴ装置２と、画像処理装置３とを備えている。対向型ＰＥＴ装置２は、互いに対向する一対の検出プレート１０Ａ，１０Ｂを有しており、これら検出プレート１０Ａ，１０Ｂ間に挟み込まれた被測定物２０から放出された光子対を検出する。一方、画像処理装置３は、光子対の検出結果に基づいて被測定物２０の内部構造を表す２次元画像または３次元画像を推定（再構成）する機能を有している。
【００１６】
対向型ＰＥＴ装置２は、更に、検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１を有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００１７】
一方の検出プレート１０Ａは、格子状に配列された検出器要素１１Ａ，…，１１Ａを有しており、これら検出器要素１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。他方の検出プレート１０Ｂも、格子状に配列された検出器要素１１Ｂ，…，１１Ｂを有しており、これら検出器要素１１Ｂ，…，１１Ｂは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。また、検出器要素群１１Ａ，…，１１Ａの配列と検出器要素１１Ｂ，…，１１Ｂの配列は、検出プレート１０Ａと検出プレート１０Ｂとの重心Ｏに関して点対称の関係を有している。
【００１８】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に応じて、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を通じて制御部４０に指示を入力することにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００１９】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対の検出器要素１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対の検出器要素１１Ａ，１１Ｂを結ぶ線（投影線）上に、陽電子が消滅した地点が存在する。
【００２０】
検出器要素１１Ａ，１１Ｂとしては、入射光に応じた蛍光出力を与えるシンチレータが使用される。シンチレータは、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯ（Ｂｉ４Ｇｅ３Ｏ１２）やＬＳＯ（Ｌｕ２（ＳｉＯ４）Ｏ：Ｃｅ）などの公知の結晶体を使用すればよい。検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、光子を検出した検出器要素１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００２１】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、検出器要素１１Ａ，１１Ｂが検出した光子対を計数し、この計数結果に基づいて投影線毎の実測投影データの値（投影値）を算出する。実測投影データは、メモリ３７内のデータ記憶部３８に格納される。図２は、投影線を説明するための模式図である。図２に示されるように、被測定物２０内の或る投影線ＰＬｉ上の陽電子消滅地点（線源）から放出された光子対が、検出プレート１０Ａ，１０Ｂの一対の検出器要素１１Ａ，１１Ｂで検出される。
【００２２】
投影データ生成部３１は、投影線ＰＬｉ上で放出された光子対を計数し、この計数結果に基づいて投影線ＰＬｉの投影値ｙｉを算出する。ここで、添え字ｉは、投影線を一意に表す番号である。投影線ＰＬｉは、一方の検出プレート１０Ａの検出器要素１１Ａと他方の検出プレート１０Ｂの検出器要素１１Ｂとを結ぶラインであるから、検出器要素１１Ａ，１１Ｂの組み合わせの数だけ投影線が存在することとなる。本実施形態では、ｎ個（ｎは正整数）の投影線ＰＬｉ（すなわち、検出器要素１１Ａ，１１Ｂのｎ個の組み合わせ）が存在するものとする（ｉ＝１〜ｎ）。たとえば、癌細胞のように活性の高い部位には線源が高密度で集中する傾向にあるので、投影線ＰＬｉ上に癌細胞が存在すると、この投影線ＰＬｉに対応する投影値ｙｉは比較的大きな値を示すこととなる。
【００２３】
濃度推定部３２は、被測定物２０が配置されている領域内に複数のボクセルＢｊを立体的に配列した解析モデルを使用して画像再構成処理を実行する。ここで、ボクセルＢｊは、立方形状または直方形状の微少な六面体要素からなる。本実施形態では、ｍ個のボクセルＢｊが配列されているものとする（ｊ＝１〜ｍ）。
【００２４】
濃度推定部３２は、データ記憶部３８から投影値の集合｛ｙｉ｝を読み出し、ＯＳ（Ordered Subset）アルゴリズムに従ってこの投影値の集合｛ｙｉ｝を複数のサブセット（部分集合）Ｓ１〜ＳΜ（Ｍは２以上の整数）に分割する。そして、濃度推定部３２は、サブセット毎に、画像再構成法に従って投影値｛ｙｉ｝を使用した反復演算を実行することにより、ボクセルＢｊ毎の放射能濃度（以下、単に「濃度」と呼ぶ。）ｘｊを推定する。
【００２５】
本実施形態では、投影値の集合｛ｙｉ｝は、投影線ＰＬｉの方向（以下、投影方向と呼ぶ。）に基づいて複数のサブセットＳ１〜ＳＭに分割される。ここで、各サブセットＳμは、自己に属する投影値｛ｙｉ｝に対応する投影線の方位が所定の方位範囲内で均等に分布するように構成されていることに特徴がある。対向型ＰＥＴ装置２の場合、投影線ＰＬの平均方位は、検出プレート１０Ａ，１０Ｂの検出面の法線方向と略一致する。このため、各サブセットＳμに属する投影値｛ｙｉ｝に対応する投影線の方位は、その法線方向を中心にした広い範囲に亘って均等に分布する。
【００２６】
このようなサブセットＳμの構成方法を以下に説明する。投影線の方位としてＰ個の方位（Ｐは、十分大きな正整数）が存在するものとする。これらＰ個の方位の中から、当該方位の傾きの小さい順または大きい順にｐ個置き（ｐは、Ｐよりも小さな正整数）に選択された方位の部分集合に対応する投影値の集合を１つのサブセットＳμとすることができる。ここでの「傾き」は、或る基準方向（たとえば、検出プレート１０Ａ，１０Ｂの検出面の法線方向）に対する所定方向（Ｘ軸方向やＹ軸方向）への傾きをいう。
【００２７】
たとえば、Ｐ個の方位の中から、Ｘ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ１とし、Ｐ個の方位のうち残りの奇数番目の方位を部分集合ＳＤＸ２とすることができる。更に、部分集合ＳＤＸ１の中から、Ｙ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ１Ｙ１とし、部分集合ＳＤＸ１のうち残りの奇数番目の方位を部分集合ＳＤＸ１Ｙ２とすることができる。一方、部分集合ＳＤＸ２の中から、Ｙ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ２Ｙ１とし、部分集合ＳＤＸ２のうち残りの奇数番目の方位を部分集合ＳＤＸ２Ｙ２とすることができる。この結果、Ｐ個の方位の集合は、４つの部分集合ＳＤＸ１Ｙ１，ＳＤＸ１Ｙ２，ＳＤＸ２Ｙ１，ＳＤＸ２Ｙ２に分割される。この場合、投影値の集合｛ｙｉ｝は、投影線の方位の部分集合ＳＤＸ１Ｙ１，ＳＤＸ１Ｙ２，ＳＤＸ２Ｙ１，ＳＤＸ２Ｙ２にそれぞれ対応する４つのサブセットＳ１，Ｓ２，Ｓ３，Ｓ４に分割される。
【００２８】
図３（Ａ）〜（Ｃ）は、Ｘ軸方向への傾きを有する投影方向の３つのパターンを概略的に例示する図である。図３（Ａ）〜（Ｃ）において、一方の検出プレート１０Ａの検出器要素１１ＡのＸ軸方向の位置は変数ｘ１で表され、他方の検出プレート１０Ｂの検出器要素１１ＢのＸ軸方向の位置は変数ｘ２で表される。このように投影線ＰＬの傾き（すなわち投影方向の傾き）は、検出器要素１１Ａ，１１Ｂの組み合わせにより定まる。
【００２９】
図４は、投影線ＰＬのＸ軸方向への傾きをｄｘ＝ｘ２−ｘ１で表した場合のサブセットの構成方法を説明するための図である。図４に示されるように、ハッチングが付されていないサブセットＳ１は、投影線ＰＬのＸ軸方向への傾きｄｘ，−ｄｘの大きさがともに偶数の場合の投影値の部分集合である。また、ハッチングが付されたサブセットＳ２は、投影線ＰＬのＸ軸方向への傾きｄｘ，−ｄｘの大きさがともに奇数の場合の投影値の部分集合である。たとえば、図４において傾きｄｘがゼロの場合の投影線ＰＬは、図３（Ａ）に示すようなパターンを構成し、図４において傾きｄｘ，−ｄｘの大きさがともに１の場合の投影線ＰＬは、図３（Ｂ）に示すようなパターンを構成する。また、図４において傾きｄｘと傾き−ｄｘの大きさとが略等しい場合の投影線ＰＬは、図３（Ｃ）に示すようなパターンを構成することとなる。
【００３０】
Ｘ軸方向の傾きに加えてＹ軸方向への傾きも考慮すれば、Ｘ軸方向の傾きとＹ軸方向の傾きとの組み合わせで、合計４つのサブセットを構成することができる。このようにサブセットを構成することで、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができる。
【００３１】
図１に示されるように、濃度推定部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これら機能ブロック３３〜３５は、ＯＳ−ＥＭ法に従った反復演算を実行する。ＯＳ−ＥＭ法では、μ番目のサブセットＳμの反復演算には、以下の推定式（１）が使用される。
【００３２】
【数１】
【００３３】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００３４】
【数２】
【００３５】
上記推定式（１）中、ｋは、反復回数であり、ｘｊ（ｋ−１）は、（ｋ−１）回目の処理におけるｊ番目のボクセルＢｊの放射能濃度であり、ｙｉは、ｉ番目の投影線ＰＬｉ（検出器要素１１Ａ，１１Ｂのｉ番目の組み合わせ）についての投影値である。サブセットＳμは、投影値ｙｉの部分集合であるが、上式（１），（１Ａ）では、説明の便宜上、サブセットＳμは、投影値ｙｉの添え字ｉの部分集合とされている。また、ｘｊ（ｋ）は、ｋ回目の処理におけるｊ番目のボクセルＢｊの放射能濃度である。よって、（ｋ−１）回目の濃度ｘｊ（ｋ−１）から、ｋ回目の濃度ｘｊ（ｋ）が算出される。初期値ｘｊ（ｋ＝０）は正の実数であれば、どのような値に設定されてもよい。
【００３６】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｊから飛来した光子対がｉ番目の検出器要素対（検出器要素１１Ａ，１１Ｂのｉ番目の組み合わせ）で検出される確率（検出確率）である。ボクセルＢｊから飛来した光子対がｉ番目の検出器要素対に入射する数は、検出確率Ｃｉｊに依存する。検出確率ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００３７】
検出確率Ｃｉｊは、行列（「システムマトリクス」と呼ばれている。）のｉ行ｊ列目の要素とみなすことができる。検出確率Ｃｉｊの値には、放射線の散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、検出確率Ｃｉｊの値は、ｉ番目の検出器要素対に対するｊ番目のボクセルＢｊの幾何学的位置のみにより予め定められるものとする。
【００３８】
次に、図５および図６を参照しつつ、画像処理装置３の動作を以下に説明する。図５は、画像処理装置３の動作手順の一例を概略的に示すフローチャートであり、図６は、μ番目のサブセットＳμについての画像再構成処理の手順の一例を概略的に示すフローチャートである。
【００３９】
図５を参照すると、先ず、ステップＳ９では、濃度推定部３２は、サブセットＳμの添え字（サブセット番号）μの値を１に初期化する。次のステップＳ１０では、濃度推定部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝（ｉ∈Ｓμ）と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から投影値｛ｙｉ｝（ｉ∈Ｓμ）を読み出す（ステップＳ１１）。なお、本実施形態では、検出確率｛Ｃｉｊ｝（ｉ∈Ｓμ）、係数｛Ｃｊ｝および投影値｛ｙｉ｝（ｉ∈Ｓμ）が一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００４０】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２）、その後、画像再構成処理（ステップＳ１３）を実行する。
【００４１】
図６は、図５のステップＳ１３の画像再構成処理の手順を示すフローチャートである。図６を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化される（ステップＳ１３１）。その後、図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影値｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００４２】
【数３】
【００４３】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影値｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００４４】
【数４】
【００４５】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００４６】
その後、濃度算出部３５は、濃度｛ｘｊ｝が収束したか否かを判定する（ステップＳ１３５）。具体的には、対数尤度関数（または尤度関数）を用いて濃度｛ｘｊ｝の収束の有無を判定することが可能である。非特許文献１（L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, v. 1, pp. 113-122.）によれば、上式（１）は、ポアッソン分布の式により与えられる対数尤度関数Ｌ（｛ｘｊ｝）に基づいて導出される。濃度算出部３５は、ステップＳ１３４で算出された濃度｛ｘｊ｝を入力としたときの対数尤度関数Ｌ（｛ｘｊ｝）の出力値が飽和し略一定に達した場合に、濃度｛ｘｊ｝が収束したと判定することができる（ステップＳ１３５のＹＥＳ）。このとき、濃度｛ｘｊ｝は、当該サブセットＳμについての濃度｛ｘｊ｝μとしてメモリに記憶される（ステップＳ１３９）。その後、処理は図５のステップＳ１４に移行する。
【００４７】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、更に、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされ（ステップＳ１３７）、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）は、濃度算出部３５は、濃度｛ｘｊ｝を当該サブセットＳμに関する濃度｛ｘｊ｝μとしてメモリに記憶する（ステップＳ１３９）。その後、処理は図５のステップＳ１４に移行する。
【００４８】
その後、濃度推定部３２は、サブセット番号μが最大値Ｍであるか否かを判定する（図５のステップＳ１４）。サブセット番号μが最大値Ｍでないと判定された場合（ステップＳ１４のＮＯ）、濃度推定部３２は、サブセット番号μをインクリメントして（ステップＳ１５）、ステップＳ１０を実行する。
【００４９】
サブセット番号μが最大値Ｍに達したと判定された場合（ステップＳ１４のＹＥＳ）、濃度算出部３５は、全てのサブセットＳ１〜ＳＭに関するＭ組の濃度｛ｘｊ｝μ＝１，｛ｘｊ｝μ＝２，…，｛ｘｊ｝μ＝Ｍに基づいて最終的に一組の推定放射能濃度｛ｘｊ｝を算出する（ステップＳ１６）。たとえば、各ボクセルについて、Ｍ個の濃度ｘｊ（μ＝１），ｘｊ（μ＝２），…，ｘｊ（μ＝Ｍ）を平均化することで最終的な推定濃度ｘｊを算出することができる。あるいは、濃度｛ｘｊ｝μ＝１，｛ｘｊ｝μ＝２，…，｛ｘｊ｝μ＝Ｍのうち最も収束の良好な一組の濃度を推定放射能濃度として算出してもよい。
【００５０】
その後、図１の表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００５１】
次に、図７および図８を参照しつつ、再構成画像の他の生成方法について説明する。図７は、画像処理装置３の動作手順の他の例を概略的に示すフローチャートであり、図８は、画像再構成処理の手順の他の例を概略的に示すフローチャートである。
【００５２】
図７を参照すると、先ず、ステップＳ１０Ｔでは、濃度推定部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から投影値｛ｙｉ｝を読み出す（ステップＳ１１Ｔ）。なお、本実施形態では、検出確率｛Ｃｉｊ｝、係数｛Ｃｊ｝および投影値｛ｙｉ｝が一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００５３】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２Ｔ）、その後、画像再構成処理（ステップＳ１３Ｔ）を実行する。
【００５４】
図８は、図７のステップＳ１３Ｔの画像再構成処理の手順を示すフローチャートである。図８を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化されるとともに、サブセット番号μの値が１に初期化される（ステップＳ１４１）。その後、図１の前方投影算出部３３は、上式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影値｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００５５】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、上式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影値｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００５６】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００５７】
その後、濃度推定部３２は、サブセット番号μが最大値Ｍであるか否かを判定する（ステップＳ１４２）。サブセット番号μが最大値Ｍでないと判定された場合（ステップＳ１４２のＮＯ）、濃度推定部３２は、サブセット番号μをインクリメントするとともに、反復回数ｋをインクリメントして（ステップＳ１４３）、ステップＳ１３２に処理を戻す。その後、以前のサブセットＳμ−１に関してステップＳ１３４で算出された推定濃度｛ｘｊ（ｋ−１）｝を用いて、新たなサブセットＳμに関する推定濃度｛ｘｊ（ｋ）｝が算出される（ステップＳ１３２〜Ｓ１３４）。
【００５８】
ステップＳ１４２でサブセット番号μが最大値Ｍであると判定されたとき（ステップＳ１４２のＹＥＳ）、濃度算出部３５は、濃度｛ｘｊ｝（＝｛ｘｊ（ｋ）｝）が収束したか否かを判定する（ステップＳ１３５）。濃度｛ｘｊ｝が収束したと判定されたとき（ステップＳ１３５のＹＥＳ）、処理は図５のステップＳ１４に移行する。
【００５９】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、更に、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされるとともにサブセット番号μが初期値（＝１）に設定される（ステップＳ１３７Ｔ）。次いで、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）、以上の処理は終了する。
【００６０】
その後、表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００６１】
上述した通り、本実施形態の画像再構成システム１は、投影線ＰＬｉの方向に基づいて投影値の集合｛ｙｉ｝を複数のサブセットＳ１〜ＳＭに分割し、サブセットＳμ毎に、投影値｛ｙｉ｝（ｉ∈Ｓμ）を使用した反復演算を実行することによりボクセル毎の濃度｛ｘｊ｝を推定する。ここで、各サブセットＳμは、投影線ＰＬｉの方位が所定の方位範囲内で均等に分布するように構成されているので、異なるサブセット間の投影方向の統計誤差の偏りを極力抑えることができる。したがって、各ボクセルの濃度を正確に算出することができ、良好な再構成画像を得ることが可能である。
【００６２】
以上、図面を参照して本発明の実施形態について述べた。図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。
【００６３】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【００６４】
上記実施形態は本発明の例示であり、上記以外の様々な構成を採用することもできる。たとえば、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１の対向型ＰＥＴ装置２に代えて、連続的に分布する単一の検出面あるいは３個以上の検出面を有する検出装置を使用してもよい。また、検出面の形状は平面に限定されるものではなく曲面であってもよい。
【００６５】
また、上記実施形態では、投影値の集合｛ｙｉ｝は、投影線ＰＬｉの方位（投影方向）に基づいて複数のサブセットＳ１〜ＳＭに分割されたが、これに限定されるものではない。たとえば、投影値ｙｉがランダムにサブセットＳ１〜ＳＭに振り分けられてもよい。この場合でも、各サブセットにおいて投影線ＰＬｉの方位が略均等に分布するので、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができるので、良好な再構成画像を得ることが可能である。
【００６６】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置２に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
【００６７】
上記実施形態では、画像再構成法としてＯＳ−ＥＭ法（すなわち、上記実施形態のＯＳ（Ordered Subset）アルゴリムが適用されたＭＬ−ＥＭ法）が使用されているが、これに限定されるものではない。ＭＬ−ＥＭ法に限らず、ＭＡＰ−ＥＭ（Maximum A Posteriori - Expectation Maximization）法などの他の統計的画像再構成にＯＳアルゴリズムを適用することも可能である。
【符号の説明】
【００６８】
１ 画像再構成システム
２ 対向型ＰＥＴ装置
３ 画像処理装置
１０Ａ，１０Ｂ 検出プレート
１１Ａ，１１Ｂ 検出器要素
２０ 被測定物
３１ 投影データ生成部
３２ 濃度推定部
３３ 前方投影算出部
３４ 後方投影算出部
３５ 濃度算出部
３６ 表示制御部
３７ メモリ
３８ データ記憶部
３９ パラメータ格納部
４０ 制御部
４１ 駆動機構
５０ ディスプレイ
５１ 入力操作部
【技術分野】
【０００１】
本発明は、被測定物から放出された光を検出し、その検出結果に基づいて被測定物の内部構造を表す２次元画像または３次元画像を再構成するための画像処理技術に関する。
【背景技術】
【０００２】
人体などの被測定物を透過したＸ線、あるいは被測定物から放出されたγ線といった放射線を検出し、その検出結果に基づいて被測定物の内部構造の２次元画像または３次元画像を再構成するための画像処理技術の研究開発が進められている。このような画像処理技術は、Ｘ線ＣＴ装置（Ｘ線コンピュータ断層撮影装置）や陽電子放射断層撮像装置（ＰＥＴ：Positron Emission Tomography）装置で使用されている。
【０００３】
画像再構成のアルゴリズムとしては、ＭＬ−ＥＭ（最尤推定−期待値最大化：Maximum Likelihood - Expectation Maximization）法やＯＳ−ＥＭ（Ordered Subsets - Expectation Maximization）法などの逐次近似的な画像再構成法が広く知られている。ＭＬ−ＥＭ法では、後に詳述するように、ｋ回目（ｋは２以上の整数）の濃度値を使用してｋ＋１回目の濃度値を算出するという逐次演算が実行される。ＭＬ−ＥＭ法は、各回で全ての投影データを使用した逐次演算を実行するので膨大な計算量を必要とし、濃度値が収束するまでに非常に時間がかかるという問題がある。ＯＳ−ＥＭ法は、投影データをいくつかのサブセット（部分集合）に分割し、サブセットごとに逐次演算を実行することで計算時間の短縮を図るアルゴリズムである。よって、ＭＬ−ＥＭ法は、サブセットの数を１つとした場合のＯＳ−ＥＭ法と同じである。また、サブセットの数が多いほどに、濃度値の収束に要する時間も短くなる。
【０００４】
ＭＬ−ＥＭ法は、たとえば、非特許文献１や非特許文献２に開示されており、ＯＳ−ＥＭ法は、たとえば、非特許文献３に開示されている。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, Vol. MI-1, No. 2, pp. 113-122.
【非特許文献２】K. Lange and R. Carson, "EM reconstruction algorithms for emission and transmission tomography," J. Comput. Assist. Tomogra., 1984, Vol. 8, No. 2, pp. 306-316.
【非特許文献３】H. M. Hudson and R. S. Larkin, "Accelerated image reconstruction using ordered subsets of projection data," IEEE Trans. Med. Imaging, 1994, Vol. 13, pp. 601-609.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
ＯＳ−ＥＭ法に基づいた逐次演算により濃度値が一様にかつ精度良く収束するためには、異なるサブセット間で投影方向の統計誤差の偏りが低いことが要求される。特に、互いに対向する一対の検出プレートを有する対向型ＰＥＴ装置では、後述するように、一対の検出プレートの検出器要素の対で投影方向が決まるので、逐次演算で考慮しなければならない投影方向の範囲（バラツキ）が広い。それ故、サブセットの構成法が悪ければ、或るサブセットで投影方向の偏りが生じてしまい、濃度値の精度や収束速度が低くなることがある。
【０００７】
上記に鑑みて本発明の目的は、サブセット間で投影方向の統計誤差の偏りを極力小さくし得、これにより良好な再構成画像を得ることを可能にする画像処理装置、画像処理システム、画像処理方法およびプログラムを提供する点にある。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理装置が提供される。この画像処理装置は、光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【０００９】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理システムが提供される。この画像処理システムは、互いに対向する一対の検出プレートを有する光検出器と、前記一対の検出プレートを結ぶ投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【００１０】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理方法が提供される。この画像処理方法は、（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出するステップと、（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定するステップと、を備えており、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【００１１】
本発明によれば、被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成するために、コンピュータ読み取り可能な記録媒体から読み出されて画像再構成処理をコンピュータに実行させるプログラムが提供される。この画像再構成処理は、（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する処理と、（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する処理と、を含み、前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されている。
【発明の効果】
【００１２】
上述の通り、本発明による画像処理装置、画像処理システム、画像処理方法およびプログラムでは、投影値の集合を複数のサブセットに分割し、サブセット毎に、投影値を使用した反復演算を実行することによりボクセル毎の濃度が推定される。各サブセットは、投影線の方位が均等に分布するように構成されているので、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができる。したがって、各ボクセルの濃度を正確に算出することができ、良好な再構成画像を得ることが可能である。
【図面の簡単な説明】
【００１３】
【図１】本発明に係る一実施形態の画像再構成システム（画像処理システム）の概略構成を示す機能ブロック図である。
【図２】投影線を説明するための模式図である。
【図３】（Ａ）〜（Ｃ）は、Ｘ軸方向への傾きを有する投影方向の３つのパターンを概略的に例示する図である。
【図４】サブセットの構造を模式的に示す図である。
【図５】画像処理装置の動作手順の一例を概略的に示すフローチャートである。
【図６】画像再構成処理の手順の一例を概略的に示すフローチャートである。
【図７】画像処理装置の動作手順の他の例を概略的に示すフローチャートである。
【図８】画像再構成処理の手順の他の例を概略的に示すフローチャートである。
【発明を実施するための形態】
【００１４】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
【００１５】
図１は、本発明に係る一実施形態の画像再構成システム（画像処理システム）１の概略構成を示す機能ブロック図である。この画像再構成システム１は、光検出器である対向型ＰＥＴ装置２と、画像処理装置３とを備えている。対向型ＰＥＴ装置２は、互いに対向する一対の検出プレート１０Ａ，１０Ｂを有しており、これら検出プレート１０Ａ，１０Ｂ間に挟み込まれた被測定物２０から放出された光子対を検出する。一方、画像処理装置３は、光子対の検出結果に基づいて被測定物２０の内部構造を表す２次元画像または３次元画像を推定（再構成）する機能を有している。
【００１６】
対向型ＰＥＴ装置２は、更に、検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１を有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００１７】
一方の検出プレート１０Ａは、格子状に配列された検出器要素１１Ａ，…，１１Ａを有しており、これら検出器要素１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。他方の検出プレート１０Ｂも、格子状に配列された検出器要素１１Ｂ，…，１１Ｂを有しており、これら検出器要素１１Ｂ，…，１１Ｂは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。また、検出器要素群１１Ａ，…，１１Ａの配列と検出器要素１１Ｂ，…，１１Ｂの配列は、検出プレート１０Ａと検出プレート１０Ｂとの重心Ｏに関して点対称の関係を有している。
【００１８】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に応じて、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を通じて制御部４０に指示を入力することにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００１９】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対の検出器要素１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対の検出器要素１１Ａ，１１Ｂを結ぶ線（投影線）上に、陽電子が消滅した地点が存在する。
【００２０】
検出器要素１１Ａ，１１Ｂとしては、入射光に応じた蛍光出力を与えるシンチレータが使用される。シンチレータは、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯ（Ｂｉ４Ｇｅ３Ｏ１２）やＬＳＯ（Ｌｕ２（ＳｉＯ４）Ｏ：Ｃｅ）などの公知の結晶体を使用すればよい。検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、光子を検出した検出器要素１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００２１】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、検出器要素１１Ａ，１１Ｂが検出した光子対を計数し、この計数結果に基づいて投影線毎の実測投影データの値（投影値）を算出する。実測投影データは、メモリ３７内のデータ記憶部３８に格納される。図２は、投影線を説明するための模式図である。図２に示されるように、被測定物２０内の或る投影線ＰＬｉ上の陽電子消滅地点（線源）から放出された光子対が、検出プレート１０Ａ，１０Ｂの一対の検出器要素１１Ａ，１１Ｂで検出される。
【００２２】
投影データ生成部３１は、投影線ＰＬｉ上で放出された光子対を計数し、この計数結果に基づいて投影線ＰＬｉの投影値ｙｉを算出する。ここで、添え字ｉは、投影線を一意に表す番号である。投影線ＰＬｉは、一方の検出プレート１０Ａの検出器要素１１Ａと他方の検出プレート１０Ｂの検出器要素１１Ｂとを結ぶラインであるから、検出器要素１１Ａ，１１Ｂの組み合わせの数だけ投影線が存在することとなる。本実施形態では、ｎ個（ｎは正整数）の投影線ＰＬｉ（すなわち、検出器要素１１Ａ，１１Ｂのｎ個の組み合わせ）が存在するものとする（ｉ＝１〜ｎ）。たとえば、癌細胞のように活性の高い部位には線源が高密度で集中する傾向にあるので、投影線ＰＬｉ上に癌細胞が存在すると、この投影線ＰＬｉに対応する投影値ｙｉは比較的大きな値を示すこととなる。
【００２３】
濃度推定部３２は、被測定物２０が配置されている領域内に複数のボクセルＢｊを立体的に配列した解析モデルを使用して画像再構成処理を実行する。ここで、ボクセルＢｊは、立方形状または直方形状の微少な六面体要素からなる。本実施形態では、ｍ個のボクセルＢｊが配列されているものとする（ｊ＝１〜ｍ）。
【００２４】
濃度推定部３２は、データ記憶部３８から投影値の集合｛ｙｉ｝を読み出し、ＯＳ（Ordered Subset）アルゴリズムに従ってこの投影値の集合｛ｙｉ｝を複数のサブセット（部分集合）Ｓ１〜ＳΜ（Ｍは２以上の整数）に分割する。そして、濃度推定部３２は、サブセット毎に、画像再構成法に従って投影値｛ｙｉ｝を使用した反復演算を実行することにより、ボクセルＢｊ毎の放射能濃度（以下、単に「濃度」と呼ぶ。）ｘｊを推定する。
【００２５】
本実施形態では、投影値の集合｛ｙｉ｝は、投影線ＰＬｉの方向（以下、投影方向と呼ぶ。）に基づいて複数のサブセットＳ１〜ＳＭに分割される。ここで、各サブセットＳμは、自己に属する投影値｛ｙｉ｝に対応する投影線の方位が所定の方位範囲内で均等に分布するように構成されていることに特徴がある。対向型ＰＥＴ装置２の場合、投影線ＰＬの平均方位は、検出プレート１０Ａ，１０Ｂの検出面の法線方向と略一致する。このため、各サブセットＳμに属する投影値｛ｙｉ｝に対応する投影線の方位は、その法線方向を中心にした広い範囲に亘って均等に分布する。
【００２６】
このようなサブセットＳμの構成方法を以下に説明する。投影線の方位としてＰ個の方位（Ｐは、十分大きな正整数）が存在するものとする。これらＰ個の方位の中から、当該方位の傾きの小さい順または大きい順にｐ個置き（ｐは、Ｐよりも小さな正整数）に選択された方位の部分集合に対応する投影値の集合を１つのサブセットＳμとすることができる。ここでの「傾き」は、或る基準方向（たとえば、検出プレート１０Ａ，１０Ｂの検出面の法線方向）に対する所定方向（Ｘ軸方向やＹ軸方向）への傾きをいう。
【００２７】
たとえば、Ｐ個の方位の中から、Ｘ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ１とし、Ｐ個の方位のうち残りの奇数番目の方位を部分集合ＳＤＸ２とすることができる。更に、部分集合ＳＤＸ１の中から、Ｙ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ１Ｙ１とし、部分集合ＳＤＸ１のうち残りの奇数番目の方位を部分集合ＳＤＸ１Ｙ２とすることができる。一方、部分集合ＳＤＸ２の中から、Ｙ軸方向への傾きの小さい順に１個置きに選択された偶数番目の方位を部分集合ＳＤＸ２Ｙ１とし、部分集合ＳＤＸ２のうち残りの奇数番目の方位を部分集合ＳＤＸ２Ｙ２とすることができる。この結果、Ｐ個の方位の集合は、４つの部分集合ＳＤＸ１Ｙ１，ＳＤＸ１Ｙ２，ＳＤＸ２Ｙ１，ＳＤＸ２Ｙ２に分割される。この場合、投影値の集合｛ｙｉ｝は、投影線の方位の部分集合ＳＤＸ１Ｙ１，ＳＤＸ１Ｙ２，ＳＤＸ２Ｙ１，ＳＤＸ２Ｙ２にそれぞれ対応する４つのサブセットＳ１，Ｓ２，Ｓ３，Ｓ４に分割される。
【００２８】
図３（Ａ）〜（Ｃ）は、Ｘ軸方向への傾きを有する投影方向の３つのパターンを概略的に例示する図である。図３（Ａ）〜（Ｃ）において、一方の検出プレート１０Ａの検出器要素１１ＡのＸ軸方向の位置は変数ｘ１で表され、他方の検出プレート１０Ｂの検出器要素１１ＢのＸ軸方向の位置は変数ｘ２で表される。このように投影線ＰＬの傾き（すなわち投影方向の傾き）は、検出器要素１１Ａ，１１Ｂの組み合わせにより定まる。
【００２９】
図４は、投影線ＰＬのＸ軸方向への傾きをｄｘ＝ｘ２−ｘ１で表した場合のサブセットの構成方法を説明するための図である。図４に示されるように、ハッチングが付されていないサブセットＳ１は、投影線ＰＬのＸ軸方向への傾きｄｘ，−ｄｘの大きさがともに偶数の場合の投影値の部分集合である。また、ハッチングが付されたサブセットＳ２は、投影線ＰＬのＸ軸方向への傾きｄｘ，−ｄｘの大きさがともに奇数の場合の投影値の部分集合である。たとえば、図４において傾きｄｘがゼロの場合の投影線ＰＬは、図３（Ａ）に示すようなパターンを構成し、図４において傾きｄｘ，−ｄｘの大きさがともに１の場合の投影線ＰＬは、図３（Ｂ）に示すようなパターンを構成する。また、図４において傾きｄｘと傾き−ｄｘの大きさとが略等しい場合の投影線ＰＬは、図３（Ｃ）に示すようなパターンを構成することとなる。
【００３０】
Ｘ軸方向の傾きに加えてＹ軸方向への傾きも考慮すれば、Ｘ軸方向の傾きとＹ軸方向の傾きとの組み合わせで、合計４つのサブセットを構成することができる。このようにサブセットを構成することで、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができる。
【００３１】
図１に示されるように、濃度推定部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これら機能ブロック３３〜３５は、ＯＳ−ＥＭ法に従った反復演算を実行する。ＯＳ−ＥＭ法では、μ番目のサブセットＳμの反復演算には、以下の推定式（１）が使用される。
【００３２】
【数１】
【００３３】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００３４】
【数２】
【００３５】
上記推定式（１）中、ｋは、反復回数であり、ｘｊ（ｋ−１）は、（ｋ−１）回目の処理におけるｊ番目のボクセルＢｊの放射能濃度であり、ｙｉは、ｉ番目の投影線ＰＬｉ（検出器要素１１Ａ，１１Ｂのｉ番目の組み合わせ）についての投影値である。サブセットＳμは、投影値ｙｉの部分集合であるが、上式（１），（１Ａ）では、説明の便宜上、サブセットＳμは、投影値ｙｉの添え字ｉの部分集合とされている。また、ｘｊ（ｋ）は、ｋ回目の処理におけるｊ番目のボクセルＢｊの放射能濃度である。よって、（ｋ−１）回目の濃度ｘｊ（ｋ−１）から、ｋ回目の濃度ｘｊ（ｋ）が算出される。初期値ｘｊ（ｋ＝０）は正の実数であれば、どのような値に設定されてもよい。
【００３６】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｊから飛来した光子対がｉ番目の検出器要素対（検出器要素１１Ａ，１１Ｂのｉ番目の組み合わせ）で検出される確率（検出確率）である。ボクセルＢｊから飛来した光子対がｉ番目の検出器要素対に入射する数は、検出確率Ｃｉｊに依存する。検出確率ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００３７】
検出確率Ｃｉｊは、行列（「システムマトリクス」と呼ばれている。）のｉ行ｊ列目の要素とみなすことができる。検出確率Ｃｉｊの値には、放射線の散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、検出確率Ｃｉｊの値は、ｉ番目の検出器要素対に対するｊ番目のボクセルＢｊの幾何学的位置のみにより予め定められるものとする。
【００３８】
次に、図５および図６を参照しつつ、画像処理装置３の動作を以下に説明する。図５は、画像処理装置３の動作手順の一例を概略的に示すフローチャートであり、図６は、μ番目のサブセットＳμについての画像再構成処理の手順の一例を概略的に示すフローチャートである。
【００３９】
図５を参照すると、先ず、ステップＳ９では、濃度推定部３２は、サブセットＳμの添え字（サブセット番号）μの値を１に初期化する。次のステップＳ１０では、濃度推定部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝（ｉ∈Ｓμ）と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から投影値｛ｙｉ｝（ｉ∈Ｓμ）を読み出す（ステップＳ１１）。なお、本実施形態では、検出確率｛Ｃｉｊ｝（ｉ∈Ｓμ）、係数｛Ｃｊ｝および投影値｛ｙｉ｝（ｉ∈Ｓμ）が一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００４０】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２）、その後、画像再構成処理（ステップＳ１３）を実行する。
【００４１】
図６は、図５のステップＳ１３の画像再構成処理の手順を示すフローチャートである。図６を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化される（ステップＳ１３１）。その後、図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影値｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００４２】
【数３】
【００４３】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影値｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００４４】
【数４】
【００４５】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００４６】
その後、濃度算出部３５は、濃度｛ｘｊ｝が収束したか否かを判定する（ステップＳ１３５）。具体的には、対数尤度関数（または尤度関数）を用いて濃度｛ｘｊ｝の収束の有無を判定することが可能である。非特許文献１（L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, v. 1, pp. 113-122.）によれば、上式（１）は、ポアッソン分布の式により与えられる対数尤度関数Ｌ（｛ｘｊ｝）に基づいて導出される。濃度算出部３５は、ステップＳ１３４で算出された濃度｛ｘｊ｝を入力としたときの対数尤度関数Ｌ（｛ｘｊ｝）の出力値が飽和し略一定に達した場合に、濃度｛ｘｊ｝が収束したと判定することができる（ステップＳ１３５のＹＥＳ）。このとき、濃度｛ｘｊ｝は、当該サブセットＳμについての濃度｛ｘｊ｝μとしてメモリに記憶される（ステップＳ１３９）。その後、処理は図５のステップＳ１４に移行する。
【００４７】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、更に、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされ（ステップＳ１３７）、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）は、濃度算出部３５は、濃度｛ｘｊ｝を当該サブセットＳμに関する濃度｛ｘｊ｝μとしてメモリに記憶する（ステップＳ１３９）。その後、処理は図５のステップＳ１４に移行する。
【００４８】
その後、濃度推定部３２は、サブセット番号μが最大値Ｍであるか否かを判定する（図５のステップＳ１４）。サブセット番号μが最大値Ｍでないと判定された場合（ステップＳ１４のＮＯ）、濃度推定部３２は、サブセット番号μをインクリメントして（ステップＳ１５）、ステップＳ１０を実行する。
【００４９】
サブセット番号μが最大値Ｍに達したと判定された場合（ステップＳ１４のＹＥＳ）、濃度算出部３５は、全てのサブセットＳ１〜ＳＭに関するＭ組の濃度｛ｘｊ｝μ＝１，｛ｘｊ｝μ＝２，…，｛ｘｊ｝μ＝Ｍに基づいて最終的に一組の推定放射能濃度｛ｘｊ｝を算出する（ステップＳ１６）。たとえば、各ボクセルについて、Ｍ個の濃度ｘｊ（μ＝１），ｘｊ（μ＝２），…，ｘｊ（μ＝Ｍ）を平均化することで最終的な推定濃度ｘｊを算出することができる。あるいは、濃度｛ｘｊ｝μ＝１，｛ｘｊ｝μ＝２，…，｛ｘｊ｝μ＝Ｍのうち最も収束の良好な一組の濃度を推定放射能濃度として算出してもよい。
【００５０】
その後、図１の表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００５１】
次に、図７および図８を参照しつつ、再構成画像の他の生成方法について説明する。図７は、画像処理装置３の動作手順の他の例を概略的に示すフローチャートであり、図８は、画像再構成処理の手順の他の例を概略的に示すフローチャートである。
【００５２】
図７を参照すると、先ず、ステップＳ１０Ｔでは、濃度推定部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から投影値｛ｙｉ｝を読み出す（ステップＳ１１Ｔ）。なお、本実施形態では、検出確率｛Ｃｉｊ｝、係数｛Ｃｊ｝および投影値｛ｙｉ｝が一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００５３】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２Ｔ）、その後、画像再構成処理（ステップＳ１３Ｔ）を実行する。
【００５４】
図８は、図７のステップＳ１３Ｔの画像再構成処理の手順を示すフローチャートである。図８を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化されるとともに、サブセット番号μの値が１に初期化される（ステップＳ１４１）。その後、図１の前方投影算出部３３は、上式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影値｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００５５】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、上式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影値｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００５６】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００５７】
その後、濃度推定部３２は、サブセット番号μが最大値Ｍであるか否かを判定する（ステップＳ１４２）。サブセット番号μが最大値Ｍでないと判定された場合（ステップＳ１４２のＮＯ）、濃度推定部３２は、サブセット番号μをインクリメントするとともに、反復回数ｋをインクリメントして（ステップＳ１４３）、ステップＳ１３２に処理を戻す。その後、以前のサブセットＳμ−１に関してステップＳ１３４で算出された推定濃度｛ｘｊ（ｋ−１）｝を用いて、新たなサブセットＳμに関する推定濃度｛ｘｊ（ｋ）｝が算出される（ステップＳ１３２〜Ｓ１３４）。
【００５８】
ステップＳ１４２でサブセット番号μが最大値Ｍであると判定されたとき（ステップＳ１４２のＹＥＳ）、濃度算出部３５は、濃度｛ｘｊ｝（＝｛ｘｊ（ｋ）｝）が収束したか否かを判定する（ステップＳ１３５）。濃度｛ｘｊ｝が収束したと判定されたとき（ステップＳ１３５のＹＥＳ）、処理は図５のステップＳ１４に移行する。
【００５９】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、更に、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされるとともにサブセット番号μが初期値（＝１）に設定される（ステップＳ１３７Ｔ）。次いで、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）、以上の処理は終了する。
【００６０】
その後、表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００６１】
上述した通り、本実施形態の画像再構成システム１は、投影線ＰＬｉの方向に基づいて投影値の集合｛ｙｉ｝を複数のサブセットＳ１〜ＳＭに分割し、サブセットＳμ毎に、投影値｛ｙｉ｝（ｉ∈Ｓμ）を使用した反復演算を実行することによりボクセル毎の濃度｛ｘｊ｝を推定する。ここで、各サブセットＳμは、投影線ＰＬｉの方位が所定の方位範囲内で均等に分布するように構成されているので、異なるサブセット間の投影方向の統計誤差の偏りを極力抑えることができる。したがって、各ボクセルの濃度を正確に算出することができ、良好な再構成画像を得ることが可能である。
【００６２】
以上、図面を参照して本発明の実施形態について述べた。図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。
【００６３】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【００６４】
上記実施形態は本発明の例示であり、上記以外の様々な構成を採用することもできる。たとえば、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１の対向型ＰＥＴ装置２に代えて、連続的に分布する単一の検出面あるいは３個以上の検出面を有する検出装置を使用してもよい。また、検出面の形状は平面に限定されるものではなく曲面であってもよい。
【００６５】
また、上記実施形態では、投影値の集合｛ｙｉ｝は、投影線ＰＬｉの方位（投影方向）に基づいて複数のサブセットＳ１〜ＳＭに分割されたが、これに限定されるものではない。たとえば、投影値ｙｉがランダムにサブセットＳ１〜ＳＭに振り分けられてもよい。この場合でも、各サブセットにおいて投影線ＰＬｉの方位が略均等に分布するので、異なるサブセット間の投影方向の統計誤差の偏りを低く抑えることができるので、良好な再構成画像を得ることが可能である。
【００６６】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置２に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
【００６７】
上記実施形態では、画像再構成法としてＯＳ−ＥＭ法（すなわち、上記実施形態のＯＳ（Ordered Subset）アルゴリムが適用されたＭＬ−ＥＭ法）が使用されているが、これに限定されるものではない。ＭＬ−ＥＭ法に限らず、ＭＡＰ−ＥＭ（Maximum A Posteriori - Expectation Maximization）法などの他の統計的画像再構成にＯＳアルゴリズムを適用することも可能である。
【符号の説明】
【００６８】
１ 画像再構成システム
２ 対向型ＰＥＴ装置
３ 画像処理装置
１０Ａ，１０Ｂ 検出プレート
１１Ａ，１１Ｂ 検出器要素
２０ 被測定物
３１ 投影データ生成部
３２ 濃度推定部
３３ 前方投影算出部
３４ 後方投影算出部
３５ 濃度算出部
３６ 表示制御部
３７ メモリ
３８ データ記憶部
３９ パラメータ格納部
４０ 制御部
４１ 駆動機構
５０ ディスプレイ
５１ 入力操作部
【特許請求の範囲】
【請求項１】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理装置であって、
光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、
前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理装置。
【請求項２】
請求項１に記載の画像処理装置であって、
前記投影線の方位にはＰ個の方位（Ｐは正整数）が存在し、
前記各サブセットに属する投影値は、前記Ｐ個の方位の中から、前記投影線の傾きの順にｐ個置きに（ｐ＜Ｐ）抽出された方向に対応した値である、画像処理装置。
【請求項３】
請求項２に記載の画像処理装置であって、
前記複数のサブセットは、
前記Ｐ個の方向のうち偶数番目の方向に対応する投影値からなる第１のサブセットと、
前記Ｐ個の方向のうち奇数番目の方向に対応する投影値からなる第２のサブセットと、
を含む、画像処理装置。
【請求項４】
請求項１から３のうちのいずれか１項に記載の画像処理装置であって、
前記光検出器は、互いに対向する一対の検出プレートを有し、
前記投影線は、前記一対の検出プレートのうち一方の検出器要素と、前記一対の検出プレートのうち他方の検出器要素とを結ぶ線である、画像処理装置。
【請求項５】
請求項４に記載の画像処理装置であって、前記検出器要素は、入射光に応じた蛍光出力を与えるシンチレータ結晶体を含む、画像処理装置。
【請求項６】
請求項１から５のうちのいずれか１項に記載の画像処理装置であって、前記画像再構成法は、ＯＳ−ＥＭ（Ordered Subsets - Expectation Maximization）法である、画像処理装置。
【請求項７】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理システムであって、
互いに対向する一対の検出プレートを有する光検出器と、
前記一対の検出プレートを結ぶ投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、
前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理システム。
【請求項８】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理方法であって、
（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出するステップと、
（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定するステップと、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理方法。
【請求項９】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成するために、コンピュータ読み取り可能な記録媒体から読み出されて画像再構成処理をコンピュータに実行させるプログラムであって、
前記画像再構成処理は、
（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する処理と、
（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する処理と、
を含み、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とするプログラム。
【請求項１】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理装置であって、
光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、
前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理装置。
【請求項２】
請求項１に記載の画像処理装置であって、
前記投影線の方位にはＰ個の方位（Ｐは正整数）が存在し、
前記各サブセットに属する投影値は、前記Ｐ個の方位の中から、前記投影線の傾きの順にｐ個置きに（ｐ＜Ｐ）抽出された方向に対応した値である、画像処理装置。
【請求項３】
請求項２に記載の画像処理装置であって、
前記複数のサブセットは、
前記Ｐ個の方向のうち偶数番目の方向に対応する投影値からなる第１のサブセットと、
前記Ｐ個の方向のうち奇数番目の方向に対応する投影値からなる第２のサブセットと、
を含む、画像処理装置。
【請求項４】
請求項１から３のうちのいずれか１項に記載の画像処理装置であって、
前記光検出器は、互いに対向する一対の検出プレートを有し、
前記投影線は、前記一対の検出プレートのうち一方の検出器要素と、前記一対の検出プレートのうち他方の検出器要素とを結ぶ線である、画像処理装置。
【請求項５】
請求項４に記載の画像処理装置であって、前記検出器要素は、入射光に応じた蛍光出力を与えるシンチレータ結晶体を含む、画像処理装置。
【請求項６】
請求項１から５のうちのいずれか１項に記載の画像処理装置であって、前記画像再構成法は、ＯＳ−ＥＭ（Ordered Subsets - Expectation Maximization）法である、画像処理装置。
【請求項７】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理システムであって、
互いに対向する一対の検出プレートを有する光検出器と、
前記一対の検出プレートを結ぶ投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する投影データ生成部と、
前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する濃度推定部と、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理システム。
【請求項８】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して、前記被測定物の内部構造を表す画像を生成する画像処理方法であって、
（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出するステップと、
（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定するステップと、
を備え、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とする画像処理方法。
【請求項９】
被測定物が置かれている領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成するために、コンピュータ読み取り可能な記録媒体から読み出されて画像再構成処理をコンピュータに実行させるプログラムであって、
前記画像再構成処理は、
（ａ）光検出器への投影線上で前記ボクセルから放出され当該光検出器で検出された光子の計数結果に基づいて前記投影線毎に投影値を算出する処理と、
（ｂ）前記投影値の集合を複数のサブセットに分割し、前記サブセット毎に、画像再構成法に従って前記投影値を使用した反復演算を実行することにより前記ボクセル毎の濃度を推定する処理と、
を含み、
前記各サブセットは、自己に属する投影値に対応する投影線の方位が均等に分布するように構成されていることを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【公開番号】特開２０１０−２０３８０７（Ｐ２０１０−２０３８０７Ａ）
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願番号】特願２００９−４７１７７（Ｐ２００９−４７１７７）
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
[ Back to top ]