(57)【要約】【課題】画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる画像再構成処理方法および画像処理方法を提供することを目的とする。【解決手段】従来のような画像の微分空間を正則化する手法と相違して、直線状空間に沿った、重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成する。つまり、投影空間を正則化する処理（正則化処理）を再構成処理に組み込んで行う。その際に、直線状空間に沿って生じるストリークアーティファクトＳＡに沿って重み付き加算値に相当する線積分値を求める。上述した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成するので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。【選択図】図７
2016018245.tif 000002
【特許請求の範囲】
【請求項１】
画像信号計測装置によって得られた被検体の計測データに基づいて、再構成処理を行って複数次元のデジタル画像を作成する画像再構成処理方法であって、
数学的最適化における正則化において、複数次元の直線状空間または曲線状空間に沿った、前記複数次元のデジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いて前記再構成処理を行って前記再構成処理後の複数次元のデジタル画像を作成する再構成処理工程
を備えることを特徴とする画像再構成処理方法。
【請求項２】
請求項１に記載の画像再構成処理方法において、
前記計測データは、投影データ，計数データ，時系列イベントデータまたは周波数空間データのいずれかであることを特徴とする画像再構成処理方法。
【請求項３】
請求項１または請求項２に記載の画像再構成処理方法において、
前記直線状空間は線分であり、前記重み付き加算値は推定対象である再構成画像の前記線分上の画素値の積分値である投影値であることを特徴とする画像再構成処理方法。
【請求項４】
請求項３に記載の画像再構成処理方法において、
前記投影値は、前記線分に平行に沿って加算した、前記複数次元のデジタル画像の画素値の積分値であることを特徴とする画像再構成処理方法。
【請求項５】
請求項１から請求項４のいずれかに記載の画像再構成処理方法において、
前記画像信号計測装置は陽電子放射断層撮影装置であり、
前記再構成処理工程では、前記正則化関数を用いて、当該陽電子放射断層撮影装置で得られた前記計測データに基づいて逐次近似法による前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項６】
請求項５に記載の画像再構成処理方法において、
前記陽電子放射断層撮影装置は、一部が開口して構成された検出器ユニットを備え、
当該検出器ユニットで検出された前記計測データに基づいて前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項７】
請求項６に記載の画像再構成処理方法において、
前記直線状空間は、前記検出器ユニットの開口領域を通過する直線状空間を含むことを特徴とする画像再構成方法。
【請求項８】
請求項５から請求項７のいずれかに記載の画像再構成処理方法において、
前記陽電子放射断層撮影装置は、平均同時計数時間分解能が１ナノ秒以下である放射線検出器を備え、
当該放射線検出器で検出された前記計測データに基づいて前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項９】
請求項１から請求項４のいずれかに記載の画像再構成処理方法において、
前記画像信号計測装置は単一光子放射断層撮影装置であり、
前記再構成処理工程では、前記正則化関数を用いて、当該単一光子放射断層撮影装置で得られた前記計測データに基づいて逐次近似法による前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項１０】
請求項１から請求項４のいずれかに記載の画像再構成処理方法において、
前記画像信号計測装置はＸ線コンピュータ断層撮影装置であり、
前記再構成処理工程では、前記正則化関数を用いて、当該Ｘ線コンピュータ断層撮影装置で得られた前記計測データに基づいて逐次近似法による前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項１１】
請求項１から請求項４のいずれかに記載の画像再構成処理方法において、
前記画像信号計測装置は核磁気共鳴断層撮影装置であり、
前記再構成処理工程では、前記正則化関数を用いて、当該核磁気共鳴断層撮影装置で得られた前記計測データに基づいて逐次近似法による前記再構成処理を行うことを特徴とする画像再構成処理方法。
【請求項１２】
入力されたデジタル画像に対して画像処理を行う画像処理方法であって、
数学的最適化における正則化において、複数次元の直線状空間または曲線状空間に沿った、推定対象であるデジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いることを特徴とする画像処理方法。
【請求項１３】
請求項１２に記載の画像処理方法において、
前記直線状空間は線分であり、前記重み付き加算値は推定対象である再構成画像の前記線分上の画素値の積分値である投影値であることを特徴とする画像処理方法。
【請求項１４】
請求項１３に記載の画像処理方法において、
前記投影値は、前記線分に平行に沿って加算した前記画素値の積分値であることを特徴とする画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、再構成処理を行って複数次元のデジタル画像を作成する画像再構成処理方法、および入力されたデジタル画像に対して画像処理を行う画像処理方法に係り、特に、正則化処理によって統計誤差からなるノイズやそれ以外からなるノイズを抑制する技術に関する。
【背景技術】
【０００２】
この画像再構成処理方法は、断層画像撮影装置（ＣＴ(Computed Tomography)装置）全般の画像再構成技術に用いられ、この画像処理方法は、上述した断層画像撮影装置における画像再構成以外の画像処理技術やデジタル画像撮影装置全般の画像処理技術に用いられる。断層画像撮影装置としては、例えば、核医学診断装置やＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）や核磁気共鳴断層撮影装置（ＭＲＩ(Magnetic Resonance Imaging)装置）があり、画像信号計測装置で構成されており、画像信号計測装置によって得られた被検体の計測データに基づいて、再構成処理を行って複数次元のデジタル画像（断層画像や３次元再構成画像など）を作成する。デジタル画像撮影装置としては、顕微鏡カメラやデジタルカメラやデジタルビデオカメラなどの断層画像撮影装置以外の汎用カメラ装置があり、入力されたデジタル画像に対して画像処理を行う。
【０００３】
また、核医学診断装置としては、陽電子放射断層撮影装置（ＰＥＴ(Positron Emission Tomography)装置）や単一光子放射断層撮影装置（ＳＰＥＣＴ（Single Photon Emission CT）装置）がある。ＰＥＴ装置は、陽電子（ポジトロン）の消滅によって発生する複数本の放射線（γ線）を検出して複数個の検出器で放射線（γ線）を同時に検出したときのみ（つまり同時計数したときのみ）検出情報を記録する。そして、記録された検出情報に対して再構成処理を行って被検体の断層画像を作成する。ＳＰＥＣＴ装置は、単一の放射線（γ線）を検出して検出情報を記録する。そして、記録された検出情報に対して再構成処理を行って被検体の断層画像を作成する。
【０００４】
かかる処理によって得られる各種の画像には統計誤差からなるノイズやそれ以外からなるノイズが必ず含まれており、ノイズ低減処理が行われる。従来、この種のノイズ低減処理として、数学的最適化における正則化関数を用いた正則化処理が知られ、正則化処理によってノイズを抑制する技術がある（例えば、特許文献１および非特許文献１参照）。
【０００５】
特許文献１は、Ｘ線ＣＴ装置の画像生成処理において、画像再構成の前処理として、実測された投影データに対して正則化処理（ノイズ低減処理）を施し、画像再構成処理では、当該正則化処理後の投影データを利用するという内容である。Ｆｉｇ．３は当該特許文献１に記載されている画像生成処理のフローチャートである。当該特許文献１の効果は、再構成画像の画質改善（例えばストリークアーティファクトの抑制）である。
【０００６】
非特許文献１は、ＰＥＴ装置の統計的画像再構成処理において、画像の微分変換（互いに隣接する、あるいは互いに近傍の画素値の差分）に対する正則化関数（下記（１）式を参照）を目的関数に加えることにより、再構成画像に生じる誤差（主に統計誤差）を抑制するという内容である。なお、統計的画像再構成処理とは、画像信号計測装置（ここではＰＥＴ装置）によって得られた被検体の計測データ（測定データ）の統計的性質から導かれる目的関数（尤度関数）の最適化計算に基づく画像再構成処理である。正則化関数は下記（１）式で表され、画像の微分変換は下記（２）式で定義される。
【０００７】
【数１】
2016018245.tif 000003
【０００８】
ここで、ｘは画像ベクトルであり、ｘｊはｊ番目の画素の画素値であり、Ｎｊはｊ番目の画素に対する近傍画素領域（隣接画素領域）であり、ｋはｋ番目の画素であり近傍画素領域Ｎｊ内の画素であり、κｊｋは、ｊ番目の画素の画素値ｘｊとｋ番目の画素の画素値ｘｋとを近づけるように設定するための係数（平滑化係数）である。上記（２）式中の「otherwise」は当該（２）式における条件文（｜ｘｊ−ｘｋ｜≦δｈ）以外の場合、すなわち｜ｘｊ−ｘｋ｜＞δｈの条件の場合を示す。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】米国特許第７，２０６，４４０号明細書（Ｆｉｇ．３）
【非特許文献】
【００１０】
【非特許文献１】J. Qi et al., “Fully 3D Bayesian Image Reconstruction for the ECAT EXACT HR+”, IEEE Transaction on Nuclear Science, vol.2, 953-957, 1997.
【発明の概要】
【発明が解決しようとする課題】
【００１１】
しかしながら、このような構成を有する従来例の場合には、次のような問題がある。
すなわち、上述した特許文献１の画像処理技術は、実測値である投影データ（計測データ）に含まれるノイズが再構成画像に悪影響を与える場合に、その影響を抑制する手法として有効である。しかし、再構成画像に生じる悪影響が実測されていない投影データに由来する場合（例えば投影データの一部分が欠損していることに由来する場合）には有効でないという問題がある。
【００１２】
また、上述した非特許文献１の画像再構成技術は、再構成画像に生じる誤差の空間的相関性が弱い場合に、その誤差を抑制する手法として有効である。しかし、再構成画像に生じる誤差の空間的相関性が強い場合（例えばストリークアーティファクトが生じる場合）には有効でないという問題がある。すなわち、ストリークアーティファクトのように空間的相関性が強い場合には、ストリークアーティファクトから外れた箇所の誤差を弱めることができるが、ストリークアーティファクト自体はノイズでない画素値と認識され、ストリークアーティファクト自体を消すことができない。
【００１３】
この発明は、このような事情に鑑みてなされたものであって、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる画像再構成処理方法および画像処理方法を提供することを目的とする。
【課題を解決するための手段】
【００１４】
この発明は、このような目的を達成するために、次のような構成をとる。
すなわち、この発明に係る画像再構成処理方法は、画像信号計測装置によって得られた被検体の計測データに基づいて、再構成処理を行って複数次元のデジタル画像を作成する画像再構成処理方法であって、数学的最適化における正則化において、複数次元の直線状空間または曲線状空間に沿った、前記複数次元のデジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いて前記再構成処理を行って前記再構成処理後の複数次元のデジタル画像を作成する再構成処理工程を備えることを特徴とするものである。
【００１５】
［作用・効果］この発明に係る画像再構成処理方法によれば、数学的最適化における正則化において、正則化処理が組み込まれた再構成処理を行う。すなわち、複数次元の直線状空間または曲線状空間に沿った、複数次元のデジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成する。ここで、重み付き加算値を計算する場合には、重み係数を“１”にして単に加算することも含まれることに留意されたい。
【００１６】
上述した非特許文献１では画像の微分空間を正則化する。したがって、従来のような画像の微分空間を正則化する手法では、誤差の空間的相関性が強い箇所ではノイズでない画素値と認識され、誤差の空間的相関性が強い箇所でのノイズ自体を抑制することができない。それに対して、この発明に係る画像再構成処理方法では、複数の重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成する。つまり、投影空間を正則化する処理（正則化処理）を再構成処理に組み込んで行う。したがって、従来のような画像の微分空間を正則化する手法と比較すると誤差の空間的相関性が強い箇所でのノイズ自体を抑制することができる。また、従来の正則化手法（例えば上述した特許文献１の手法）では計測データが投影データの場合には、正則化処理後のデータ（つまり復元されるデータ）も投影データであるが、この発明に係る画像再構成処理方法では正則化処理を再構成処理に組み込んで行うので、正則化処理後のデータ（つまり復元されるデータ）は再構成処理後の複数次元のデジタル画像（例えば断層画像や３次元再構成画像）である。また、複数の重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成するので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
【００１７】
上述した計測データは、投影データ，計数データ，時系列イベントデータまたは周波数空間データのいずれかである。投影データとしては、Ｘ線ＣＴ装置で得られた測定対象物（減弱係数分布）の投影データやＰＥＴ装置で得られた測定対象物（放射能分布）のサイノグラムなどがある。計数データとしては、ＰＥＴ装置の放射線検出器で同時計数されたカウント値などがある。時系列イベントデータとしては、ＰＥＴ装置の放射線検出器で得られた検出イベント情報（検出器番号、検出時間、γ線のエネルギ等）を時系列で保存したリストモードデータなどがある。周波数空間データとしては、ＭＲＩ装置で得られた測定対象物（プロトン密度分布など）の周波数空間データなどがある。
【００１８】
また、直線状空間に沿って重み付き加算値を求める場合には、直線状空間は線分であり、重み付き加算値は再構成画像の線分上の、複数次元のデジタル画像の画素値（例えば、逐次近似計算における推定途中の画像の画素値）の積分値である投影値である。このように投影値を求めることにより、画像に生じる誤差の空間的相関性を正確に抽出することができて、ノイズをより効果的に低減させることができる。
【００１９】
また、投影値は、線分に平行に沿って加算した上述した（複数次元のデジタル画像の）画素値の積分値であるのが好ましい。もし、図８（ａ）に示すようにストリークアーティファクトＳＡに直交した線分Ｌ１に平行に沿って加算（すなわちストリークアーティファクトＳＡに直交して加算）すると、主に投影値の低周波成分にストリークアーティファクトの影響が表れてしまう。つまり、投影値にはストリークアーティファクトの影響が現れにくい（すなわち、ストリークアーティファクトの影響を正確に抽出することができない）ので、ストリークアーティファクトが低減されにくくなる。なぜなら、投影空間の正則化では投影空間の高周波成分を抑制することで再構成画像のノイズ低減を図るためである。そこで、図８（ｂ）に示すようにストリークアーティファクトＳＡに沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２に平行に沿って加算することにより、得られた投影値は、ストリークアーティファクトの影響が現れた（すなわち、ストリークアーティファクトの影響を正確に抽出した）投影値となる。もちろん、投影値については、線分に平行に沿って加算した画素値の積分値に限定されず、投影線の空間的な拡がりを考慮した投影値であってもよい。つまり、投影線上の厳密な積分値（線積分値）でなくてもよい。
【００２０】
なお、上述した画像信号計測装置の一例は陽電子放射断層撮影装置（ＰＥＴ装置）である。この場合には、上述した再構成処理工程では、上述した正則化関数を用いて、当該陽電子放射断層撮影装置で得られた計測データに基づいて逐次近似法による再構成処理を行う。
【００２１】
画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、一部が開口して構成された検出器ユニットを備え、当該検出器ユニットで検出された計測データに基づいて再構成処理を行うのが有効である。すなわち、リング状の検出器ユニットとは相違して、一部が開口して構成された検出器ユニットの場合には、開口部分（抜け部分）を通過する放射線は検出されない。よって、投影データの一部分が欠損することに起因して、再構成画像などに代表される再構成処理後の複数次元のデジタル画像に空間的相関性の強いノイズ（例えばストリークアーティファクト）が生じる。この発明に係る画像再構成処理方法を、（一部が開口して構成された）検出器ユニットで検出された計測データに基づく再構成処理に適用することにより、空間的相関性の強いノイズを低減させることができる。
【００２２】
また、（一部が開口して構成された）検出器ユニットを備えた場合において、直線状空間は、検出器ユニットの開口領域を通過する直線状空間を含むのが好ましい。投影データの一部分が欠損することに起因したストリークアーティファクトＳＡ（図８（ｂ）を参照）に沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２（図８（ｂ）を参照）に平行に沿って加算することにより、ストリークアーティファクトをより一層抑制することができる。
【００２３】
また、画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、陽電子放射断層撮影装置は、平均同時計数時間分解能が１ナノ秒以下である放射線検出器を備え、当該放射線検出器で検出された計測データに基づいて再構成処理を行うのが好ましい。平均同時計数時間分解能が１ナノ秒を超えるとポジトロンの消滅位置や放射線による発光位置を正確に推定することができない。平均同時計数時間分解能が１ナノ秒以下である放射線検出器を用いることで、ポジトロンの消滅位置や放射線による発光位置を正確に推定することができ、より正確な断層画像を取得することができる。
【００２４】
また、画像信号計測装置の他の一例は単一光子放射断層撮影装置（ＳＰＥＣＴ装置）である。この場合には、上述した再構成処理工程では、上述した正則化関数を用いて、当該単一光子放射断層撮影装置で得られた計測データに基づいて逐次近似法による再構成処理を行う。
【００２５】
また、画像信号計測装置のさらなる他の一例はＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）である。上述した再構成処理工程では、上述した正則化関数を用いて、当該Ｘ線コンピュータ断層撮影装置で得られた計測データに基づいて逐次近似法による再構成処理を行う。
【００２６】
また、画像信号計測装置のさらなる他の一例は核磁気共鳴断層撮影装置（ＭＲＩ装置）である。上述した再構成処理工程では、上述した正則化関数を用いて、当該核磁気共鳴断層撮影装置で得られた計測データに基づいて逐次近似法による再構成処理を行う。
【００２７】
また、この発明に係る画像処理方法は、入力されたデジタル画像に対して画像処理を行う画像処理方法であって、数学的最適化における正則化において、複数次元の直線状空間または曲線状空間に沿った、前記デジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いることを特徴とするものである。
【００２８】
［作用・効果］この発明に係る画像処理方法によれば、画像再構成処理方法と同様に、数学的最適化における正則化において、正則化処理が組み込まれた画像処理を行う。ここでの画像処理は、再構成処理以外の画像処理も含まれる。すなわち、複数次元の直線状空間または曲線状空間に沿った、デジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いることにより入力されたデジタル画像に対して画像処理を行う。画像再構成処理方法と同様に、重み付き加算値を計算する場合には、重み係数を“１”にして単に加算することも含まれることに留意されたい。
【００２９】
上述した特許文献１では実測された投影データに対して正則化処理を行う。それに対して、この発明に係る画像処理方法では当該直線状空間または曲線状空間に沿った、デジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いて正則化処理を行う。つまり、未知の求めたいデジタル画像の投影空間が正則である（すなわち滑らかである）ことを直接的に利用することにより、デジタル画像の投影空間を正則化する。また、複数の重み付き加算値の相互類似度を評価した正則化関数を用いて正則化処理を行うので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
【００３０】
なお、上述した入力されたデジタル画像は、画像に出力される構造物が実際の形態で表されるものであれば良く、例えば、上述したＰＥＴ画像，ＳＰＥＣＴ画像，Ｘ線ＣＴ画像，ＭＲＩ画像，顕微鏡カメラ画像，デジタルカメラ画像，デジタル画像などであり、カラー画像も含まれる。
【００３１】
画像再構成処理方法と同様に、直線状空間に沿って重み付き加算値を求める場合には、直線状空間は線分であり、重み付き加算値は再構成画像（例えば、逐次近似計算における推定途中の画像）の線分上の画素値の積分値である投影値である。このように投影値を求めることにより、画像に生じる誤差の空間的相関性を正確に抽出することができて、ノイズをより効果的に低減させることができる。
【００３２】
また、画像再構成処理方法と同様に、投影値は、線分に平行に沿って加算した画素値の積分値であるのが好ましい。図８（ｂ）に示すようにストリークアーティファクトＳＡに沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２に平行に沿って加算することにより、得られた投影値は、ストリークアーティファクトの影響が現れた（すなわち、ストリークアーティファクトの影響を正確に抽出した）投影値となる。もちろん、画像再構成処理方法と同様に、投影値については、線分に平行に沿って加算した画素値の積分値に限定されず、投影線の空間的な拡がりを考慮した投影値であってもよい。つまり、投影線上の厳密な積分値（線積分値）でなくてもよい。
【発明の効果】
【００３３】
この発明に係る画像再構成処理方法によれば、数学的最適化における正則化において、正則化処理が組み込まれた再構成処理を行う。すなわち、複数次元の直線状空間または曲線状空間に沿った、複数次元のデジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成する。複数の重み付き加算値の相互類似度を評価した正則化関数を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成するので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
また、この発明に係る画像処理方法によれば、数学的最適化における正則化において、正則化処理が組み込まれた画像処理を行う。すなわち、複数次元の直線状空間または曲線状空間に沿った、デジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数を用いることにより入力されたデジタル画像に対して画像処理を行う。複数の重み付き加算値の相互類似度を評価した正則化関数を用いて正則化処理を行うので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
【図面の簡単な説明】
【００３４】
【図１】各実施例に係る部分リング型ＰＥＴ装置のγ線検出器配置の一実施態様を示す概略斜視図およびブロック図である。
【図２】γ線検出器の概略斜視図である。
【図３】実施例１に係る画像再構成処理のフローチャートである。
【図４】実施例１に係る画像再構成処理のフローチャートである。
【図５】（ａ）は部分リング型ＰＥＴ装置のジオメトリの一例を示す模式図、（ｂ）は（ａ）に対応する２次元投影空間の欠損領域および非欠損領域である。
【図６】投影関数の設定の説明に供する模式図であり、（ａ）は投影線が有限な体積を有さない場合の模式図、（ｂ）は投影線が有限な体積を有する場合の模式図である。
【図７】ストリークアーティファクトと投影値との関係を示した模式図である。
【図８】（ａ）はストリークアーティファクトに直交して直線状空間を設定したときの模式図、（ｂ）はストリークアーティファクトに沿って直線状空間を設定したときの模式図である。
【図９】実施例２に係る画像処理のフローチャートである。
【図１０】実施例２に係る画像処理のフローチャートである。
【図１１】（ａ）および（ｂ）は実施例３に係る曲線状空間を設定したときの模式図である。
【図１２】変形例に係るマンモグラフィ装置の側面図およびブロック図である。
【実施例１】
【００３５】
以下、図面を参照してこの発明の実施例１を説明する。
図１は、各実施例に係る部分リング型ＰＥＴ装置のγ線検出器配置の一実施態様を示す概略斜視図およびブロック図であり、図２は、γ線検出器の概略斜視図である。後述する実施例２、３も含めて本実施例１では、画像信号計測装置として、陽電子放射断層撮影装置（ＰＥＴ装置）を例に採って説明する。また、図１および図２は各実施例とも共通の構成である。
【００３６】
図１に示すように、部分リング型ＰＥＴ装置１は、検出器ユニット２Ａ，２Ｂを備えている。検出器ユニット２Ａ，２Ｂ内には複数のγ線検出器３が埋設されている。部分リング型ＰＥＴ装置１は、この発明における画像信号計測装置に相当し、この発明における陽電子放射断層撮影装置にも相当する。また、検出器ユニット２Ａ，２Ｂは、検出器ユニットに相当し、γ線検出器３は、この発明における放射線検出器に相当する。
【００３７】
なお、本実施例１では、検出器ユニット２Ａ，２Ｂは、平均同時計数時間分解能が５００ピコ秒（ｐｓ）あるいは３００ピコ秒（ｐｓ）程度の性能を有している。平均同時計数時間分解能については上述の値（５００ｐｓあるいは３００ｐｓ）に限定されないが、平均同時計数時間分解能が１ナノ秒（ｎｓ）以下であるのが好ましい。
【００３８】
また、検出器ユニット２Ａ，２Ｂは、一部が開口して構成されている。すなわち、検出器ユニット２Ａ，２Ｂ間には、γ線検出器３が存在しない開放領域（開口領域）が存在する。図１の場合には、開放領域（開口領域）はｘｚ平面（図５（ａ）および図６（ａ）を参照）に沿って存在しているので、検出器ユニット２Ａ，２Ｂは左右配置型のジオメトリとなる。もちろん、開放領域（開口領域）はｘｚ平面方向に限定されず、ｙｚ平面方向に沿って開放領域（開口領域）が存在するようにγ線検出器３を配置してもよい（このときは検出器ユニット２Ａ，２Ｂは上下配置型のジオメトリとなる）。また、ｙｚ平面，ｘｚ平面以外の平面に沿って開放領域（開口領域）が存在するようにγ線検出器３を配置してもよい。
【００３９】
その他にも、部分リング型ＰＥＴ装置１は、同時計数回路４と演算回路５とを備えている。図１では、γ線検出器３から同時計数回路４への結線を２つのみ図示しているが、実際には、γ線検出器３の光電子増倍管(PMT: Photo Multiplier Tube)３３（図２を参照）の総チャンネル数分、同時計数回路４に接続されている。
【００４０】
放射性薬剤が投与された被検体（図示省略）から発生したγ線をγ線検出器３のシンチレータブロック３１（図２を参照）が光に変換して、変換されたその光をγ線検出器３の光電子増倍管（ＰＭＴ）３３（図２を参照）は増倍させて電気信号に変換する。その電気信号を同時計数回路４に送り込み、カウント値のデータを生成する。カウント値は、この発明における計数データに相当し、この発明における計測データにも相当する。
【００４１】
具体的には、被検体（図示省略）に放射性薬剤を投与すると、ポジトロン放出型のＲＩのポジトロンが消滅することにより、２本のγ線が発生する。同時計数回路４は、シンチレータブロック３１（図２を参照）の位置とγ線の入射タイミングとをチェックし、被検体の両側にある２つのシンチレータブロック３１でγ線が同時に入射したときのみ、送り込まれた電気信号を適正なデータと判定する。一方のシンチレータブロック３１のみにγ線が入射したときには、同時計数回路４は棄却する。つまり、同時計数回路４は、上述した電気信号に基づいて、２つのγ線検出器３においてγ線が同時観測（すなわち同時計数）されたことを検出する。
【００４２】
同時計数回路４により同時計数と判定された適正なデータにより構成された計測データ（ここではカウント値）を、演算回路５に送り込む。演算回路５は、後述するステップＳ１〜Ｓ９（図３および図４を参照）を行って、部分リング型ＰＥＴ装置１によって得られた被検体（図示省略）の計測データ（ここではカウント値）に基づいて、再構成処理を行って複数次元のデジタル画像（ここでは再構成画像）を作成する。演算回路５の具体的な機能については後述する。
【００４３】
γ線検出器３は、図２に示すようにシンチレータブロック３１と、そのシンチレータブロック３１に対して光学的に結合されたライトガイド３２と、そのライトガイド３２に対して光学的に結合された光電子増倍管（以下、単に「ＰＭＴ」と略記する）３３とを備えている。シンチレータブロック３１を構成する各シンチレータ素子は、γ線の入射に伴って発光することでγ線から光に変換する。この変換によってシンチレータ素子はγ線を検出する。シンチレータ素子において発光した光がシンチレータブロック３１で十分に拡散されて、ライトガイド３２を介してＰＭＴ３３に入力される。ＰＭＴ３３は、シンチレータブロック３１で変換された光を増倍させて電気信号に変換する。その電気信号は同時計数回路４（図１を参照）に送り込まれる。
【００４４】
また、γ線検出器３は、図２に示すように、３次元的に配置されたシンチレータ素子からなり、深さ方向に複数の層からなるＤＯＩ検出器である。図２では、４層のＤＯＩ検出器を図示しているが、層の数については、複数であれば特に限定されない。
【００４５】
ここで、ＤＯＩ検出器は、各々のシンチレータ素子を放射線の深さ方向に積層して構成されたものであり、相互作用を起こした深さ(DOI: Depth of Interaction)方向と横方向（入射面に平行な方向）との座標情報を重心演算により求める。ＤＯＩ検出器を用いることにより深さ方向の空間分解能をより一層向上させることができる。よって、ＤＯＩ検出器の層の数は、深さ方向に積層されたシンチレータ素子の層の数である。
【００４６】
次に、演算回路５の具体的な機能について、図３〜図８を参照して説明する。図３および図４は、実施例１に係る画像再構成処理のフローチャートであり、図５（ａ）は、部分リング型ＰＥＴ装置のジオメトリの一例を示す模式図であり、図５（ｂ）は、図５（ａ）に対応する２次元投影空間の欠損領域および非欠損領域であり、図６は、投影関数の設定の説明に供する模式図であり、図６（ａ）は投影線が有限な体積を有さない場合の模式図であり、図６（ｂ）は投影線が有限な体積を有する場合の模式図であり、図７は、ストリークアーティファクトと投影値との関係を示した模式図であり、図８（ａ）は、ストリークアーティファクトに直交して直線状空間を設定したときの模式図であり、図８（ｂ）はストリークアーティファクトに沿って直線状空間を設定したときの模式図である。
【００４７】
（ステップＳ１）投影空間の設定
図３に示すように、先ずＰＥＴ装置のジオメトリ情報Ｇを用いて投影空間の設定を行う。投影空間は、動径ｒ（図５（ａ）を参照），投影角度φ（図５（ａ）を参照），体軸座標ｚ（図５（ａ）を参照），傾斜角度θ（図示省略）およびＴＯＦ情報ｔ（図３では「ＰＥＴ装置の時間分解能」で表記）を座標（変数）とする５次元空間である。
【００４８】
ここで、図５（ａ）に示すように、水平軸をｘ軸とし、鉛直軸をｙ軸とし、奥行き方向の水平軸（被検体の体軸）をｚ軸とする（図６（ａ）も参照）。γ線検出器が存在しない開放領域（開口領域）を含ませると円筒（図５（ａ）では２次元で図示しているので円）の視野領域Ｖ（図５（ｂ）に示すように直径はＤ）となり、当該視野領域Ｖの中心を通る軸が動径ｒとなる。動径ｒと直交する線を投影線としたときに、ｙ軸と当該投影線とのなす角度を投影角度φとする。体軸ｚに直交する平面と投影線とがなす角度を傾斜角度（ティルト角度）θとする。消滅放射線の検出時間差 (TOF:Time Of Flight)をＴＯＦ情報ｔとする。
【００４９】
上述した投影空間は、一般に「平行投影空間」と呼ばれる。したがって、ステップＳ１の投影空間の設定とは、それぞれの座標軸の定義域および標本化間隔を規定することである。ここで、投影空間は、５つの座標軸のうち１つ以上の軸を除いた低次元空間（４次元以下）であってもよい。また、投影空間は、円錐（コーンビーム）投影空間などの非平行投影空間であってもよい。
【００５０】
（ステップＳ２）投影空間の欠損領域の設定
次に、ＰＥＴ装置のジオメトリ情報Ｇを用いて投影空間の欠損領域の設定を行う。具体的には、図５（ａ）に示すようにＰＥＴ装置によって測定される投影線を非欠損投影線１１とし、ＰＥＴ装置によって測定されない投影線を欠損投影線１２とし、ＰＥＴ装置のジオメトリ情報Ｇから欠損投影線１２を算出する。このようにして投影空間の欠損領域Ｓｎｕｌｌを設定する。
【００５１】
ここで、ＰＥＴ装置の場合には同時計数されたときのみ有効なデータとなるので、一方のγ線検出器３（図１を参照）のみから検出された（すなわち他方のγ線検出器３から検出されない）ときの投影線も、欠損投影線１２に含ませる。例えば、図５（ａ）に示す左右配置型の装置ジオメトリの場合、２次元投影空間（ここでは動径ｒを横軸，投影角度φを縦軸としたサイノグラム）の欠損領域および非欠損領域は図５（ｂ）のようになる。符号１３は非欠損領域であり、符号１４は欠損領域である。
【００５２】
（ステップＳ３）投影関数の設定
次に、投影関数の設定を行う。投影関数をｒｍ（・）（∀ｍ∈Ｓｎｕｌｌ）とする。なお、∀は任意の集合に属する全ての要素を意味する。例えば、再構成画像をｘ（Ｊ次元ベクトル）としたときに、図６（ａ）に示すように、Ｊ次元ベクトルの再構成画像ｘのｊ番目（１≦ｊ≦Ｊ）の画素（ボクセル）２１（太枠で図示）とｍ番目の欠損投影線２２との交差長をｒｍｊとする。そして、交差長ｒｍｊを第ｊ要素とするＪ次元ベクトルｒｍと再構成画像ｘとの内積ｒｍ・ｘを投影関数ｒｍ（ｘ）（＝ｒｍ・ｘ）とする。ここで、図６（ｂ）に示すように、投影線を有限な体積を有するｍ番目のチューブ２４と考え、ｊ番目の画素（ボクセル）２３と当該チューブ２４との交差体積をｒｍｊとしてもよい。図６（ａ）および図６（ｂ）では、図示の便宜上、２次元で図示しているが実際の空間は３次元であることに留意されたい。したがって、画素は３次元のボクセルとなり、図６（ｂ）の投影線はチューブとなる。符号２５は全画素の集合からなる画像空間である。
【００５３】
このようにして、内積ｒｍ・ｘを表す投影関数ｒｍ（ｘ）は、複数次元の直線状空間または曲線状空間（本実施例１では欠損領域Ｓｎｕｌｌに属する投影線）に沿った、再構成画像（例えば、逐次近似再構成計算の途中の画像）の画素値に関する重み付き加算値を意味する。
【００５４】
また、本実施例１では、図５（ａ）および図６（ａ）に示すように直線状空間（平行投影空間）を採用しており、直線状空間は線分であり、重み付き加算値に相当する投影関数ｒｍ（ｘ），ｒｎ（ｘ）は、線分ｍ，ｎ上の画素値の積分値である投影値（つまり線積分値）と見なすことができる。部分リング型ＰＥＴ装置１（図１を参照）において、図１や図５（ａ）や図６（ａ）に示す検出器配置では、正則化処理を伴わない画像再構成処理を行った場合、検出器ユニット２Ａ，２Ｂの開口領域を通過する直線状空間に沿って、図７に示すようにストリークアーティファクトＳＡが生じる。よって、ストリークアーティファクトＳＡが生じる部分での投影値は、他の箇所での投影値と比較すると、山３１あるいは谷３２となって急峻に変化する。
【００５５】
このことを鑑みると、図８（ａ）に示すようにストリークアーティファクトＳＡに直交した線分Ｌ１に平行に沿って加算すると、主に投影値の低周波成分にストリークアーティファクトの影響が表れる。そこで、図８（ｂ）に示すようにストリークアーティファクトＳＡに沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２に平行に沿って加算することにより、得られた投影値は、ストリークアーティファクトの影響が現れた（すなわち、ストリークアーティファクトの影響を正確に抽出した）投影値となる。このステップＳ３は、後述する図４におけるステップＳ９（目的関数の最適化計算）で行われる。
【００５６】
（ステップＳ４）距離関数の設定
次に、距離関数の設定を行う。互いに隣接する、あるいは互いに近傍の２つの欠損投影線をｍ，ｎ（∀ｍ，∀ｎ∈Ｓｎｕｌｌ）とし、２つの欠損投影線ｍ，ｎの投影値の距離を表す距離関数をφｍｎ（ｕ，ｖ）とする。ただし、距離関数φｍｎ（ｕ，ｖ）は、ｕ＝ｖのときに“０”、ｕ≠ｖのときに０以上の値をとる任意の実数値関数である。例えば、距離関数φｍｎ（ｕ，ｖ）を２次関数ｗｍｎ（ｕ−ｖ）２／２とする（ただし、ｗｍｎは任意の正の値である）。ここで、距離関数φｍｎ（ｕ，ｖ）は２次関数以外の関数（例えば、絶対値関数ｗｍｎ｜ｕ−ｖ｜、Ｌｏｇ−ｃｏｓｈ関数ｗｍｎｌｏｇ［ｃｏｓｈ｛（ｕ−ｖ）／δ｝］（ただしδ＞０）など）であってもよい。
【００５７】
（ステップＳ５）投影値の期待値の設定？
次に、投影値の期待値を設定するのか否か？の設定を行う。投影値の期待値をｐｍ（∀ｍ∈Ｓｎｕｌｌ）とする。
【００５８】
（ステップＳ６）ｐｍ＝１（∀ｍ∈Ｓｎｕｌｌ）
設定しない場合（図３では「Ｎｏ」で表記）にはｐｍ＝１（∀ｍ∈Ｓｎｕｌｌ）とする。
【００５９】
（ステップＳ７）ｐｍ（∀ｍ∈Ｓｎｕｌｌ）の設定
設定する場合（図３では「Ｙｅｓ」で表記）には、以下の手順により投影値の期待値ｐｍを設定する。
【００６０】
ＴＯＦ情報ｔ（図３では「ＰＥＴ装置の時間分解能」で表記）を用いて、図３に示すようにＰＥＴ装置の測定データｄを周知の画像再構成手法（ＭＬ−ＥＭ法、ＯＳ−ＥＭ法、ＤＲＡＭＡ法などの逐次近似法）で再構成する。当該ＰＥＴ装置の測定データｄは、実際の被検体を用いて得られたデータであり、後述する図４におけるステップＳ８（画像再構成問題の目的関数の設定）以降での実際の再構成処理と同一のデータが用いられる。よって、２回に分けて撮影を行う必要はない。再構成によって得られた仮の再構成画像をｘ´とし、そのときの投影データをｒｍ（ｘ´）とする。再構成によって得られた仮の再構成画像ｘ´の投影データｒｍ（ｘ´）を計算し、その投影データｒｍ（ｘ´）に対してノイズ低減処理（例えば移動平均フィルタ）を施した後の投影値をｐｍとして、投影値の期待値（ノイズ低減処理後の投影値）ｐｍを求める。
【００６１】
（ステップＳ８）画像再構成問題の目的関数の設定
次に、図４に示すように、ＰＥＴ装置のジオメトリ情報Ｇ，ＴＯＦ情報ｔ（図４では「ＰＥＴ装置の時間分解能」で表記）およびＰＥＴ装置の測定データｄに基づいて、画像再構成問題（すなわち統計的画像再構成）の目的関数の設定を行う。画像再構成問題（統計的画像再構成）の目的関数をＦ（ｘ）とする。正則化関数をＵ（ｘ）とし、画像再構成問題（統計的画像再構成）の目的関数Ｆ（ｘ）に、未知のＰＥＴ画像（再構成画像）ｘの投影データ（ラドン変換）に基づいた正則化関数Ｕ（ｘ）を加えることにより、部分リング型ＰＥＴ装置１（ここでは部分リング型ＴＯＦ−ＰＥＴ装置）の再構成画像に生じるストリークアーティファクトを抑制することを目的としている。
【００６２】
この正則化関数Ｕ（ｘ）は、「ＰＥＴ画像の投影データの値は局所的相関性が強い」という先験的な知識（事前情報）を鑑みてデザインされた関数である。つまり、正則化関数Ｕ（ｘ）は、複数の重み付き加算値（投影関数ｒｍ（ｘ），ｒｎ（ｘ））の相互類似度を評価した関数である。
【００６３】
本実施例１では、具体的には、再構成画像ｘの投影データを計算したときに、ＰＥＴ装置のジオメトリ情報（検出器配置）から定まる投影データの欠損領域Ｓｎｕｌｌに属する、互いに隣接する、あるいは互いに近傍の２つの欠損投影線ｍ，ｎの線積分値（すなわち投影関数）ｒｍ（ｘ），ｒｎ（ｘ）の相互相関性（相互類似度）が高まるように、下記（３）式に示す正則化関数Ｕ（ｘ）を用いる。
【００６４】
【数２】
2016018245.tif 000004
【００６５】
ステップＳ４で設定された距離関数φｍｎ（ｕ，ｖ）において、ｕ＝ｒｍ（ｘ）／ｐｍを代入するとともに、ｖ＝ｒｎ（ｘ）／ｐｎを代入して、上記（３）式により正則化関数Ｕ（ｘ）を設定する。上述したように、距離関数φｍｎ（ｕ，ｖ）は、ｕ＝ｖのときに“０”、ｕ≠ｖのときに０以上の値をとる任意の実数値関数である。また、ｐｍはステップＳ６あるいはＳ７で設定されたｍ番目の投影値の期待値であり、ｐｎはステップＳ６あるいはＳ７で設定されたｎ番目の投影値の期待値である。ここで、ｐｍやｐｎで除算したのは規格化のためで、図７に示すように端部の場合には投影値において互いに隣接する画素値の差分が大きくなり、その部分までが平滑化されるのを防止するためである。本明細書では、上記（３）式の正則化関数Ｕ（ｘ）を「投影空間の正則化関数」と呼ぶ。
【００６６】
上記（３）式で設定された正則化関数Ｕ（ｘ）を用いて目的関数Ｆ（ｘ）を設定する。例えば、ポアソン分布の対数尤度関数をＬ（ｘ）とすると、ポアソン分布の対数尤度関数Ｌ（ｘ）と投影空間の正則化関数Ｕ（ｘ）とを用いて、Ｌ（ｘ）＋βＵ（ｘ）を目的関数Ｆ（ｘ）とする。ただし、βは正の定数であり、正則化の強さを調整するパラメータである。なお、β＝０のときには投影空間の正則化関数Ｕ（ｘ）の項が消えて、Ｆ（ｘ）＝Ｌ（ｘ）となり通常の従来の再構成処理となる。ポアソン分布の対数尤度関数Ｌ（ｘ）は、下記（４）式で定義される。
【００６７】
【数３】
2016018245.tif 000005
【００６８】
ただし、ＩはＰＥＴ装置の測定データｄのデータ点数であり、ａｉｊはｊ番目の画素から放射された消滅放射線がｉ番目のデータ点で検出される確率であり、ｙｉはｉ番目のデータ点での計数値（カウント値）であり、ｓｉはｙｉに含まれる偶発同時計数および散乱同時計数の計数値（カウント値）の期待値である。ここで、対数尤度関数Ｌ（ｘ）は、リストモードデータの対数尤度関数であってもよい。
【００６９】
（ステップＳ９）目的関数の最適化計算
次に、目的関数Ｆ（ｘ）の最適化計算を行う。「最適化計算」とは目的関数Ｆ（ｘ）を最小にする画像ベクトル（Ｊ次元ベクトルの再構成画像）ｘを求めることであり、反復アルゴリズム（「逐次近似法」とも呼ばれる）を用いて実行する。つまり、目的関数Ｆ（ｘ）と投影空間の正則化関数Ｕ（ｘ）との和が最小になるように逆問題を解くことになる。
【００７０】
反復アルゴリズム（逐次近似法）としては、例えば、De PierroのMAP-EM法(A.R. De Pierro. A modified expectation maximization algorithm for penalized likelihood estimation tomography. IEEE Trans. Med. Imaging, Vol. 14, pp. 132-137, 1995.)、GreenのOSL-MAP-EM法(P.J. Green. Bayesian reconstructions from emission tomography data using a modified EM algorithm. IEEE Trans. Med. Imaging, Vol. 9, pp. 84-93, 1990.)がある。下記（５）式は、De PierroのMAP-EM法を参考にして導出した反復アルゴリズムである。
【００７１】
【数４】
2016018245.tif 000006
【００７２】
ただし、ｋは反復回数（図４では「反復回数入力Ｉｔ」で表記）であり、ｃｊは正則化関数Ｕ（ｘ）の曲率の近似値であり、εは小さな正の定数であり、上記（５）式中のｍａｘ（ａ，ｂ）はａ，ｂのうち大きいの方の値を返す関数である。上記（５）式の反復アルゴリズムは、ordered-subsets法を導入することで、反復計算の収束性を高めることも可能である。反復アルゴリズムの終了条件は、ユーザが指定（入力）した反復回数（反復回数入力Ｉｔ）に従う。
【００７３】
先ず、初期画像ｘ（０）を設定する。初期画像ｘ（０）については、例えば一様な画素値を有する再構成画像であればよく、初期画像ｘ（０）＞０とする。
【００７４】
上記（５）式にしたがって、反復回数ｋの値を繰り上げて再構成画像の画素値ｘ（ｋ）を更新する。反復アルゴリズムが終了する反復回数ｋに達したら、そのときに更新された画素値ｘ（ｋ）を、それに対応する画素ｊごとに並べることで再構成画像（図４では「再構成画像Ｉ０」で表記）が得られる。
【００７５】
このようにして得られた再構成画像Ｉ０は、正則化によってストリークアーティファクトが抑制された画像となる。このステップＳ９は、この発明における（正則化処理が組み込まれた）再構成処理工程に相当する。
【００７６】
以上の図３および図４に示すステップＳ１〜Ｓ９により再構成画像Ｉ０が生成され、一連の処理が終了する。
【００７７】
本実施例１に係る画像再構成処理方法によれば、数学的最適化における正則化において、正則化処理が組み込まれた再構成処理を行う。すなわち、複数次元の直線状空間または曲線状空間（本実施例１ではストリークアーティファクトに沿った直線状空間を示す欠損領域Ｓｎｕｌｌ）に沿った、複数次元のデジタル画像（本実施例１では再構成画像）の画素値に関する重み付き加算値（本実施例１ではステップＳ３（投影関数の設定）で設定された投影関数ｒｍ（ｘ），ｒｎ（ｘ））の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて再構成処理（ステップＳ９：目的関数の最適化計算）を行って再構成処理後の複数次元のデジタル画像を作成する。「課題を解決するための手段」の欄でも述べたように、重み付き加算値を計算する場合には、重み係数を“１”にして単に加算することも含まれることに留意されたい。
【００７８】
上述した非特許文献１では画像の微分空間を正則化する。したがって、従来のような画像の微分空間を正則化する手法では、誤差の空間的相関性が強い箇所ではノイズでない画素値と認識され、誤差の空間的相関性が強い箇所でのノイズ自体を抑制することができない。それに対して、本実施例１に係る画像再構成処理方法では複数の重み付き加算値（投影関数ｒｍ（ｘ），ｒｎ（ｘ））の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像（再構成画像）を作成する。つまり、投影空間を正則化する処理（正則化処理）を再構成処理に組み込んで行う。したがって、従来のような画像の微分空間を正則化する手法と比較すると誤差の空間的相関性が強い箇所でのノイズ（本実施例１ではストリークアーティファクト）自体を抑制することができる。また、従来の正則化手法（例えば上述した特許文献１の手法）では計測データが投影データの場合には、正則化処理後のデータ（つまり復元されるデータ）も投影データであるが、本実施例１に係る画像再構成処理方法では正則化処理を再構成処理に組み込んで行うので、正則化処理後のデータ（つまり復元されるデータ）は再構成処理後の複数次元のデジタル画像（例えば断層画像や３次元再構成画像）である。また、複数の重み付き加算値（投影関数ｒｍ（ｘ），ｒｎ（ｘ））の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて再構成処理を行って再構成処理後の複数次元のデジタル画像を作成するので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
【００７９】
本実施例１では、計測データはカウント値（計数値）である。当該カウント値は、本実施例１のようにＰＥＴ装置のγ線検出器３で同時計数された計測データである。
【００８０】
本実施例１では、直線状空間に沿って重み付き加算値を求める場合には、直線状空間は線分であり、重み付き加算値に相当する投影関数ｒｍ（ｘ），ｒｎ（ｘ）は再構成画像の線分上の、複数次元のデジタル画像の画素値（例えば、逐次近似計算における推定途中の画像の画素値）の積分値である投影値（つまり線積分値）である。このように投影値を求めることにより、画像に生じる誤差の空間的相関性を正確に抽出することができて、ノイズをより効果的に低減させることができる。
【００８１】
本実施例１のように、投影値は、線分に平行に沿って加算した上述した（複数次元のデジタル画像の）画素値の積分値であるのが好ましい。「課題を解決するための手段」の欄でも述べたように、もし、図８（ａ）に示すようにストリークアーティファクトＳＡに直交した線分Ｌ１に平行に沿って加算（すなわちストリークアーティファクトＳＡに直交して加算）すると、主に投影値の低周波成分にストリークアーティファクトの影響が表れてしまう。つまり、投影値にはストリークアーティファクトの影響が現れにくい（すなわち、ストリークアーティファクトの影響を正確に抽出することができない）ので、ストリークアーティファクトが低減されにくくなる。なぜなら、投影空間の正則化では投影空間の高周波成分を抑制することで再構成画像のノイズ低減を図るためである。そこで、図８（ｂ）に示すようにストリークアーティファクトＳＡに沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２に平行に沿って加算することにより、得られた投影値は、ストリークアーティファクトの影響が現れた（すなわち、ストリークアーティファクトの影響を正確に抽出した）投影値となる。
【００８２】
本実施例１では、画像信号計測装置として、陽電子放射断層撮影装置（ＰＥＴ装置）に例に採って説明している。この場合には、上述した再構成処理工程（ステップＳ９：目的関数の最適化計算）では、上述した正則化関数Ｕ（ｘ）を用いて、当該陽電子放射断層撮影装置で得られた計測データに基づいて逐次近似法による再構成処理を行う。
【００８３】
本実施例１のように、画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、一部が開口して構成された検出器ユニット２Ａ，２Ｂを備え、当該検出器ユニット２Ａ，２Ｂで検出された計測データに基づいて再構成処理を行うのが有効である。すなわち、リング状の検出器ユニットとは相違して、一部が開口して構成された検出器ユニット２Ａ，２Ｂの場合には、開口部分（抜け部分）を通過する放射線は検出されない。よって、投影データの一部分が欠損することに起因して、再構成画像などに代表される再構成処理後の複数次元のデジタル画像に空間的相関性の強いノイズ（例えばストリークアーティファクト）が生じる。本実施例１に係る画像再構成処理方法を、（一部が開口して構成された）検出器ユニット２Ａ，２Ｂで検出された計測データに基づく再構成処理に適用することにより、空間的相関性の強いノイズを低減させることができる。
【００８４】
本実施例１のように、（一部が開口して構成された）検出器ユニット２Ａ，２Ｂを備えた場合において、直線状空間は、検出器ユニット２Ａ，２Ｂの開口領域（すなわち、γ線検出器３が存在しない開放領域）を通過する直線状空間を含むのが好ましい。投影データの一部分が欠損することに起因したストリークアーティファクトＳＡ（図８（ｂ）を参照）に沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２（図８（ｂ）を参照）に平行に沿って加算することにより、ストリークアーティファクトをより一層抑制することができる。
【００８５】
本実施例１のように、画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、陽電子放射断層撮影装置は、平均同時計数時間分解能が１ナノ秒（ｎｓ）以下（本実施例１では５００ピコ秒（ｐｓ）あるいは３００ピコ秒（ｐｓ）程度）であるγ線検出器３を備え、当該γ線検出器３で検出された計測データ（本実施例１ではカウント値）に基づいて再構成処理を行うのが好ましい。平均同時計数時間分解能が１ナノ秒を超えるとポジトロンの消滅位置や放射線による発光位置を正確に推定することができない。平均同時計数時間分解能が１ナノ秒以下である放射線検出器を用いることで、ポジトロンの消滅位置や放射線による発光位置を正確に推定することができ、より正確な断層画像を取得することができる。
【実施例２】
【００８６】
次に、図面を参照してこの発明の実施例２を説明する。
図９および図１０は、実施例２に係る画像処理のフローチャートである。本実施例２では、入力されたデジタル画像に対して画像処理を行う実施形態である。上述した実施例１では、投影空間の正則化関数Ｕ（ｘ）をＰＥＴ画像の画像再構成処理において利用して述べたが、本実施例２のように、再構成処理を除いた画像処理（画質改善処理）に利用することも可能である。
【００８７】
具体的には、入力されたデジタル画像（例えばＰＥＴ画像）にストリークアーティファクトなどのノイズが生じており、その画像のストリークアーティファクトを低減させる処理に転用することができる。そこで、本実施例２において、以下では、投影空間の正則化を利用したストリークアーティファクトの低減処理方法について説明する。すなわち、図９および図１０は、投影空間の正則化関数Ｕ（ｘ）（上記（３）式）を利用したストリークアーティファクトの低減処理方法のフローチャートでもある。
【００８８】
（ステップＴ１）投影空間の設定
図９に示すように、先ず、アーティファクトが生じている入力画像（入力されたデジタル画像）のサイズに基づいて投影空間の設定を行う。投影空間をＳとする。上述した実施例１でも述べたように、投影空間Ｓの設定とは、多次元投影空間の各座標軸の定義域および標本化間隔を規定することである。
【００８９】
上述した実施例１では、ＰＥＴ装置のジオメトリ情報などを用いてステップＳ２：投影空間の欠損領域の設定（図３を参照）を行ったが、本実施例２では投影空間の欠損領域Ｓｎｕｌｌの替わりに正則化領域Ｓｒｅｇを用いる。正則化領域Ｓｒｅｇについては詳しく後述する。
【００９０】
（ステップＴ２）投影関数の設定
次に、投影関数ｒｍ（・）（∀ｍ∈Ｓ）の設定を行う。投影関数ｒｍ（・）（∀ｍ∈Ｓ）の設定については、上述した実施例１と同様なので、詳しい説明を省略する。
【００９１】
（ステップＴ３）入力画像の投影データの計算
次に、投影関数ｒｍ（・）（∀ｍ∈Ｓ）を用いて、入力されたデジタル画像（以下、「入力画像Ｉ１０」と略記する）の投影データの計算を行う。入力画像の投影データをＰｒｏ０とする。
【００９２】
（ステップＴ４）投影空間の正則化領域の設定
次に、入力画像の投影データＰｒｏ０を解析することにより、投影空間の正則化領域Ｓｒｅｇの設定を行う。この際、入力画像Ｉ１０の生成条件に関する情報（図９では「入力画像の生成条件Ｉｎｆ」で表記）がある場合には、その情報を用いて正則化領域Ｓｒｅｇを設定してもよい。入力画像Ｉ１０がＰＥＴ画像であれば、ＰＥＴ画像のジオメトリ条件に関する情報（図３および図４の「ＰＥＴ画像のジオメトリ情報Ｇ」を参照）などが入力画像の生成条件Ｉｎｆとなり得る。
【００９３】
入力画像の生成条件Ｉｎｆが全くの未知であり、正則化領域Ｓｒｅｇを入力画像の投影データＰｒｏ０自体から判断する場合には、以下の方法が考えられる。
（方法１）ユーザが投影データを観察し、ノイズが有意に大きいと思われる領域を視覚的に（主観的に）判断し、その領域を正則化領域Ｓｒｅｇとする。
（方法２）ＣＰＵなどの計算機（例えば図１の演算回路５）によって、投影データの局所的な正規化標準偏差(normalized standard deviation: NSD)をノイズ指標として計算し、ＮＳＤ値が閾値以上である領域を正則化領域Ｓｒｅｇとする。ＮＳＤ計算の標本値の数（局所領域の大きさ）についてはユーザが指定する。閾値については、ユーザが指定するか、あるいはＮＳＤ値の分析解析（例えばクラスタリング解析）により自動算出する。
【００９４】
（ステップＴ５）距離関数の設定
次に、距離関数φｍｎ（ｕ，ｖ）（∀ｍ，∀ｎ∈Ｓｒｅｇ）の設定を行う。距離関数φｍｎ（ｕ，ｖ）（∀ｍ，∀ｎ∈Ｓｒｅｇ）の設定については、上述した実施例１と同様なので、詳しい説明を省略する。
【００９５】
（ステップＴ６）投影値の期待値の設定？
次に、図１０に示すように、上述した実施例１と同様に投影値の期待値ｐｍ（∀ｍ∈Ｓｒｅｇ）を設定するのか否か？の設定を行う。
【００９６】
（ステップＴ７）ｐｍ＝１（∀ｍ∈Ｓｒｅｇ）
上述した実施例１と同様に、設定しない場合（図１０では「Ｎｏ」で表記）にはｐｍ＝１（∀ｍ∈Ｓｒｅｇ）とする。
【００９７】
（ステップＴ８）ｐｍ（∀ｍ∈Ｓｒｅｇ）の設定
設定する場合（図１０では「Ｙｅｓ」で表記）には、以下の手順により投影値の期待値ｐｍを設定する。
【００９８】
（ステップＴ９）ノイズ低減処理
入力画像の投影データＰｒｏ０（図９を参照）に対してノイズ低減処理（例えば移動平均フィルタ）を施した後の投影データをｑｍとして、ノイズ低減処理後の投影データｑｍ（∀ｍ∈Ｓ）を取得する（図１０では「ノイズ低減された投影データＰｒｏ１」で表記）。そして、ノイズ低減処理後の投影データｑｍ（∀ｍ∈Ｓ）での正則化領域Ｓｒｅｇの投影データｑｍ（∀ｍ∈ｒｅｇ）を投影値の期待値ｐｍ（∀ｍ∈ｒｅｇ）とする。つまり、投影値の期待値ｐｍ（∀ｍ∈ｒｅｇ）＝ノイズ低減処理後の投影データｑｍ（∀ｍ∈ｒｅｇ）となる。
【００９９】
（ステップＴ１０）アーティファクト除去問題の目的関数の設定
次に、アーティファクト除去問題の目的関数Ｆ（ｘ）の設定を行う。ただし、ｘは出力画像（変数ベクトル）（図１０では「出力画像Ｉ１１」で表記）である。上記（３）式で設定された正則化関数Ｕ（ｘ）を用いて目的関数Ｆ（ｘ）を設定する。例えば、入力画像をｙとし、入力画像ｙと出力画像ｘとの距離を測るデータフィット関数をＤ（ｘ，ｙ）としたときに、データフィット関数Ｄ（ｘ，ｙ）と投影空間の正則化関数Ｕ（ｘ）とを用いて、Ｄ（ｘ，ｙ）＋βＵ（ｘ）を目的関数Ｆ（ｘ）とする。ここで、βは正の定数であり、上述した実施例１でも述べたように正則化の強さを調整するパラメータである。
【０１００】
データフィット関数Ｄ（ｘ，ｙ）の具体例としては、ベクトルｘとｙとの二乗距離｜｜ｘ−ｙ｜｜２が考えられる。ここで、入力画像ｘの分散を表す共分散行列Σが既知である場合には、重み付き二乗距離（二乗マハラノビス距離）（ｘ−ｙ）ＴΣ−１（ｘ−ｙ）を用いるのが好ましい。ただし、上付きの“−１”は逆行列、上付きの“Ｔ”はベクトルの転置を表す。
【０１０１】
（ステップＴ１１）目的関数の最適化計算
次に、目的関数Ｆ（ｘ）の最適化計算を行う。目的関数Ｆ（ｘ）の最適化計算については、上述した実施例１と同様なので、詳しい説明を省略する。
【０１０２】
以上の図９および図１０に示すステップＴ１〜Ｔ１１により出力画像Ｉ１１が生成され、一連の処理が終了する。
【０１０３】
本実施例２に係る画像処理方法によれば、上述した実施例１に係る画像再構成処理方法と同様に、数学的最適化における正則化において、正則化処理が組み込まれた画像処理を行う。「課題を解決するための手段」の欄でも述べたように、本実施例２での画像処理は、再構成処理以外の画像処理も含まれる。すなわち、複数次元の直線状空間または曲線状空間（本実施例２ではストリークアーティファクトに沿った直線状空間）に沿った、デジタル画像（例えばＰＥＴ画像）の画素値に関する重み付き加算値（本実施例２ではステップＴ２（投影関数の設定）で設定された投影関数ｒｍ（・）（∀ｍ∈Ｓ））の相互類似度を評価した正則化関数Ｕ（ｘ）を用いることにより入力されたデジタル画像に対して画像処理を行う。上述した実施例１に係る画像再構成処理方法と同様に、重み付き加算値を計算する場合には、重み係数を“１”にして単に加算することも含まれることに留意されたい。つまり、ステップＴ１１（目的関数の最適化計算）では、複数の重み付き加算値の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて正則化処理を行う。
【０１０４】
上述した特許文献１では実測された投影データに対して正則化処理を行う。それに対して、本実施例２に係る画像処理方法では当該直線状空間に沿った、デジタル画像の画素値に関する重み付き加算値の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて正則化処理を行う。つまり、未知の求めたいデジタル画像の投影空間が正則である（すなわち滑らかである）ことを直接的に利用することにより、デジタル画像の投影空間を正則化する。また、複数の重み付き加算値の相互類似度を評価した正則化関数Ｕ（ｘ）を用いて正則化処理を行うので、画像に生じる誤差の空間的相関性が強い場合においてもノイズを低減させることができる。
【０１０５】
なお、入力されたデジタル画像（入力画像）は、画像に出力される構造物が実際の形態で表されるものであれば良く、例えば、上述したようなＰＥＴ画像，ＳＰＥＣＴ画像，Ｘ線ＣＴ画像，ＭＲＩ画像，顕微鏡カメラ画像，デジタルカメラ画像，デジタル画像などであり、カラー画像も含まれる。
【０１０６】
本実施例２では、上述した実施例１に係る画像再構成処理方法と同様に、直線状空間は線分であり、重み付き加算値に相当する投影関数ｒｍ（・）は線分上の画素値の積分値である投影値（つまり線積分値）である。このように投影値を求めることにより、画像に生じる誤差の空間的相関性を正確に抽出することができて、ノイズをより効果的に低減させることができる。
【０１０７】
本実施例２では、上述した実施例１に係る画像再構成処理方法と同様に、投影値は、線分に平行に沿って加算した画素値の積分値であるのが好ましい。図８（ｂ）に示すようにストリークアーティファクトＳＡに沿って直線状空間を設定し、当該直線状空間からなる線分Ｌ２に平行に沿って加算することにより、得られた投影値は、ストリークアーティファクトの影響が現れた（すなわち、ストリークアーティファクトの影響を正確に抽出した）投影値となる。
【実施例３】
【０１０８】
次に、図面を参照してこの発明の実施例３を説明する。
図１１は、実施例３に係る曲線状空間を設定したときの模式図である。上述した実施例１、２では、直線状空間に沿った、重み付き加算値の相互類似度を評価した正則化関数を用いて画像再構成処理あるいは画像処理を行う場合について述べた。これに対して、本実施例３の図１１に示すように、曲線状空間に沿った、重み付き加算値の相互類似度を評価した正則化関数を用いて画像再構成処理あるいは画像処理を行うことも可能である。
【０１０９】
つまり、上述した実施例１、２では、画像に生じ得るストリークアーティファクト、または画像に既に生じているストリークアーティファクトの低減方法について述べた。本実施例３では、投影関数ｒｍ（・）の定義を変更することで、本発明の特徴部分である上記（３）式で設定された正則化関数Ｕ（ｘ）を用いて、ストリークアーティファクト（つまり直線状のアーティファクト）以外の曲線状のアーティファクトを低減させることも可能である。
【０１１０】
例えば、図１１（ａ）に示すように、同心円状のアーティファクト（リングアーティファクト）Ａが画像に生じ得る場合には、投影線を円にして、投影関数ｒｍ（・）によって円上の画素値の重み付き加算値（積分値）を計算すれば、上記（３）式の正則化関数Ｕ（ｘ）によって、同心円状のアーティファクト（リングアーティファクト）Ａを低減させることができる。
【０１１１】
また、例えば、図１１（ｂ）に示すように、波線状のアーティファクトＡが画像に生じ得る場合には、投影線を波線（に近似した関数）にして、投影関数ｒｍ（・）によって波線上の画素値の重み付き加算値（積分値）を計算すれば、上記（３）式の正則化関数Ｕ（ｘ）によって、波線状のアーティファクトＡを低減させることができる。
【０１１２】
この発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
【０１１３】
（１）上述した各実施例では、画像信号計測装置として、陽電子放射断層撮影装置（ＰＥＴ装置）を例に採って説明したが、被検体の計測データを取得する装置であれば、特に限定されない。単一光子放射断層撮影装置（ＳＰＥＣＴ装置）やＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）や核磁気共鳴断層撮影装置（ＭＲＩ装置）などに適用してもよい。特に、上述した実施例１のような画像再構成処理方法に適用する場合には、正則化関数を用いて、ＳＰＥＣＴ装置やＸ線ＣＴ装置やＭＲＩ装置などで得られた計測データに基づいて逐次近似法による再構成処理を行う。
【０１１４】
（２）上述した各実施例において撮影対象については、特に限定されない。特に、上述した実施例１の場合には、被検体の全身を撮影する装置や、被検体の頭部を撮影する装置や、被検体の乳房を撮影する装置に適用してもよいし、上述した実施例２の場合には、デジタル画像に対して画像処理を行えばよく、撮影対象については人体などの被検体に限定されない。
【０１１５】
（３）上述した各実施例では、図１のような部分リング型ＰＥＴ装置１であったが、人体（符号Ｍを参照）の乳房を撮影するマンモグラフィ装置のように、放射線検出器を対向配置した装置にも適用してもよい。図１の検出器ユニット２Ａ，２Ｂが、図１２に示すように、乳房検査部２Ｃに置き換わるのを除けば、図１と同じ構成を有する。なお、図１２の場合には、乳房検査部２Ｃは切り欠きとなっており、この切り欠きに脇で挟むことで乳房を検査する。また、γ線検出器３（図１２では図示省略）は、この切り欠きに合わせて乳房検査部２Ｃ内に複数に並設されている。
【０１１６】
（４）上述した各実施例では、ＤＯＩ検出器であったが、深さ方向を弁別しない放射線検出器に適用してもよい。
【０１１７】
（５）上述した実施例１では、計測データはカウント値（計数値）であったが、画像信号計測装置によって得られる計測データであれば、特に限定されない。カウント値などに代表される計数データ以外には、計測データとしては、投影データや時系列イベントデータや周波数空間データなどがある。投影データとしては、Ｘ線ＣＴ装置で得られた測定対象物（減弱係数分布）の投影データやＰＥＴ装置で得られた測定対象物（放射能分布）のサイノグラムなどがある。時系列イベントデータとしては、ＰＥＴ装置の放射線検出器で得られた検出イベント情報（検出器番号、検出時間、γ線のエネルギ等）を時系列で保存したリストモードデータなどがある。周波数空間データとしては、ＭＲＩ装置で得られた測定対象物（プロトン密度分布など）の周波数空間データなどがある。
【０１１８】
（６）上述した実施例１、２では、直線状空間に沿って重み付き加算値を求める場合には、直線状空間は線分であり、重み付き加算値は再構成画像の線分上の、複数次元のデジタル画像の画素値の積分値である投影値であったが、投影値については、線分に平行に沿って加算した（複数次元のデジタル画像の）画素値の積分値に限定されない。投影線の空間的な拡がりを考慮した投影値であってもよい。つまり、投影線上の厳密な積分値（線積分値）でなくてもよい。
【０１１９】
（７）上述した各実施例では、画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、一部が開口して構成された検出器ユニットを備えたＰＥＴ装置（図１では部分リング型ＰＥＴ装置１）であったが、必ずしも一部が開口して構成された検出器ユニットでなくとも、通常のフルリング型ＰＥＴ装置に適用してもよい。
【０１２０】
（８）上述した各実施例では、画像信号計測装置が陽電子放射断層撮影装置（ＰＥＴ装置）である場合において、平均同時計数時間分解能が１ナノ秒以下である放射線検出器（図１および図２ではγ線検出器３）を備えたが、必ずしも平均同時計数時間分解能が１ナノ秒以下である必要はない。特に、上述した実施例２の場合には、ＴＯＦ情報は必ずしも必要でないので、平均同時計数時間分解能が１ナノ秒を超える放射線検出器を備えたＰＥＴ装置に適用してもよい。
【０１２１】
（９）上述した実施例２では、入力されたデジタル画像はＰＥＴ装置で得られたＰＥＴ画像であったが、ＰＥＴ画像に限定されない。実施例２でも述べたように、画像に出力される構造物が実際の形態で表されるものであれば、ＰＥＴ画像以外にも、ＳＰＥＣＴ画像やＸ線ＣＴ画像やＭＲＩ画像や顕微鏡カメラ画像やデジタルカメラ画像やデジタル画像などにも適用することができる。また、グレースケール画像に限定されず、カラー画像にも適用することができる。
【０１２２】
（１０）上述した実施例１では、投影空間の正則化領域Ｓｎｕｌｌは欠損投影線のみからなる集合であったが、Ｓｎｕｌｌは実測投影線と欠損投影線とを含んだ集合（つまり、投影空間全体）でもよい。
【０１２３】
（１１）上述した実施例２では、投影空間の正則化領域Ｓｒｅｇは、上述した（方法１）または（方法２）により設定するとしたが、正則化領域Ｓｒｅｇは投影空間全体でもよい。
【産業上の利用可能性】
【０１２４】
以上のように、この発明は、陽電子放射断層撮影装置（ＰＥＴ装置）や単一光子放射断層撮影装置（ＳＰＥＣＴ装置）やＸ線コンピュータ断層撮影装置（Ｘ線ＣＴ装置）や核磁気共鳴断層撮影装置（ＭＲＩ装置）などの断層画像撮影装置の画像再構成技術や、デジタル画像撮影装置の画像処理技術などに適している。
【符号の説明】
【０１２５】
１  …  部分リング型ＰＥＴ装置
２Ａ，２Ｂ  …  検出器ユニット
３  …  γ線検出器
Ｕ（ｘ）  …  正則化関数
ｒｍ（・），ｒｍ（ｘ），ｒｎ（ｘ）  …  投影関数
Ｉ０  …  再構成画像
Ｉ１０  …  入力画像
Ｉ１１  …  出力画像
【図１】
2016018245.tif 000007
【図２】
2016018245.tif 000008
【図３】
2016018245.tif 000009
【図４】
2016018245.tif 000010
【図５】
2016018245.tif 000011
【図６】
2016018245.tif 000012
【図７】
2016018245.tif 000013
【図８】
2016018245.tif 000014
【図９】
2016018245.tif 000015
【図１０】
2016018245.tif 000016
【図１１】
2016018245.tif 000017
【図１２】
2016018245.tif 000018