(57)【要約】【課題】放射線断層撮影装置により収集されたリストデータに基づいて高性能にノイズ除去された断層画像を作成することができる装置および方法を提供する。【解決手段】画像処理装置１０は、放射線断層撮影装置２により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する装置であって、順投影計算部１４、逆投影計算部１５、画像更新計算部１６等を備える。画像更新計算部１６は、複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行って、当該処理後の画像guided(λ,Ｉ)を求める。画像更新計算部１６は、複数のフレームそれぞれについて、画像guided(λ,Ｉ)に基づく正則化項Ｕ(λ)を用いて、逆投影計算部１５による逆投影計算結果に基づいて画像更新計算を行う。【選択図】図３
2019120530.tif 000002
【特許請求の範囲】
【請求項１】
放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する装置であって、
前記リストデータを収集順に複数のフレームに分割し、前記リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成部と、
初期画像を作成する初期画像作成部と、
前記複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算部と、
前記複数のフレームそれぞれについて前記順投影計算部により作成されたサイノグラムと前記サイノグラム作成部により作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算部と、
前記複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、前記逆投影計算部による逆投影計算結果に基づいて画像更新計算を行う画像更新計算部と、
を備え、
前記順投影計算部による順投影計算、前記逆投影計算部による逆投影計算および前記画像更新計算部による画像更新計算を繰り返し行い、
前記順投影計算部は、その繰り返し処理の初回においては前記初期画像を前記入力画像とし、以降においては前記画像更新計算部による画像更新計算の後の画像を前記入力画像とする、
画像処理装置。
【請求項２】
前記画像更新計算部は、前記ガイディド・フィルタ処理の際のガイダンス画像として、前記リストデータのうち前記サイノグラム作成部において各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いる、
請求項１に記載の画像処理装置。
【請求項３】
前記画像更新計算部は、前記ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いる、
請求項１に記載の画像処理装置。
【請求項４】
被検体の断層画像を再構成するためのリストデータを収集する放射線断層撮影装置と、
前記放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する請求項１〜３の何れか１項に記載の画像処理装置と、
を備える放射線断層撮影システム。
【請求項５】
放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する方法であって、
前記リストデータを収集順に複数のフレームに分割し、前記リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成ステップと、
初期画像を作成する初期画像作成ステップと、
前記複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算ステップと、
前記複数のフレームそれぞれについて前記順投影計算ステップにおいて作成されたサイノグラムと前記サイノグラム作成ステップにおいて作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算ステップと、
前記複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、前記逆投影計算ステップにおける逆投影計算結果に基づいて画像更新計算を行う画像更新計算ステップと、
を備え、
前記順投影計算ステップ、前記逆投影計算ステップおよび前記画像更新計算ステップを繰り返し行い、
前記順投影計算ステップにおいて、その繰り返し処理の初回においては前記初期画像を前記入力画像とし、以降においては前記画像更新計算ステップにおける画像更新計算の後の画像を前記入力画像とする、
画像処理方法。
【請求項６】
前記画像更新計算ステップは、前記ガイディド・フィルタ処理の際のガイダンス画像として、前記リストデータのうち前記サイノグラム作成ステップにおいて各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いる、
請求項５に記載の画像処理方法。
【請求項７】
前記画像更新計算ステップは、前記ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いる、
請求項５に記載の画像処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、放射線断層撮影装置により収集されたリストデータに基づいて断層画像を作成する装置および方法に関するものである。
【背景技術】
【０００２】
被検体（生体）の断層画像を取得することができる放射線断層撮影装置として、ＰＥＴ（Positron  Emission  Tomography）装置およびＳＰＥＣＴ（Single  Photon  Emission  Computed  Tomography）装置が挙げられる。
【０００３】
ＰＥＴ装置は、被検体が置かれる測定空間の周囲に配列された多数の小型の放射線検出器を有する検出部を備えている。ＰＥＴ装置は、陽電子放出アイソトープ（ＲＩ線源）が投入された被検体内における電子・陽電子の対消滅に伴って発生するエネルギ５１１ｋｅＶの光子対を検出部により同時計数法で検出し、この同時計数情報を収集する。そして、この収集した多数の同時計数情報に基づいて、測定空間における光子対の発生頻度の空間分布（すなわち、ＲＩ線源の空間分布）を表す断層画像を再構成することができる。このとき、ＰＥＴ装置により収集された同時計数情報を時系列に並べたリストデータを収集順に複数のフレームに分割し、そのリストデータのうち各フレームに含まれるデータ群を用いて画像再構成処理を行うことで、複数のフレームの断層画像からなる動態ＰＥＴ画像を得ることができる。このＰＥＴ装置は核医学分野等で重要な役割を果たしており、これを用いて例えば生体機能や脳の高次機能の研究を行うことができる。
【０００４】
このようにして再構成された断層画像は、ノイズを多く含んでいるので、画像フィルタによるノイズ除去処理が必要である。ノイズ除去の為に用いられる画像フィルタとしては、ガウシアン・フィルタ（Gaussian  Filter）およびガイディド・フィルタ（Guided  Filter）が挙げられる。従来からガウシアン・フィルタが用いられている。これに対して、ガイディド・フィルタは、近年になって開発されたものであって、ガウシアン・フィルタと比べると画像中の濃淡の境界をよく保存することができるという特徴がある。
【０００５】
特許文献１および非特許文献１には、ガイディド・フィルタにより動態ＰＥＴ画像のノイズを除去する技術が記載されている。特許文献１および非特許文献１に記載された技術は、ガイディド・フィルタによるノイズ除去処理に際して、複数のフレームの断層画像からなる動態ＰＥＴ画像を積算して得られる画像をガイダンス画像として用いている。
【０００６】
非特許文献２〜４には、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理において、画像更新計算の際にノイズ除去を図る技術が記載されている。非特許文献２には、ＭＡＰ-ＥＭ（Maximum  a  posteriori  Expectation  Maximization）法が記載されている。また、非特許文献３，４には、ＭＡＰ-ＥＭ法の一種であるＭＲＰ-ＥＭ（Median  Root  Prior  Expectation  Maximization）法が記載されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】中国特許出願公開第１０３９５５８９９号明細書
【非特許文献】
【０００８】
【非特許文献１】Lijun Lu, et al, "DynamicPET Denoising Incorporating a Composite Image Guided Filter," IEEE NuclearScience Symposium and Medical Imaging Conference (NSS/MIC), 2014.
【非特許文献２】小林哲哉、「エミッションＣＴにおける統計的画像再構成の基礎と展開」、MEDICAL IMAGING TECHNOLOGY, Vol.31, No.1, pp.21-25, January 2013.
【非特許文献３】E. Karali, et al, "MRPapproaches to PET image reconstruction," 2011 4th International Conferenceon Biomedical Engineering and Informatics (BMEI), pp.399-403 (2011).
【非特許文献４】篠原広行、他、「ML-EM法とMaximum a Posterior (Map-EM)法」、断層映像研究会雑誌、第４０巻、第２号、pp.11-20 (2013).
【発明の概要】
【発明が解決しようとする課題】
【０００９】
上記のノイズ除去処理技術を含む従来のノイズ除去処理技術は、断層画像に含まれるノイズを或る程度は除去することができるものの、その性能は十分ではない。特に、動態ＰＥＴ画像の場合のように、断層画像の生成に用いるデータの数が少なくノイズが多い場合には、従来のノイズ除去処理の性能は不十分である。ＰＥＴ画像およびＳＰＥＣＴ画像に対して更なるノイズ除去性能の向上が望まれている。
【００１０】
本発明は、上記問題点を解消する為になされたものであり、放射線断層撮影装置により収集されたリストデータに基づいて高性能にノイズ除去された断層画像を作成することができる装置および方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明の画像処理装置は、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する装置であって、(1) リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成部と、(2) 初期画像を作成する初期画像作成部と、(3) 複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算部と、(4) 複数のフレームそれぞれについて順投影計算部により作成されたサイノグラムとサイノグラム作成部により作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算部と、(5) 複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、逆投影計算部による逆投影計算結果に基づいて画像更新計算を行う画像更新計算部と、を備える。そして、この画像処理装置は、順投影計算部による順投影計算、逆投影計算部による逆投影計算および画像更新計算部による画像更新計算を繰り返し行い、順投影計算部は、その繰り返し処理の初回においては初期画像を入力画像とし、以降においては画像更新計算部による画像更新計算の後の画像を入力画像とする。
【００１２】
本発明の画像処理装置において、画像更新計算部は、ガイディド・フィルタ処理の際のガイダンス画像として、リストデータのうちサイノグラム作成部において各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いるのが好適である。また、画像更新計算部は、ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いるのも好適である。
【００１３】
本発明の放射線断層撮影システムは、被検体の断層画像を再構成するためのリストデータを収集する放射線断層撮影装置と、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する上記の本発明の画像処理装置と、を備える。
【００１４】
本発明の画像処理方法は、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する方法であって、(1) リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成ステップと、(2) 初期画像を作成する初期画像作成ステップと、(3) 複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算ステップと、(4) 複数のフレームそれぞれについて順投影計算ステップにおいて作成されたサイノグラムとサイノグラム作成ステップにおいて作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算ステップと、(5) 複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、逆投影計算ステップにおける逆投影計算結果に基づいて画像更新計算を行う画像更新計算ステップと、を備える。そして、この画像処理方法は、順投影計算ステップ、逆投影計算ステップおよび画像更新計算ステップを繰り返し行い、順投影計算ステップにおいて、その繰り返し処理の初回においては初期画像を入力画像とし、以降においては画像更新計算ステップにおける画像更新計算の後の画像を入力画像とする。
【００１５】
本発明の画像処理方法において、画像更新計算ステップは、ガイディド・フィルタ処理の際のガイダンス画像として、リストデータのうちサイノグラム作成ステップにおいて各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いるのが好適である。また、画像更新計算ステップは、ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いるのも好適である。
【発明の効果】
【００１６】
本発明によれば、放射線断層撮影装置により収集されたリストデータに基づいて高性能にノイズ除去された断層画像を作成することができる。
【図面の簡単な説明】
【００１７】
【図１】図１は、ＲＩ線源が投入された被検体における同時計数情報の取得頻度の時間変化を示すグラフである。
【図２】図２は、ＭＡＰ-ＥＭ法の反復式を説明する図である。
【図３】図３は、本実施形態の放射線断層撮影システム１の構成を示す図である。
【図４】図４は、本実施形態の画像処理方法を説明するフローチャートである。
【図５】図５は、シミュレーションにおいで用いた数値ファントムの画像を示す図である。
【図６】図６は、シミュレーションにおいて用いた白質部（ＷＭ）および灰白質部（ＧＭ）それぞれのアクティビティの時間変化のモデルを示すグラフである。
【図７】図７は、シミュレーションにおいてガイダンス画像として用いたＭＲＩ画像を示す図である。
【図８】図８は、比較例１で再構成された断層画像を示す図である。
【図９】図９は、比較例２で再構成された断層画像を示す図である。
【図１０】図１０は、実施例１で再構成された断層画像を示す図である。
【図１１】図１１は、実施例２で再構成された断層画像を示す図である。
【図１２】図１２は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＭＳＥの時間変化を示すグラフである。
【図１３】図１３は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＰＳＮＲの時間変化を示すグラフである。
【図１４】図１４は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＳＳＩＭの時間変化を示すグラフである。
【発明を実施するための形態】
【００１８】
以下、添付図面を参照して、本発明を実施するための形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。本発明は、これらの例示に限定されるものではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。
【００１９】
本実施形態の画像処理装置および画像処理方法は、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する。特に、本実施形態の画像処理装置および画像処理方法は、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理において、画像更新計算の際に、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いてノイズ除去を図る。以下では、本実施形態の画像処理装置および画像処理方法の説明に先だって、ガイディド・フィルタ、ＭＬ-ＥＭ法、ＭＡＰ-ＥＭ法およびＭＲＰ-ＥＭ法それぞれの概要について説明する。なお、ＭＬ-ＥＭ（Maximum  Likelihood  Expectation  Maximization）法は、逐次近似的画像再構成処理の基本的なものである（非特許文献２，４を参照）。
【００２０】
初めに、ガイディド・フィルタの概要について説明する。ガイディド・フィルタは、下記（１）式で表される線形変換モデルを仮定して、入力画像Ｐおよびガイダンス画像Ｉを用いてノイズ除去処理を行い、ノイズ除去処理後の出力画像Ｑを作成する。Ｑｉは、出力画像Ｑにおける画素ｉの画素値である。Ｉｉは、ガイダンス画像Ｉにおける画素ｉの画素値である。ωｋは、入力画像Ｐにおいて画素ｋを中心とする複数の画素からなる小領域（パッチ領域）を表す。ａｋ，ｂｋは係数である。
【００２１】
【数１】
2019120530.tif 000003
【００２２】
誤差関数Ｅ（ａｋ，ｂｋ）は下記（２）式で表される。この誤差関数Ｅ（ａｋ，ｂｋ）を最小化する係数ａｋ，ｂｋを求める。εは、ａｋをペナルティとする正則化パラメータであり、下記（３）式で表される。σは、入力画像Ｐの全体におけるノイズの推定標準偏差である。αは定数である。
【００２３】
【数２】
2019120530.tif 000004
【００２４】
【数３】
2019120530.tif 000005
【００２５】
この（２）式を正則化最小二乗法（線形リッジ回帰）により解くと、下記（４）式および（５）式が得られる。<Ｐｋ> は、入力画像Ｐの小領域ωｋ内の平均画素値である。μｋは、ガイダンス画像Ｉの小領域ωｋ内の平均画素値である。σｋは、ガイダンス画像Ｉの小領域ωｋ内の画素値の分散値である。|ωｋ|は、小領域ωｋ内の画素数である。
【００２６】
【数４】
2019120530.tif 000006
【００２７】
【数５】
2019120530.tif 000007
【００２８】
ノイズの影響を低減する為に、下記（６）式および（７）式で表されるように、小領域ωｉ内で係数ａｋ，ｂｋの平均値<ａｋ>，<ｂｋ> を求める。そして、この <ａｋ>，<ｂｋ> を用いて、下記（８）式により、ノイズ除去処理後の出力画像Ｑを得ることができる。
【００２９】
【数６】
2019120530.tif 000008
【００３０】
【数７】
2019120530.tif 000009
【００３１】
【数８】
2019120530.tif 000010
【００３２】
一般に、ガイディド・フィルタでは、ガイダンス画像Ｉは入力画像Ｐと同じであってもよい。すなわち、動態ＰＥＴ画像がＭ個のフレームの断層画像からなるとして、これらＭ個のフレームのうち第ｍフレームの断層画像を入力画像Ｐｍおよびガイダンス画像Ｉｍの双方として用いて、ガイディド・フィルタにより、第ｍフレームのノイズ除去処理後の出力画像Ｑｍを作成してもよい。
【００３３】
特許文献１および非特許文献１に記載されたガイディド・フィルタによる動態ＰＥＴ画像に対するノイズ除去処理技術は、動態ＰＥＴ画像のＭ個のフレームの断層画像（入力画像Ｐｍ）の総和をガイダンス画像Ｉとして用いる（下記（９）式）。このガイダンス画像Ｉは、各入力画像Ｐｍよりノイズが低減されたものとなる。したがって、特許文献１および非特許文献１では、このガイダンス画像Ｉを用いることで、各入力画像Ｐｍに対して、より高性能のノイズ除去処理が可能であるとされている。
【００３４】
【数９】
2019120530.tif 000011
【００３５】
一般に、ＲＩ線源が投入された被検体内における電子・陽電子の対消滅の頻度（同時計数情報の取得頻度）は、図１に示されるように、ＲＩ線源投入直後から時間の経過とともに次第に高くなっていき、やがて或る時刻Ｔｐでピークに達し、そのピークの後は緩やかに次第に低くなっていく。図１は、ＲＩ線源が投入された被検体における同時計数情報の取得頻度の時間変化を示すグラフである。
【００３６】
このことから、ＲＩ線源投入直後から同時計数情報の取得の頻度が高い期間においては各フレームの期間は相対的に短く設定され、その後の同時計数情報の取得の頻度が緩やかに低くなっていく期間においては各フレームの期間は相対的に長く設定される。また、ＲＩ線源投入直後から同時計数情報の取得の頻度がピークに達するまでの期間（時刻Ｔｐ以前）と、そのピークの後の期間（時刻Ｔｐ以降）とでは、断層画像が大きく相違していることが知られている。
【００３７】
特許文献１および非特許文献１に記載されたノイズ除去処理技術は、各フレームの期間が一定でない点を考慮することなく、また、ピーク前後の断層画像が大きく相違する点をも考慮することなく、Ｍ個のフレームの断層画像（入力画像Ｐｍ）の単純な総和をガイダンス画像Ｉとして用いる。このことから、各フレームの断層画像に対するノイズ除去の性能には限界がある。
【００３８】
次に、ＭＬ-ＥＭ法の概要について説明する。投影データ（サイノグラム）ｙの観測モデルは、下記（１０）式のように、ＲＩ分布の線積分と統計ノイズとの和として表される。投影データ（サイノグラム）ｙは、放射線断層撮影装置の放射線検出器の対ごとに同時計数情報をヒストグラム化したものである。λは、測定空間における光子対の発生頻度の空間分布（すなわち、ＲＩ線源の空間分布）である。このλは、逐次近似的画像再構成処理により得たい断層画像である。Ｃは、線形変換の作用素である。Ｃλは、或る線分の線積分すなわち投影変換（順投影）を表す。ｎは、統計ノイズ（ポアソンノイズ）を表す。
【００３９】
【数１０】
2019120530.tif 000012
【００４０】
投影データはポアソンノイズによって統計的な変動を含むので、尤度とよばれる生起確率分布を求め、それを最大化するパラメータを投影データの推定値とする。このとき、投影データ推定値は、統計的な意味において最も確からしい。ポアソン分布に従う投影データｙの尤度Ｐ(ｙ|λ)は、下記（１１）式で表される。Ｓは、互いに独立した投影データｙの投影数である。下記（１２）式で示すパラメータは、推定した投影データである。
【００４１】
【数１１】
2019120530.tif 000013
【００４２】
【数１２】
2019120530.tif 000014
【００４３】
この尤度Ｐ(ｙ|λ)を最大化するλを解くことで、断層画像を求めることができる。尤度Ｐ(ｙ|λ)の除算を回避するとともに、最小化問題として捉えるために、−log Ｐ(ｙ|λ)を最小化する問題として下記（１３）式で定式化する。Ｌ(λ)は、-log Ｐ(ｙ|λ)から定数項を除いた関数である。これを、λが非負である制約条件の下で解く。
【００４４】
【数１３】
2019120530.tif 000015
【００４５】
Ｌ(λ)を最小化する問題は、ＥＭ法により、下記（１４）式の反復式として表される。ここで、ｊは断層画像の画素番号を表す。ｋは繰り返し回数を表す。この反復式は、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理を表している。すなわち、順投影計算は、（１４）式右辺の第２因子（総和を求める因子）の分母で記述されている。逆投影計算は、（１４）式右辺の第２因子で総和を求める演算として記述されている。また、画像更新計算は、（１４）式右辺において第１因子と第２因子との積として記述されている。
【００４６】
【数１４】
2019120530.tif 000016
【００４７】
ＭＬ-ＥＭ法では必ずしもノイズに頑健ではないことが知られている。そこで、断層画像の画像品質向上を目的として、ＭＡＰ（maximum  a  posteriori）推定に基づいた画像再構成法であるＭＡＰ-ＥＭ法が提案されている。
【００４８】
次に、ＭＡＰ-ＥＭ法の概要について説明する。ＭＡＰ-ＥＭ法は下記（１５）式で定式化される。ここで、Ｌ(λ)は、前述した対数尤度（（１３）式）である。Ｕ(λ)は、正則化項である。正則化項Ｕ(λ)は、ＲＩ分布の尤もらしさを評価する関数であり、ＲＩ分布に対する先験情報（事前情報ともいう。）を基に設計される。βは、正値をとる正則化パラメータである。これを、λが非負である制約条件の下で解く。
【００４９】
【数１５】
2019120530.tif 000017
【００５０】
Ｆ(λ)を最小化する問題は、ＯＳＬ（one  step  late）法により下記（１６）式の反復式として表される。図２は、ＭＡＰ-ＥＭ法の反復式（（１６）式）を説明する図である。この図では、反復式（（１６）式）において順投影計算、逆投影計算および画像更新計算それぞれの部分を矩形枠で囲っている。ＭＬ-ＥＭ法の反復式（（１４）式）と比較すると、ＭＡＰ-ＥＭ法の反復式（（１６）式）は、順投影計算および逆投影計算については同じであるが、画像更新計算については相違している。
【００５１】
【数１６】
2019120530.tif 000018
【００５２】
次に、ＭＲＰ-ＥＭ法の概要について説明する。ＭＲＰ-ＥＭ法は、ＭＡＰ-ＥＭ法において、下記（１７）式で表される正則化項Ｕ(λ)を用いる。この正則化項Ｕ(λ)の式の右辺にあるmed(λ)は、更新前の画像λに対してメディアン・フィルタ処理を行った後の画像を表す。メディアン・フィルタ処理は、注目画素およびこれの周囲にある画素の各画素値のうちの中央値を注目画素の画素値とする処理である。この正則化項Ｕ(λ)は、平均値および分散の双方がmed(λ)であるガウス関数として捉えることができる。ＭＲＰ-ＥＭ法の反復式は下記（１８）式で表される。
【００５３】
【数１７】
2019120530.tif 000019
【００５４】
【数１８】
2019120530.tif 000020
【００５５】
次に、本実施形態の画像処理方法について説明する。本実施形態の画像処理方法は、ＭＡＰ-ＥＭ法において、下記（１９）式で表される正則化項Ｕ(λ)を用いる。この正則化項Ｕ(λ)の式の右辺にあるguided(λ,Ｉ)は、更新前の画像λに対してガイダンス画像Ｉを用いてガイディド・フィルタ処理を行った後の画像を表す。この正則化項Ｕ(λ)は、平均値および分散の双方がguided(λ,Ｉ)であるガウス関数として捉えることができる。本実施形態の画像処理方法の反復式は下記（２０）式で表される。
【００５６】
【数１９】
2019120530.tif 000021
【００５７】
【数２０】
2019120530.tif 000022
【００５８】
図３は、本実施形態の放射線断層撮影システム１の構成を示す図である。放射線断層撮影システム１は、放射線断層撮影装置２および画像処理装置１０を備える。画像処理装置１０は、サイノグラム作成部１１、ガイダンス画像作成部１２、初期画像作成部１３、順投影計算部１４、逆投影計算部１５、画像更新計算部１６および記憶部１７を備える。画像処理装置１０としてコンピュータが用いられる。
【００５９】
放射線断層撮影装置２は、被検体の断層画像を再構成するためのリストデータを収集する装置である。放射線断層撮影装置２として、ＰＥＴ装置およびＳＰＥＣＴ装置が挙げられる。以下では、放射線断層撮影装置２がＰＥＴ装置であるとして説明をする。
【００６０】
放射線断層撮影装置２は、被検体が置かれる測定空間の周囲に配列された多数の小型の放射線検出器を有する検出部を備えている。放射線断層撮影装置２は、陽電子放出アイソトープ（ＲＩ線源）が投入された被検体内における電子・陽電子の対消滅に伴って発生するエネルギ５１１ｋｅＶの光子対を検出部により同時計数法で検出し、この同時計数情報を蓄積する。そして、放射線断層撮影装置２は、この蓄積した多数の同時計数情報を時系列に並べたものをリストデータとして画像処理装置１０へ出力する。
【００６１】
リストデータは、光子対を同時計数した１対の放射線検出器の識別情報および検出時刻情報を含む。リストデータは、さらに、各放射線検出器が検出した光子のエネルギ情報、および、１対の放射線検出器の検出時間差情報をも含んでいてもよい。画像処理装置１０は、放射線断層撮影装置２により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する。
【００６２】
サイノグラム作成部１１は、リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成する。ガイダンス画像作成部１２は、画像更新計算部１６による画像更新計算の際に用いるガイダンス画像Ｉｍを作成する。初期画像作成部１３は、初期画像を作成する。初期画像は、任意であるが、例えば、各画素値が正の一定値（例えば値１）である画像である。
【００６３】
順投影計算部１４は、複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する。逆投影計算部１５は、複数のフレームそれぞれについて順投影計算部１４により作成されたサイノグラムとサイノグラム作成部１１により作成されたサイノグラムとの比較に基づいて逆投影計算を行う。画像更新計算部１６は、複数のフレームそれぞれについて、ガイダンス画像Ｉｍを用いて更新前の画像（処理対象画像Ｐｍ）に対してガイディド・フィルタ処理を行って、当該処理後の画像guided(λ,Ｉ)を求める。そして、画像更新計算部１６は、複数のフレームそれぞれについて、画像guided(λ,Ｉ)に基づく正則化項Ｕ(λ)を用いて、逆投影計算部１５による逆投影計算結果に基づいて画像更新計算を行う。なお、Ｍ個のフレームのうちの第ｍフレームの処理対象画像をＰｍと記し、これに対応するガイダンス画像をＩｍと記す。
【００６４】
記憶部１７は、リストデータを記憶し、各フレームのサイノグラム、処理対象画像Ｐｍおよびガイダンス画像Ｉｍ等を記憶する。また、記憶部１７は、各計算の途中の画像を記憶する。
【００６５】
順投影計算部１４による順投影計算、逆投影計算部１５による逆投影計算および画像更新計算部１６による画像更新計算は、繰り返し回数が所定の反復回数に達するまで、繰り返し行われる。順投影計算部１４は、その繰り返し処理の初回においては、初期画像作成部１３により作成された初期画像を入力画像とする。順投影計算部１４は、以降においては、画像更新計算部１６による画像更新計算の後の画像を入力画像とする。繰り返し回数が所定の反復回数に達したときの画像更新計算部１６による画像更新計算後の画像が、最終的な断層画像として出力される。
【００６６】
図４は、本実施形態の画像処理方法を説明するフローチャートである。本実施形態の画像処理方法は、放射線断層撮影装置２により収集されたリストデータを取得するリストデータ取得ステップＳ１０と、サイノグラム作成部１１が行うサイノグラム作成ステップＳ１１と、ガイダンス画像作成部１２が行うガイダンス画像作成ステップＳ１２と、初期画像作成部１３が行う初期画像作成ステップＳ１３と、順投影計算部１４が行う順投影計算ステップＳ１４と、逆投影計算部１５が行う逆投影計算ステップＳ１５と、画像更新計算部１６が行う画像更新計算ステップＳ１６と、繰り返し回数ｋの値を１増する繰り返し回数更新ステップＳ１７と、繰り返し回数ｋが所定の反復回数に達したか否かを判定する判定ステップＳ１８と、を備える。
【００６７】
リストデータ取得ステップＳ１０では、放射線断層撮影装置２により収集されたリストデータを取得する。サイノグラム作成ステップＳ１１では、リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成する。ガイダンス画像作成ステップＳ１２では、画像更新計算ステップＳ１６における画像更新計算の際に用いるガイダンス画像Ｉｍを作成する。初期画像作成ステップＳ１３では、初期画像を作成する。
【００６８】
順投影計算ステップＳ１４では、複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する。逆投影計算ステップＳ１５では、複数のフレームそれぞれについて順投影計算ステップＳ１４において作成されたサイノグラムとサイノグラム作成ステップＳ１１により作成されたサイノグラムとの比較に基づいて逆投影計算を行う。画像更新計算ステップＳ１６では、複数のフレームそれぞれについて、ガイダンス画像Ｉｍを用いて更新前の画像（処理対象画像Ｐｍ）に対してガイディド・フィルタ処理を行って、当該処理後の画像guided(λ,Ｉ)を求める。そして、画像更新計算ステップＳ１６では、複数のフレームそれぞれについて、画像guided(λ,Ｉ)に基づく正則化項Ｕ(λ)を用いて、逆投影計算ステップＳ１５における逆投影計算結果に基づいて画像更新計算を行う。
【００６９】
繰り返し回数更新ステップＳ１７では、繰り返し回数ｋの値を１増する。判定ステップＳ１８では、繰り返し回数ｋが所定の反復回数に達したか否かを判定する。判定ステップＳ１８において繰り返し回数ｋが所定の反復回数に達したと判定されるまで、順投影計算ステップＳ１４、逆投影計算ステップＳ１５、画像更新計算ステップＳ１６および繰り返し回数更新ステップＳ１７を、繰り返し行う。
【００７０】
順投影計算ステップＳ１４において、その繰り返し処理の初回（ｋ＝０）においては、初期画像作成ステップＳ１３において作成された初期画像を入力画像とする。順投影計算ステップＳ１４において、以降（ｋ≧１）においては、画像更新計算ステップＳ１６における画像更新計算の後の画像を入力画像とする。
【００７１】
判定ステップＳ１８において繰り返し回数ｋが所定の反復回数に達したと判定されると、画像更新計算ステップＳ１６における画像更新計算後の画像が、最終的な断層画像として出力される。
【００７２】
画像更新計算部１６（画像更新計算ステップＳ１６）においてガイディド・フィルタ処理の際に用いられるガイダンス画像Ｉｍは、放射線断層撮影装置２により取得されたリストデータに基づいて再構成された画像であってもよい。第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、リストデータのうちサイノグラム作成部１１において第ｍフレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成された画像であるのが好適である。第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、第ｍフレームのデータ群の２倍以上のデータ数のデータ群を用いて作成されるのが好適である。第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、第ｍフレームのデータ群を含むデータ群を用いて作成されてもよいし、第ｍフレームのデータ群の前後のデータ群を用いて作成されてもよい。
【００７３】
第ｍフレームが図１において時刻Ｔｐ以前（血流依存領域）である場合には、第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、リストデータのうち時刻Ｔｐ以前のデータ群を用いて作成されるのが好適である。第ｍフレームが図１において時刻Ｔｐ以降（ＰＥＴリガンド動態依存領域）である場合には、第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、リストデータのうち時刻Ｔｐ以降のデータ群を用いて作成されるのが好適である。第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、リストデータの全体を用いて作成されてもよい。
【００７４】
また、ガイダンス画像作成部１２（ガイダンス画像作成ステップＳ１２）において、第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、ガイダンス画像Ｉｍの最大画素値が処理対象画像Ｐｍの最大画素値と等しくなるよう正規化を行って作成されるのが好適である。すなわち、下記（２１）式で表される正規化ガイダンス画像Ｉ'ｍを用いるのが好適である。MaxＰｍは、第ｍフレームの処理対象画像Ｐｍの最大画素値である。MaxＩｍは、第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍの最大画素値である。
【００７５】
【数２１】
2019120530.tif 000023
【００７６】
第ｍフレームの処理対象画像Ｐｍに対応するガイダンス画像Ｉｍは、放射線断層撮影装置により取得されたリストデータに基づいて作成されたものに限らず、他の装置により取得されたデータに基づいて作成されたものであってもよい。例えば、ガイダンス画像Ｉｍは、ＭＲＩ（Magnetic  Resonance  Imaging）装置により得られたＭＲＩ画像であってもよい。この場合、各フレームにおいて共通のＭＲＩ画像がガイダンス画像として用いられる。
【００７７】
次に、シミュレーション結果について説明する。本シミュレーションでは、先ず数値ファントム（図５）を用意した。この数値ファントムは、18Ｆ-ＦＤＧ（フルオロデオキシグルコース）を投入したサルの脳の断層画像を模擬したものである。この数値ファントムには、脳の白質部（white  matter、ＷＭ）および灰白質部（gray  matter、ＧＭ）が含まれている。図６は、シミュレーションにおいて用いた白質部（ＷＭ）および灰白質部（ＧＭ）それぞれのアクティビティの時間変化のモデルを示すグラフである。アクティビティの時間変化（Time-Activity  Curve、ＴＡＣ）は、一般に、ＲＩ線源投入直後から時間の経過とともに次第に高くなっていき、やがて或る時刻でピークに達し、そのピークの後は緩やかに次第に低くなっていく。総カウント数（同時計数情報の数）を１５００万とした。
【００７８】
この数値ファントムに基づいて各フレームのサイノグラムを作成した。第１〜第４のフレームの各期間を２０秒（小計８０秒）とし、第５〜第８のフレームの各期間を４０秒（小計１６０秒）とし、第９〜第１２のフレームの各期間を６０秒（小計２４０秒）とし、第１３〜第１６のフレームの各期間を１８０秒（小計７２０秒）とし、第１７〜第３０のフレームの各期間を３００秒（小計４２００秒）として、全体のリストデータの収集時間を５４００秒とした。各フレームのサイノグラムに対して、該フレームのカウント数に応じたポアソンノイズを与えて、ノイズ付加サイノグラムを作成した。このノイズ付加サイノグラムに基づいて画像再構成処理を行った。
【００７９】
比較例１では、ＭＬ-ＥＭ法により画像再構成処理を行った。比較例２では、ＭＲＰ-ＥＭ法により画像再構成処理を行った。実施例１，２では、本実施形態の画像処理方法により画像再構成処理を行った。実施例１では、画像更新計算におけるガイディド・フィルタ処理の際に、リストデータの全体を用いて作成された正規化ガイダンス画像を用いた。実施例２では、画像更新計算におけるガイディド・フィルタ処理の際に、ＭＲＩ画像を用いた。反復回数は４０回であった。
【００８０】
図７は、ガイダンス画像として用いたＭＲＩ画像を示す図である。図８は、比較例１で再構成された断層画像を示す図である。図９は、比較例２で再構成された断層画像を示す図である。図１０は、実施例１で再構成された断層画像を示す図である。図１１は、実施例２で再構成された断層画像を示す図である。図８〜図１２は、第２６フレーム（３９００〜４２００秒の期間）の断層画像を示す。これらの図の対比から分かるように、比較例１，２と比べて、実施例１，２では、画像中の濃淡の境界がよく保存された上で、ノイズがよく除去されている。
【００８１】
図１２は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＭＳＥの時間変化を示すグラフである。ＭＳＥ（Mean  Square  Error）は、真の画像（ここでは図５（ａ）の数値ファントム画像）との誤差を表しており、値が低いほど真の画像に近いことを意味する。図１３は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＰＳＮＲの時間変化を示すグラフである。ＰＳＮＲ（Peak  Signal  to  Noise  Ratio）は、画像の品質をデシベル（ｄＢ）で表したものであり、値が高いほど良好な画質であることを意味する。図１４は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＳＳＩＭの時間変化を示すグラフである。ＳＳＩＭ（Structural  Similarity  Index）は、画像の輝度およびコントラストならびに構造の変化を定量化する指標であり、値が高いほど良好な画質であることを意味する。ＭＳＥ、ＰＳＮＲおよびＳＳＩＭの何れの指標も、比較例２と比べて実施例１，２の方が性能が優れていることを示している。
【００８２】
以上のとおり、本実施形態では、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する。特に、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理において、画像更新計算の際に、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いてノイズ除去を図る。これにより、放射線断層撮影装置により収集されたリストデータに基づいて高性能にノイズ除去された断層画像を作成することができる。
【００８３】
なお、放射線断層撮影装置２は、上記の実施形態ではＰＥＴ装置であるとしたが、ＳＰＥＣＴ装置であってもよい。
【符号の説明】
【００８４】
１…放射線断層撮影システム、２…放射線断層撮影装置、１０…画像処理装置、１１…サイノグラム作成部、１２…ガイダンス画像作成部、１３…初期画像作成部、１４…順投影計算部、１５…逆投影計算部、１６…画像更新計算部、１７…記憶部。
【図１】
2019120530.tif 000024
【図２】
2019120530.tif 000025
【図３】
2019120530.tif 000026
【図４】
2019120530.tif 000027
【図５】
2019120530.tif 000028
【図６】
2019120530.tif 000029
【図７】
2019120530.tif 000030
【図８】
2019120530.tif 000031
【図９】
2019120530.tif 000032
【図１０】
2019120530.tif 000033
【図１１】
2019120530.tif 000034
【図１２】
2019120530.tif 000035
【図１３】
2019120530.tif 000036
【図１４】
2019120530.tif 000037