画像処理装置、画像再構成システム、画像処理方法およびプログラム
【課題】画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置を提供する。
【解決手段】画像処理装置３において、反復演算部３２は、パラメータ格納部３９から読み出された検出確率Ｃｉｊとｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理と、検出確率Ｃｉｊとｉ番目の検出器要素にかかる投影データとの第二積和演算処理と、を含む演算処理を行う。そして、反復演算部３２は、投影データが所定値以下であるシンチレーション検出器１１に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【解決手段】画像処理装置３において、反復演算部３２は、パラメータ格納部３９から読み出された検出確率Ｃｉｊとｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理と、検出確率Ｃｉｊとｉ番目の検出器要素にかかる投影データとの第二積和演算処理と、を含む演算処理を行う。そして、反復演算部３２は、投影データが所定値以下であるシンチレーション検出器１１に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被測定物から放出され、または被測定物を透過した光の量に基づいて、被測定物の内部構造の三次元画像を再構成するための画像処理技術に関する。
【背景技術】
【０００２】
人体などの被測定物を透過したＸ線、あるいは被測定物から放出されたγ線といった放射線を検出し、その検出結果に基づいて被測定物の内部構造の二次元画像または三次元画像を再構成するための画像処理技術の研究開発が進められている。このような画像処理技術は、Ｘ線ＣＴ装置（Ｘ線コンピュータ断層撮影装置）や陽電子放射断層撮像装置（ＰＥＴ：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.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、たとえ従来の高速演算アルゴリズムを用いたとしても、計算量が膨大であるので、その計算速度は臨床上実用的であるとは言い難いという問題がある。たとえば、画像を再構成するまでの計算時間が長大であると、診断に要する時間が長くなるので、臨床上実用的であるとは言い難い。
【０００７】
本発明は上記課題に鑑みてなされたものであり、画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置、画像再構成システム、画像処理方法およびプログラムを提供するものである。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データと、前記被測定物が配置された検出領域の内部に複数のボクセルが立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて前記投影データより前記ボクセルごとのボクセル値を推定する演算処理を反復的に実行して前記被測定物の内部構造の三次元画像を再構成する反復演算部と、
ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊを示すデータを格納しているパラメータ格納部と、
を備える画像処理装置が提供される。
この画像処理装置は、前記反復演算部が、
（Ａ）前記パラメータ格納部から読み出された検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理と、
前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理と、
を含む前記演算処理を行うとともに、
（Ｂ）前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする。
【０００９】
また、本発明の画像処理装置においては、より具体的な態様として、前記反復演算部が、前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルのボクセル値をゼロに設定してもよい。
【００１０】
また、本発明の画像処理装置においては、より具体的な態様として、前記投影データが前記所定値以下である検出器要素に対応するボクセルの集合であるゼロ領域と、前記投影データが前記所定値よりも大きい前記検出器要素に対応するボクセルの集合である非ゼロ領域と、を含む、ボクセル値の二次元画像を生成する二次元画像生成部と、
前記二次元画像を内包する三次元初期画像を生成し、前記三次元初期画像における前記ゼロ領域に対応するボクセルを決定する三次元初期画像生成部と、
をさらに備え、
前記反復演算部が、前記三次元初期画像における前記ゼロ領域に対応するボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップしてもよい。
【００１１】
また、本発明の画像処理装置においては、より具体的な態様として、前記ボクセルのサイズを複数段階に微細化して前記演算処理を反復的に実行し、
微細化前に前記三次元初期画像における前記ゼロ領域に対応していたボクセルに内包される微細化後のボクセルについて、前記反復演算部が、前記第一または第二積和演算処理の少なくとも一方の実行をスキップしてもよい。
【００１２】
また、本発明によれば、対向して配置されたシンチレーション検出器の対を前記検出器要素として含む上記記載の画像処理装置を備えるとともに、前記シンチレーション検出器が、前記被測定物から所定方向の両側に放射された光子対を検出する画像再構成システムが提供される。
【００１３】
また、本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像処理方法が提供される。
この画像処理方法は、（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影ステップと、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前記前方投影ステップおよび後方投影ステップを反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算ステップと、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００１４】
また、本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像再構成処理を情報処理装置に実行させるプログラムが提供される。
このプログラムでは、前記画像再構成処理が、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成処理と、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影処理と、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影処理と、
（ｄ）前記前方投影処理および後方投影処理を反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算処理と、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００１５】
上記発明において、被測定物から検出器要素に光子が飛来するとは、光子と検出器要素との相互作用により光子が検出されることを意味する。また、被測定物から光子が飛来するとは、被測定物の後方から照射された光子が被測定物を透過して検出器要素に飛来する場合と、被測定物を光源として光子が放射される場合とを含む。以下、被測定物から光子が放出される表現した場合、上記の両態様を含む意味である。
【００１６】
なお、本発明の各種の構成要素は、その機能を実現するように形成されていればよく、たとえば、所定の機能を発揮する専用のハードウェア、所定の機能がコンピュータプログラムにより付与された画像処理装置、コンピュータプログラムにより画像処理装置に実現された所定の機能、これらの任意の組み合わせ、等として実現することができる。
また、本発明の各種の構成要素は、個々に独立した存在である必要はなく、複数の構成要素が一個の部材として形成されていること、一つの構成要素が複数の部材で形成されていること、ある構成要素が他の構成要素の一部であること、ある構成要素の一部と他の構成要素の一部とが重複していること、等を許容する。
【００１７】
また、本発明の画像処理方法は、複数の工程を順番に記載してあるが、その記載の順番は複数の工程を実行する順番を限定するものではない。このため、本発明の画像処理方法を実施するときには、その複数の工程の順番は内容的に支障しない範囲で変更することができる。
さらに、本発明の画像処理方法は、複数の工程が個々に相違するタイミングで実行されることに限定されない。このため、ある工程の実行中に他の工程が発生すること、ある工程の実行タイミングと他の工程の実行タイミングとの一部ないし全部が重複していること、等でもよい。
【００１８】
また、本発明で云う反復的な演算処理は、コンピュータプログラムを読み取って対応するデータ処理を実行できるように、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、Ｉ／Ｆ（Interface）ユニット、等の汎用デバイスで構築されたハードウェア、所定のデータ処理を実行するように構築された専用の論理回路、これらの組み合わせ、等を用いて実施することができる。
なお、本発明でプログラムに対応した各種動作を情報処理装置に実行させることは、各種デバイスを情報処理装置に動作制御させることを含む意味である。
【００１９】
さらに、本発明において、データを記憶または格納しているとは、本発明の装置が、少なくとも演算処理を実行するときに、当該データを保持している状態となる機能を有することを意味している。したがって、たとえば、光量を示す投影データや検出確率Ｃｉｊを示すデータに関しては、本発明の装置が出荷されるときに既に登録されていることの他、出荷されるときには登録されていないこれらのデータが、演算処理の時点で、または演算処理の最中に登録されることも許容する。
【発明の効果】
【００２０】
本発明による画像処理技術によれば、飛来した光子の光量を示す投影データが所定値以下である検出器要素に対応するボクセルについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。これは、光子の飛来しない領域に存在するボクセルに関してはボクセル値を一律に所定値（たとえばゼロ）とみなすことが可能であることによる。そして、第一および第二積和演算処理は統計的画像再構成法による反復的な演算処理における計算量の多い処理である。したがって、上記のボクセルに関して、かかる積和演算処理をスキップすることにより、演算処理の総時間を大幅に短縮することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明に係る一実施形態の画像再構成システムの概略構成を示す機能ブロック図である。
【図２】実測投影データを説明するための模式図である。
【図３】ＭＬ−ＥＭ法による処理手順を概略的に示すフローチャートである。
【図４】二次元画像の作成処理の具体的な手順を示すフローチャートである。
【図５】三次元初期画像の作成出力の具体的な手順を示すフローチャートである。
【図６】画像再構成処理の具体的な手順を示すフローチャートである。
【図７】（ａ）は解析モデルおよび検出領域、（ｂ）は解析モデルおよび被測定物、（ｃ）は二次元画像の一例を示す斜視図である。
【図８】三次元初期画像を示す斜視図である。
【図９】（ａ）は概略再構成処理に用いる概略モデル、（ｂ）は精細再構成処理に用いる精細モデルの斜視図である。
【発明を実施するための形態】
【００２２】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
図１は、本発明に係る一実施形態の画像再構成システム１および画像処理装置３の概略構成を示す機能ブロック図である。
【００２３】
本実施形態の画像再構成システム１は、画像処理装置３と検出部２とを備えている。
画像処理装置３は、反復演算部３２とパラメータ格納部３９とを備えている。
反復演算部３２は、被測定物２０から複数の検出器要素（結晶対１４：図２を参照）の各々に飛来した光子の光量を示す投影データと、被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘ（図７を参照）が立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて投影データよりボクセルＢｘごとのボクセル値ｘを推定する演算処理を反復的に実行して被測定物２０の内部構造の三次元画像を再構成する。
また、パラメータ格納部３９は、ｊ番目のボクセルＢｘｊを通過した光子がｉ番目のシンチレーション検出器１１ｉで検出される検出確率Ｃｉｊを示すデータを格納している。
【００２４】
そして、本実施形態の反復演算部３２は、パラメータ格納部３９から読み出された検出確率Ｃｉｊと、ｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理と、検出確率Ｃｉｊと、ｉ番目の検出器要素にかかる投影データとの第二積和演算処理と、を含む演算処理を行う。
また、本実施形態の反復演算部３２は、投影データが所定値以下であるシンチレーション検出器１１に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００２５】
画像再構成システム１は、互いに対向する一対の検出プレート１０（１０Ａ，１０Ｂ）を有する対向型ＰＥＴ装置を構成し、被測定物２０から放射された放射線の線量を測定して、被測定物２０の外形とともにその内部構造の三次元画像を推定（再構成）する。
【００２６】
検出部２は、互いに対向する検出面１２（１２Ａ、１２Ｂ）を有する一対の検出プレート１０Ａ，１０Ｂと、これら検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１と、を有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００２７】
検出プレート１０には、多数のシンチレーション検出器１１（１１Ａ，１１Ｂ）が二次元平面内に配列されている。
シンチレーション検出器１１の配列の態様は特に限定されず、格子状、千鳥状、またはそれ以外でもよい。本実施形態では、シンチレーション検出器１１がＸ−Ｙ平面に沿って格子状（マトリクス状）に配列された検出プレート１０を例示する。
各方向のシンチレーション検出器１１の配列数は任意であるが、本実施形態では、Ｘ，Ｙ軸方向に沿って、それぞれｍ個ずつのシンチレーション検出器１１が配列されて検出プレート１０を構成しているものとする。
【００２８】
シンチレーション検出器１１は、シンチレータ結晶（図示せず）と、シンチレータ結晶の内部で放射線が消失することにより発生する蛍光を受光する受光素子（図示せず）とを有している。
検出プレート１０Ａ，１０Ｂの対向する検出面１２は、シンチレータ結晶の端面の集合により構成されている。各シンチレータ結晶は面一に配置され、検出面１２は平坦に形成されている。
本実施形態のシンチレータ結晶は、検出面１２を構成する端面が正方形をなす四角柱をなしている。ただし、シンチレータ結晶の端面の形状および寸法はこれに限られるものではない。
そして、本実施形態は、対向するシンチレーション検出器１１Ａと１１Ｂとの組み合わせごとに、被測定物２０からの放射線の飛来を検出する。
【００２９】
すなわち、本実施形態の画像再構成システム１は、対向して配置されたシンチレーション検出器１１の対（結晶対１４）を、検出器要素として備えている。そして、シンチレーション検出器１１では、被測定物２０から所定方向の両側に放射された光子対を検出する。
【００３０】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に従って、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。
したがって、本実施形態の画像再構成システム１は、駆動機構４１の駆動により、一対の検出面１２Ａ、１２Ｂの間隔が可変である。したがって、検出プレート１０Ａ，１０Ｂに挟み込まれた被測定物２０は検出面１２Ａ、１２Ｂに押圧されて拡幅する。画像再構成システム１では、かかる押圧状態の被測定物２０から飛来する放射線を受光する。
ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を通じて制御部４０に指示を出すことにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００３１】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対のシンチレーション検出器１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対のシンチレーション検出器１１Ａ，１１Ｂを結ぶ線上に、陽電子が消滅した地点が存在する。
【００３２】
シンチレーション検出器１１Ａ，１１Ｂに用いられるシンチレータ結晶は、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯやＬＳＯなどの公知の結晶体を使用すればよい。
検出プレート１０Ａ，１０Ｂが備える受光素子は、シンチレーション検出器１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、ケーブルなどの転送手段（図示せず）を介して、光子を検出したシンチレーション検出器１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００３３】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、結晶対１４ごとに検出した光子を計数し、この計数結果に基づいて、結晶対１４（検出器要素）毎の投影データ（実測投影データ）を生成する。この実測投影データは、メモリ３７内のデータ記憶部３８に格納される。
【００３４】
図２は、実測投影データｙを説明するための模式図である。被測定物２０内の投影線上の陽電子消滅地点（線源）２１から所定方向の両側に放射された光子は、検出プレート１０Ａ，１０Ｂのシンチレーション検出器１１Ａ，１１Ｂでそれぞれほぼ同時に検出される。
線源２１は被測定物２０のほぼ全体に分布しており、光子は被測定物２０のほぼ全体から、さまざまな方向に対して放射される。そして、癌細胞のように活性の高い部位に線源２１は高密度で集中しているため、癌細胞を挟む結晶対１４では光子が多く検出される。
図２に示すように、特定のシンチレーション検出器１１Ｂに着目した場合、検出プレート１０Ａを構成する各シンチレーション検出器１１Ａのうち、線源２１を挟むシンチレーション検出器１１Ａ１と、線源２１を挟まないシンチレーション検出器１１Ａ２とが存在する場合がある。
他のシンチレーション検出器１１Ｂに関しても同様である。
【００３５】
投影データ生成部３１（図１を参照）は、対向するシンチレーション検出器１１Ａ，１１Ｂで所定の時間差内に光子を検出したことを検知すると、かかるシンチレーション検出器１１Ａ，１１Ｂの組み合わせに相当する結晶対１４ｉ（ｉは結晶対の番号）について、当該光子の光量を積算してカウントする。かかる計測を所定時間に亘っておこなうことにより、結晶対の番号ごとの実測投影データｙ（実測投影データｙｉ）を求めることができる。
そして、積算された実測投影データｙｉの分布に基づいて、統計的画像再構成法を用いることにより線源２１の分布密度が再構成され、もって癌細胞２３の有無およびその位置が推定される。
【００３６】
図１に示す反復演算部３２は、被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘ（立方形状または長方形状などの微小な立体要素、図７を参照）が立体的に配列された解析モデル２５を利用して画像再構成処理を実行する。
反復演算部３２は、データ記憶部３８から実測投影データｙｉを読み出し、統計的画像再構成法に従って、この実測投影データｙｉを用いてボクセルＢｘ毎の放射能濃度に相当するボクセル値ｘを推定する演算処理を反復的に実行する。
これにより、被測定物２０の内部構造の断層画像または三次元画像を再構成する機能を有している。本実施形態では、ボクセル値の三次元画像に基づいて、二次元画像である断層画像を出力することを含めて、三次元画像の再構成という。
【００３７】
本実施形態では、統計的画像再構成法として、典型的なＭＬ−ＥＭ法を使用することができる。ただし、使用可能な統計的画像再構成法はこれに限定されるものではない。
【００３８】
反復演算部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これらの機能ブロックは、ＭＬ−ＥＭ法に従って処理を反復的に実行する。ＭＬ−ＥＭ法では、以下の推定式（１）が使用される。
【００３９】
【数１】
【００４０】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００４１】
【数２】
【００４２】
上記推定式（１）中、ｋは反復回数であり、ｘｊ（ｋ）はｋ回目の処理におけるｊ番目のボクセルＢｘｊの放射能濃度（ボクセル値ｘ）である。ｙｉはｉ番目の結晶対（結晶対１４ｉ）による検出結果に基づいて得られた、上述の実測投影データの値（投影値）である。
本実施形態の画像再構成システム１の場合、検出プレート１０は、シンチレーション検出器１１がＸ，Ｙ軸方向に沿ってそれぞれｍ個ずつ配列されて構成されている。したがって、検出プレート１０はそれぞれｍ２個のシンチレーション検出器１１を有している。そして、シンチレーション検出器１１の組み合わせである結晶対１４はｍ４組だけ存在し、ｉの最大値はｍ４となる。
【００４３】
また、ｘｊ（ｋ＋１）は、（ｋ＋１）回目の処理におけるｊ番目のボクセルＢｘｊにおける放射能濃度である。よって、ｋ回目のボクセル値ｘｊ（ｋ）から、（ｋ＋１）回目のボクセル値ｘｊ（ｋ＋１）が算出される。
ここで、ｘｊ（ｋ）の初期値ｘｊ（０）は、正の実数であるかぎりどのような値に設定されても、反復回数ｋを十分に大きく設定して演算処理を反復的に行うことで、所定の確度の三次元画像を再構成することができる。
【００４４】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｘｊから飛来した光子対がｉ番目の結晶対１４ｉで検出される確率（検出確率）である。ボクセルＢｘｊから飛来した光子対が結晶対１４ｉに入射する数は、検出確率Ｃｉｊに依存することが知られている。検出確率ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００４５】
検出確率Ｃｉｊは、行列（「システムマトリクス」と呼ばれる。）のｉ行ｊ列目の要素とみなすことができる。検出確率Ｃｉｊの値には、散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、検出確率Ｃｉｊの値は、ｉ番目の結晶対１４ｉに対するｊ番目のボクセルＢｘｊの幾何学的位置のみにより定めることができる。
検出確率Ｃｉｊの値には、レイトレース法（Ray-tracing method）や、モンテカルロシミュレーション法（Monte Carlo simulation method）等によって求まる値を設定することができる。
【００４６】
次に、画像処理装置３の動作として実現される本実施形態の画像処理方法（以下、本方法という）を説明する。
【００４７】
まず、本方法の概要について説明する。
本方法は、被測定物２０から複数の検出器要素（結晶対１４）の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて被測定物２０の内部構造の三次元画像を再構成する方法である。
この本方法は、
（ａ）被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルＢｘｊを通過した光子がｉ番目の結晶対１４ｉで検出される検出確率Ｃｉｊと、ｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理をおこなう前方投影ステップと、
（ｃ）検出確率Ｃｉｊと、ｉ番目の結晶対１４ｉにかかる投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前方投影ステップおよび後方投影ステップを反復的に実行してボクセルＢｘごとのボクセル値ｘを推定する反復演算ステップと、
を含む。
そして、本方法は、投影データが所定値以下である結晶対１４に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００４８】
また、本方法は、（ｅ）被測定物２０を押圧して拡幅させた状態で光子の光量を測定する実測ステップをさらに含む。
【００４９】
つぎに、図面を参照して本方法をより詳細に説明する。
図３は、ＭＬ−ＥＭ法による処理手順を概略的に示すフローチャートである。
図４は、二次元画像の作成処理（ステップＳ１２）の具体的な手順を示すフローチャートである。
図５は、三次元初期画像の作成出力（ステップＳ１３）の具体的な手順を示すフローチャートである。
図６は、画像再構成処理（ステップＳ１４およびＳ１５）の具体的な手順を示すフローチャートである。
【００５０】
ステップＳ１０では、反復演算部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。
続けて、反復演算部３２は、データ記憶部３８から実測投影データ｛ｙｉ｝を読み出す（ステップＳ１１）。なお、本実施形態では、検出確率｛Ｃｉｊ｝、係数｛Ｃｊ｝および実測投影データ｛ｙｉ｝が最初に一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００５１】
続けて、反復演算部３２は、放射能濃度に相当するボクセル値｛ｘｊ｝の値を初期化するため、実測投影データｙｉに基づいて作成される二次元画像、および二次元画像に基づいて作成される三次元初期画像を作成する（ステップＳ１３、Ｓ１４）。
反復演算部３２は、その後、ボクセルサイズの大きな概略モデルを用いた画像再構成処理（ステップＳ１５）と、ボクセルサイズを微細化した精細モデルを用いた画像再構成処理（ステップＳ１６）とを実行する。
なお、本方法では、ボクセルサイズを二段階に微細化して画像再構成処理をおこなう態様を例示的に説明するが、微細化の段階はこれに限られず三段階以上としてもよい。
また、ボクセルサイズの微細化は、本発明においては任意であり、単一の解析モデルのみを利用して画像再構成処理をおこなってもよい。
【００５２】
まず、画像再構成処理（ステップＳ１５、Ｓ１６）について、その具体的な手段を、図６を参照して説明する。かかる画像再構成処理に好適に用いられる初期画像として本実施形態で作成される三次元初期画像については後述する。
【００５３】
図６を参照すると、先ず、反復演算部３２は、反復回数の上限を適切な値に設定する（ステップＳ１４０）。
概略モデルによる画像再構成処理（以下、「概略再構成処理」という場合がある。）で設定される反復回数ｋ１と、精細モデルによる画像再構成処理（以下、「精細再構成処理」という場合がある。）で設定される反復回数ｋ２とは、互いに同一でも相違してもよい。
また、反復回数ｋ１とｋ２の値、およびその大小関係も任意であるが、演算処理の総所要時間を短縮して同レベルの精度の画像再構成をおこなう観点からは、ｋ１をｋ２よりも大きく設定するとよい。
【００５４】
図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｊと、推定濃度に相当するボクセル値ｘｊ（ｋ）との積和演算（第一積和演算）を実行して前方投影値ｐｉを算出する（ステップＳ１４１）。
【００５５】
【数３】
【００５６】
推定式（１）中で演算量の多大な部分が、式（２）により与えられる前方投影データの値（前方投影値）ｐｉの計算である。
前方投影算出部３３は、上式（２）の第一積和演算を、複数のスレッド（グループ）に分解して並列演算により実行してもよい。
【００５７】
上記ステップＳ１４１の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、実測投影データｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算（第二積和演算）を実行して後方投影値Ｄｊを算出する（ステップＳ１４２）。
【００５８】
【数４】
【００５９】
この後方投影値Ｄｊの計算も、前方投影値の算出と同様に非常に演算量の多い部分である。
【００６０】
そして、濃度算出部３５は、後方投影値Ｄｊを用いて、（ｋ＋１）回目の推定濃度に相当するボクセル値ｘｊ（ｋ＋１）（＝Ｄｊ×ｘｊ（ｋ））を算出する（ステップＳ１４３）。その後、反復回数ｋ（ｋ１、ｋ２）が上限に達したか否かを判定し（ステップＳ１４４）、反復回数ｋが上限に達するまで、上記ステップＳ１４１〜Ｓ１４４の手順を繰り返し実行する。最終的に反復回数ｋが上限に達したとき（ステップＳ１４４のＹＥＳ）は、図３のフローチャートに戻って、概略再構成処理（ステップＳ１４）または精細再構成処理（ステップＳ１５）を終了する。
【００６１】
そして、濃度算出部３５は、ボクセル値｛ｘｊ｝を出力する（ステップＳ１６）。
【００６２】
表示制御部３６は、濃度算出部３５から出力されたボクセル値｛ｘｊ｝から三次元画像を生成する。生成画像は、二次元断層画像または三次元画像としてディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力されたボクセル値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００６３】
＜三次元初期画像作成処理＞
ここで、従来の統計的画像再構成法では、ボクセル値の初期値として、すべてのボクセルに対して正の値を設定し、ボクセル値ｘｊ（ｋ）の反復演算をすべてのボクセルＢｘｊに関しておこなっていた。
しかし、上述のように、（ｋ＋１）回目の推定濃度に相当するボクセル値ｘｊ（ｋ＋１）は、ｋ回目の推定濃度に相当するボクセル値ｘｊ（ｋ）に後方投影値Ｄｊを乗じて算出される。したがって、ひとたびｘｊ（ｋ）がゼロになったボクセルでは、かりに以降の反復演算をおこなったとしてもボクセル値はゼロが繰り返される。
【００６４】
このため、本方法においては、反復演算部３２は、実測投影データｙｉが所定値以下である結晶対１４に対応するボクセルＢｘのボクセル値ｘｊをゼロに設定する。そして、ボクセル値ｘｊがゼロに設定されたボクセルＢｘｊに関しては、前方投影ステップ（ステップＳ１４１）における前方投影算出部３３による第一積和演算処理（上式（２））、または後方投影ステップ（ステップＳ１４２）における後方投影算出部３４による第二積和演算処理（上式（３））の少なくとも一方の実行をスキップする。
これにより、反復演算部３２による演算処理の総所要時間の短縮を図っている。
【００６５】
図１に示すように、本実施形態の画像処理装置３は、実測投影データｙｉが所定値以下である結晶対１４に対応するボクセルＢｘに対して、ボクセル値ｘｊの初期値ｘｊ（０）としてゼロを与える初期値設定部６０を有する。
【００６６】
初期値設定部６０は、二次元画像生成部６１と三次元初期画像生成部６２とを備えている。
図７を参照して、初期値設定部６０の機能を説明する。
【００６７】
初期値設定部６０は、上記モデル生成ステップ（ａ）として、検出プレート１０Ａ，１０Ｂに挟まれた検出領域２２の内部に複数のボクセルＢｘが立体的に配列された解析モデル２５を作成する。三次元初期画像生成部６２を作成するにあたり、本方法では、ボクセルＢｘのセルサイズとシンチレーション検出器１１のサイズとを一致させている。
なお、本方法においては、検出プレート１０Ａ，１０Ｂに挟まれた全領域を検出領域２２とし、かつ検出領域２２の全体をボクセルＢｘで分割して解析モデル２５を構成しているが、本発明はこれに限られない。検出プレート１０Ａ，１０Ｂで挟まれた領域の一部を検出領域２２としてもよく、解析モデル２５を検出領域２２の一部領域に対応して生成してもよい。
【００６８】
図７（ａ）は解析モデル２５および検出領域２２の斜視図である。同図に示すように、本実施形態ではボクセルＢｘの形状を立方体とする。検出プレート１０Ａ，１０ＢをＸ−Ｙ平面内に設置し、シンチレーション検出器１１Ａ，１１Ｂの配列方向をＸ軸およびＹ軸方向とする。すなわち、検出プレート１０Ａ，１０Ｂは、Ｚ軸方向に対向している。
また、本方法では、検出領域２２に対してＸ，Ｙ，Ｚ軸方向にｍ個ずつのボクセルＢｘを配列するものとする。
【００６９】
図７（ｂ）は解析モデル２５および被測定物２０の斜視図である。同図に示すように、被測定物２０は検出領域２２の内部にセットされる。同図では、検出プレート１０Ａ，１０Ｂは図示を省略している。
被測定物２０としては、被験者の乳房を例示する。乳房は、ほぼ回転対称形をなすことを特徴とする部位である。具体的には、半球状または回転楕円体（半体）をなしている。
ただし、本実施形態の画像再構成システム１では、被験者のリンパ節を含む領域や、胃や肝臓などの内臓を含む領域などの人体、または他の物体を被測定物２０とすることもできる。
被測定物２０は、その回転軸ＡＸをＸ軸またはＹ軸方向に一致させて、検出プレート１０Ａ，１０Ｂにて挟持する。同図では、回転軸ＡＸをＸ軸方向に一致させている。
【００７０】
ここで、図２を参照して上述したように、投影データ生成部３１は、被測定物２０から放射される光子の光量を結晶対１４ｉごとに積算して実測投影データｙｉを計測する（実測ステップ（ｅ））。かかる計測は、駆動機構４１を駆動して検出プレート１０Ａ，１０Ｂの対向間隔を調整して、被測定物２０を検出プレート１０Ａ，１０Ｂで押圧して拡幅させた状態でおこなう。
これにより、自然状態における被測定物２０のＺ軸方向投影図よりも広い領域において光子が検出されることとなる。一方、癌細胞２３は被測定物２０の内部に存在する。このため、被測定物２０を押圧することにより、被測定物２０における癌細胞２３の周囲領域を拡大した状態で実測投影データｙｉを測定することとなる。よって、後述する三次元初期画像２６の非ゼロ領域２６ＮＺは癌細胞２３を十分に内包することとなり、被測定物２０のうち線源２１（図２を参照）が含まれる可能性のある領域に対してボクセル値ｘｊの初期値をゼロと誤設定してしまうことがない。
【００７１】
二次元画像生成部６１は、実測投影データｙｉが所定値以下である結晶対１４ｉに対応するボクセルＢｘの集合であるゼロ領域２４ＺＡと、実測投影データｙｉが所定値よりも大きい結晶対１４ｉに対応するボクセルＢｘの集合である非ゼロ領域２４ＮＺと、を含む、ボクセル値ｘｊの二次元画像２４を作成する。
【００７２】
図４を参照して、二次元画像生成部６１による二次元画像の作成処理（図３におけるステップＳ１２）を具体的に説明する。
【００７３】
二次元画像生成部６１は結晶対１４を順番に指定する（ステップＳ１２０）。本実施形態の検出プレート１０はシンチレーション検出器１１をそれぞれｍ２個ずつ備えており、結晶対１４はｍ４組存在する。
次に、二次元画像生成部６１は、指定された結晶対１４ｉに関する実測投影データｙｉをデータ記憶部３８から読み出す（ステップＳ１２１）。
【００７４】
二次元画像生成部６１は、実測投影データｙｉのゼロ判定をおこない（ステップＳ１２２）、ゼロの場合は（ステップＳ１２２：Ｎｏ）次の結晶対１４ｉに関する処理をおこなう。
実測投影データｙｉのゼロ判定においては、所定の微小な閾値を設定し、実測投影データｙｉが閾値以下の場合には、再構成される画像に与える影響が無視できる程度にゼロに近接する（実質的にゼロである）ものとして、実測投影データｙｉをゼロと判定してもよい。以下、実測投影データｙｉまたはボクセル値ｘｊの値がゼロであるとは、当該値がゼロの場合と、実質的にゼロの場合とを含むものとする。
【００７５】
そして、実測投影データｙｉが非ゼロの場合（ステップＳ１２２：Ｙｅｓ）、当該結晶対１４ｉを構成するシンチレーション検出器１１Ａ，１１Ｂの間に線源２１が存在することが分かる。二次元画像生成部６１は、結晶対１４ｉに挟まれるボクセルＢｘｊのうち、Ｚ軸方向の中間位置に存在するものに対して、当該実測投影データｙｉの値を加算する（ステップＳ１２３）。
本実施形態の検出プレート１０は、シンチレーション検出器１１が平面的に平坦に配置されているため、結晶対１４ｉの中間位置は、検出プレート１０ＡのＺ高さと，検出プレート１０ＢのＺ高さの中間に存在する二次元平面２７を描く。
これにより、当該二次元平面２７を構成するボクセルＢｘｊに対して非ゼロのボクセル値ｘｊが設定される。
【００７６】
二次元画像生成部６１は、ｉがｍ４に至るまで、上記処理を繰り返す（ステップＳ１２４）。
これにより、シンチレーション検出器１１Ａ，１１Ｂの中間高さに配置されたボクセルＢｘｊのうち、結晶対１４での検出確率の高いものに対して、大きなボクセル値ｘｊが設定される。
【００７７】
そして、二次元画像生成部６１は、積算されたボクセル値ｘｊを係数｛Ｃｊ｝で除するなどの適宜演算をおこない、二次元画像２４の濃度値を算出する（ステップＳ１２５）。
【００７８】
図７（ｃ）は、かかる二次元作成処理により生成された二次元画像２４の一例を示す斜視図である。
二次元画像２４には、濃度値がゼロであるゼロ領域２４ＺＡと、濃度値が非ゼロである非ゼロ領域２４ＮＺとが存在する。
そして、二次元画像２４の非ゼロ領域２４ＮＺにおける複数のボクセル値ｘｊは、実測投影データｙｉの値に応じて互いに相違する。すなわち、本方法では、二次元平面２７上のボクセルＢｘに対して、ステップＳ１２３において実測投影データｙｉの値を積算していくため、実測投影データｙｉの大小が二次元画像２４の画素濃度値として反映される。
そして、二次元画像２４の非ゼロ領域２４ＮＺのうち、画素濃度値の大きい領域（高濃度領域２８）は、線源２１が密集した癌細胞２３（図２を参照）を二次元平面２７に投影した形状に相当する。
【００７９】
つぎに、三次元初期画像生成部６２は、二次元画像２４を内包する三次元初期画像２６を作成し、三次元初期画像２６におけるゼロ領域２６ＺＡに対応するボクセルを決定する。
【００８０】
三次元初期画像２６の形状は特に限定されるものではない。本方法では、第一の例として、二次元画像２４を軸回転してなる三次元初期画像２６を三次元初期画像生成部６２が作成する場合を説明する。図８（ａ）は、本方法により作成された三次元初期画像２６を示す斜視図である。
【００８１】
三次元初期画像２６は、二次元平面２７上において非ゼロ領域２４ＮＺおよびゼロ領域２４ＺＡにより構成された二次元画像２４のボクセルＢｘを回転軸ＡＸまわりに１８０度または３６０度回転させたものである。
検出領域２２に立体的に存在するボクセルＢｘのうち、二次元平面２７におけるゼロ領域２４ＺＡを構成するボクセルの回転位置に対応するものについては、いずれもボクセル値ｘｊをゼロに設定する。
そして、二次元平面２７における非ゼロ領域２４ＮＺを構成するボクセルの回転位置に対応するものについては、いずれもボクセル値ｘｊを非ゼロの値に設定する。かかる非ゼロの設定値については限定されないが、二次元画像２４における画素濃度値の濃淡を三次元初期画像２６に対しても反映することにより、統計的画像再構成法に基づく画像再構成演算の収束を早めることができる。
【００８２】
図５を参照して、三次元初期画像生成部６２による三次元初期画像の作成処理（図３におけるステップＳ１３）を具体的に説明する。
【００８３】
三次元初期画像生成部６２は、Ｘ軸方向およびＹ軸方向ごとに、二次元画像２４の各画素に関し、当該画素の座標値と画素濃度値との積算値の総和を算出する（ステップＳ１３０）。
つぎに、三次元初期画像生成部６２は、上記総和を、Ｘ軸方向またはＹ軸方向の実測投影データｙｉの総量によって除することにより、各軸方向における画素濃度値の重心座標（Ｇｘ，Ｇｙ）を求める（ステップＳ１３１）。
そして、三次元初期画像生成部６２は、二次元画像２４の回転軸ＡＸとして、重心座標（Ｇｘ，Ｇｙ）を通り、Ｘ軸に平行な直線を設定する。具体的には、ｙ＝Ｇｙ、ｚ＝ｍ／２の直線が回転軸ＡＸとなる。
【００８４】
三次元初期画像生成部６２は、ボクセルＢｘを順番に指定して、当該ボクセルＢｘｊと回転軸ＡＸとのＹ−Ｚ平面内における距離を算出する（ステップＳ１３２）。
本方法では、検出領域２２に対してＸ，Ｙ，Ｚ軸方向にｍ個ずつのボクセルＢｘが配列されているため、ボクセルＢｘは全部でｍ３個存在することとなる。したがって、たとえばｊ＝１からｍ３までの引数によってボクセルＢｘを順に指定する。
【００８５】
そして、三次元初期画像生成部６２は、当該ボクセルＢｘｊのボクセル値ｘｊとして、回転軸ＡＸからの距離がこれと等しい二次元画像２４上のボクセルＢｘのボクセル値ｘを代入して与える（ステップＳ１３３）。これにより、二次元画像２４におけるゼロ領域２４ＺＡの回転位置に存在するボクセルＢｘには初期値ｘｊ（０）としてゼロが与えられる。そして、二次元画像２４における非ゼロ領域２４ＮＺの回転位置に存在するボクセルＢｘには初期値ｘｊ（０）として非ゼロの値が与えられる。
かかる処理により、二次元画像２４における画素濃度値が回転軸ＡＸまわりに与えられて三次元初期画像２６が作成される。
なお、解析モデル２５内の任意のボクセルＢｘｊに関して、回転軸ＡＸからの距離が等しい二次元平面２７上のボクセルＢｘは、回転軸ＡＸから±Ｙ軸方向に２つ存在する場合がある。かかる場合、当該ボクセルＢｘｊに対して、二次元平面２７上のどちらのボクセルＢｘのボクセル値ｘｊを与えてもよい。
【００８６】
ここで、二次元画像２４を回転軸ＡＸまわりに回転してなる回転立体は円柱形状となるため、立方体状の解析モデル２５のうちＹ軸およびＺ軸成分の絶対値の大きいボクセルＢｘは当該円柱形状よりも外側に位置することとなる。かかるボクセルＢｘに関しても、被測定物２０の外部領域にあたるとして、本方法では初期値ｘｊ（０）にゼロを与える。
【００８７】
三次元初期画像生成部６２は、ボクセル値ｘｊの初期値ｘｊ（０）の設定を、すべてのボクセルＢｘｊに対しておこなう（ステップＳ１３４）。
【００８８】
そして、図３のステップＳ１４に戻り、反復演算部３２は、三次元初期画像２６におけるゼロ領域２６ＺＡに対応するボクセルについて、前方投影ステップ（ステップＳ１４１：図６を参照）または後方投影ステップ（ステップＳ１４２）における第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００８９】
本方法においては、三次元初期画像２６のゼロ領域２６ＺＡと非ゼロ領域２６ＮＺとの境界を、被測定物２０の外表面と一致させる、または被測定物２０の外表面の僅かに外側に設定することが好ましい。これにより、被測定物２０の外表面を含むボクセルに関しては積和演算処理の反復演算を行い、被測定物２０の外部領域に関する積和演算処理をスキップすることで、画像再構成の品質と演算時間をバランスして改善することができる。
ここで、被測定物２０が軸対称形状の場合、三次元初期画像２６の各ボクセルＢｘｊの初期値ｘｊ（０）として、被測定物２０を内包する最小の領域に対して非ゼロの値を設定することにより、積和演算処理のスキップ回数が最大化されるため好適である。
【００９０】
これに対し、被測定物２０の形状が軸対称形状ではない場合、または被測定物２０の形状が不特定である場合については、三次元初期画像生成部６２は、二次元画像２４を二次元平面２７に直交する方向（Ｚ方向）に掃引してなる三次元初期画像２６を作成するとよい。これにより、被測定物２０を含むボクセルＢｘｊの初期値ｘｊ（０）として、誤ってゼロを設定してしまうことがない。
図８（ｂ）は、三次元初期画像２６の第二の例を示す斜視図である。本例では、矢印にて示すように、二次元平面２７の非ゼロ領域２４ＮＺをＺ軸方向に引き出して、柱状の三次元初期画像２６を作成している。これにより、一般的形状の被測定物２０を内包する立体形状が特定され、かかる立体形状を一部または全部に有するボクセルＢｘを非ゼロ領域２６ＮＺに設定することができる。
【００９１】
また、本方法においては、三次元初期画像生成処理（ステップＳ１３）で作成される三次元初期画像２６として、図８（ａ）に示す軸対称形状と、同図（ｂ）に示す柱状とを選択可能としてもよい。
すなわち、本実施形態の三次元初期画像生成部６２は、作成する三次元初期画像２６のタイプを指定する入力情報を受け付ける指定入力部（入力操作部５１）と、入力操作部５１が受け付けた入力情報に基づいて、二次元画像２４を軸回転してなる三次元初期画像２６、または二次元画像２４を当該二次元平面に直交する方向に掃引してなる三次元初期画像２６、のいずれかを作成してもよい。
これにより、被測定物２０が軸対称形状に近似できるときは三次元初期画像２６を軸対称形状とし、被測定物２０がその他の形状である場合には三次元初期画像２６を柱状とすることができる。
【００９２】
本方法の二次元画像作成処理（ステップＳ１２）の場合、二次元画像２４を構成するボクセル値ｘｊがゼロ領域２４ＺＡに設定されるのは、当該ボクセルＢｘを間に挟んで対向する結晶対１４のすべてについて実測投影データｙｉの計測値がゼロの場合である。すなわち、図２に示すように、シンチレーション検出器１１Ａ１と１１Ｂとの間には線源２１が存在するため、当該結晶対１４ｉに関する実測投影データｙｉは非ゼロとなり、当該結晶対１４ｉの中間位置に存在するボクセルＢｘは二次元画像２４において非ゼロ領域２４ＮＺに属することとなる。
【００９３】
一方、図２に示すシンチレーション検出器１１Ａ２と１１Ｂとの間には線源２１が存在しないため、当該結晶対１４ｉに関する実測投影データｙｉはゼロとなり、当該結晶対１４ｉの中間位置に存在するボクセルＢｘｊは二次元画像２４においてゼロ領域２４ＺＡに属することとなる。
そして、本方法の場合、二次元画像２４を構成する任意のボクセルＢｘｊの初期値ｘｊ（０）にゼロが設定されるのは、当該ボクセルＢｘｊを挟むすべての結晶対１４にかかる実測投影データｙｉがゼロの場合である。
すなわち、ボクセルＢｘに対応する複数の異なる結晶対１４ｉにかかる実測投影データｙｉが、いずれも所定の閾値以下である場合に、ボクセルＢｘｊについて第一または第二積和演算処理の少なくとも一方の実行がスキップされる。
【００９４】
したがって、任意のボクセルＢｘを挟む結晶対１４の一つにおいて、確率的な揺らぎによって光子がたまたま検出されなかったとしても、当該ボクセルＢｘｊを挟む他の結晶対１４ｉのいずれかにおいて所定の閾値以上の光量が検出された場合は、当該ボクセルＢｘの初期値ｘｊ（０）には非ゼロの値が設定される。このため、前方投影ステップ（ステップＳ１４１）や後方投影ステップ（ステップＳ１４２）をおこなうべきボクセルＢｘに関する演算処理が誤ってスキップされることはない。
【００９５】
＜画像再構成処理＞
初期値ｘｊ（０）が設定されたボクセルＢｘｊは、検出確率Ｃｉｊ、係数Ｃｊおよび実測投影データｙｉを用いた統計的画像再構成法により、ボクセル値ｘｊが推定される。
初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘｊでは、上式（２）で算出される前方投影値ｐｉは常にゼロとなる。
さらに、初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘｊは、対応する実測投影データｙｉがゼロであったことに起因するため、上式（３）で算出される後方投影値Ｄｊもゼロとなる。
したがって、第一積和演算にかかるステップＳ１４１の冒頭において、前方投影算出部３３はボクセルＢｘｊのゼロ判定をおこなう。そして、ボクセルＢｘｊがゼロである場合には、図６における第一積和演算（ステップＳ１４１）および第二積和演算（ステップＳ１４２）をともにスキップする。
【００９６】
ここで、高解像度の三次元画像を再構成するためには、ボクセルサイズを微細化することが求められる。
本方法では、画像再構成処理（ステップＳ１４およびＳ１５）において、反復演算部３２は、ボクセルＢｘのサイズを複数段階に微細化して演算処理を反復的に実行する。
そして、微細化前に三次元初期画像２６におけるゼロ領域２６ＺＡに対応していたボクセルＢｘに内包される微細化後のボクセルＢｘｊについて、反復演算部３２は、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００９７】
すなわち、本方法では、概略再構成処理（ステップＳ１４）の反復演算と、精細再構成処理（ステップＳ１５）の反復演算とを、異なる解析モデル２５を用いて行う。
図９（ａ）は、概略再構成処理に用いる解析モデル２５（概略モデル２５１）のボクセル配置を示す斜視図である。そして、同図（ｂ）は、精細再構成処理に用いる解析モデル２５（精細モデル２５２）のボクセル配置を示す斜視図である。
精細モデル２５２におけるボクセルＢｘ（ボクセルＢｘ２）は、概略モデル２５１におけるボクセルＢｘ（ボクセルＢｘ１）の辺長を各軸方向に二分の一としたものである。
したがって、精細モデル２５２は、概略モデル２５１に対して八倍のボクセル数で構成され、各軸方向に関して二倍の解像度にて被測定物２０の内部構造の三次元画像を再構成することができる。
精細モデル２５２を構成するボクセルＢｘ２は、各軸方向に隣接する他のボクセルとあわせた合計８個が、概略モデル２５１のボクセルＢｘ１の個々に対応する。すなわち、ボクセルＢｘ１とボクセルＢｘ２とは、１：Ｎ（Ｎは２以上の整数）で対応する。
【００９８】
また、精細モデル２５２においては、検出領域２２の内部が各軸方向に２ｍ個のボクセルで分割されることになる。そして、個々のシンチレーション検出器１１に対して４個のボクセルが対応する。
【００９９】
本方法では、概略モデル２５１により高速で第一および第二積和演算を実行して得られたボクセルＢｘ１のボクセル値ｘｊを、精細モデル２５２による第一および第二積和演算の初期値としてボクセルＢｘ２に与える。
したがって、概略再構成処理の初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘ１に包含されるボクセルＢｘ２に関しては、精細再構成処理の初期値として引き続きゼロが与えられる。
これにより、実測投影データｙｉがゼロである結晶対１４ｉに対応する三次元初期画像２６のゼロ領域２６ＺＡに含まれるボクセルＢｘ２に対しては、精細再構成処理における第一および第二積和演算がスキップされるため、精細再構成処理の総所要時間を抑制することができる。
【０１００】
なお、精細モデル２５２においては、一つのシンチレーション検出器１１に対して、４個のボクセルが対応する。ここで、上式（３）に示す実測投影データｙｉに関しては、投影データ生成部３１により生成されて概略再構成処理に供された投影値をそのまま用いてもよく、またはボクセルＢｘ２に対応して投影値を内挿補間した補間データを用いてもよい。
また、検出確率Ｃｉｊおよび係数Ｃｊに関しては、精細モデル２５２のボクセルＢｘ２に対応して予め算出された値を用いる。
【０１０１】
なお、本方法においては、概略モデル２５１のボクセルＢｘ１と精細モデル２５２のボクセルＢｘ２とが１：Ｎで対応する場合を例示したが、本発明はこれに限られない。すなわち、概略モデル２５１のボクセルＢｘ１と精細モデル２５２のボクセルＢｘ２とをＭ：Ｎ（ＭはＮよりも小さな２以上の整数）で対応させてもよい。
この場合、概略再構成処理で求められたボクセルＢｘ１ごとのボクセル値ｘｊを、ボクセルＢｘ１およびボクセルＢｘ２の位置関係に基づいて適宜内挿補間して、ボクセルＢｘ２の初期値を求めるとよい。
【０１０２】
ところで、図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。
【０１０３】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【０１０４】
以上、図面を参照して本発明の実施形態について述べたが、これらは本発明の例示であり、上記以外の様々な構成を採用することもできる。
【０１０５】
たとえば、上記実施形態では、画像再構成処理の初期値ｘｊ（０）を設定するために三次元初期画像２６を作成しているが、これに限られない。たとえば、画像処理装置３が、実測投影データｙｉが閾値以下である結晶対１４ｉに挟まれるボクセルＢｘの番号ｉをすべて記憶しておき、当該ボクセルＢｘｊに関する第一または第二積和演算処理をスキップしてもよい。
【０１０６】
また、上記実施形態では、ボクセル値ｘｊの初期値ｘｊ（０）としてゼロを設定する態様を例示したが、これに限られない。所定のボクセルＢｘに対して、予め定めた特定のフラグ値を初期値ｘｊ（０）として与え、前方投影算出部３３または後方投影算出部３４にてボクセル値ｘｊとフラグ値との一致判定をおこなってもよい。
【０１０７】
また、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１に示す検出部２に代えて、連続的に分布する単一の検出面１２あるいは３個以上の検出面１２を有する検出装置を使用してもよい。また、検出面１２の形状は平面に限定されるものではなく曲面であってもよい。
【０１０８】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
Ｘ線ＣＴ装置の場合、検出部２は必ずしも対向する一対の検出プレート１０を備えてはいない。この場合、光源より照射されて被測定物２０を透過するＸ線の光量を、複数の位置にて検出プレート１０で測定する。そして、被測定物２０に対する既知の入射強度と、検出プレート１０における検出強度が実質的に等しい検出器要素を検出することにより、当該検出プレート１０の位置と光源との間に被測定物２０が存在しないことが知得される。そして、かかる不存在領域に配置されたボクセルＢｘの初期値ｘｊ（０）をゼロに設定することにより、Ｘ線ＣＴ装置に関しても上記実施形態と同様の効果を得ることができる。
【０１０９】
また、上記実施形態では、統計的画像再構成法としてＭＬ−ＥＭ法が使用されているが、検出確率を用いて前方投影データまたは後方投影データが使用される統計的画像再構成法であればよく、ＭＬ−ＥＭ法に限定されるものではない。たとえば、ＯＳ−ＥＭ法やＭＡＰ−ＥＭ（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.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、たとえ従来の高速演算アルゴリズムを用いたとしても、計算量が膨大であるので、その計算速度は臨床上実用的であるとは言い難いという問題がある。たとえば、画像を再構成するまでの計算時間が長大であると、診断に要する時間が長くなるので、臨床上実用的であるとは言い難い。
【０００７】
本発明は上記課題に鑑みてなされたものであり、画像再構成処理を効率的に実行し得、計算速度を大幅に向上させ得る画像処理装置、画像再構成システム、画像処理方法およびプログラムを提供するものである。
【課題を解決するための手段】
【０００８】
本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データと、前記被測定物が配置された検出領域の内部に複数のボクセルが立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて前記投影データより前記ボクセルごとのボクセル値を推定する演算処理を反復的に実行して前記被測定物の内部構造の三次元画像を再構成する反復演算部と、
ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊを示すデータを格納しているパラメータ格納部と、
を備える画像処理装置が提供される。
この画像処理装置は、前記反復演算部が、
（Ａ）前記パラメータ格納部から読み出された検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理と、
前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理と、
を含む前記演算処理を行うとともに、
（Ｂ）前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする。
【０００９】
また、本発明の画像処理装置においては、より具体的な態様として、前記反復演算部が、前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルのボクセル値をゼロに設定してもよい。
【００１０】
また、本発明の画像処理装置においては、より具体的な態様として、前記投影データが前記所定値以下である検出器要素に対応するボクセルの集合であるゼロ領域と、前記投影データが前記所定値よりも大きい前記検出器要素に対応するボクセルの集合である非ゼロ領域と、を含む、ボクセル値の二次元画像を生成する二次元画像生成部と、
前記二次元画像を内包する三次元初期画像を生成し、前記三次元初期画像における前記ゼロ領域に対応するボクセルを決定する三次元初期画像生成部と、
をさらに備え、
前記反復演算部が、前記三次元初期画像における前記ゼロ領域に対応するボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップしてもよい。
【００１１】
また、本発明の画像処理装置においては、より具体的な態様として、前記ボクセルのサイズを複数段階に微細化して前記演算処理を反復的に実行し、
微細化前に前記三次元初期画像における前記ゼロ領域に対応していたボクセルに内包される微細化後のボクセルについて、前記反復演算部が、前記第一または第二積和演算処理の少なくとも一方の実行をスキップしてもよい。
【００１２】
また、本発明によれば、対向して配置されたシンチレーション検出器の対を前記検出器要素として含む上記記載の画像処理装置を備えるとともに、前記シンチレーション検出器が、前記被測定物から所定方向の両側に放射された光子対を検出する画像再構成システムが提供される。
【００１３】
また、本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像処理方法が提供される。
この画像処理方法は、（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影ステップと、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前記前方投影ステップおよび後方投影ステップを反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算ステップと、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００１４】
また、本発明によれば、被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像再構成処理を情報処理装置に実行させるプログラムが提供される。
このプログラムでは、前記画像再構成処理が、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成処理と、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影処理と、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影処理と、
（ｄ）前記前方投影処理および後方投影処理を反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算処理と、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００１５】
上記発明において、被測定物から検出器要素に光子が飛来するとは、光子と検出器要素との相互作用により光子が検出されることを意味する。また、被測定物から光子が飛来するとは、被測定物の後方から照射された光子が被測定物を透過して検出器要素に飛来する場合と、被測定物を光源として光子が放射される場合とを含む。以下、被測定物から光子が放出される表現した場合、上記の両態様を含む意味である。
【００１６】
なお、本発明の各種の構成要素は、その機能を実現するように形成されていればよく、たとえば、所定の機能を発揮する専用のハードウェア、所定の機能がコンピュータプログラムにより付与された画像処理装置、コンピュータプログラムにより画像処理装置に実現された所定の機能、これらの任意の組み合わせ、等として実現することができる。
また、本発明の各種の構成要素は、個々に独立した存在である必要はなく、複数の構成要素が一個の部材として形成されていること、一つの構成要素が複数の部材で形成されていること、ある構成要素が他の構成要素の一部であること、ある構成要素の一部と他の構成要素の一部とが重複していること、等を許容する。
【００１７】
また、本発明の画像処理方法は、複数の工程を順番に記載してあるが、その記載の順番は複数の工程を実行する順番を限定するものではない。このため、本発明の画像処理方法を実施するときには、その複数の工程の順番は内容的に支障しない範囲で変更することができる。
さらに、本発明の画像処理方法は、複数の工程が個々に相違するタイミングで実行されることに限定されない。このため、ある工程の実行中に他の工程が発生すること、ある工程の実行タイミングと他の工程の実行タイミングとの一部ないし全部が重複していること、等でもよい。
【００１８】
また、本発明で云う反復的な演算処理は、コンピュータプログラムを読み取って対応するデータ処理を実行できるように、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、Ｉ／Ｆ（Interface）ユニット、等の汎用デバイスで構築されたハードウェア、所定のデータ処理を実行するように構築された専用の論理回路、これらの組み合わせ、等を用いて実施することができる。
なお、本発明でプログラムに対応した各種動作を情報処理装置に実行させることは、各種デバイスを情報処理装置に動作制御させることを含む意味である。
【００１９】
さらに、本発明において、データを記憶または格納しているとは、本発明の装置が、少なくとも演算処理を実行するときに、当該データを保持している状態となる機能を有することを意味している。したがって、たとえば、光量を示す投影データや検出確率Ｃｉｊを示すデータに関しては、本発明の装置が出荷されるときに既に登録されていることの他、出荷されるときには登録されていないこれらのデータが、演算処理の時点で、または演算処理の最中に登録されることも許容する。
【発明の効果】
【００２０】
本発明による画像処理技術によれば、飛来した光子の光量を示す投影データが所定値以下である検出器要素に対応するボクセルについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。これは、光子の飛来しない領域に存在するボクセルに関してはボクセル値を一律に所定値（たとえばゼロ）とみなすことが可能であることによる。そして、第一および第二積和演算処理は統計的画像再構成法による反復的な演算処理における計算量の多い処理である。したがって、上記のボクセルに関して、かかる積和演算処理をスキップすることにより、演算処理の総時間を大幅に短縮することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明に係る一実施形態の画像再構成システムの概略構成を示す機能ブロック図である。
【図２】実測投影データを説明するための模式図である。
【図３】ＭＬ−ＥＭ法による処理手順を概略的に示すフローチャートである。
【図４】二次元画像の作成処理の具体的な手順を示すフローチャートである。
【図５】三次元初期画像の作成出力の具体的な手順を示すフローチャートである。
【図６】画像再構成処理の具体的な手順を示すフローチャートである。
【図７】（ａ）は解析モデルおよび検出領域、（ｂ）は解析モデルおよび被測定物、（ｃ）は二次元画像の一例を示す斜視図である。
【図８】三次元初期画像を示す斜視図である。
【図９】（ａ）は概略再構成処理に用いる概略モデル、（ｂ）は精細再構成処理に用いる精細モデルの斜視図である。
【発明を実施するための形態】
【００２２】
以下、本発明に係る実施の形態について図面を参照しつつ説明する。なお、すべての図面において、同様な構成要素には同一符号を付し、その詳細な説明は重複しないように適宜省略される。
図１は、本発明に係る一実施形態の画像再構成システム１および画像処理装置３の概略構成を示す機能ブロック図である。
【００２３】
本実施形態の画像再構成システム１は、画像処理装置３と検出部２とを備えている。
画像処理装置３は、反復演算部３２とパラメータ格納部３９とを備えている。
反復演算部３２は、被測定物２０から複数の検出器要素（結晶対１４：図２を参照）の各々に飛来した光子の光量を示す投影データと、被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘ（図７を参照）が立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて投影データよりボクセルＢｘごとのボクセル値ｘを推定する演算処理を反復的に実行して被測定物２０の内部構造の三次元画像を再構成する。
また、パラメータ格納部３９は、ｊ番目のボクセルＢｘｊを通過した光子がｉ番目のシンチレーション検出器１１ｉで検出される検出確率Ｃｉｊを示すデータを格納している。
【００２４】
そして、本実施形態の反復演算部３２は、パラメータ格納部３９から読み出された検出確率Ｃｉｊと、ｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理と、検出確率Ｃｉｊと、ｉ番目の検出器要素にかかる投影データとの第二積和演算処理と、を含む演算処理を行う。
また、本実施形態の反復演算部３２は、投影データが所定値以下であるシンチレーション検出器１１に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００２５】
画像再構成システム１は、互いに対向する一対の検出プレート１０（１０Ａ，１０Ｂ）を有する対向型ＰＥＴ装置を構成し、被測定物２０から放射された放射線の線量を測定して、被測定物２０の外形とともにその内部構造の三次元画像を推定（再構成）する。
【００２６】
検出部２は、互いに対向する検出面１２（１２Ａ、１２Ｂ）を有する一対の検出プレート１０Ａ，１０Ｂと、これら検出プレート１０Ａ，１０Ｂを並進移動または回転移動させる駆動機構４１と、を有している。人体などの被測定物２０は、これら検出プレート１０Ａ，１０Ｂによって挟み込まれる。
【００２７】
検出プレート１０には、多数のシンチレーション検出器１１（１１Ａ，１１Ｂ）が二次元平面内に配列されている。
シンチレーション検出器１１の配列の態様は特に限定されず、格子状、千鳥状、またはそれ以外でもよい。本実施形態では、シンチレーション検出器１１がＸ−Ｙ平面に沿って格子状（マトリクス状）に配列された検出プレート１０を例示する。
各方向のシンチレーション検出器１１の配列数は任意であるが、本実施形態では、Ｘ，Ｙ軸方向に沿って、それぞれｍ個ずつのシンチレーション検出器１１が配列されて検出プレート１０を構成しているものとする。
【００２８】
シンチレーション検出器１１は、シンチレータ結晶（図示せず）と、シンチレータ結晶の内部で放射線が消失することにより発生する蛍光を受光する受光素子（図示せず）とを有している。
検出プレート１０Ａ，１０Ｂの対向する検出面１２は、シンチレータ結晶の端面の集合により構成されている。各シンチレータ結晶は面一に配置され、検出面１２は平坦に形成されている。
本実施形態のシンチレータ結晶は、検出面１２を構成する端面が正方形をなす四角柱をなしている。ただし、シンチレータ結晶の端面の形状および寸法はこれに限られるものではない。
そして、本実施形態は、対向するシンチレーション検出器１１Ａと１１Ｂとの組み合わせごとに、被測定物２０からの放射線の飛来を検出する。
【００２９】
すなわち、本実施形態の画像再構成システム１は、対向して配置されたシンチレーション検出器１１の対（結晶対１４）を、検出器要素として備えている。そして、シンチレーション検出器１１では、被測定物２０から所定方向の両側に放射された光子対を検出する。
【００３０】
駆動機構４１は、画像処理装置３に組み込まれた制御部４０からの制御信号に従って、一対の検出プレート１０Ａ，１０Ｂを駆動する機能を有する。本実施形態では、駆動機構４１は、検出プレート１０Ａ，１０Ｂの間の間隔を狭めたり、一対の検出プレート１０Ａ，１０ＢをＸ軸の周りに回転させたり、一対の検出プレート１０Ａ，１０ＢをＸ軸方向に移動させたりすることができる。
したがって、本実施形態の画像再構成システム１は、駆動機構４１の駆動により、一対の検出面１２Ａ、１２Ｂの間隔が可変である。したがって、検出プレート１０Ａ，１０Ｂに挟み込まれた被測定物２０は検出面１２Ａ、１２Ｂに押圧されて拡幅する。画像再構成システム１では、かかる押圧状態の被測定物２０から飛来する放射線を受光する。
ユーザは、キーボードやポインティングデバイスなどの入力操作部５１を通じて制御部４０に指示を出すことにより、駆動機構４１に検出プレート１０Ａ，１０Ｂを駆動させることができる。
【００３１】
被測定物２０の内部には、陽電子を発生させる放射性同位元素（ＲＩ：Radioactive Isotope）が分布している。たとえば、被測定物２０が人体や動物などの被検体であれば、この被検体に放射性同位元素を投与することで放射線同位元素を被測定物２０内に分布させることができる。陽電子は、放射性同位元素の崩壊過程で発生する。陽電子が被測定物２０内の電子と相互作用して消滅すると、約５１１ｋｅＶのエネルギーを持つ２本のガンマ線（光子対）が発生し、これらガンマ線は互いに反対方向に放出される。よって、陽電子が消滅した地点に関して点対称性を有する一対のシンチレーション検出器１１Ａ，１１Ｂが、それぞれ、当該２本のガンマ線をほぼ同時に検出することとなる。よって、光子対をそれぞれ検出した一対のシンチレーション検出器１１Ａ，１１Ｂを結ぶ線上に、陽電子が消滅した地点が存在する。
【００３２】
シンチレーション検出器１１Ａ，１１Ｂに用いられるシンチレータ結晶は、Ｘ線やガンマ線などの高エネルギー放射線の光子を、よりエネルギーの低い多数の光子（蛍光）に変換する結晶体である。この種の結晶体としては、たとえば、ＢＧＯやＬＳＯなどの公知の結晶体を使用すればよい。
検出プレート１０Ａ，１０Ｂが備える受光素子は、シンチレーション検出器１１Ａ，１１Ｂで生成された蛍光出力を増幅し、電気信号に変換する。その電気信号は、ケーブルなどの転送手段（図示せず）を介して、光子を検出したシンチレーション検出器１１Ａ，１１Ｂの番地情報とともに、画像処理装置３の投影データ生成部３１に転送される。
【００３３】
投影データ生成部３１は、検出プレート１０Ａ，１０Ｂから転送された電気信号に基づいて、結晶対１４ごとに検出した光子を計数し、この計数結果に基づいて、結晶対１４（検出器要素）毎の投影データ（実測投影データ）を生成する。この実測投影データは、メモリ３７内のデータ記憶部３８に格納される。
【００３４】
図２は、実測投影データｙを説明するための模式図である。被測定物２０内の投影線上の陽電子消滅地点（線源）２１から所定方向の両側に放射された光子は、検出プレート１０Ａ，１０Ｂのシンチレーション検出器１１Ａ，１１Ｂでそれぞれほぼ同時に検出される。
線源２１は被測定物２０のほぼ全体に分布しており、光子は被測定物２０のほぼ全体から、さまざまな方向に対して放射される。そして、癌細胞のように活性の高い部位に線源２１は高密度で集中しているため、癌細胞を挟む結晶対１４では光子が多く検出される。
図２に示すように、特定のシンチレーション検出器１１Ｂに着目した場合、検出プレート１０Ａを構成する各シンチレーション検出器１１Ａのうち、線源２１を挟むシンチレーション検出器１１Ａ１と、線源２１を挟まないシンチレーション検出器１１Ａ２とが存在する場合がある。
他のシンチレーション検出器１１Ｂに関しても同様である。
【００３５】
投影データ生成部３１（図１を参照）は、対向するシンチレーション検出器１１Ａ，１１Ｂで所定の時間差内に光子を検出したことを検知すると、かかるシンチレーション検出器１１Ａ，１１Ｂの組み合わせに相当する結晶対１４ｉ（ｉは結晶対の番号）について、当該光子の光量を積算してカウントする。かかる計測を所定時間に亘っておこなうことにより、結晶対の番号ごとの実測投影データｙ（実測投影データｙｉ）を求めることができる。
そして、積算された実測投影データｙｉの分布に基づいて、統計的画像再構成法を用いることにより線源２１の分布密度が再構成され、もって癌細胞２３の有無およびその位置が推定される。
【００３６】
図１に示す反復演算部３２は、被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘ（立方形状または長方形状などの微小な立体要素、図７を参照）が立体的に配列された解析モデル２５を利用して画像再構成処理を実行する。
反復演算部３２は、データ記憶部３８から実測投影データｙｉを読み出し、統計的画像再構成法に従って、この実測投影データｙｉを用いてボクセルＢｘ毎の放射能濃度に相当するボクセル値ｘを推定する演算処理を反復的に実行する。
これにより、被測定物２０の内部構造の断層画像または三次元画像を再構成する機能を有している。本実施形態では、ボクセル値の三次元画像に基づいて、二次元画像である断層画像を出力することを含めて、三次元画像の再構成という。
【００３７】
本実施形態では、統計的画像再構成法として、典型的なＭＬ−ＥＭ法を使用することができる。ただし、使用可能な統計的画像再構成法はこれに限定されるものではない。
【００３８】
反復演算部３２は、前方投影算出部３３、後方投影算出部３４および濃度算出部３５という機能ブロックを有する。これらの機能ブロックは、ＭＬ−ＥＭ法に従って処理を反復的に実行する。ＭＬ−ＥＭ法では、以下の推定式（１）が使用される。
【００３９】
【数１】
【００４０】
ただし、推定式（１）中のＣｊは、次式（１Ａ）で与えられる。
【００４１】
【数２】
【００４２】
上記推定式（１）中、ｋは反復回数であり、ｘｊ（ｋ）はｋ回目の処理におけるｊ番目のボクセルＢｘｊの放射能濃度（ボクセル値ｘ）である。ｙｉはｉ番目の結晶対（結晶対１４ｉ）による検出結果に基づいて得られた、上述の実測投影データの値（投影値）である。
本実施形態の画像再構成システム１の場合、検出プレート１０は、シンチレーション検出器１１がＸ，Ｙ軸方向に沿ってそれぞれｍ個ずつ配列されて構成されている。したがって、検出プレート１０はそれぞれｍ２個のシンチレーション検出器１１を有している。そして、シンチレーション検出器１１の組み合わせである結晶対１４はｍ４組だけ存在し、ｉの最大値はｍ４となる。
【００４３】
また、ｘｊ（ｋ＋１）は、（ｋ＋１）回目の処理におけるｊ番目のボクセルＢｘｊにおける放射能濃度である。よって、ｋ回目のボクセル値ｘｊ（ｋ）から、（ｋ＋１）回目のボクセル値ｘｊ（ｋ＋１）が算出される。
ここで、ｘｊ（ｋ）の初期値ｘｊ（０）は、正の実数であるかぎりどのような値に設定されても、反復回数ｋを十分に大きく設定して演算処理を反復的に行うことで、所定の確度の三次元画像を再構成することができる。
【００４４】
推定式（１）において、Ｃｉｊは、ｊ番目のボクセルＢｘｊから飛来した光子対がｉ番目の結晶対１４ｉで検出される確率（検出確率）である。ボクセルＢｘｊから飛来した光子対が結晶対１４ｉに入射する数は、検出確率Ｃｉｊに依存することが知られている。検出確率ＣｉｊとＣｊの値は、予め算出されたうえでパラメータ格納部３９に格納されている。
【００４５】
検出確率Ｃｉｊは、行列（「システムマトリクス」と呼ばれる。）のｉ行ｊ列目の要素とみなすことができる。検出確率Ｃｉｊの値には、散乱や減衰などの物理的要因を組み込むことが可能であるが、本実施形態では、検出確率Ｃｉｊの値は、ｉ番目の結晶対１４ｉに対するｊ番目のボクセルＢｘｊの幾何学的位置のみにより定めることができる。
検出確率Ｃｉｊの値には、レイトレース法（Ray-tracing method）や、モンテカルロシミュレーション法（Monte Carlo simulation method）等によって求まる値を設定することができる。
【００４６】
次に、画像処理装置３の動作として実現される本実施形態の画像処理方法（以下、本方法という）を説明する。
【００４７】
まず、本方法の概要について説明する。
本方法は、被測定物２０から複数の検出器要素（結晶対１４）の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて被測定物２０の内部構造の三次元画像を再構成する方法である。
この本方法は、
（ａ）被測定物２０が配置された検出領域２２の内部に複数のボクセルＢｘを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルＢｘｊを通過した光子がｉ番目の結晶対１４ｉで検出される検出確率Ｃｉｊと、ｊ番目のボクセルＢｘｊのボクセル値ｘｊとの第一積和演算処理をおこなう前方投影ステップと、
（ｃ）検出確率Ｃｉｊと、ｉ番目の結晶対１４ｉにかかる投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前方投影ステップおよび後方投影ステップを反復的に実行してボクセルＢｘごとのボクセル値ｘを推定する反復演算ステップと、
を含む。
そして、本方法は、投影データが所定値以下である結晶対１４に対応するボクセルＢｘについて、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００４８】
また、本方法は、（ｅ）被測定物２０を押圧して拡幅させた状態で光子の光量を測定する実測ステップをさらに含む。
【００４９】
つぎに、図面を参照して本方法をより詳細に説明する。
図３は、ＭＬ−ＥＭ法による処理手順を概略的に示すフローチャートである。
図４は、二次元画像の作成処理（ステップＳ１２）の具体的な手順を示すフローチャートである。
図５は、三次元初期画像の作成出力（ステップＳ１３）の具体的な手順を示すフローチャートである。
図６は、画像再構成処理（ステップＳ１４およびＳ１５）の具体的な手順を示すフローチャートである。
【００５０】
ステップＳ１０では、反復演算部３２は、パラメータ格納部３９から、検出確率｛Ｃｉｊ｝と係数｛Ｃｊ｝とを読み出す。
続けて、反復演算部３２は、データ記憶部３８から実測投影データ｛ｙｉ｝を読み出す（ステップＳ１１）。なお、本実施形態では、検出確率｛Ｃｉｊ｝、係数｛Ｃｊ｝および実測投影データ｛ｙｉ｝が最初に一括して読み出されるが、これに限定されず、必要なときに読み出されてもよい。
【００５１】
続けて、反復演算部３２は、放射能濃度に相当するボクセル値｛ｘｊ｝の値を初期化するため、実測投影データｙｉに基づいて作成される二次元画像、および二次元画像に基づいて作成される三次元初期画像を作成する（ステップＳ１３、Ｓ１４）。
反復演算部３２は、その後、ボクセルサイズの大きな概略モデルを用いた画像再構成処理（ステップＳ１５）と、ボクセルサイズを微細化した精細モデルを用いた画像再構成処理（ステップＳ１６）とを実行する。
なお、本方法では、ボクセルサイズを二段階に微細化して画像再構成処理をおこなう態様を例示的に説明するが、微細化の段階はこれに限られず三段階以上としてもよい。
また、ボクセルサイズの微細化は、本発明においては任意であり、単一の解析モデルのみを利用して画像再構成処理をおこなってもよい。
【００５２】
まず、画像再構成処理（ステップＳ１５、Ｓ１６）について、その具体的な手段を、図６を参照して説明する。かかる画像再構成処理に好適に用いられる初期画像として本実施形態で作成される三次元初期画像については後述する。
【００５３】
図６を参照すると、先ず、反復演算部３２は、反復回数の上限を適切な値に設定する（ステップＳ１４０）。
概略モデルによる画像再構成処理（以下、「概略再構成処理」という場合がある。）で設定される反復回数ｋ１と、精細モデルによる画像再構成処理（以下、「精細再構成処理」という場合がある。）で設定される反復回数ｋ２とは、互いに同一でも相違してもよい。
また、反復回数ｋ１とｋ２の値、およびその大小関係も任意であるが、演算処理の総所要時間を短縮して同レベルの精度の画像再構成をおこなう観点からは、ｋ１をｋ２よりも大きく設定するとよい。
【００５４】
図１の前方投影算出部３３は、次式（２）に従い、検出確率Ｃｉｊと、推定濃度に相当するボクセル値ｘｊ（ｋ）との積和演算（第一積和演算）を実行して前方投影値ｐｉを算出する（ステップＳ１４１）。
【００５５】
【数３】
【００５６】
推定式（１）中で演算量の多大な部分が、式（２）により与えられる前方投影データの値（前方投影値）ｐｉの計算である。
前方投影算出部３３は、上式（２）の第一積和演算を、複数のスレッド（グループ）に分解して並列演算により実行してもよい。
【００５７】
上記ステップＳ１４１の処理が終了した後、図１の後方投影算出部３４は、次式（３）に従い、実測投影データｙｉと前方投影値ｐｉの逆数との積ｑｉを算出するとともに、検出確率Ｃｉｊと当該積ｑｉとの第２の積和演算（第二積和演算）を実行して後方投影値Ｄｊを算出する（ステップＳ１４２）。
【００５８】
【数４】
【００５９】
この後方投影値Ｄｊの計算も、前方投影値の算出と同様に非常に演算量の多い部分である。
【００６０】
そして、濃度算出部３５は、後方投影値Ｄｊを用いて、（ｋ＋１）回目の推定濃度に相当するボクセル値ｘｊ（ｋ＋１）（＝Ｄｊ×ｘｊ（ｋ））を算出する（ステップＳ１４３）。その後、反復回数ｋ（ｋ１、ｋ２）が上限に達したか否かを判定し（ステップＳ１４４）、反復回数ｋが上限に達するまで、上記ステップＳ１４１〜Ｓ１４４の手順を繰り返し実行する。最終的に反復回数ｋが上限に達したとき（ステップＳ１４４のＹＥＳ）は、図３のフローチャートに戻って、概略再構成処理（ステップＳ１４）または精細再構成処理（ステップＳ１５）を終了する。
【００６１】
そして、濃度算出部３５は、ボクセル値｛ｘｊ｝を出力する（ステップＳ１６）。
【００６２】
表示制御部３６は、濃度算出部３５から出力されたボクセル値｛ｘｊ｝から三次元画像を生成する。生成画像は、二次元断層画像または三次元画像としてディスプレイ５０に表示させることができる。あるいは、濃度算出部３５から出力されたボクセル値｛ｘｊ｝は、メモリ３７に記憶されてもよい。
【００６３】
＜三次元初期画像作成処理＞
ここで、従来の統計的画像再構成法では、ボクセル値の初期値として、すべてのボクセルに対して正の値を設定し、ボクセル値ｘｊ（ｋ）の反復演算をすべてのボクセルＢｘｊに関しておこなっていた。
しかし、上述のように、（ｋ＋１）回目の推定濃度に相当するボクセル値ｘｊ（ｋ＋１）は、ｋ回目の推定濃度に相当するボクセル値ｘｊ（ｋ）に後方投影値Ｄｊを乗じて算出される。したがって、ひとたびｘｊ（ｋ）がゼロになったボクセルでは、かりに以降の反復演算をおこなったとしてもボクセル値はゼロが繰り返される。
【００６４】
このため、本方法においては、反復演算部３２は、実測投影データｙｉが所定値以下である結晶対１４に対応するボクセルＢｘのボクセル値ｘｊをゼロに設定する。そして、ボクセル値ｘｊがゼロに設定されたボクセルＢｘｊに関しては、前方投影ステップ（ステップＳ１４１）における前方投影算出部３３による第一積和演算処理（上式（２））、または後方投影ステップ（ステップＳ１４２）における後方投影算出部３４による第二積和演算処理（上式（３））の少なくとも一方の実行をスキップする。
これにより、反復演算部３２による演算処理の総所要時間の短縮を図っている。
【００６５】
図１に示すように、本実施形態の画像処理装置３は、実測投影データｙｉが所定値以下である結晶対１４に対応するボクセルＢｘに対して、ボクセル値ｘｊの初期値ｘｊ（０）としてゼロを与える初期値設定部６０を有する。
【００６６】
初期値設定部６０は、二次元画像生成部６１と三次元初期画像生成部６２とを備えている。
図７を参照して、初期値設定部６０の機能を説明する。
【００６７】
初期値設定部６０は、上記モデル生成ステップ（ａ）として、検出プレート１０Ａ，１０Ｂに挟まれた検出領域２２の内部に複数のボクセルＢｘが立体的に配列された解析モデル２５を作成する。三次元初期画像生成部６２を作成するにあたり、本方法では、ボクセルＢｘのセルサイズとシンチレーション検出器１１のサイズとを一致させている。
なお、本方法においては、検出プレート１０Ａ，１０Ｂに挟まれた全領域を検出領域２２とし、かつ検出領域２２の全体をボクセルＢｘで分割して解析モデル２５を構成しているが、本発明はこれに限られない。検出プレート１０Ａ，１０Ｂで挟まれた領域の一部を検出領域２２としてもよく、解析モデル２５を検出領域２２の一部領域に対応して生成してもよい。
【００６８】
図７（ａ）は解析モデル２５および検出領域２２の斜視図である。同図に示すように、本実施形態ではボクセルＢｘの形状を立方体とする。検出プレート１０Ａ，１０ＢをＸ−Ｙ平面内に設置し、シンチレーション検出器１１Ａ，１１Ｂの配列方向をＸ軸およびＹ軸方向とする。すなわち、検出プレート１０Ａ，１０Ｂは、Ｚ軸方向に対向している。
また、本方法では、検出領域２２に対してＸ，Ｙ，Ｚ軸方向にｍ個ずつのボクセルＢｘを配列するものとする。
【００６９】
図７（ｂ）は解析モデル２５および被測定物２０の斜視図である。同図に示すように、被測定物２０は検出領域２２の内部にセットされる。同図では、検出プレート１０Ａ，１０Ｂは図示を省略している。
被測定物２０としては、被験者の乳房を例示する。乳房は、ほぼ回転対称形をなすことを特徴とする部位である。具体的には、半球状または回転楕円体（半体）をなしている。
ただし、本実施形態の画像再構成システム１では、被験者のリンパ節を含む領域や、胃や肝臓などの内臓を含む領域などの人体、または他の物体を被測定物２０とすることもできる。
被測定物２０は、その回転軸ＡＸをＸ軸またはＹ軸方向に一致させて、検出プレート１０Ａ，１０Ｂにて挟持する。同図では、回転軸ＡＸをＸ軸方向に一致させている。
【００７０】
ここで、図２を参照して上述したように、投影データ生成部３１は、被測定物２０から放射される光子の光量を結晶対１４ｉごとに積算して実測投影データｙｉを計測する（実測ステップ（ｅ））。かかる計測は、駆動機構４１を駆動して検出プレート１０Ａ，１０Ｂの対向間隔を調整して、被測定物２０を検出プレート１０Ａ，１０Ｂで押圧して拡幅させた状態でおこなう。
これにより、自然状態における被測定物２０のＺ軸方向投影図よりも広い領域において光子が検出されることとなる。一方、癌細胞２３は被測定物２０の内部に存在する。このため、被測定物２０を押圧することにより、被測定物２０における癌細胞２３の周囲領域を拡大した状態で実測投影データｙｉを測定することとなる。よって、後述する三次元初期画像２６の非ゼロ領域２６ＮＺは癌細胞２３を十分に内包することとなり、被測定物２０のうち線源２１（図２を参照）が含まれる可能性のある領域に対してボクセル値ｘｊの初期値をゼロと誤設定してしまうことがない。
【００７１】
二次元画像生成部６１は、実測投影データｙｉが所定値以下である結晶対１４ｉに対応するボクセルＢｘの集合であるゼロ領域２４ＺＡと、実測投影データｙｉが所定値よりも大きい結晶対１４ｉに対応するボクセルＢｘの集合である非ゼロ領域２４ＮＺと、を含む、ボクセル値ｘｊの二次元画像２４を作成する。
【００７２】
図４を参照して、二次元画像生成部６１による二次元画像の作成処理（図３におけるステップＳ１２）を具体的に説明する。
【００７３】
二次元画像生成部６１は結晶対１４を順番に指定する（ステップＳ１２０）。本実施形態の検出プレート１０はシンチレーション検出器１１をそれぞれｍ２個ずつ備えており、結晶対１４はｍ４組存在する。
次に、二次元画像生成部６１は、指定された結晶対１４ｉに関する実測投影データｙｉをデータ記憶部３８から読み出す（ステップＳ１２１）。
【００７４】
二次元画像生成部６１は、実測投影データｙｉのゼロ判定をおこない（ステップＳ１２２）、ゼロの場合は（ステップＳ１２２：Ｎｏ）次の結晶対１４ｉに関する処理をおこなう。
実測投影データｙｉのゼロ判定においては、所定の微小な閾値を設定し、実測投影データｙｉが閾値以下の場合には、再構成される画像に与える影響が無視できる程度にゼロに近接する（実質的にゼロである）ものとして、実測投影データｙｉをゼロと判定してもよい。以下、実測投影データｙｉまたはボクセル値ｘｊの値がゼロであるとは、当該値がゼロの場合と、実質的にゼロの場合とを含むものとする。
【００７５】
そして、実測投影データｙｉが非ゼロの場合（ステップＳ１２２：Ｙｅｓ）、当該結晶対１４ｉを構成するシンチレーション検出器１１Ａ，１１Ｂの間に線源２１が存在することが分かる。二次元画像生成部６１は、結晶対１４ｉに挟まれるボクセルＢｘｊのうち、Ｚ軸方向の中間位置に存在するものに対して、当該実測投影データｙｉの値を加算する（ステップＳ１２３）。
本実施形態の検出プレート１０は、シンチレーション検出器１１が平面的に平坦に配置されているため、結晶対１４ｉの中間位置は、検出プレート１０ＡのＺ高さと，検出プレート１０ＢのＺ高さの中間に存在する二次元平面２７を描く。
これにより、当該二次元平面２７を構成するボクセルＢｘｊに対して非ゼロのボクセル値ｘｊが設定される。
【００７６】
二次元画像生成部６１は、ｉがｍ４に至るまで、上記処理を繰り返す（ステップＳ１２４）。
これにより、シンチレーション検出器１１Ａ，１１Ｂの中間高さに配置されたボクセルＢｘｊのうち、結晶対１４での検出確率の高いものに対して、大きなボクセル値ｘｊが設定される。
【００７７】
そして、二次元画像生成部６１は、積算されたボクセル値ｘｊを係数｛Ｃｊ｝で除するなどの適宜演算をおこない、二次元画像２４の濃度値を算出する（ステップＳ１２５）。
【００７８】
図７（ｃ）は、かかる二次元作成処理により生成された二次元画像２４の一例を示す斜視図である。
二次元画像２４には、濃度値がゼロであるゼロ領域２４ＺＡと、濃度値が非ゼロである非ゼロ領域２４ＮＺとが存在する。
そして、二次元画像２４の非ゼロ領域２４ＮＺにおける複数のボクセル値ｘｊは、実測投影データｙｉの値に応じて互いに相違する。すなわち、本方法では、二次元平面２７上のボクセルＢｘに対して、ステップＳ１２３において実測投影データｙｉの値を積算していくため、実測投影データｙｉの大小が二次元画像２４の画素濃度値として反映される。
そして、二次元画像２４の非ゼロ領域２４ＮＺのうち、画素濃度値の大きい領域（高濃度領域２８）は、線源２１が密集した癌細胞２３（図２を参照）を二次元平面２７に投影した形状に相当する。
【００７９】
つぎに、三次元初期画像生成部６２は、二次元画像２４を内包する三次元初期画像２６を作成し、三次元初期画像２６におけるゼロ領域２６ＺＡに対応するボクセルを決定する。
【００８０】
三次元初期画像２６の形状は特に限定されるものではない。本方法では、第一の例として、二次元画像２４を軸回転してなる三次元初期画像２６を三次元初期画像生成部６２が作成する場合を説明する。図８（ａ）は、本方法により作成された三次元初期画像２６を示す斜視図である。
【００８１】
三次元初期画像２６は、二次元平面２７上において非ゼロ領域２４ＮＺおよびゼロ領域２４ＺＡにより構成された二次元画像２４のボクセルＢｘを回転軸ＡＸまわりに１８０度または３６０度回転させたものである。
検出領域２２に立体的に存在するボクセルＢｘのうち、二次元平面２７におけるゼロ領域２４ＺＡを構成するボクセルの回転位置に対応するものについては、いずれもボクセル値ｘｊをゼロに設定する。
そして、二次元平面２７における非ゼロ領域２４ＮＺを構成するボクセルの回転位置に対応するものについては、いずれもボクセル値ｘｊを非ゼロの値に設定する。かかる非ゼロの設定値については限定されないが、二次元画像２４における画素濃度値の濃淡を三次元初期画像２６に対しても反映することにより、統計的画像再構成法に基づく画像再構成演算の収束を早めることができる。
【００８２】
図５を参照して、三次元初期画像生成部６２による三次元初期画像の作成処理（図３におけるステップＳ１３）を具体的に説明する。
【００８３】
三次元初期画像生成部６２は、Ｘ軸方向およびＹ軸方向ごとに、二次元画像２４の各画素に関し、当該画素の座標値と画素濃度値との積算値の総和を算出する（ステップＳ１３０）。
つぎに、三次元初期画像生成部６２は、上記総和を、Ｘ軸方向またはＹ軸方向の実測投影データｙｉの総量によって除することにより、各軸方向における画素濃度値の重心座標（Ｇｘ，Ｇｙ）を求める（ステップＳ１３１）。
そして、三次元初期画像生成部６２は、二次元画像２４の回転軸ＡＸとして、重心座標（Ｇｘ，Ｇｙ）を通り、Ｘ軸に平行な直線を設定する。具体的には、ｙ＝Ｇｙ、ｚ＝ｍ／２の直線が回転軸ＡＸとなる。
【００８４】
三次元初期画像生成部６２は、ボクセルＢｘを順番に指定して、当該ボクセルＢｘｊと回転軸ＡＸとのＹ−Ｚ平面内における距離を算出する（ステップＳ１３２）。
本方法では、検出領域２２に対してＸ，Ｙ，Ｚ軸方向にｍ個ずつのボクセルＢｘが配列されているため、ボクセルＢｘは全部でｍ３個存在することとなる。したがって、たとえばｊ＝１からｍ３までの引数によってボクセルＢｘを順に指定する。
【００８５】
そして、三次元初期画像生成部６２は、当該ボクセルＢｘｊのボクセル値ｘｊとして、回転軸ＡＸからの距離がこれと等しい二次元画像２４上のボクセルＢｘのボクセル値ｘを代入して与える（ステップＳ１３３）。これにより、二次元画像２４におけるゼロ領域２４ＺＡの回転位置に存在するボクセルＢｘには初期値ｘｊ（０）としてゼロが与えられる。そして、二次元画像２４における非ゼロ領域２４ＮＺの回転位置に存在するボクセルＢｘには初期値ｘｊ（０）として非ゼロの値が与えられる。
かかる処理により、二次元画像２４における画素濃度値が回転軸ＡＸまわりに与えられて三次元初期画像２６が作成される。
なお、解析モデル２５内の任意のボクセルＢｘｊに関して、回転軸ＡＸからの距離が等しい二次元平面２７上のボクセルＢｘは、回転軸ＡＸから±Ｙ軸方向に２つ存在する場合がある。かかる場合、当該ボクセルＢｘｊに対して、二次元平面２７上のどちらのボクセルＢｘのボクセル値ｘｊを与えてもよい。
【００８６】
ここで、二次元画像２４を回転軸ＡＸまわりに回転してなる回転立体は円柱形状となるため、立方体状の解析モデル２５のうちＹ軸およびＺ軸成分の絶対値の大きいボクセルＢｘは当該円柱形状よりも外側に位置することとなる。かかるボクセルＢｘに関しても、被測定物２０の外部領域にあたるとして、本方法では初期値ｘｊ（０）にゼロを与える。
【００８７】
三次元初期画像生成部６２は、ボクセル値ｘｊの初期値ｘｊ（０）の設定を、すべてのボクセルＢｘｊに対しておこなう（ステップＳ１３４）。
【００８８】
そして、図３のステップＳ１４に戻り、反復演算部３２は、三次元初期画像２６におけるゼロ領域２６ＺＡに対応するボクセルについて、前方投影ステップ（ステップＳ１４１：図６を参照）または後方投影ステップ（ステップＳ１４２）における第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００８９】
本方法においては、三次元初期画像２６のゼロ領域２６ＺＡと非ゼロ領域２６ＮＺとの境界を、被測定物２０の外表面と一致させる、または被測定物２０の外表面の僅かに外側に設定することが好ましい。これにより、被測定物２０の外表面を含むボクセルに関しては積和演算処理の反復演算を行い、被測定物２０の外部領域に関する積和演算処理をスキップすることで、画像再構成の品質と演算時間をバランスして改善することができる。
ここで、被測定物２０が軸対称形状の場合、三次元初期画像２６の各ボクセルＢｘｊの初期値ｘｊ（０）として、被測定物２０を内包する最小の領域に対して非ゼロの値を設定することにより、積和演算処理のスキップ回数が最大化されるため好適である。
【００９０】
これに対し、被測定物２０の形状が軸対称形状ではない場合、または被測定物２０の形状が不特定である場合については、三次元初期画像生成部６２は、二次元画像２４を二次元平面２７に直交する方向（Ｚ方向）に掃引してなる三次元初期画像２６を作成するとよい。これにより、被測定物２０を含むボクセルＢｘｊの初期値ｘｊ（０）として、誤ってゼロを設定してしまうことがない。
図８（ｂ）は、三次元初期画像２６の第二の例を示す斜視図である。本例では、矢印にて示すように、二次元平面２７の非ゼロ領域２４ＮＺをＺ軸方向に引き出して、柱状の三次元初期画像２６を作成している。これにより、一般的形状の被測定物２０を内包する立体形状が特定され、かかる立体形状を一部または全部に有するボクセルＢｘを非ゼロ領域２６ＮＺに設定することができる。
【００９１】
また、本方法においては、三次元初期画像生成処理（ステップＳ１３）で作成される三次元初期画像２６として、図８（ａ）に示す軸対称形状と、同図（ｂ）に示す柱状とを選択可能としてもよい。
すなわち、本実施形態の三次元初期画像生成部６２は、作成する三次元初期画像２６のタイプを指定する入力情報を受け付ける指定入力部（入力操作部５１）と、入力操作部５１が受け付けた入力情報に基づいて、二次元画像２４を軸回転してなる三次元初期画像２６、または二次元画像２４を当該二次元平面に直交する方向に掃引してなる三次元初期画像２６、のいずれかを作成してもよい。
これにより、被測定物２０が軸対称形状に近似できるときは三次元初期画像２６を軸対称形状とし、被測定物２０がその他の形状である場合には三次元初期画像２６を柱状とすることができる。
【００９２】
本方法の二次元画像作成処理（ステップＳ１２）の場合、二次元画像２４を構成するボクセル値ｘｊがゼロ領域２４ＺＡに設定されるのは、当該ボクセルＢｘを間に挟んで対向する結晶対１４のすべてについて実測投影データｙｉの計測値がゼロの場合である。すなわち、図２に示すように、シンチレーション検出器１１Ａ１と１１Ｂとの間には線源２１が存在するため、当該結晶対１４ｉに関する実測投影データｙｉは非ゼロとなり、当該結晶対１４ｉの中間位置に存在するボクセルＢｘは二次元画像２４において非ゼロ領域２４ＮＺに属することとなる。
【００９３】
一方、図２に示すシンチレーション検出器１１Ａ２と１１Ｂとの間には線源２１が存在しないため、当該結晶対１４ｉに関する実測投影データｙｉはゼロとなり、当該結晶対１４ｉの中間位置に存在するボクセルＢｘｊは二次元画像２４においてゼロ領域２４ＺＡに属することとなる。
そして、本方法の場合、二次元画像２４を構成する任意のボクセルＢｘｊの初期値ｘｊ（０）にゼロが設定されるのは、当該ボクセルＢｘｊを挟むすべての結晶対１４にかかる実測投影データｙｉがゼロの場合である。
すなわち、ボクセルＢｘに対応する複数の異なる結晶対１４ｉにかかる実測投影データｙｉが、いずれも所定の閾値以下である場合に、ボクセルＢｘｊについて第一または第二積和演算処理の少なくとも一方の実行がスキップされる。
【００９４】
したがって、任意のボクセルＢｘを挟む結晶対１４の一つにおいて、確率的な揺らぎによって光子がたまたま検出されなかったとしても、当該ボクセルＢｘｊを挟む他の結晶対１４ｉのいずれかにおいて所定の閾値以上の光量が検出された場合は、当該ボクセルＢｘの初期値ｘｊ（０）には非ゼロの値が設定される。このため、前方投影ステップ（ステップＳ１４１）や後方投影ステップ（ステップＳ１４２）をおこなうべきボクセルＢｘに関する演算処理が誤ってスキップされることはない。
【００９５】
＜画像再構成処理＞
初期値ｘｊ（０）が設定されたボクセルＢｘｊは、検出確率Ｃｉｊ、係数Ｃｊおよび実測投影データｙｉを用いた統計的画像再構成法により、ボクセル値ｘｊが推定される。
初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘｊでは、上式（２）で算出される前方投影値ｐｉは常にゼロとなる。
さらに、初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘｊは、対応する実測投影データｙｉがゼロであったことに起因するため、上式（３）で算出される後方投影値Ｄｊもゼロとなる。
したがって、第一積和演算にかかるステップＳ１４１の冒頭において、前方投影算出部３３はボクセルＢｘｊのゼロ判定をおこなう。そして、ボクセルＢｘｊがゼロである場合には、図６における第一積和演算（ステップＳ１４１）および第二積和演算（ステップＳ１４２）をともにスキップする。
【００９６】
ここで、高解像度の三次元画像を再構成するためには、ボクセルサイズを微細化することが求められる。
本方法では、画像再構成処理（ステップＳ１４およびＳ１５）において、反復演算部３２は、ボクセルＢｘのサイズを複数段階に微細化して演算処理を反復的に実行する。
そして、微細化前に三次元初期画像２６におけるゼロ領域２６ＺＡに対応していたボクセルＢｘに内包される微細化後のボクセルＢｘｊについて、反復演算部３２は、第一または第二積和演算処理の少なくとも一方の実行をスキップする。
【００９７】
すなわち、本方法では、概略再構成処理（ステップＳ１４）の反復演算と、精細再構成処理（ステップＳ１５）の反復演算とを、異なる解析モデル２５を用いて行う。
図９（ａ）は、概略再構成処理に用いる解析モデル２５（概略モデル２５１）のボクセル配置を示す斜視図である。そして、同図（ｂ）は、精細再構成処理に用いる解析モデル２５（精細モデル２５２）のボクセル配置を示す斜視図である。
精細モデル２５２におけるボクセルＢｘ（ボクセルＢｘ２）は、概略モデル２５１におけるボクセルＢｘ（ボクセルＢｘ１）の辺長を各軸方向に二分の一としたものである。
したがって、精細モデル２５２は、概略モデル２５１に対して八倍のボクセル数で構成され、各軸方向に関して二倍の解像度にて被測定物２０の内部構造の三次元画像を再構成することができる。
精細モデル２５２を構成するボクセルＢｘ２は、各軸方向に隣接する他のボクセルとあわせた合計８個が、概略モデル２５１のボクセルＢｘ１の個々に対応する。すなわち、ボクセルＢｘ１とボクセルＢｘ２とは、１：Ｎ（Ｎは２以上の整数）で対応する。
【００９８】
また、精細モデル２５２においては、検出領域２２の内部が各軸方向に２ｍ個のボクセルで分割されることになる。そして、個々のシンチレーション検出器１１に対して４個のボクセルが対応する。
【００９９】
本方法では、概略モデル２５１により高速で第一および第二積和演算を実行して得られたボクセルＢｘ１のボクセル値ｘｊを、精細モデル２５２による第一および第二積和演算の初期値としてボクセルＢｘ２に与える。
したがって、概略再構成処理の初期値ｘｊ（０）としてゼロが与えられたボクセルＢｘ１に包含されるボクセルＢｘ２に関しては、精細再構成処理の初期値として引き続きゼロが与えられる。
これにより、実測投影データｙｉがゼロである結晶対１４ｉに対応する三次元初期画像２６のゼロ領域２６ＺＡに含まれるボクセルＢｘ２に対しては、精細再構成処理における第一および第二積和演算がスキップされるため、精細再構成処理の総所要時間を抑制することができる。
【０１００】
なお、精細モデル２５２においては、一つのシンチレーション検出器１１に対して、４個のボクセルが対応する。ここで、上式（３）に示す実測投影データｙｉに関しては、投影データ生成部３１により生成されて概略再構成処理に供された投影値をそのまま用いてもよく、またはボクセルＢｘ２に対応して投影値を内挿補間した補間データを用いてもよい。
また、検出確率Ｃｉｊおよび係数Ｃｊに関しては、精細モデル２５２のボクセルＢｘ２に対応して予め算出された値を用いる。
【０１０１】
なお、本方法においては、概略モデル２５１のボクセルＢｘ１と精細モデル２５２のボクセルＢｘ２とが１：Ｎで対応する場合を例示したが、本発明はこれに限られない。すなわち、概略モデル２５１のボクセルＢｘ１と精細モデル２５２のボクセルＢｘ２とをＭ：Ｎ（ＭはＮよりも小さな２以上の整数）で対応させてもよい。
この場合、概略再構成処理で求められたボクセルＢｘ１ごとのボクセル値ｘｊを、ボクセルＢｘ１およびボクセルＢｘ２の位置関係に基づいて適宜内挿補間して、ボクセルＢｘ２の初期値を求めるとよい。
【０１０２】
ところで、図１の画像処理装置３の機能ブロック３１〜３６の全部または一部は、半導体集積回路などのハードウェアで実現されてもよいし、あるいは、不揮発性メモリや光ディスクなどの記録媒体に記録されたプログラムまたはプログラムコードで実現されてもよい。このようなプログラムまたはプログラムコードは、機能ブロック３１〜３６の全部または一部の処理を、ＣＰＵなどのプロセッサを有するコンピュータに実行させるものである。
【０１０３】
また、データ記憶部３８やパラメータ格納部３９は、揮発性メモリまたは不揮発性メモリなどの記録媒体（たとえば、半導体メモリや磁気記録媒体）と、この記録媒体に対してデータの書込と読出を行うための回路やプログラムとで構成することができる。データ記憶部３８やパラメータ格納部３９は、予め記録媒体の所定の記憶領域上に構成されていてもよいし、あるいは、画像処理装置３の動作時に割り当てられた適当な記憶領域上に構成されてもよい。
【０１０４】
以上、図面を参照して本発明の実施形態について述べたが、これらは本発明の例示であり、上記以外の様々な構成を採用することもできる。
【０１０５】
たとえば、上記実施形態では、画像再構成処理の初期値ｘｊ（０）を設定するために三次元初期画像２６を作成しているが、これに限られない。たとえば、画像処理装置３が、実測投影データｙｉが閾値以下である結晶対１４ｉに挟まれるボクセルＢｘの番号ｉをすべて記憶しておき、当該ボクセルＢｘｊに関する第一または第二積和演算処理をスキップしてもよい。
【０１０６】
また、上記実施形態では、ボクセル値ｘｊの初期値ｘｊ（０）としてゼロを設定する態様を例示したが、これに限られない。所定のボクセルＢｘに対して、予め定めた特定のフラグ値を初期値ｘｊ（０）として与え、前方投影算出部３３または後方投影算出部３４にてボクセル値ｘｊとフラグ値との一致判定をおこなってもよい。
【０１０７】
また、上記実施形態では、一対の検出プレート１０Ａ，１０Ｂを使用し、これら検出プレート１０Ａ，１０Ｂの検出面は平面であったが、これらに限定されるものではない。図１に示す検出部２に代えて、連続的に分布する単一の検出面１２あるいは３個以上の検出面１２を有する検出装置を使用してもよい。また、検出面１２の形状は平面に限定されるものではなく曲面であってもよい。
【０１０８】
上記実施形態では、画像処理装置３は対向型ＰＥＴ装置に適用されたが、これに限定されるものではない。画像処理装置３は、Ｘ線ＣＴ装置などの他の検出装置にも適用することが可能である。
Ｘ線ＣＴ装置の場合、検出部２は必ずしも対向する一対の検出プレート１０を備えてはいない。この場合、光源より照射されて被測定物２０を透過するＸ線の光量を、複数の位置にて検出プレート１０で測定する。そして、被測定物２０に対する既知の入射強度と、検出プレート１０における検出強度が実質的に等しい検出器要素を検出することにより、当該検出プレート１０の位置と光源との間に被測定物２０が存在しないことが知得される。そして、かかる不存在領域に配置されたボクセルＢｘの初期値ｘｊ（０）をゼロに設定することにより、Ｘ線ＣＴ装置に関しても上記実施形態と同様の効果を得ることができる。
【０１０９】
また、上記実施形態では、統計的画像再構成法としてＭＬ−ＥＭ法が使用されているが、検出確率を用いて前方投影データまたは後方投影データが使用される統計的画像再構成法であればよく、ＭＬ−ＥＭ法に限定されるものではない。たとえば、ＯＳ−ＥＭ法やＭＡＰ−ＥＭ（Maximum A Posteriori - Expectation Maximization）法などのＭＬ−ＥＭ法の改良版を使用することが可能である。
【０１１０】
さらに、画像再構成処理に用いる実測投影データｙｉは、上記実施形態のように検出部２で計測して取得してもよく、または記憶媒体から画像処理装置３のデータ記憶部３８に取り込んで記憶してもよい。
【符号の説明】
【０１１１】
１ 画像再構成システム
２ 検出部
３ 画像処理装置
１０ 検出プレート
１１ シンチレーション検出器
１２ 検出面
１４ 結晶対
２０ 被測定物
２１ 線源
２２ 検出領域
２３ 癌細胞
２４ 二次元画像
２４ＮＺ 非ゼロ領域
２４ＺＡ ゼロ領域
２５ 解析モデル
２５１ 概略モデル
２５２ 精細モデル
２６ 三次元初期画像
２６ＮＺ 非ゼロ領域
２６ＺＡ ゼロ領域
２７ 二次元平面
２８ 高濃度領域
３１ 投影データ生成部
３２ 反復演算部
３３ 前方投影算出部
３４ 後方投影算出部
３５ 濃度算出部
３６ 表示制御部
３７ メモリ
３８ データ記憶部
３９ パラメータ格納部
４０ 制御部
４１ 駆動機構
５０ ディスプレイ
５１ 入力操作部
６０ 初期値設定部
６１ 二次元画像生成部
６２ 三次元初期画像生成部
ＡＸ 回転軸
Ｂｘ ボクセル
ｘ ボクセル値
ｙ 実測投影データ
【特許請求の範囲】
【請求項１】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データと、前記被測定物が配置された検出領域の内部に複数のボクセルが立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて前記投影データより前記ボクセルごとのボクセル値を推定する演算処理を反復的に実行して前記被測定物の内部構造の三次元画像を再構成する反復演算部と、
ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊを示すデータを格納しているパラメータ格納部と、
を備える画像処理装置であって、
前記反復演算部が、
（Ａ）前記パラメータ格納部から読み出された検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理と、
前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理と、
を含む前記演算処理を行うとともに、
（Ｂ）前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理装置。
【請求項２】
前記反復演算部が、前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルのボクセル値をゼロに設定する請求項１に記載の画像処理装置。
【請求項３】
前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルに対して、前記ボクセル値の初期値としてゼロを与える初期値設定部を有する請求項２に記載の画像処理装置。
【請求項４】
ボクセルに対応する複数の異なる検出器要素にかかる前記投影データが、いずれも前記所定値以下である場合に、当該ボクセルについて前記第一または第二積和演算処理の少なくとも一方の実行をスキップする請求項１から３のいずれかに記載の画像処理装置。
【請求項５】
前記投影データが前記所定値以下である前記検出器要素に対応するボクセルの集合であるゼロ領域と、前記投影データが前記所定値よりも大きい前記検出器要素に対応するボクセルの集合である非ゼロ領域と、を含む、ボクセル値の二次元画像を生成する二次元画像生成部と、
前記二次元画像を内包する三次元初期画像を生成し、前記三次元初期画像における前記ゼロ領域に対応するボクセルを決定する三次元初期画像生成部と、
をさらに備え、
前記反復演算部が、前記三次元初期画像における前記ゼロ領域に対応するボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする請求項１から４のいずれかに記載の画像処理装置。
【請求項６】
前記三次元初期画像生成部が、前記二次元画像を軸回転してなる前記三次元初期画像を生成する請求項５に記載の画像処理装置。
【請求項７】
前記三次元初期画像生成部が、前記二次元画像を当該二次元平面に直交する方向に掃引してなる前記三次元初期画像を生成する請求項５に記載の画像処理装置。
【請求項８】
前記三次元初期画像生成部が、
生成する前記三次元初期画像のタイプを指定する入力情報を受け付ける指定入力部と、
前記指定入力部が受け付けた前記入力情報に基づいて、前記二次元画像を軸回転してなる前記三次元初期画像、または前記二次元画像を当該二次元平面に直交する方向に掃引してなる前記三次元初期画像、のいずれかを生成することを特徴とする請求項５に記載の画像処理装置。
【請求項９】
前記二次元画像の前記非ゼロ領域における複数のボクセル値が、前記投影データの値に応じて互いに相違することを特徴とする請求項５から８のいずれかに記載の画像処理装置。
【請求項１０】
前記ボクセルのサイズを複数段階に微細化して前記演算処理を反復的に実行する請求項５から９のいずれかに記載の画像処理装置であって、
微細化前に前記三次元初期画像における前記ゼロ領域に対応していたボクセルに内包される微細化後のボクセルについて、前記反復演算部が、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理装置。
【請求項１１】
対向して配置されたシンチレーション検出器の対を前記検出器要素として含む請求項１から１０のいずれかに記載の画像処理装置を備えるとともに、
前記シンチレーション検出器が、前記被測定物から所定方向の両側に放射された光子対を検出することを特徴とする画像再構成システム。
【請求項１２】
複数のシンチレーション検出器がそれぞれ配列された一対の対向する検出面を備えるとともに、前記検出面同士の対向間隔が可変であることを特徴とする請求項１１に記載の画像再構成システム。
【請求項１３】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像処理方法であって、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影ステップと、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前記前方投影ステップおよび後方投影ステップを反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算ステップと、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理方法。
【請求項１４】
前記被測定物を押圧して拡幅させた状態で前記光子の光量を測定する実測ステップをさらに含む請求項１３に記載の画像処理方法。
【請求項１５】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像再構成処理を情報処理装置に実行させるプログラムであって、
前記画像再構成処理が、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成処理と、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影処理と、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影処理と、
（ｄ）前記前方投影処理および後方投影処理を反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算処理と、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とするプログラム。
【請求項１】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データと、前記被測定物が配置された検出領域の内部に複数のボクセルが立体的に配列された解析モデルと、に基づいて、統計的画像再構成法を用いて前記投影データより前記ボクセルごとのボクセル値を推定する演算処理を反復的に実行して前記被測定物の内部構造の三次元画像を再構成する反復演算部と、
ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊを示すデータを格納しているパラメータ格納部と、
を備える画像処理装置であって、
前記反復演算部が、
（Ａ）前記パラメータ格納部から読み出された検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理と、
前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理と、
を含む前記演算処理を行うとともに、
（Ｂ）前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理装置。
【請求項２】
前記反復演算部が、前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルのボクセル値をゼロに設定する請求項１に記載の画像処理装置。
【請求項３】
前記投影データが前記所定値以下である検出器要素に対応する前記ボクセルに対して、前記ボクセル値の初期値としてゼロを与える初期値設定部を有する請求項２に記載の画像処理装置。
【請求項４】
ボクセルに対応する複数の異なる検出器要素にかかる前記投影データが、いずれも前記所定値以下である場合に、当該ボクセルについて前記第一または第二積和演算処理の少なくとも一方の実行をスキップする請求項１から３のいずれかに記載の画像処理装置。
【請求項５】
前記投影データが前記所定値以下である前記検出器要素に対応するボクセルの集合であるゼロ領域と、前記投影データが前記所定値よりも大きい前記検出器要素に対応するボクセルの集合である非ゼロ領域と、を含む、ボクセル値の二次元画像を生成する二次元画像生成部と、
前記二次元画像を内包する三次元初期画像を生成し、前記三次元初期画像における前記ゼロ領域に対応するボクセルを決定する三次元初期画像生成部と、
をさらに備え、
前記反復演算部が、前記三次元初期画像における前記ゼロ領域に対応するボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップする請求項１から４のいずれかに記載の画像処理装置。
【請求項６】
前記三次元初期画像生成部が、前記二次元画像を軸回転してなる前記三次元初期画像を生成する請求項５に記載の画像処理装置。
【請求項７】
前記三次元初期画像生成部が、前記二次元画像を当該二次元平面に直交する方向に掃引してなる前記三次元初期画像を生成する請求項５に記載の画像処理装置。
【請求項８】
前記三次元初期画像生成部が、
生成する前記三次元初期画像のタイプを指定する入力情報を受け付ける指定入力部と、
前記指定入力部が受け付けた前記入力情報に基づいて、前記二次元画像を軸回転してなる前記三次元初期画像、または前記二次元画像を当該二次元平面に直交する方向に掃引してなる前記三次元初期画像、のいずれかを生成することを特徴とする請求項５に記載の画像処理装置。
【請求項９】
前記二次元画像の前記非ゼロ領域における複数のボクセル値が、前記投影データの値に応じて互いに相違することを特徴とする請求項５から８のいずれかに記載の画像処理装置。
【請求項１０】
前記ボクセルのサイズを複数段階に微細化して前記演算処理を反復的に実行する請求項５から９のいずれかに記載の画像処理装置であって、
微細化前に前記三次元初期画像における前記ゼロ領域に対応していたボクセルに内包される微細化後のボクセルについて、前記反復演算部が、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理装置。
【請求項１１】
対向して配置されたシンチレーション検出器の対を前記検出器要素として含む請求項１から１０のいずれかに記載の画像処理装置を備えるとともに、
前記シンチレーション検出器が、前記被測定物から所定方向の両側に放射された光子対を検出することを特徴とする画像再構成システム。
【請求項１２】
複数のシンチレーション検出器がそれぞれ配列された一対の対向する検出面を備えるとともに、前記検出面同士の対向間隔が可変であることを特徴とする請求項１１に記載の画像再構成システム。
【請求項１３】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像処理方法であって、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成ステップと、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影ステップと、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影ステップと、
（ｄ）前記前方投影ステップおよび後方投影ステップを反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算ステップと、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とする画像処理方法。
【請求項１４】
前記被測定物を押圧して拡幅させた状態で前記光子の光量を測定する実測ステップをさらに含む請求項１３に記載の画像処理方法。
【請求項１５】
被測定物から複数の検出器要素の各々に飛来した光子の光量を示す投影データに基づいて、統計的画像再構成法を用いて前記被測定物の内部構造の三次元画像を再構成する画像再構成処理を情報処理装置に実行させるプログラムであって、
前記画像再構成処理が、
（ａ）前記被測定物が配置された検出領域の内部に複数のボクセルを立体的に配列して解析モデルを生成するモデル生成処理と、
（ｂ）ｊ番目のボクセルを通過した前記光子がｉ番目の検出器要素で検出される検出確率Ｃｉｊと、前記ｊ番目のボクセルのボクセル値との第一積和演算処理をおこなう前方投影処理と、
（ｃ）前記検出確率Ｃｉｊと、前記ｉ番目の検出器要素にかかる前記投影データとの第二積和演算処理をおこなう後方投影処理と、
（ｄ）前記前方投影処理および後方投影処理を反復的に実行して前記ボクセルごとのボクセル値を推定する反復演算処理と、
を含むとともに、
前記投影データが所定値以下である前記検出器要素に対応する前記ボクセルについて、前記第一または第二積和演算処理の少なくとも一方の実行をスキップすることを特徴とするプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１０−２０４７５５（Ｐ２０１０−２０４７５５Ａ）
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願番号】特願２００９−４７１７９（Ｐ２００９−４７１７９）
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
【公開日】平成２２年９月１６日（２０１０．９．１６）
【国際特許分類】
【出願日】平成２１年２月２７日（２００９．２．２７）
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【出願人】（５０９０２３１３７）医療法人みやぎクリニック (4)
【Ｆターム（参考）】
[ Back to top ]