画像処理装置、画像再構成システム、画像処理方法およびプログラム
【課題】画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置を提供する。
【解決手段】画像処理装置３では、前方投影算出部３３は、システムマトリクス要素と当該システムマトリクス要素に対応する濃度との積和演算を実行して投影線毎に前方投影値を算出する。濃度算出部３５は、前方投影値と実測投影値を用いて新たな濃度を算出する。投影線の配列とシステムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類される。前方投影算出部３３は、積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行する。
【解決手段】画像処理装置３では、前方投影算出部３３は、システムマトリクス要素と当該システムマトリクス要素に対応する濃度との積和演算を実行して投影線毎に前方投影値を算出する。濃度算出部３５は、前方投影値と実測投影値を用いて新たな濃度を算出する。投影線の配列とシステムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類される。前方投影算出部３３は、積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被測定物から放出された光を検出し、その検出結果に基づいて被測定物の内部の２次元画像または３次元画像を再構成するための画像処理技術に関する。
【背景技術】
【０００２】
人体などの被測定物を透過したＸ線、あるいは被測定物から放出されたγ線といった放射線を検出し、その検出結果に基づいて被測定物の内部構造の２次元画像または３次元画像を再構成するための画像処理技術の研究開発が進められている。このような画像処理技術は、Ｘ線ＣＴ装置（Ｘ線コンピュータ断層撮影装置）や陽電子放射断層撮像装置（ＰＥＴ：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.
【非特許文献４】A. Motta et al., "Fast 3D-EM reconstruction using Planograms for stationary planar positron emission mammography camera," Computerized Medical Imaging and Graphics, 2005, Vol. 29, pp. 587-596.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、たとえ従来の高速演算アルゴリズムを用いたとしても、計算量が膨大であるので、その計算速度は臨床上実用的であるとはいい難いという問題がある。たとえば、画像を再構成するまでの計算時間が長大であると、診断に要する時間が長くなるので、臨床上実用的であるとはいい難い。
【０００７】
上記に鑑みて本発明の目的は、画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置、画像処理方法およびプログラムを提供する点にある。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理装置が提供される。この画像処理装置は、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、を備える。前記濃度推定部は、前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する。
【０００９】
また、本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成する画像再構成システムが提供される。この画像再構成システムは、前記領域を挟み込む複数の検出器要素対を有する検出装置と、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、を備える。また、この画像再構成システムの濃度推定部は、前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する。
【００１０】
また、本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理方法が提供される。この画像処理方法は、（ａ）前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出するステップと、（ｂ）前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行するステップと、を備え、前記ステップ（ｂ）は、（ｂ−１）前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出するステップと、（ｂ−２）前記前方投影値および前記実測投影値を用いて新たな濃度を算出するステップと、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記ステップ（ｂ−１）では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される。
【００１１】
そして、本発明によれば、コンピュータ読み取り可能な記録媒体から読み出されてプロセッサに画像再構成処理を実行させるプログラムが提供される。前記画像再構成処理は、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成処理と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定処理と、を備える。前記濃度推定処理は、前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出処理と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出処理と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出処理では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される。
【発明の効果】
【００１２】
上述の通り、本発明による画像処理装置、画像再構成システム、画像処理方法および画像再構成処理では、いずれも、投影線は、複数種の幾何学的対称性に基づいた複数のグループに分類される。第１の積和演算は、当該第１の積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することにより実行される。これにより、前方投影値を短時間で算出し、各ボクセルの濃度を短時間で算出することが可能となる。
【図面の簡単な説明】
【００１３】
【図１】本発明に係る一実施形態の画像再構成システムの概略構成を示す機能ブロック図である。
【図２】投影線を説明するための模式図である。
【図３】画像処理装置による画像生成処理の手順を概略的に示すフローチャートである。
【図４】図３の画像再構成処理の手順を概略的に示すフローチャートである。
【図５】図３の前方投影算出処理の手順を概略的に示すフローチャートである。
【図６】図３の後方投影算出処理の手順を概略的に示すフローチャートである。
【図７】ボクセルの配列パターンを概略的に例示する図である。
【図８】投影線の配列の幾何学的対称性を説明するための概略図である。
【図９】（Ａ）は、傾きｄｘ＝３を有する検出器要素対（結晶対）と投影線（ハッチ部分）とを概略的に示す図であり、（Ｂ）は、傾きｄｘ＝−３を有する検出器要素対（結晶対）と投影線（ハッチ部分）とを概略的に示す図である。
【図１０】回転対称性を説明するための図である。
【図１１】ｄｘ−ｄｙ平面を概略的に示す図である。
【図１２】並進対称性を有する複数のシステムマトリクス要素を例示する図である。
【発明を実施するための形態】
【００１４】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
【００１５】
図１は、本発明に係る一実施形態の画像再構成システム１の概略構成を示す機能ブロック図である。
【００１６】
この画像再構成システム１は、図１に示すように被測定物２０が配置されるべき領域内に複数のボクセル（図示せず）を立体的に配列した解析モデルを使用して被測定物２０の内部構造を表す断層画像または３次元画像を生成（再構成）するものである。画像再構成システム１は、対向型ＰＥＴ装置（光検出装置）２および画像処理装置３を有する。
【００１７】
対向型ＰＥＴ装置２は、被測定物２０が配置された領域を挟み込む一対の検出プレート１０Ａ，１０Ｂを有し、これら検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂを含む。後述するように、被測定物２０から放出された光子対の一方は検出プレート１０Ａで検出され、当該光子対の他方は検出プレート１０Ｂで検出される。よって、被測定物２０から放出された光子対は、検出器要素１１Ａ，１１Ｂの対（検出器要素対）で検出される。検出プレート１０Ａ，１０Ｂは、光子対の検出結果を電気信号として画像処理装置３に供給する。画像処理装置３は、その検出結果に基づいて、被測定物２０の内部構造を表す断層画像または３次元画像を再構成する機能を有する。
【００１８】
図１に示されるように、対向型ＰＥＴ装置２は、互いに対向する検出面を有する一対の検出プレート１０Ａ，１０Ｂと、これら検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１とを有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００１９】
検出プレート１０Ａ，１０Ｂは、被測定物２０から放出された光子対をそれぞれ受光する検出面を有している。一方の検出プレート１０Ａは、格子状に配列された検出器要素１１Ａ，…，１１Ａを有しており、これら検出器要素１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。他方の検出プレート１０Ｂも、格子状に配列された検出器要素１１Ｂ，…，１１Ｂを有しており、これら検出器要素１１Ｂ，…，１１Ｂは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。検出器要素群１１Ａ，…，１１Ａの配列と検出器要素１１Ｂ，…，１１Ｂの配列とは、Ｘ−Ｙ平面に関して互いに対称である。また、検出器要素群１１Ａ，…，１１Ａの配列は、Ｘ軸方向とＹ軸方向に沿ってそれぞれ並進対称性を有する。すなわち、検出器要素群１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿って格子状に配列されているので、或る検出器要素１１ＡをＸ軸方向またはＹ軸方向に格子間距離の整数倍だけ移動させたとき、この検出器要素１１Ａと同じ構造を有する検出器要素１１Ａがその移動先に存在する。
【００２０】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に従って、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を操作して制御部４０に指示を入力することにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００２１】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、運動量保存則とエネルギー保存則とにより、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対の検出器要素１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対の検出器要素１１Ａ，１１Ｂを結ぶ線上に、陽電子が消滅した地点が存在する。
【００２２】
検出器要素１１Ａ，１１Ｂとしては、入射光に応じた蛍光出力を与えるシンチレータが使用される。シンチレータは、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯ（Ｂｉ４Ｇｅ３Ｏ１２）やＬＳＯ（Ｌｕ２（ＳｉＯ４）Ｏ：Ｃｅ）などの公知の結晶体を使用すればよい。検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、ケーブルなどの転送手段（図示せず）を介して、光子を検出した検出器要素１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００２３】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、検出器要素１１Ａ，１１Ｂが検出した光子対を計数し、この計数結果に基づいて投影線毎の実測投影データの値（実測投影値）を算出する。実測投影データは、メモリ３７内のデータ記憶部３８に格納される。図２は、投影線を説明するための模式図である。図２に示されるように、被測定物２０内の或る投影線ＰＬ上の陽電子消滅地点（線源）から放出された光子対が、検出プレート１０Ａ，１０Ｂの一対の検出器要素１１Ａ，１１Ｂで検出される。
【００２４】
投影データ生成部３１は、投影線ＰＬ上で放出された光子対を計数し、この計数結果に基づいて投影線ＰＬの投影値ｙｉを算出する。ここで、添え字ｉは、投影線を一意に表す番号である。投影線ＰＬは、一方の検出プレート１０Ａの検出器要素１１Ａと他方の検出プレート１０Ｂの検出器要素１１Ｂとを結ぶラインであるから、検出器要素１１Ａ，１１Ｂの組み合わせの数だけ投影線が存在することとなる。本実施形態では、ｎ個（ｎは正整数）の投影線ＰＬ（すなわち、検出器要素１１Ａ，１１Ｂのｎ個の組み合わせ）が存在するものとする（ｉ＝１〜ｎ）。たとえば、癌細胞のように活性の高い部位には線源が高密度で集中する傾向にあるので、投影線ＰＬ上に癌細胞が存在すると、この投影線ＰＬに対応する投影値ｙｉは比較的大きな値を示すこととなる。
【００２５】
濃度推定部３２は、被測定物２０が配置されている領域内に複数のボクセルＢ１〜Ｂｍ（ｍは正整数）が立体的に配列された解析モデルを利用して画像再構成処理を実行する。ここで、ボクセルＢｊは、立方形状または直方形状を有する微少な六面体要素からなる。
【００２６】
濃度推定部３２は、データ記憶部３８から実測投影値の集合｛ｙｉ｝を読み出し、統計的画像再構成法に従ってこれら実測投影値｛ｙｉ｝を用いてボクセル毎の放射能濃度（以下、単に「濃度」と呼ぶ。）を推定する処理を反復的に実行することにより、被測定物２０の内部構造の断層画像または３次元画像を再構成する機能を有している。本実施形態では、統計的画像再構成法として典型的なＭＬ−ＥＭ法が使用されるが、これに限定されるものではない。
【００２７】
図１に示されるように、濃度推定部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これら機能ブロック３３〜３５は、ＭＬ−ＥＭ法に従って、ボクセルの濃度を推定する処理を反復的に実行する。ＭＬ−ＥＭ法では、以下の推定式（１）が使用される。
【００２８】
【数１】
【００２９】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００３０】
【数２】
【００３１】
上記推定式（１）中、ｋは、反復回数であり、ｘｊ（ｋ−１）は、ｋ−１回目の処理におけるｊ番目のボクセルＢｊの濃度であり、ｙｉは、ｉ番目の検出器要素対（検出器要素１１Ａ，１１Ｂの対のうちｉ番目の対）による検出結果に基づいて得られた実測投影値である。また、ｘｊ（ｋ）は、（ｋ）回目の処理におけるｊ番目のボクセルＢｊの濃度である。上式（１）によれば、ｋ−１回目の濃度ｘｊ（ｋ−１）から、ｋ回目の濃度ｘｊ（ｋ）が算出される。初期値ｘｊ（ｋ＝０）は正の実数であれば、どのような値に設定されてもよい。
【００３２】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｊから放出された光子対がｉ番目の検出器要素対で検出される確率である。この確率Ｃｉｊは、検出確率と呼ばれており、この確率Ｃｉｊの集合はシステムマトリクスと呼ばれている行列である。確率Ｃｉｊは、システムマトリクスのｉ行ｊ列目の要素（以下、「システムマトリクス要素」と呼ぶ。）である。ｊ番目のボクセルＢｊから放出された光子対がｉ番目の検出器要素対に入射する数は、システムマトリクス要素Ｃｉｊに依存することが知られている。システムマトリクス要素ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００３３】
なお、システムマトリクス要素Ｃｉｊの値には、散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、システムマトリクス要素Ｃｉｊの値は、ｉ番目の検出器要素対に対するｊ番目のボクセルＢｊの幾何学的位置のみにより定められる。たとえば、ｉ番目の検出器要素対の一方から他方を結ぶ投影線ＰＬと、ｊ番目のボクセルＢｊとが重複する領域の体積またはこれに比例する値に、システムマトリクス要素Ｃｉｊの値を設定することができる。
【００３４】
次に、図３〜図６を参照しつつ、画像処理装置３の動作を以下に説明する。図３は、画像処理装置３による画像生成処理の手順を概略的に示すフローチャートであり、図４は、図３の画像再構成処理の手順を概略的に示すフローチャートである。また、図５は、図３の前方投影算出処理の手順を概略的に示すフローチャートであり、図６は、図３の後方投影算出処理の手順を概略的に示すフローチャートである。
【００３５】
図３を参照すると、先ず、ステップＳ１０では、濃度推定部３２は、パラメータ格納部３９から、システムマトリクス｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から実測投影値｛ｙｉ｝を読み出す（ステップＳ１１）。なお、本実施形態では、システムマトリクス｛Ｃｉｊ｝、係数｛Ｃｊ｝および実測投影値｛ｙｉ｝が最初に一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００３６】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２）、その後、画像再構成処理（ステップＳ１３）を実行する。
【００３７】
図４は、ステップＳ１３の画像再構成処理の手順を示すフローチャートである。図４を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化される（ステップＳ１３１）。その後、図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影データの値（前方投影値）｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００３８】
【数３】
【００３９】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影データの値（後方投影値）｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００４０】
【数４】
【００４１】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００４２】
その後、濃度算出部３５は、濃度｛ｘｊ｝が収束したか否かを判定する（ステップＳ１３５）。具体的には、対数尤度関数（または尤度関数）を用いて濃度｛ｘｊ｝の収束の有無を判定することが可能である。非特許文献１（L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, v. 1, pp. 113-122.）によれば、上式（１）は、ポアッソン分布の式により与えられる対数尤度関数Ｌ（｛ｘｊ｝）に基づいて導出される。濃度算出部３５は、ステップＳ１３４で算出された濃度｛ｘｊ｝を入力としたときの対数尤度関数Ｌ（｛ｘｊ｝）の出力値が飽和し略一定に達した場合に、濃度｛ｘｊ｝が収束したと判定することができる（ステップＳ１３５のＹＥＳ）。その後、処理は図３のステップＳ１４に移行する。そして、濃度｛ｘｊ｝が出力される（ステップＳ１４）。
【００４３】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、さらに、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされ（ステップＳ１３７）、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）は、処理は図３のステップＳ１４に移行する。そして、濃度｛ｘｊ｝が出力される（ステップＳ１４）。
【００４４】
その後、表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００４５】
上記推定式（１）中で演算量の多大な部分が、前方投影値ｐｉの計算（式（２））と後方投影値Ｄｊの計算（式（３））である。以下、前方投影値ｐｉの計算方法と後方投影値Ｄｊの計算方法について説明する。
【００４６】
本実施形態では、図７に例示されるように、複数のボクセルＢｉが検出器要素対１１Ａ，１１Ｂの間の領域に設けられている。また、検出器要素対１１Ａ，１１Ｂの一方から他方を結ぶ領域である投影線ＰＬが設けられている。検出器要素対１１Ａ，１１Ｂの組み合わせに応じて無数の投影線ＰＬが存在する。なお、これらボクセルＢｊや投影線ＰＬは現実に当該領域に設けられているのではなく、画像再構成処理のために仮想的に設けられたものである。このような投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、複数種の幾何学的対称性を共有している。本実施形態では、後述するように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、反転対称性および回転対称性を共有する。また、投影線ＰＬやシステムマトリクス要素｛Ｃｉｊ｝は、これら幾何学的対称性に基づいて複数のグループに分類される。
【００４７】
図８（Ａ），（Ｂ）および図９（Ａ），（Ｂ）は、投影線ＰＬの配列の幾何学的対称性を説明するための概略図である。図８（Ａ）には、Ｘ軸方向への傾きを有する投影線ＰＬが示されており、図８（Ｂ）には、Ｙ軸方向への傾きを有する投影線ＰＬが示されている。図８（Ａ）に示されるように、投影線ＰＬのＸ軸方向への傾きは、検出器要素１１ＡのＸ軸上の位置ｘ１と、検出器要素１１ＢのＸ軸上の位置ｘ２との差分ｄｘ＝ｘ１−ｘ２で表現することができる。また、図８（Ｂ）に示されるように、投影線ＰＬのＹ軸方向への傾きは、検出器要素１１ＢのＹ軸上の位置ｙ１と、検出器要素１１ＢのＹ軸上の位置ｙ２との差分ｄｙ＝ｙ１−ｙ２で表現することができる。したがって、投影線ＰＬの傾きは、ｄｘ−ｄｙ平面上の点（ｄｘ，ｄｙ）で表現することが可能である。言い換えれば、ｄｘ−ｄｙ平面上の点（ｄｘ，ｄｙ）は、当該点に対応する傾きを有する投影線ＰＬの集合を示すものである。
【００４８】
投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、複数種の反転対称性を共有する。たとえば、図９（Ａ）には、傾きｄｘ＝＋３を有する投影線ＰＬ，ＰＬ，…が示されており、図９（Ｂ）には、傾きｄｘ＝−３を有する投影線ＰＬ，ＰＬ，…が示されている。図９（Ａ）の投影線ＰＬと図９（Ｂ）の投影線ＰＬとは、Ｙ−Ｚ平面に対する反転対称の関係にある。ボクセルＢｉと投影線ＰＬとが重複する領域の体積またはこれに比例する値にシステムマトリクス要素Ｃｉｊの値を設定することができる。この場合、システムマトリクス要素Ｃｉｊの値は、投影線ＰＬの配列と共通する反転対称性を有する。また、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、Ｘ軸方向の並進対称性を共有するとともに、Ｙ軸方向の並進対称性を共有することが分かる。
【００４９】
さらに、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、図１の原点Ｏを基準としたときのＺ軸に関する回転対称性をも共有している。たとえば、図１０は、傾きｄｘ＝＋３とｄｙ＝０を有する投影線ＰＬと、傾きｄｘ＝０とｄｙ＝＋３を有する投影線ＰＬとを示す図である。図１０に示される一方の投影線ＰＬは、Ｚ軸の周りに回転されれば、他方の投影線ＰＬと一致する。
【００５０】
図１１は、ｄｘ−ｄｙ平面を概略的に示す図である。前述したように、図１１のｄｘ−ｄｙ平面上の各点が投影線ＰＬの傾きに対応する。図１１に示されるように、２点Ｐ１，Ｐ２は、ｄｙ軸に関して反転対称の関係にあり、２点Ｐ１，Ｐ３は、ｄｘ軸に関して反転対称の関係にあり、２点Ｐ１，Ｐ４は、ｄｘ−ｄｙ平面の原点Ｏに関して反転対称の関係にあり、２点Ｐ１，Ｐ５は、π／４（４５度）の回転対称の関係にあり、２点Ｐ２，Ｐ６は、π／４（４５度）の回転対称の関係にあり、２点Ｐ３，Ｐ７は、π／４（４５度）の回転対称の関係にあり、２点Ｐ４，Ｐ８は、π／４（４５度）の回転対称の関係にある。
【００５１】
よって、反転対称性と回転対称性との組み合わせにより、図１１に示されるように、ｄｘ−ｄｙ平面は、以下の８つの領域Ａ１〜Ａ８に区分けすることができる。領域Ａ１に対応するシステムマトリクス要素｛Ｃｉｊ｝が独立な成分であり、他の領域Ａ２〜Ａ８に対応するシステムマトリクス要素は、領域Ａ１に対応するシステムマトリクス要素｛Ｃｉｊ｝のコピーとみなすことができる。
【００５２】
・領域Ａ１：ｄｘ≧０、ｄｙ≧０、かつ、ｄｘ≧ｄｙ
・領域Ａ２：ｄｘ＜０、ｄｙ≧０、かつ、｜ｄｘ｜≧ｄｙ（領域Ａ１のｘ反転）
・領域Ａ３：ｄｘ≧０、ｄｙ＜０、かつ、ｄｘ≧｜ｄｙ｜（領域Ａ１のｙ反転）
・領域Ａ４：ｄｘ＜０、ｄｙ＜０、かつ、｜ｄｘ｜≧｜ｄｙ｜（領域Ａ１のｘ反転およびｙ反転）
・領域Ａ５：ｄｘ≧０、ｄｙ≧０、かつ、ｄｘ≦ｄｙ（領域Ａ１のπ／４回転）
・領域Ａ６：ｄｘ＜０、ｄｙ≧０、かつ、｜ｄｘ｜≦ｄｙ（領域Ａ２のπ／４回転）
・領域Ａ７：ｄｘ≧０、ｄｙ＜０、かつ、ｄｘ≦｜ｄｙ｜（領域Ａ３のπ／４回転）
・領域Ａ８：ｄｘ＜０、ｄｙ＜０、かつ、｜ｄｘ｜≦｜ｄｙ｜（領域Ａ４のπ／４回転）
【００５３】
したがって、投影線ＰＬは、領域Ａ１〜Ａ８にそれぞれ対応する８つのグループＧ１〜Ｇ８に分類できる。前方投影算出部３３は、上式（２）の第１の積和演算を投影線のグループＧ１〜Ｇ８にそれぞれ対応する８つのスレッドＴＰ１〜ＴＰ８に分割し、これらスレッドＴＰ１〜ＴＰ８を並列に実行することができる。
【００５４】
図５は、図３の前方投影算出処理の手順を概略的に示すフローチャートである。図５に示されるように、前方投影算出処理の手順は、８つのスレッドに分割される。ステップＳ２００Ａ〜Ｓ２００Ｂのループ処理はスレッドＴＰ１に相当し、ステップＳ２１０Ａ〜Ｓ２１０Ｂのループ処理はスレッドＴＰ２に相当し、ステップＳ２２０Ａ〜Ｓ２２０Ｂのループ処理はスレッドＴＰ３に相当し、ステップＳ２３０Ａ〜Ｓ２３０Ｂのループ処理はスレッドＴＰ４に相当し、ステップＳ２４０Ａ〜Ｓ２４０Ｂのループ処理はスレッドＴＰ５に相当し、ステップＳ２５０Ａ〜Ｓ２５０Ｂのループ処理はスレッドＴＰ６に相当し、 ステップＳ２６０Ａ〜Ｓ２６０Ｂのループ処理はスレッドＴＰ７に相当し、ステップＳ２７０Ａ〜Ｓ２７０Ｂのループ処理はスレッドＴＰ８に相当する。
【００５５】
スレッドＴＰ１〜ＴＰ７では、ｄｘ，ｄｙの値がそれぞれインクリメントにより指定される（ステップＳ２００Ａ，Ｓ２１０Ａ，Ｓ２２０Ａ，Ｓ２３０Ａ，Ｓ２４０Ａ，Ｓ２５０Ａ，Ｓ２６０Ａ，Ｓ２７０Ａ）。ただし、当該ステップの最初の処理では、インクリメントが実行される代わりに、ｄｘ，ｄｙの値として初期値が与えられる。
【００５６】
次いで、ｄｘ，ｄｙの組により、システムマトリクス要素の値Ｃｘ＝Ｃｉｊ（ｄｘ，ｄｙ）が指定され、読み出される（ステップＳ２０１，Ｓ２１１，Ｓ２２１，Ｓ２３１，Ｓ２４１，Ｓ２５１，Ｓ２６１，Ｓ２７１）。
【００５７】
次に、前方投影算出部３３は、システムマトリクス要素の共通の値Ｃｘを用いたベクトル演算により第１の積和演算を実行する（ステップＳ２０２，Ｓ２１２，Ｓ２２２，Ｓ２３２，Ｓ２４２，Ｓ２５２，Ｓ２６２，Ｓ２７２）。上述したように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、Ｘ軸方向の並進対称性とＹ軸方向の並進対称性とを共有している。よって、これら並進対称性に基づいて、同一値Ｃｘを共有する複数のシステムマトリクス要素が必ず存在する。前方投影算出部３３は、この共通の値Ｃｘを用いたベクトル演算を実行することができる。
【００５８】
たとえば、図１２に示されるように、Ｘ軸方向またはＹ軸方向への同一の傾きを有する複数の投影線ＰＬ，ＰＬ，…上におけるボクセルＢ１，Ｂ３，Ｂ５，Ｂ７，Ｂ９，Ｂ１１，Ｂ１３，Ｂ１５に対応するシステムマトリクス要素は、Ｘ軸方向またはＹ軸方向の並進対称性を有する。よって、これらボクセルＢ１，Ｂ３，Ｂ５，Ｂ７，Ｂ９，Ｂ１１，Ｂ１３，Ｂ１５に対応するシステムマトリクス要素は共通の値を有している。同様に、複数の投影線ＰＬ，ＰＬ，…上におけるボクセルＢ２，Ｂ４，Ｂ６，Ｂ８，Ｂ１０，Ｂ１２，Ｂ１４，Ｂ１６に対応するシステムマトリクス要素も、共通の値を有している。
【００５９】
たとえば、前方投影値ｐ１，ｐ２，ｐ３，ｐ４が、システムマトリクス要素の共通の値Ｃ１，Ｃ２，…，ＣＮを用いて以下のように計算されるものとする。
【００６０】
ｐ１＝Ｃ１・ｘ１＋Ｃ２・ｘｊ＋…
ｐ２＝Ｃ１・ｘ２＋Ｃ２・ｘｊ＋１＋…
ｐ３＝Ｃ１・ｘ３＋Ｃ２・ｘｊ＋２＋…
ｐ４＝Ｃ１・ｘ４＋Ｃ２・ｘｊ＋３＋…
【００６１】
このとき、前方投影算出部３３は、上式の右辺第１項の計算を、（ｘ１、ｘ２、ｘ３、ｘ４）をベクトルとし、このベクトルに係数Ｃ１を乗算することで行うことができる。続けて、前方投影算出部３３は、上式の右辺第２項の計算を、（ｘｊ、ｘｊ＋１、ｘｊ＋２、ｘｊ＋３）をベクトルとし、このベクトルに係数Ｃ２を乗算することで行うことができる。なお、ベクトルの要素の個数は、４個に限らず、ｎ個（ｎは５以上）であってよい。
【００６２】
スレッドＴＰ１〜ＴＰ８のループ処理は、ｄｘ，ｄｙの値の全ての組み合わせについて実行された後に終了する。その後、図４のステップＳ１３３に処理が移行する。
【００６３】
また、後方投影算出部３４も、上式（３）の第２の積和演算を８つのグループＧ１〜Ｇ８にそれぞれ対応する８つのスレッドＴＢ１〜ＴＢ８に分割し、これらスレッドＴＢ１〜ＴＢ８を並列に実行することができる。
【００６４】
図６は、図３の後方投影算出処理の手順を概略的に示すフローチャートである。図６に示されるように、後方投影算出処理の手順は、８つのスレッドに分割される。ステップＳ３００Ａ〜Ｓ３００Ｂのループ処理はスレッドＴＢ１に相当し、ステップＳ３１０Ａ〜Ｓ３１０Ｂのループ処理はスレッドＴＢ２に相当し、ステップＳ３２０Ａ〜Ｓ３２０Ｂのループ処理はスレッドＴＢ３に相当し、ステップＳ３３０Ａ〜Ｓ３３０Ｂのループ処理はスレッドＴＢ４に相当し、ステップＳ３４０Ａ〜Ｓ３４０Ｂのループ処理はスレッドＴＢ５に相当し、ステップＳ３５０Ａ〜Ｓ３５０Ｂのループ処理はスレッドＴＢ６に相当し、ステップＳ３６０Ａ〜Ｓ３６０Ｂのループ処理はスレッドＴＢ７に相当し、ステップＳ３７０Ａ〜Ｓ３７０Ｂのループ処理はスレッドＴＢ８に相当する。
【００６５】
スレッドＴＢ１〜ＴＢ８では、ｄｘ，ｄｙの値がそれぞれインクリメントにより指定される（ステップＳ３００Ａ，Ｓ３１０Ａ，Ｓ３２０Ａ，Ｓ３３０Ａ，Ｓ３４０Ａ，Ｓ３５０Ａ，Ｓ３６０Ａ，Ｓ３７０Ａ）。ただし、当該ステップの最初の処理では、インクリメントが実行される代わりに、ｄｘ，ｄｙの値として初期値が与えられる。
【００６６】
次いで、ｄｘ，ｄｙの組により、システムマトリクス要素の値Ｃｘ＝Ｃｉｊ（ｄｘ，ｄｙ）が読み出される（ステップＳ３０１，Ｓ３１１，Ｓ３２１，Ｓ３３１，Ｓ３４１，Ｓ３５１，Ｓ３６１，Ｓ３７１）。
【００６７】
次に、後方投影算出部３４は、システムマトリクス要素の共通の値Ｃｘを用いたベクトル演算により第１の積和演算を実行する（ステップＳ３０２，Ｓ３１２，Ｓ３２２，Ｓ３３２，Ｓ３４２，Ｓ３５２，Ｓ３６２，Ｓ３７２）。前述したように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とが共有する並進対称性に基づいて、システムマトリクス要素の値Ｃｘを共有する複数のシステムマトリクス要素が必ず存在する。後方投影算出部３４は、この共通の値Ｃｘを用いたベクトル演算を実行することができる。
【００６８】
たとえば、後方投影値Ｄ１，Ｄ２，Ｄ３，Ｄ４が、システムマトリクス要素の共通の値Ｃ１，Ｃ２，…，ＣＮを用いて以下のように計算されるものとする。
【００６９】
Ｄ１＝Ｃ１・ｑ１＋Ｃ２・ｑｉ＋…
Ｄ２＝Ｃ１・ｑ２＋Ｃ２・ｑｉ＋１＋…
Ｄ３＝Ｃ１・ｑ３＋Ｃ２・ｑｉ＋２＋…
Ｄ４＝Ｃ１・ｑ４＋Ｃ２・ｑｉ＋３＋…
【００７０】
このとき、後方投影算出部３４は、上式の右辺第１項の計算を、（ｑ１、ｑ２、ｑ３、ｑ４）をベクトルとし、このベクトルに係数Ｃ１を乗算することで行うことができる。続けて、後方投影算出部３４は、上式の右辺第２項の計算を、（ｑｉ、ｑｉ＋１、ｑｉ＋２、ｑｉ＋３）をベクトルとし、このベクトルに係数Ｃ２を乗算することで行うことができる。なお、ベクトルの要素の個数は、４個に限らず、ｍ個（ｍは５以上）であってもよい。
【００７１】
スレッドＴＢ１〜ＴＢ８のループ処理は、ｄｘ，ｄｙの値の全ての組み合わせについて実行された後に終了する。その後、図４のステップＳ１３４に処理が移行する。
【００７２】
上述した通り、画像再構成システム１は、投影線ＰＬを、投影線の配列とシステムマトリクス要素の値とが共有する複数種の幾何学的対称性に基づいて複数のグループに分類する。前方投影算出部３３は、第１の積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することができる。これにより、前方投影値を短時間で算出することができるので、各ボクセルの濃度を短時間で算出することが可能となる。同様に、後方投影算出部３４は、第２の積和演算を複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することができる。これにより、後方投影値を短時間で算出することができるので、各ボクセルの濃度を短時間で算出することが可能となる。
【００７３】
さらに、投影線ＰＬの配列とシステムマトリクス要素の値とは同一方向の並進対称性を共有するので、この並進対称性に基づいた同一値を有する複数のシステムマトリクス要素が存在する。前方投影算出部３３と後方投影算出部３４は、このシステムマトリクス要素を用いたベクトル演算を実行して前方投影値と後方投影値を短時間で算出することができる。
【００７４】
ところで、図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。ここで、プロセッサは、上記ベクトル演算を実行するためのベクトル演算機能を有する。
【００７５】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【００７６】
以上、図面を参照して本発明の実施形態について述べたが、これらは本発明の例示であり、上記以外の様々な構成を採用することもできる。
【００７７】
たとえば、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１の対向型ＰＥＴ装置２に代えて、連続的に分布する単一の検出面あるいは３個以上の検出面を有する検出装置を使用してもよい。また、検出面の形状は平面に限定されるものではなく曲面であってもよい。
【００７８】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置２に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
【００７９】
上記実施形態では、統計的画像再構成法としてＭＬ−ＥＭ法が使用されているが、システムマトリクス要素（検出確率）を用いた画像再構成法を使用すればよく、ＭＬ−ＥＭ法に限定されるものではない。たとえば、ＯＳ−ＥＭ法やＭＡＰ−ＥＭ（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.
【非特許文献４】A. Motta et al., "Fast 3D-EM reconstruction using Planograms for stationary planar positron emission mammography camera," Computerized Medical Imaging and Graphics, 2005, Vol. 29, pp. 587-596.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、たとえ従来の高速演算アルゴリズムを用いたとしても、計算量が膨大であるので、その計算速度は臨床上実用的であるとはいい難いという問題がある。たとえば、画像を再構成するまでの計算時間が長大であると、診断に要する時間が長くなるので、臨床上実用的であるとはいい難い。
【０００７】
上記に鑑みて本発明の目的は、画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置、画像処理方法およびプログラムを提供する点にある。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理装置が提供される。この画像処理装置は、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、を備える。前記濃度推定部は、前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する。
【０００９】
また、本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成する画像再構成システムが提供される。この画像再構成システムは、前記領域を挟み込む複数の検出器要素対を有する検出装置と、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、を備える。また、この画像再構成システムの濃度推定部は、前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する。
【００１０】
また、本発明によれば、被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理方法が提供される。この画像処理方法は、（ａ）前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出するステップと、（ｂ）前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行するステップと、を備え、前記ステップ（ｂ）は、（ｂ−１）前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出するステップと、（ｂ−２）前記前方投影値および前記実測投影値を用いて新たな濃度を算出するステップと、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記ステップ（ｂ−１）では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される。
【００１１】
そして、本発明によれば、コンピュータ読み取り可能な記録媒体から読み出されてプロセッサに画像再構成処理を実行させるプログラムが提供される。前記画像再構成処理は、前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成処理と、前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定処理と、を備える。前記濃度推定処理は、前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出処理と、前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出処理と、を含み、前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、前記前方投影算出処理では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される。
【発明の効果】
【００１２】
上述の通り、本発明による画像処理装置、画像再構成システム、画像処理方法および画像再構成処理では、いずれも、投影線は、複数種の幾何学的対称性に基づいた複数のグループに分類される。第１の積和演算は、当該第１の積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することにより実行される。これにより、前方投影値を短時間で算出し、各ボクセルの濃度を短時間で算出することが可能となる。
【図面の簡単な説明】
【００１３】
【図１】本発明に係る一実施形態の画像再構成システムの概略構成を示す機能ブロック図である。
【図２】投影線を説明するための模式図である。
【図３】画像処理装置による画像生成処理の手順を概略的に示すフローチャートである。
【図４】図３の画像再構成処理の手順を概略的に示すフローチャートである。
【図５】図３の前方投影算出処理の手順を概略的に示すフローチャートである。
【図６】図３の後方投影算出処理の手順を概略的に示すフローチャートである。
【図７】ボクセルの配列パターンを概略的に例示する図である。
【図８】投影線の配列の幾何学的対称性を説明するための概略図である。
【図９】（Ａ）は、傾きｄｘ＝３を有する検出器要素対（結晶対）と投影線（ハッチ部分）とを概略的に示す図であり、（Ｂ）は、傾きｄｘ＝−３を有する検出器要素対（結晶対）と投影線（ハッチ部分）とを概略的に示す図である。
【図１０】回転対称性を説明するための図である。
【図１１】ｄｘ−ｄｙ平面を概略的に示す図である。
【図１２】並進対称性を有する複数のシステムマトリクス要素を例示する図である。
【発明を実施するための形態】
【００１４】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
【００１５】
図１は、本発明に係る一実施形態の画像再構成システム１の概略構成を示す機能ブロック図である。
【００１６】
この画像再構成システム１は、図１に示すように被測定物２０が配置されるべき領域内に複数のボクセル（図示せず）を立体的に配列した解析モデルを使用して被測定物２０の内部構造を表す断層画像または３次元画像を生成（再構成）するものである。画像再構成システム１は、対向型ＰＥＴ装置（光検出装置）２および画像処理装置３を有する。
【００１７】
対向型ＰＥＴ装置２は、被測定物２０が配置された領域を挟み込む一対の検出プレート１０Ａ，１０Ｂを有し、これら検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂを含む。後述するように、被測定物２０から放出された光子対の一方は検出プレート１０Ａで検出され、当該光子対の他方は検出プレート１０Ｂで検出される。よって、被測定物２０から放出された光子対は、検出器要素１１Ａ，１１Ｂの対（検出器要素対）で検出される。検出プレート１０Ａ，１０Ｂは、光子対の検出結果を電気信号として画像処理装置３に供給する。画像処理装置３は、その検出結果に基づいて、被測定物２０の内部構造を表す断層画像または３次元画像を再構成する機能を有する。
【００１８】
図１に示されるように、対向型ＰＥＴ装置２は、互いに対向する検出面を有する一対の検出プレート１０Ａ，１０Ｂと、これら検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１とを有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００１９】
検出プレート１０Ａ，１０Ｂは、被測定物２０から放出された光子対をそれぞれ受光する検出面を有している。一方の検出プレート１０Ａは、格子状に配列された検出器要素１１Ａ，…，１１Ａを有しており、これら検出器要素１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。他方の検出プレート１０Ｂも、格子状に配列された検出器要素１１Ｂ，…，１１Ｂを有しており、これら検出器要素１１Ｂ，…，１１Ｂは、Ｘ−Ｙ平面に沿ってマトリクス状に配列されている。検出器要素群１１Ａ，…，１１Ａの配列と検出器要素１１Ｂ，…，１１Ｂの配列とは、Ｘ−Ｙ平面に関して互いに対称である。また、検出器要素群１１Ａ，…，１１Ａの配列は、Ｘ軸方向とＹ軸方向に沿ってそれぞれ並進対称性を有する。すなわち、検出器要素群１１Ａ，…，１１Ａは、Ｘ−Ｙ平面に沿って格子状に配列されているので、或る検出器要素１１ＡをＸ軸方向またはＹ軸方向に格子間距離の整数倍だけ移動させたとき、この検出器要素１１Ａと同じ構造を有する検出器要素１１Ａがその移動先に存在する。
【００２０】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に従って、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を操作して制御部４０に指示を入力することにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００２１】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、運動量保存則とエネルギー保存則とにより、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対の検出器要素１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対の検出器要素１１Ａ，１１Ｂを結ぶ線上に、陽電子が消滅した地点が存在する。
【００２２】
検出器要素１１Ａ，１１Ｂとしては、入射光に応じた蛍光出力を与えるシンチレータが使用される。シンチレータは、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯ（Ｂｉ４Ｇｅ３Ｏ１２）やＬＳＯ（Ｌｕ２（ＳｉＯ４）Ｏ：Ｃｅ）などの公知の結晶体を使用すればよい。検出プレート１０Ａ，１０Ｂは、それぞれ、検出器要素１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、ケーブルなどの転送手段（図示せず）を介して、光子を検出した検出器要素１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００２３】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、検出器要素１１Ａ，１１Ｂが検出した光子対を計数し、この計数結果に基づいて投影線毎の実測投影データの値（実測投影値）を算出する。実測投影データは、メモリ３７内のデータ記憶部３８に格納される。図２は、投影線を説明するための模式図である。図２に示されるように、被測定物２０内の或る投影線ＰＬ上の陽電子消滅地点（線源）から放出された光子対が、検出プレート１０Ａ，１０Ｂの一対の検出器要素１１Ａ，１１Ｂで検出される。
【００２４】
投影データ生成部３１は、投影線ＰＬ上で放出された光子対を計数し、この計数結果に基づいて投影線ＰＬの投影値ｙｉを算出する。ここで、添え字ｉは、投影線を一意に表す番号である。投影線ＰＬは、一方の検出プレート１０Ａの検出器要素１１Ａと他方の検出プレート１０Ｂの検出器要素１１Ｂとを結ぶラインであるから、検出器要素１１Ａ，１１Ｂの組み合わせの数だけ投影線が存在することとなる。本実施形態では、ｎ個（ｎは正整数）の投影線ＰＬ（すなわち、検出器要素１１Ａ，１１Ｂのｎ個の組み合わせ）が存在するものとする（ｉ＝１〜ｎ）。たとえば、癌細胞のように活性の高い部位には線源が高密度で集中する傾向にあるので、投影線ＰＬ上に癌細胞が存在すると、この投影線ＰＬに対応する投影値ｙｉは比較的大きな値を示すこととなる。
【００２５】
濃度推定部３２は、被測定物２０が配置されている領域内に複数のボクセルＢ１〜Ｂｍ（ｍは正整数）が立体的に配列された解析モデルを利用して画像再構成処理を実行する。ここで、ボクセルＢｊは、立方形状または直方形状を有する微少な六面体要素からなる。
【００２６】
濃度推定部３２は、データ記憶部３８から実測投影値の集合｛ｙｉ｝を読み出し、統計的画像再構成法に従ってこれら実測投影値｛ｙｉ｝を用いてボクセル毎の放射能濃度（以下、単に「濃度」と呼ぶ。）を推定する処理を反復的に実行することにより、被測定物２０の内部構造の断層画像または３次元画像を再構成する機能を有している。本実施形態では、統計的画像再構成法として典型的なＭＬ−ＥＭ法が使用されるが、これに限定されるものではない。
【００２７】
図１に示されるように、濃度推定部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これら機能ブロック３３〜３５は、ＭＬ−ＥＭ法に従って、ボクセルの濃度を推定する処理を反復的に実行する。ＭＬ−ＥＭ法では、以下の推定式（１）が使用される。
【００２８】
【数１】
【００２９】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００３０】
【数２】
【００３１】
上記推定式（１）中、ｋは、反復回数であり、ｘｊ（ｋ−１）は、ｋ−１回目の処理におけるｊ番目のボクセルＢｊの濃度であり、ｙｉは、ｉ番目の検出器要素対（検出器要素１１Ａ，１１Ｂの対のうちｉ番目の対）による検出結果に基づいて得られた実測投影値である。また、ｘｊ（ｋ）は、（ｋ）回目の処理におけるｊ番目のボクセルＢｊの濃度である。上式（１）によれば、ｋ−１回目の濃度ｘｊ（ｋ−１）から、ｋ回目の濃度ｘｊ（ｋ）が算出される。初期値ｘｊ（ｋ＝０）は正の実数であれば、どのような値に設定されてもよい。
【００３２】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｊから放出された光子対がｉ番目の検出器要素対で検出される確率である。この確率Ｃｉｊは、検出確率と呼ばれており、この確率Ｃｉｊの集合はシステムマトリクスと呼ばれている行列である。確率Ｃｉｊは、システムマトリクスのｉ行ｊ列目の要素（以下、「システムマトリクス要素」と呼ぶ。）である。ｊ番目のボクセルＢｊから放出された光子対がｉ番目の検出器要素対に入射する数は、システムマトリクス要素Ｃｉｊに依存することが知られている。システムマトリクス要素ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００３３】
なお、システムマトリクス要素Ｃｉｊの値には、散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、システムマトリクス要素Ｃｉｊの値は、ｉ番目の検出器要素対に対するｊ番目のボクセルＢｊの幾何学的位置のみにより定められる。たとえば、ｉ番目の検出器要素対の一方から他方を結ぶ投影線ＰＬと、ｊ番目のボクセルＢｊとが重複する領域の体積またはこれに比例する値に、システムマトリクス要素Ｃｉｊの値を設定することができる。
【００３４】
次に、図３〜図６を参照しつつ、画像処理装置３の動作を以下に説明する。図３は、画像処理装置３による画像生成処理の手順を概略的に示すフローチャートであり、図４は、図３の画像再構成処理の手順を概略的に示すフローチャートである。また、図５は、図３の前方投影算出処理の手順を概略的に示すフローチャートであり、図６は、図３の後方投影算出処理の手順を概略的に示すフローチャートである。
【００３５】
図３を参照すると、先ず、ステップＳ１０では、濃度推定部３２は、パラメータ格納部３９から、システムマトリクス｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。続けて、濃度推定部３２は、データ記憶部３８から実測投影値｛ｙｉ｝を読み出す（ステップＳ１１）。なお、本実施形態では、システムマトリクス｛Ｃｉｊ｝、係数｛Ｃｊ｝および実測投影値｛ｙｉ｝が最初に一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００３６】
続けて、濃度推定部３２は、放射能濃度｛ｘｊ｝の値を初期化し（ステップＳ１２）、その後、画像再構成処理（ステップＳ１３）を実行する。
【００３７】
図４は、ステップＳ１３の画像再構成処理の手順を示すフローチャートである。図４を参照すると、先ず、濃度推定部３２は、反復回数ｋの上限を適切な値に設定する（ステップＳ１３０）。次に、反復回数ｋの値が１に初期化される（ステップＳ１３１）。その後、図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｒと推定濃度ｘｒ（ｋ−１）との積和演算を実行して前方投影データの値（前方投影値）｛ｐｉ｝（＝｛ｐ１，ｐ２，…｝）を算出する（ステップＳ１３２）。
【００３８】
【数３】
【００３９】
上記ステップＳ１３２の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、投影値ｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算を実行して後方投影データの値（後方投影値）｛Ｄｊ｝＝｛Ｄ１，Ｄ２，…，Ｄｍ｝を算出する（ステップＳ１３３）。
【００４０】
【数４】
【００４１】
そして、濃度算出部３５は、ｊ＝１〜ｍについて、後方投影値Ｄｊを用いて、ｋ回目の推定濃度｛ｘｊ｝＝｛ｘｊ（ｋ）｝＝｛Ｄｊ・ｘｊ（ｋ−１）／Ｃｊ｝を算出する（ステップＳ１３４）。
【００４２】
その後、濃度算出部３５は、濃度｛ｘｊ｝が収束したか否かを判定する（ステップＳ１３５）。具体的には、対数尤度関数（または尤度関数）を用いて濃度｛ｘｊ｝の収束の有無を判定することが可能である。非特許文献１（L. A. Shepp and Y. Vardi, "Maximum likelihood reconstruction for emission tomography," IEEE Trans. Med. Imaging, 1982, v. 1, pp. 113-122.）によれば、上式（１）は、ポアッソン分布の式により与えられる対数尤度関数Ｌ（｛ｘｊ｝）に基づいて導出される。濃度算出部３５は、ステップＳ１３４で算出された濃度｛ｘｊ｝を入力としたときの対数尤度関数Ｌ（｛ｘｊ｝）の出力値が飽和し略一定に達した場合に、濃度｛ｘｊ｝が収束したと判定することができる（ステップＳ１３５のＹＥＳ）。その後、処理は図３のステップＳ１４に移行する。そして、濃度｛ｘｊ｝が出力される（ステップＳ１４）。
【００４３】
一方、濃度｛ｘｊ｝が収束しないと判定された場合（ステップＳ１３５のＮＯ）、濃度算出部３５は、さらに、反復回数ｋが上限に達したか否かを判定する（ステップＳ１３６）。反復回数ｋが上限に達しない場合（ステップＳ１３６のＮＯ）は、反復回数ｋがインクリメントされ（ステップＳ１３７）、ステップＳ１３２の処理が反復される。濃度｛ｘｊ｝が収束しないまま、最終的に反復回数ｋが上限に達したとき（ステップＳ１３６のＹＥＳ）は、処理は図３のステップＳ１４に移行する。そして、濃度｛ｘｊ｝が出力される（ステップＳ１４）。
【００４４】
その後、表示制御部３６は、濃度算出部３５から出力された推定放射能濃度｛ｘｊ｝から２次元または３次元の再構成画像を生成し、この再構成画像をディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力された推定濃度値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００４５】
上記推定式（１）中で演算量の多大な部分が、前方投影値ｐｉの計算（式（２））と後方投影値Ｄｊの計算（式（３））である。以下、前方投影値ｐｉの計算方法と後方投影値Ｄｊの計算方法について説明する。
【００４６】
本実施形態では、図７に例示されるように、複数のボクセルＢｉが検出器要素対１１Ａ，１１Ｂの間の領域に設けられている。また、検出器要素対１１Ａ，１１Ｂの一方から他方を結ぶ領域である投影線ＰＬが設けられている。検出器要素対１１Ａ，１１Ｂの組み合わせに応じて無数の投影線ＰＬが存在する。なお、これらボクセルＢｊや投影線ＰＬは現実に当該領域に設けられているのではなく、画像再構成処理のために仮想的に設けられたものである。このような投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、複数種の幾何学的対称性を共有している。本実施形態では、後述するように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、反転対称性および回転対称性を共有する。また、投影線ＰＬやシステムマトリクス要素｛Ｃｉｊ｝は、これら幾何学的対称性に基づいて複数のグループに分類される。
【００４７】
図８（Ａ），（Ｂ）および図９（Ａ），（Ｂ）は、投影線ＰＬの配列の幾何学的対称性を説明するための概略図である。図８（Ａ）には、Ｘ軸方向への傾きを有する投影線ＰＬが示されており、図８（Ｂ）には、Ｙ軸方向への傾きを有する投影線ＰＬが示されている。図８（Ａ）に示されるように、投影線ＰＬのＸ軸方向への傾きは、検出器要素１１ＡのＸ軸上の位置ｘ１と、検出器要素１１ＢのＸ軸上の位置ｘ２との差分ｄｘ＝ｘ１−ｘ２で表現することができる。また、図８（Ｂ）に示されるように、投影線ＰＬのＹ軸方向への傾きは、検出器要素１１ＢのＹ軸上の位置ｙ１と、検出器要素１１ＢのＹ軸上の位置ｙ２との差分ｄｙ＝ｙ１−ｙ２で表現することができる。したがって、投影線ＰＬの傾きは、ｄｘ−ｄｙ平面上の点（ｄｘ，ｄｙ）で表現することが可能である。言い換えれば、ｄｘ−ｄｙ平面上の点（ｄｘ，ｄｙ）は、当該点に対応する傾きを有する投影線ＰＬの集合を示すものである。
【００４８】
投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、複数種の反転対称性を共有する。たとえば、図９（Ａ）には、傾きｄｘ＝＋３を有する投影線ＰＬ，ＰＬ，…が示されており、図９（Ｂ）には、傾きｄｘ＝−３を有する投影線ＰＬ，ＰＬ，…が示されている。図９（Ａ）の投影線ＰＬと図９（Ｂ）の投影線ＰＬとは、Ｙ−Ｚ平面に対する反転対称の関係にある。ボクセルＢｉと投影線ＰＬとが重複する領域の体積またはこれに比例する値にシステムマトリクス要素Ｃｉｊの値を設定することができる。この場合、システムマトリクス要素Ｃｉｊの値は、投影線ＰＬの配列と共通する反転対称性を有する。また、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、Ｘ軸方向の並進対称性を共有するとともに、Ｙ軸方向の並進対称性を共有することが分かる。
【００４９】
さらに、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、図１の原点Ｏを基準としたときのＺ軸に関する回転対称性をも共有している。たとえば、図１０は、傾きｄｘ＝＋３とｄｙ＝０を有する投影線ＰＬと、傾きｄｘ＝０とｄｙ＝＋３を有する投影線ＰＬとを示す図である。図１０に示される一方の投影線ＰＬは、Ｚ軸の周りに回転されれば、他方の投影線ＰＬと一致する。
【００５０】
図１１は、ｄｘ−ｄｙ平面を概略的に示す図である。前述したように、図１１のｄｘ−ｄｙ平面上の各点が投影線ＰＬの傾きに対応する。図１１に示されるように、２点Ｐ１，Ｐ２は、ｄｙ軸に関して反転対称の関係にあり、２点Ｐ１，Ｐ３は、ｄｘ軸に関して反転対称の関係にあり、２点Ｐ１，Ｐ４は、ｄｘ−ｄｙ平面の原点Ｏに関して反転対称の関係にあり、２点Ｐ１，Ｐ５は、π／４（４５度）の回転対称の関係にあり、２点Ｐ２，Ｐ６は、π／４（４５度）の回転対称の関係にあり、２点Ｐ３，Ｐ７は、π／４（４５度）の回転対称の関係にあり、２点Ｐ４，Ｐ８は、π／４（４５度）の回転対称の関係にある。
【００５１】
よって、反転対称性と回転対称性との組み合わせにより、図１１に示されるように、ｄｘ−ｄｙ平面は、以下の８つの領域Ａ１〜Ａ８に区分けすることができる。領域Ａ１に対応するシステムマトリクス要素｛Ｃｉｊ｝が独立な成分であり、他の領域Ａ２〜Ａ８に対応するシステムマトリクス要素は、領域Ａ１に対応するシステムマトリクス要素｛Ｃｉｊ｝のコピーとみなすことができる。
【００５２】
・領域Ａ１：ｄｘ≧０、ｄｙ≧０、かつ、ｄｘ≧ｄｙ
・領域Ａ２：ｄｘ＜０、ｄｙ≧０、かつ、｜ｄｘ｜≧ｄｙ（領域Ａ１のｘ反転）
・領域Ａ３：ｄｘ≧０、ｄｙ＜０、かつ、ｄｘ≧｜ｄｙ｜（領域Ａ１のｙ反転）
・領域Ａ４：ｄｘ＜０、ｄｙ＜０、かつ、｜ｄｘ｜≧｜ｄｙ｜（領域Ａ１のｘ反転およびｙ反転）
・領域Ａ５：ｄｘ≧０、ｄｙ≧０、かつ、ｄｘ≦ｄｙ（領域Ａ１のπ／４回転）
・領域Ａ６：ｄｘ＜０、ｄｙ≧０、かつ、｜ｄｘ｜≦ｄｙ（領域Ａ２のπ／４回転）
・領域Ａ７：ｄｘ≧０、ｄｙ＜０、かつ、ｄｘ≦｜ｄｙ｜（領域Ａ３のπ／４回転）
・領域Ａ８：ｄｘ＜０、ｄｙ＜０、かつ、｜ｄｘ｜≦｜ｄｙ｜（領域Ａ４のπ／４回転）
【００５３】
したがって、投影線ＰＬは、領域Ａ１〜Ａ８にそれぞれ対応する８つのグループＧ１〜Ｇ８に分類できる。前方投影算出部３３は、上式（２）の第１の積和演算を投影線のグループＧ１〜Ｇ８にそれぞれ対応する８つのスレッドＴＰ１〜ＴＰ８に分割し、これらスレッドＴＰ１〜ＴＰ８を並列に実行することができる。
【００５４】
図５は、図３の前方投影算出処理の手順を概略的に示すフローチャートである。図５に示されるように、前方投影算出処理の手順は、８つのスレッドに分割される。ステップＳ２００Ａ〜Ｓ２００Ｂのループ処理はスレッドＴＰ１に相当し、ステップＳ２１０Ａ〜Ｓ２１０Ｂのループ処理はスレッドＴＰ２に相当し、ステップＳ２２０Ａ〜Ｓ２２０Ｂのループ処理はスレッドＴＰ３に相当し、ステップＳ２３０Ａ〜Ｓ２３０Ｂのループ処理はスレッドＴＰ４に相当し、ステップＳ２４０Ａ〜Ｓ２４０Ｂのループ処理はスレッドＴＰ５に相当し、ステップＳ２５０Ａ〜Ｓ２５０Ｂのループ処理はスレッドＴＰ６に相当し、 ステップＳ２６０Ａ〜Ｓ２６０Ｂのループ処理はスレッドＴＰ７に相当し、ステップＳ２７０Ａ〜Ｓ２７０Ｂのループ処理はスレッドＴＰ８に相当する。
【００５５】
スレッドＴＰ１〜ＴＰ７では、ｄｘ，ｄｙの値がそれぞれインクリメントにより指定される（ステップＳ２００Ａ，Ｓ２１０Ａ，Ｓ２２０Ａ，Ｓ２３０Ａ，Ｓ２４０Ａ，Ｓ２５０Ａ，Ｓ２６０Ａ，Ｓ２７０Ａ）。ただし、当該ステップの最初の処理では、インクリメントが実行される代わりに、ｄｘ，ｄｙの値として初期値が与えられる。
【００５６】
次いで、ｄｘ，ｄｙの組により、システムマトリクス要素の値Ｃｘ＝Ｃｉｊ（ｄｘ，ｄｙ）が指定され、読み出される（ステップＳ２０１，Ｓ２１１，Ｓ２２１，Ｓ２３１，Ｓ２４１，Ｓ２５１，Ｓ２６１，Ｓ２７１）。
【００５７】
次に、前方投影算出部３３は、システムマトリクス要素の共通の値Ｃｘを用いたベクトル演算により第１の積和演算を実行する（ステップＳ２０２，Ｓ２１２，Ｓ２２２，Ｓ２３２，Ｓ２４２，Ｓ２５２，Ｓ２６２，Ｓ２７２）。上述したように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とは、Ｘ軸方向の並進対称性とＹ軸方向の並進対称性とを共有している。よって、これら並進対称性に基づいて、同一値Ｃｘを共有する複数のシステムマトリクス要素が必ず存在する。前方投影算出部３３は、この共通の値Ｃｘを用いたベクトル演算を実行することができる。
【００５８】
たとえば、図１２に示されるように、Ｘ軸方向またはＹ軸方向への同一の傾きを有する複数の投影線ＰＬ，ＰＬ，…上におけるボクセルＢ１，Ｂ３，Ｂ５，Ｂ７，Ｂ９，Ｂ１１，Ｂ１３，Ｂ１５に対応するシステムマトリクス要素は、Ｘ軸方向またはＹ軸方向の並進対称性を有する。よって、これらボクセルＢ１，Ｂ３，Ｂ５，Ｂ７，Ｂ９，Ｂ１１，Ｂ１３，Ｂ１５に対応するシステムマトリクス要素は共通の値を有している。同様に、複数の投影線ＰＬ，ＰＬ，…上におけるボクセルＢ２，Ｂ４，Ｂ６，Ｂ８，Ｂ１０，Ｂ１２，Ｂ１４，Ｂ１６に対応するシステムマトリクス要素も、共通の値を有している。
【００５９】
たとえば、前方投影値ｐ１，ｐ２，ｐ３，ｐ４が、システムマトリクス要素の共通の値Ｃ１，Ｃ２，…，ＣＮを用いて以下のように計算されるものとする。
【００６０】
ｐ１＝Ｃ１・ｘ１＋Ｃ２・ｘｊ＋…
ｐ２＝Ｃ１・ｘ２＋Ｃ２・ｘｊ＋１＋…
ｐ３＝Ｃ１・ｘ３＋Ｃ２・ｘｊ＋２＋…
ｐ４＝Ｃ１・ｘ４＋Ｃ２・ｘｊ＋３＋…
【００６１】
このとき、前方投影算出部３３は、上式の右辺第１項の計算を、（ｘ１、ｘ２、ｘ３、ｘ４）をベクトルとし、このベクトルに係数Ｃ１を乗算することで行うことができる。続けて、前方投影算出部３３は、上式の右辺第２項の計算を、（ｘｊ、ｘｊ＋１、ｘｊ＋２、ｘｊ＋３）をベクトルとし、このベクトルに係数Ｃ２を乗算することで行うことができる。なお、ベクトルの要素の個数は、４個に限らず、ｎ個（ｎは５以上）であってよい。
【００６２】
スレッドＴＰ１〜ＴＰ８のループ処理は、ｄｘ，ｄｙの値の全ての組み合わせについて実行された後に終了する。その後、図４のステップＳ１３３に処理が移行する。
【００６３】
また、後方投影算出部３４も、上式（３）の第２の積和演算を８つのグループＧ１〜Ｇ８にそれぞれ対応する８つのスレッドＴＢ１〜ＴＢ８に分割し、これらスレッドＴＢ１〜ＴＢ８を並列に実行することができる。
【００６４】
図６は、図３の後方投影算出処理の手順を概略的に示すフローチャートである。図６に示されるように、後方投影算出処理の手順は、８つのスレッドに分割される。ステップＳ３００Ａ〜Ｓ３００Ｂのループ処理はスレッドＴＢ１に相当し、ステップＳ３１０Ａ〜Ｓ３１０Ｂのループ処理はスレッドＴＢ２に相当し、ステップＳ３２０Ａ〜Ｓ３２０Ｂのループ処理はスレッドＴＢ３に相当し、ステップＳ３３０Ａ〜Ｓ３３０Ｂのループ処理はスレッドＴＢ４に相当し、ステップＳ３４０Ａ〜Ｓ３４０Ｂのループ処理はスレッドＴＢ５に相当し、ステップＳ３５０Ａ〜Ｓ３５０Ｂのループ処理はスレッドＴＢ６に相当し、ステップＳ３６０Ａ〜Ｓ３６０Ｂのループ処理はスレッドＴＢ７に相当し、ステップＳ３７０Ａ〜Ｓ３７０Ｂのループ処理はスレッドＴＢ８に相当する。
【００６５】
スレッドＴＢ１〜ＴＢ８では、ｄｘ，ｄｙの値がそれぞれインクリメントにより指定される（ステップＳ３００Ａ，Ｓ３１０Ａ，Ｓ３２０Ａ，Ｓ３３０Ａ，Ｓ３４０Ａ，Ｓ３５０Ａ，Ｓ３６０Ａ，Ｓ３７０Ａ）。ただし、当該ステップの最初の処理では、インクリメントが実行される代わりに、ｄｘ，ｄｙの値として初期値が与えられる。
【００６６】
次いで、ｄｘ，ｄｙの組により、システムマトリクス要素の値Ｃｘ＝Ｃｉｊ（ｄｘ，ｄｙ）が読み出される（ステップＳ３０１，Ｓ３１１，Ｓ３２１，Ｓ３３１，Ｓ３４１，Ｓ３５１，Ｓ３６１，Ｓ３７１）。
【００６７】
次に、後方投影算出部３４は、システムマトリクス要素の共通の値Ｃｘを用いたベクトル演算により第１の積和演算を実行する（ステップＳ３０２，Ｓ３１２，Ｓ３２２，Ｓ３３２，Ｓ３４２，Ｓ３５２，Ｓ３６２，Ｓ３７２）。前述したように、投影線ＰＬの配列とシステムマトリクス要素｛Ｃｉｊ｝とが共有する並進対称性に基づいて、システムマトリクス要素の値Ｃｘを共有する複数のシステムマトリクス要素が必ず存在する。後方投影算出部３４は、この共通の値Ｃｘを用いたベクトル演算を実行することができる。
【００６８】
たとえば、後方投影値Ｄ１，Ｄ２，Ｄ３，Ｄ４が、システムマトリクス要素の共通の値Ｃ１，Ｃ２，…，ＣＮを用いて以下のように計算されるものとする。
【００６９】
Ｄ１＝Ｃ１・ｑ１＋Ｃ２・ｑｉ＋…
Ｄ２＝Ｃ１・ｑ２＋Ｃ２・ｑｉ＋１＋…
Ｄ３＝Ｃ１・ｑ３＋Ｃ２・ｑｉ＋２＋…
Ｄ４＝Ｃ１・ｑ４＋Ｃ２・ｑｉ＋３＋…
【００７０】
このとき、後方投影算出部３４は、上式の右辺第１項の計算を、（ｑ１、ｑ２、ｑ３、ｑ４）をベクトルとし、このベクトルに係数Ｃ１を乗算することで行うことができる。続けて、後方投影算出部３４は、上式の右辺第２項の計算を、（ｑｉ、ｑｉ＋１、ｑｉ＋２、ｑｉ＋３）をベクトルとし、このベクトルに係数Ｃ２を乗算することで行うことができる。なお、ベクトルの要素の個数は、４個に限らず、ｍ個（ｍは５以上）であってもよい。
【００７１】
スレッドＴＢ１〜ＴＢ８のループ処理は、ｄｘ，ｄｙの値の全ての組み合わせについて実行された後に終了する。その後、図４のステップＳ１３４に処理が移行する。
【００７２】
上述した通り、画像再構成システム１は、投影線ＰＬを、投影線の配列とシステムマトリクス要素の値とが共有する複数種の幾何学的対称性に基づいて複数のグループに分類する。前方投影算出部３３は、第１の積和演算をこれら複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することができる。これにより、前方投影値を短時間で算出することができるので、各ボクセルの濃度を短時間で算出することが可能となる。同様に、後方投影算出部３４は、第２の積和演算を複数のグループにそれぞれ対応する複数のスレッドに分割し、当該複数のスレッドを並列に実行することができる。これにより、後方投影値を短時間で算出することができるので、各ボクセルの濃度を短時間で算出することが可能となる。
【００７３】
さらに、投影線ＰＬの配列とシステムマトリクス要素の値とは同一方向の並進対称性を共有するので、この並進対称性に基づいた同一値を有する複数のシステムマトリクス要素が存在する。前方投影算出部３３と後方投影算出部３４は、このシステムマトリクス要素を用いたベクトル演算を実行して前方投影値と後方投影値を短時間で算出することができる。
【００７４】
ところで、図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。ここで、プロセッサは、上記ベクトル演算を実行するためのベクトル演算機能を有する。
【００７５】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【００７６】
以上、図面を参照して本発明の実施形態について述べたが、これらは本発明の例示であり、上記以外の様々な構成を採用することもできる。
【００７７】
たとえば、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１の対向型ＰＥＴ装置２に代えて、連続的に分布する単一の検出面あるいは３個以上の検出面を有する検出装置を使用してもよい。また、検出面の形状は平面に限定されるものではなく曲面であってもよい。
【００７８】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置２に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
【００７９】
上記実施形態では、統計的画像再構成法としてＭＬ−ＥＭ法が使用されているが、システムマトリクス要素（検出確率）を用いた画像再構成法を使用すればよく、ＭＬ−ＥＭ法に限定されるものではない。たとえば、ＯＳ−ＥＭ法やＭＡＰ−ＥＭ（Maximum A Posteriori - Expectation Maximization）法などのＭＬ−ＥＭ法の改良版を使用することが可能である。
【符号の説明】
【００８０】
１ 画像再構成システム
２ 対向型ＰＥＴ装置
３ 画像処理装置
１０Ａ，１０Ｂ 検出プレート
１１Ａ，１１Ｂ 検出器要素
２０ 被測定物
３１ 投影データ生成部
３２ 濃度推定部
３３ 前方投影算出部
３４ 後方投影算出部
３５ 濃度算出部
３６ 表示制御部
３７ メモリ
３８ データ記憶部
３９ パラメータ格納部
４０ 制御部
４１ 駆動機構
５０ ディスプレイ
５１ 入力操作部
【特許請求の範囲】
【請求項１】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理装置であって、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、
画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、
を備え、
前記濃度推定部は、
前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する、画像処理装置。
【請求項２】
請求項１に記載の画像処理装置であって、前記複数種の幾何学的対称性は、反転対称性と回転対称性とを含む、画像処理装置。
【請求項３】
請求項１または２に記載の画像処理装置であって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を共有しており、
前記前方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値を、当該システムマトリクス要素を用いたベクトル演算により算出する、画像処理装置。
【請求項４】
請求項１から３のうちのいずれか１項に記載の画像処理装置であって、
前記濃度推定部は、
前記実測投影値と当該実測投影値に対応する当該前方投影値の逆数との積を算出するとともに前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該積との第２の積和演算を実行して前記投影線毎に後方投影値を算出する後方投影算出部、
をさらに含み、
前記後方投影算出部は、前記第２の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行し、
前記濃度算出部は、前記後方投影値を用いて当該新たな濃度を算出する、
画像処理装置。
【請求項５】
請求項４に記載の画像処理装置であって、
前記投影線の配列と前記システムマトリクス要素の値とは、同一方向の並進対称性を有しており、
前記後方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｍ個（ｍは２以上の整数）のシステムマトリクス要素に対応するｍ個の後方投影値を、当該システムマトリクス要素を用いたベクトル演算により実行する、画像処理装置。
【請求項６】
請求項１から５のうちのいずれか１項に記載の画像処理装置であって、前記複数の検出器要素対は、前記被測定物を挟み込む一対の検出プレートを構成する、画像処理装置。
【請求項７】
請求項１から６のうちのいずれか１項に記載の画像処理装置であって、前記各検出器要素は、当該各検出器要素への入射光に応じた蛍光出力を前記検出結果として与えるシンチレータを含む、画像処理装置。
【請求項８】
請求項１から７のうちのいずれか１項に記載の画像処理装置であって、前記画像再構成法は、ＭＬ−ＥＬ（Maximum Likelihood - Expectation Maximization）法である、画像処理装置。
【請求項９】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成する画像再構成システムであって、
前記領域を挟み込む複数の検出器要素対を有する検出装置と、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、
画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、
を備え、
前記濃度推定部は、
前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する、画像再構成システム。
【請求項１０】
請求項９に記載の画像再構成システムであって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を共有しており、
前記前方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値を、当該システムマトリクス要素を用いたベクトル演算により算出する、画像再構成システム。
【請求項１１】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理方法であって、
（ａ）前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出するステップと、
（ｂ）前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行するステップと、
を備え、
前記ステップ（ｂ）は、
（ｂ−１）前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出するステップと、
（ｂ−２）前記前方投影値および前記実測投影値を用いて新たな濃度を算出するステップと、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記ステップ（ｂ−１）では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される、画像処理方法。
【請求項１２】
コンピュータ読み取り可能な記録媒体から読み出されてプロセッサに画像再構成処理を実行させるプログラムであって、
前記画像再構成処理は、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成処理と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定処理と、
を備え、
前記濃度推定処理は、
前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出処理と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出処理と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出処理では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される、プログラム。
【請求項１３】
請求項１２に記載のプログラムであって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を有しており、
前記前方投影算出処理では、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値が、当該システムマトリクス要素を用いたベクトル演算により算出される、プログラム。
【請求項１４】
請求項１３に記載のプログラムであって、前記プロセッサはベクトル演算機能を有する、プログラム。
【請求項１】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理装置であって、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、
画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、
を備え、
前記濃度推定部は、
前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する、画像処理装置。
【請求項２】
請求項１に記載の画像処理装置であって、前記複数種の幾何学的対称性は、反転対称性と回転対称性とを含む、画像処理装置。
【請求項３】
請求項１または２に記載の画像処理装置であって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を共有しており、
前記前方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値を、当該システムマトリクス要素を用いたベクトル演算により算出する、画像処理装置。
【請求項４】
請求項１から３のうちのいずれか１項に記載の画像処理装置であって、
前記濃度推定部は、
前記実測投影値と当該実測投影値に対応する当該前方投影値の逆数との積を算出するとともに前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該積との第２の積和演算を実行して前記投影線毎に後方投影値を算出する後方投影算出部、
をさらに含み、
前記後方投影算出部は、前記第２の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行し、
前記濃度算出部は、前記後方投影値を用いて当該新たな濃度を算出する、
画像処理装置。
【請求項５】
請求項４に記載の画像処理装置であって、
前記投影線の配列と前記システムマトリクス要素の値とは、同一方向の並進対称性を有しており、
前記後方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｍ個（ｍは２以上の整数）のシステムマトリクス要素に対応するｍ個の後方投影値を、当該システムマトリクス要素を用いたベクトル演算により実行する、画像処理装置。
【請求項６】
請求項１から５のうちのいずれか１項に記載の画像処理装置であって、前記複数の検出器要素対は、前記被測定物を挟み込む一対の検出プレートを構成する、画像処理装置。
【請求項７】
請求項１から６のうちのいずれか１項に記載の画像処理装置であって、前記各検出器要素は、当該各検出器要素への入射光に応じた蛍光出力を前記検出結果として与えるシンチレータを含む、画像処理装置。
【請求項８】
請求項１から７のうちのいずれか１項に記載の画像処理装置であって、前記画像再構成法は、ＭＬ−ＥＬ（Maximum Likelihood - Expectation Maximization）法である、画像処理装置。
【請求項９】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用して前記被測定物の内部構造を表す画像を生成する画像再構成システムであって、
前記領域を挟み込む複数の検出器要素対を有する検出装置と、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成部と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率がシステムマトリクス要素として記憶されているパラメータ格納部と、
画像再構成法に従って、前記システムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定部と、
を備え、
前記濃度推定部は、
前記パラメータ格納部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出部と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出部と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出部は、前記第１の積和演算を前記複数のグループにそれぞれ対応する複数のスレッドに分割するとともに当該複数のスレッドを並列に実行する、画像再構成システム。
【請求項１０】
請求項９に記載の画像再構成システムであって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を共有しており、
前記前方投影算出部は、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値を、当該システムマトリクス要素を用いたベクトル演算により算出する、画像再構成システム。
【請求項１１】
被測定物が配置されるべき領域内に複数のボクセルを立体的に配列した解析モデルを使用し、前記領域を挟み込む複数の検出器要素対による光子対の検出結果に基づいて前記被測定物の内部構造を表す画像を生成する画像処理方法であって、
（ａ）前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出するステップと、
（ｂ）前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行するステップと、
を備え、
前記ステップ（ｂ）は、
（ｂ−１）前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出するステップと、
（ｂ−２）前記前方投影値および前記実測投影値を用いて新たな濃度を算出するステップと、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記ステップ（ｂ−１）では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される、画像処理方法。
【請求項１２】
コンピュータ読み取り可能な記録媒体から読み出されてプロセッサに画像再構成処理を実行させるプログラムであって、
前記画像再構成処理は、
前記検出器要素対の一方から他方を結ぶ投影線上における前記ボクセルから放出され当該検出器要素対で検出された光子対の計数結果に基づいて、前記投影線毎に実測投影値を算出する投影データ生成処理と、
前記各ボクセルから放出された光子を前記複数の検出器要素対がそれぞれ検出する確率であるシステムマトリクス要素をデータ記憶部から読み出し、画像再構成法に従って、当該読み出されたシステムマトリクス要素と前記実測投影値とを使用して前記ボクセル毎の濃度を算出する処理を反復的に実行する濃度推定処理と、
を備え、
前記濃度推定処理は、
前記データ記憶部から読み出された当該システムマトリクス要素と当該システムマトリクス要素に対応する当該濃度との第１の積和演算を実行して前記投影線毎に前方投影値を算出する前方投影算出処理と、
前記前方投影値および前記実測投影値を用いて新たな濃度を算出する濃度算出処理と、
を含み、
前記投影線の配列と前記システムマトリクス要素の値とは、複数種の幾何学的対称性を共有し、前記投影線は、前記複数種の幾何学的対称性に基づいた複数のグループに分類されており、
前記前方投影算出処理では、前記第１の積和演算は前記複数のグループにそれぞれ対応する複数のスレッドに分割され、前記複数のスレッドが並列に実行される、プログラム。
【請求項１３】
請求項１２に記載のプログラムであって、
前記投影線の配列と前記システムマトリクス要素の値とは、さらに同一方向の並進対称性を有しており、
前記前方投影算出処理では、前記複数のシステムマトリクス要素のうち前記並進対称性に基づいた同一値を有するｎ個（ｎは２以上の整数）のシステムマトリクス要素に対応するｎ個の前方投影値が、当該システムマトリクス要素を用いたベクトル演算により算出される、プログラム。
【請求項１４】
請求項１３に記載のプログラムであって、前記プロセッサはベクトル演算機能を有する、プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【公開番号】特開２０１０−２０３８０６（Ｐ２０１０−２０３８０６Ａ）
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願番号】特願２００９−４７１７５（Ｐ２００９−４７１７５）
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
[ Back to top ]