JP6986961B2 - 画像処理装置および画像処理方法 - Google Patents

Description

本発明は、放射線断層撮影装置により収集されたリストデータに基づいて断層画像を作成する装置および方法に関するものである。

被検体（生体）の断層画像を取得することができる放射線断層撮影装置として、ＰＥＴ（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）法は、逐次近似的画像再構成処理の基本的なものである（非特許文献２，４を参照）。

初めに、ガイディド・フィルタの概要について説明する。ガイディド・フィルタは、下記（１）式で表される線形変換モデルを仮定して、入力画像Ｐおよびガイダンス画像Ｉを用いてノイズ除去処理を行い、ノイズ除去処理後の出力画像Ｑを作成する。Ｑ_ｉは、出力画像Ｑにおける画素ｉの画素値である。Ｉ_ｉは、ガイダンス画像Ｉにおける画素ｉの画素値である。ω_ｋは、入力画像Ｐにおいて画素ｋを中心とする複数の画素からなる小領域（パッチ領域）を表す。ａ_ｋ，ｂ_ｋは係数である。

誤差関数Ｅ（ａ_ｋ，ｂ_ｋ）は下記（２）式で表される。この誤差関数Ｅ（ａ_ｋ，ｂ_ｋ）を最小化する係数ａ_ｋ，ｂ_ｋを求める。εは、ａ_ｋをペナルティとする正則化パラメータであり、下記（３）式で表される。σは、入力画像Ｐの全体におけるノイズの推定標準偏差である。αは定数である。

この（２）式を正則化最小二乗法（線形リッジ回帰）により解くと、下記（４）式および（５）式が得られる。<Ｐ_ｋ> は、入力画像Ｐの小領域ω_ｋ内の平均画素値である。μ_ｋは、ガイダンス画像Ｉの小領域ω_ｋ内の平均画素値である。σ_ｋは、ガイダンス画像Ｉの小領域ω_ｋ内の画素値の分散値である。|ω_ｋ|は、小領域ω_ｋ内の画素数である。

ノイズの影響を低減する為に、下記（６）式および（７）式で表されるように、小領域ω_ｉ内で係数ａ_ｋ，ｂ_ｋの平均値<ａ_ｋ>，<ｂ_ｋ> を求める。そして、この <ａ_ｋ>，<ｂ_ｋ> を用いて、下記（８）式により、ノイズ除去処理後の出力画像Ｑを得ることができる。

一般に、ガイディド・フィルタでは、ガイダンス画像Ｉは入力画像Ｐと同じであってもよい。すなわち、動態ＰＥＴ画像がＭ個のフレームの断層画像からなるとして、これらＭ個のフレームのうち第ｍフレームの断層画像を入力画像Ｐ_ｍおよびガイダンス画像Ｉ_ｍの双方として用いて、ガイディド・フィルタにより、第ｍフレームのノイズ除去処理後の出力画像Ｑ_ｍを作成してもよい。

特許文献１および非特許文献１に記載されたガイディド・フィルタによる動態ＰＥＴ画像に対するノイズ除去処理技術は、動態ＰＥＴ画像のＭ個のフレームの断層画像（入力画像Ｐ_ｍ）の総和をガイダンス画像Ｉとして用いる（下記（９）式）。このガイダンス画像Ｉは、各入力画像Ｐ_ｍよりノイズが低減されたものとなる。したがって、特許文献１および非特許文献１では、このガイダンス画像Ｉを用いることで、各入力画像Ｐ_ｍに対して、より高性能のノイズ除去処理が可能であるとされている。

一般に、ＲＩ線源が投入された被検体内における電子・陽電子の対消滅の頻度（同時計数情報の取得頻度）は、図１に示されるように、ＲＩ線源投入直後から時間の経過とともに次第に高くなっていき、やがて或る時刻Ｔｐでピークに達し、そのピークの後は緩やかに次第に低くなっていく。図１は、ＲＩ線源が投入された被検体における同時計数情報の取得頻度の時間変化を示すグラフである。

このことから、ＲＩ線源投入直後から同時計数情報の取得の頻度が高い期間においては各フレームの期間は相対的に短く設定され、その後の同時計数情報の取得の頻度が緩やかに低くなっていく期間においては各フレームの期間は相対的に長く設定される。また、ＲＩ線源投入直後から同時計数情報の取得の頻度がピークに達するまでの期間（時刻Ｔｐ以前）と、そのピークの後の期間（時刻Ｔｐ以降）とでは、断層画像が大きく相違していることが知られている。

特許文献１および非特許文献１に記載されたノイズ除去処理技術は、各フレームの期間が一定でない点を考慮することなく、また、ピーク前後の断層画像が大きく相違する点をも考慮することなく、Ｍ個のフレームの断層画像（入力画像Ｐ_ｍ）の単純な総和をガイダンス画像Ｉとして用いる。このことから、各フレームの断層画像に対するノイズ除去の性能には限界がある。

次に、ＭＬ-ＥＭ法の概要について説明する。投影データ（サイノグラム）ｙの観測モデルは、下記（１０）式のように、ＲＩ分布の線積分と統計ノイズとの和として表される。投影データ（サイノグラム）ｙは、放射線断層撮影装置の放射線検出器の対ごとに同時計数情報をヒストグラム化したものである。λは、測定空間における光子対の発生頻度の空間分布（すなわち、ＲＩ線源の空間分布）である。このλは、逐次近似的画像再構成処理により得たい断層画像である。Ｃは、線形変換の作用素である。Ｃλは、或る線分の線積分すなわち投影変換（順投影）を表す。ｎは、統計ノイズ（ポアソンノイズ）を表す。

投影データはポアソンノイズによって統計的な変動を含むので、尤度とよばれる生起確率分布を求め、それを最大化するパラメータを投影データの推定値とする。このとき、投影データ推定値は、統計的な意味において最も確からしい。ポアソン分布に従う投影データｙの尤度Ｐ(ｙ|λ)は、下記（１１）式で表される。Ｓは、互いに独立した投影データｙの投影数である。下記（１２）式で示すパラメータは、推定した投影データである。

この尤度Ｐ(ｙ|λ)を最大化するλを解くことで、断層画像を求めることができる。尤度Ｐ(ｙ|λ)の除算を回避するとともに、最小化問題として捉えるために、−log Ｐ(ｙ|λ)を最小化する問題として下記（１３）式で定式化する。Ｌ(λ)は、-log Ｐ(ｙ|λ)から定数項を除いた関数である。これを、λが非負である制約条件の下で解く。

Ｌ(λ)を最小化する問題は、ＥＭ法により、下記（１４）式の反復式として表される。ここで、ｊは断層画像の画素番号を表す。ｋは繰り返し回数を表す。この反復式は、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理を表している。すなわち、順投影計算は、（１４）式右辺の第２因子（総和を求める因子）の分母で記述されている。逆投影計算は、（１４）式右辺の第２因子で総和を求める演算として記述されている。また、画像更新計算は、（１４）式右辺において第１因子と第２因子との積として記述されている。

ＭＬ-ＥＭ法では必ずしもノイズに頑健ではないことが知られている。そこで、断層画像の画像品質向上を目的として、ＭＡＰ（maximum a posteriori）推定に基づいた画像再構成法であるＭＡＰ-ＥＭ法が提案されている。

次に、ＭＡＰ-ＥＭ法の概要について説明する。ＭＡＰ-ＥＭ法は下記（１５）式で定式化される。ここで、Ｌ(λ)は、前述した対数尤度（（１３）式）である。Ｕ(λ)は、正則化項である。正則化項Ｕ(λ)は、ＲＩ分布の尤もらしさを評価する関数であり、ＲＩ分布に対する先験情報（事前情報ともいう。）を基に設計される。βは、正値をとる正則化パラメータである。これを、λが非負である制約条件の下で解く。

Ｆ(λ)を最小化する問題は、ＯＳＬ（one step late）法により下記（１６）式の反復式として表される。図２は、ＭＡＰ-ＥＭ法の反復式（（１６）式）を説明する図である。この図では、反復式（（１６）式）において順投影計算、逆投影計算および画像更新計算それぞれの部分を矩形枠で囲っている。ＭＬ-ＥＭ法の反復式（（１４）式）と比較すると、ＭＡＰ-ＥＭ法の反復式（（１６）式）は、順投影計算および逆投影計算については同じであるが、画像更新計算については相違している。

次に、ＭＲＰ-ＥＭ法の概要について説明する。ＭＲＰ-ＥＭ法は、ＭＡＰ-ＥＭ法において、下記（１７）式で表される正則化項Ｕ(λ)を用いる。この正則化項Ｕ(λ)の式の右辺にあるmed(λ)は、更新前の画像λに対してメディアン・フィルタ処理を行った後の画像を表す。メディアン・フィルタ処理は、注目画素およびこれの周囲にある画素の各画素値のうちの中央値を注目画素の画素値とする処理である。この正則化項Ｕ(λ)は、平均値および分散の双方がmed(λ)であるガウス関数として捉えることができる。ＭＲＰ-ＥＭ法の反復式は下記（１８）式で表される。

次に、本実施形態の画像処理方法について説明する。本実施形態の画像処理方法は、ＭＡＰ-ＥＭ法において、下記（１９）式で表される正則化項Ｕ(λ)を用いる。この正則化項Ｕ(λ)の式の右辺にあるguided(λ,Ｉ)は、更新前の画像λに対してガイダンス画像Ｉを用いてガイディド・フィルタ処理を行った後の画像を表す。この正則化項Ｕ(λ)は、平均値および分散の双方がguided(λ,Ｉ)であるガウス関数として捉えることができる。本実施形態の画像処理方法の反復式は下記（２０）式で表される。

図３は、本実施形態の放射線断層撮影システム１の構成を示す図である。放射線断層撮影システム１は、放射線断層撮影装置２および画像処理装置１０を備える。画像処理装置１０は、サイノグラム作成部１１、ガイダンス画像作成部１２、初期画像作成部１３、順投影計算部１４、逆投影計算部１５、画像更新計算部１６および記憶部１７を備える。画像処理装置１０としてコンピュータが用いられる。

放射線断層撮影装置２は、被検体の断層画像を再構成するためのリストデータを収集する装置である。放射線断層撮影装置２として、ＰＥＴ装置およびＳＰＥＣＴ装置が挙げられる。以下では、放射線断層撮影装置２がＰＥＴ装置であるとして説明をする。

放射線断層撮影装置２は、被検体が置かれる測定空間の周囲に配列された多数の小型の放射線検出器を有する検出部を備えている。放射線断層撮影装置２は、陽電子放出アイソトープ（ＲＩ線源）が投入された被検体内における電子・陽電子の対消滅に伴って発生するエネルギ５１１ｋｅＶの光子対を検出部により同時計数法で検出し、この同時計数情報を蓄積する。そして、放射線断層撮影装置２は、この蓄積した多数の同時計数情報を時系列に並べたものをリストデータとして画像処理装置１０へ出力する。

リストデータは、光子対を同時計数した１対の放射線検出器の識別情報および検出時刻情報を含む。リストデータは、さらに、各放射線検出器が検出した光子のエネルギ情報、および、１対の放射線検出器の検出時間差情報をも含んでいてもよい。画像処理装置１０は、放射線断層撮影装置２により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する。

サイノグラム作成部１１は、リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成する。ガイダンス画像作成部１２は、画像更新計算部１６による画像更新計算の際に用いるガイダンス画像Ｉ_ｍを作成する。初期画像作成部１３は、初期画像を作成する。初期画像は、任意であるが、例えば、各画素値が正の一定値（例えば値１）である画像である。

順投影計算部１４は、複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する。逆投影計算部１５は、複数のフレームそれぞれについて順投影計算部１４により作成されたサイノグラムとサイノグラム作成部１１により作成されたサイノグラムとの比較に基づいて逆投影計算を行う。画像更新計算部１６は、複数のフレームそれぞれについて、ガイダンス画像Ｉ_ｍを用いて更新前の画像（処理対象画像Ｐ_ｍ）に対してガイディド・フィルタ処理を行って、当該処理後の画像guided(λ,Ｉ)を求める。そして、画像更新計算部１６は、複数のフレームそれぞれについて、画像guided(λ,Ｉ)に基づく正則化項Ｕ(λ)を用いて、逆投影計算部１５による逆投影計算結果に基づいて画像更新計算を行う。なお、Ｍ個のフレームのうちの第ｍフレームの処理対象画像をＰ_ｍと記し、これに対応するガイダンス画像をＩ_ｍと記す。

記憶部１７は、リストデータを記憶し、各フレームのサイノグラム、処理対象画像Ｐ_ｍおよびガイダンス画像Ｉ_ｍ等を記憶する。また、記憶部１７は、各計算の途中の画像を記憶する。

順投影計算部１４による順投影計算、逆投影計算部１５による逆投影計算および画像更新計算部１６による画像更新計算は、繰り返し回数が所定の反復回数に達するまで、繰り返し行われる。順投影計算部１４は、その繰り返し処理の初回においては、初期画像作成部１３により作成された初期画像を入力画像とする。順投影計算部１４は、以降においては、画像更新計算部１６による画像更新計算の後の画像を入力画像とする。繰り返し回数が所定の反復回数に達したときの画像更新計算部１６による画像更新計算後の画像が、最終的な断層画像として出力される。

図４は、本実施形態の画像処理方法を説明するフローチャートである。本実施形態の画像処理方法は、放射線断層撮影装置２により収集されたリストデータを取得するリストデータ取得ステップＳ１０と、サイノグラム作成部１１が行うサイノグラム作成ステップＳ１１と、ガイダンス画像作成部１２が行うガイダンス画像作成ステップＳ１２と、初期画像作成部１３が行う初期画像作成ステップＳ１３と、順投影計算部１４が行う順投影計算ステップＳ１４と、逆投影計算部１５が行う逆投影計算ステップＳ１５と、画像更新計算部１６が行う画像更新計算ステップＳ１６と、繰り返し回数ｋの値を１増する繰り返し回数更新ステップＳ１７と、繰り返し回数ｋが所定の反復回数に達したか否かを判定する判定ステップＳ１８と、を備える。

リストデータ取得ステップＳ１０では、放射線断層撮影装置２により収集されたリストデータを取得する。サイノグラム作成ステップＳ１１では、リストデータを収集順に複数のフレームに分割し、リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成する。ガイダンス画像作成ステップＳ１２では、画像更新計算ステップＳ１６における画像更新計算の際に用いるガイダンス画像Ｉ_ｍを作成する。初期画像作成ステップＳ１３では、初期画像を作成する。

順投影計算ステップＳ１４では、複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する。逆投影計算ステップＳ１５では、複数のフレームそれぞれについて順投影計算ステップＳ１４において作成されたサイノグラムとサイノグラム作成ステップＳ１１により作成されたサイノグラムとの比較に基づいて逆投影計算を行う。画像更新計算ステップＳ１６では、複数のフレームそれぞれについて、ガイダンス画像Ｉ_ｍを用いて更新前の画像（処理対象画像Ｐ_ｍ）に対してガイディド・フィルタ処理を行って、当該処理後の画像guided(λ,Ｉ)を求める。そして、画像更新計算ステップＳ１６では、複数のフレームそれぞれについて、画像guided(λ,Ｉ)に基づく正則化項Ｕ(λ)を用いて、逆投影計算ステップＳ１５における逆投影計算結果に基づいて画像更新計算を行う。

繰り返し回数更新ステップＳ１７では、繰り返し回数ｋの値を１増する。判定ステップＳ１８では、繰り返し回数ｋが所定の反復回数に達したか否かを判定する。判定ステップＳ１８において繰り返し回数ｋが所定の反復回数に達したと判定されるまで、順投影計算ステップＳ１４、逆投影計算ステップＳ１５、画像更新計算ステップＳ１６および繰り返し回数更新ステップＳ１７を、繰り返し行う。

順投影計算ステップＳ１４において、その繰り返し処理の初回（ｋ＝０）においては、初期画像作成ステップＳ１３において作成された初期画像を入力画像とする。順投影計算ステップＳ１４において、以降（ｋ≧１）においては、画像更新計算ステップＳ１６における画像更新計算の後の画像を入力画像とする。

判定ステップＳ１８において繰り返し回数ｋが所定の反復回数に達したと判定されると、画像更新計算ステップＳ１６における画像更新計算後の画像が、最終的な断層画像として出力される。

画像更新計算部１６（画像更新計算ステップＳ１６）においてガイディド・フィルタ処理の際に用いられるガイダンス画像Ｉ_ｍは、放射線断層撮影装置２により取得されたリストデータに基づいて再構成された画像であってもよい。第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、リストデータのうちサイノグラム作成部１１において第ｍフレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成された画像であるのが好適である。第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、第ｍフレームのデータ群の２倍以上のデータ数のデータ群を用いて作成されるのが好適である。第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、第ｍフレームのデータ群を含むデータ群を用いて作成されてもよいし、第ｍフレームのデータ群の前後のデータ群を用いて作成されてもよい。

第ｍフレームが図１において時刻Ｔｐ以前（血流依存領域）である場合には、第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、リストデータのうち時刻Ｔｐ以前のデータ群を用いて作成されるのが好適である。第ｍフレームが図１において時刻Ｔｐ以降（ＰＥＴリガンド動態依存領域）である場合には、第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、リストデータのうち時刻Ｔｐ以降のデータ群を用いて作成されるのが好適である。第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、リストデータの全体を用いて作成されてもよい。

また、ガイダンス画像作成部１２（ガイダンス画像作成ステップＳ１２）において、第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、ガイダンス画像Ｉ_ｍの最大画素値が処理対象画像Ｐ_ｍの最大画素値と等しくなるよう正規化を行って作成されるのが好適である。すなわち、下記（２１）式で表される正規化ガイダンス画像Ｉ'_ｍを用いるのが好適である。MaxＰ_ｍは、第ｍフレームの処理対象画像Ｐ_ｍの最大画素値である。MaxＩ_ｍは、第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍの最大画素値である。

第ｍフレームの処理対象画像Ｐ_ｍに対応するガイダンス画像Ｉ_ｍは、放射線断層撮影装置により取得されたリストデータに基づいて作成されたものに限らず、他の装置により取得されたデータに基づいて作成されたものであってもよい。例えば、ガイダンス画像Ｉ_ｍは、ＭＲＩ（Magnetic Resonance Imaging）装置により得られたＭＲＩ画像であってもよい。この場合、各フレームにおいて共通のＭＲＩ画像がガイダンス画像として用いられる。

次に、シミュレーション結果について説明する。本シミュレーションでは、先ず数値ファントム（図５）を用意した。この数値ファントムは、¹⁸Ｆ-ＦＤＧ（フルオロデオキシグルコース）を投入したサルの脳の断層画像を模擬したものである。この数値ファントムには、脳の白質部（white matter、ＷＭ）および灰白質部（gray matter、ＧＭ）が含まれている。図６は、シミュレーションにおいて用いた白質部（ＷＭ）および灰白質部（ＧＭ）それぞれのアクティビティの時間変化のモデルを示すグラフである。アクティビティの時間変化（Time-Activity Curve、ＴＡＣ）は、一般に、ＲＩ線源投入直後から時間の経過とともに次第に高くなっていき、やがて或る時刻でピークに達し、そのピークの後は緩やかに次第に低くなっていく。総カウント数（同時計数情報の数）を１５００万とした。

この数値ファントムに基づいて各フレームのサイノグラムを作成した。第１〜第４のフレームの各期間を２０秒（小計８０秒）とし、第５〜第８のフレームの各期間を４０秒（小計１６０秒）とし、第９〜第１２のフレームの各期間を６０秒（小計２４０秒）とし、第１３〜第１６のフレームの各期間を１８０秒（小計７２０秒）とし、第１７〜第３０のフレームの各期間を３００秒（小計４２００秒）として、全体のリストデータの収集時間を５４００秒とした。各フレームのサイノグラムに対して、該フレームのカウント数に応じたポアソンノイズを与えて、ノイズ付加サイノグラムを作成した。このノイズ付加サイノグラムに基づいて画像再構成処理を行った。

比較例１では、ＭＬ-ＥＭ法により画像再構成処理を行った。比較例２では、ＭＲＰ-ＥＭ法により画像再構成処理を行った。実施例１，２では、本実施形態の画像処理方法により画像再構成処理を行った。実施例１では、画像更新計算におけるガイディド・フィルタ処理の際に、リストデータの全体を用いて作成された正規化ガイダンス画像を用いた。実施例２では、画像更新計算におけるガイディド・フィルタ処理の際に、ＭＲＩ画像を用いた。反復回数は４０回であった。

図７は、ガイダンス画像として用いたＭＲＩ画像を示す図である。図８は、比較例１で再構成された断層画像を示す図である。図９は、比較例２で再構成された断層画像を示す図である。図１０は、実施例１で再構成された断層画像を示す図である。図１１は、実施例２で再構成された断層画像を示す図である。図８〜図１２は、第２６フレーム（３９００〜４２００秒の期間）の断層画像を示す。これらの図の対比から分かるように、比較例１，２と比べて、実施例１，２では、画像中の濃淡の境界がよく保存された上で、ノイズがよく除去されている。

図１２は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＭＳＥの時間変化を示すグラフである。ＭＳＥ（Mean Square Error）は、真の画像（ここでは図５（ａ）の数値ファントム画像）との誤差を表しており、値が低いほど真の画像に近いことを意味する。図１３は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＰＳＮＲの時間変化を示すグラフである。ＰＳＮＲ（Peak Signal to Noise Ratio）は、画像の品質をデシベル（ｄＢ）で表したものであり、値が高いほど良好な画質であることを意味する。図１４は、シミュレーションにおいて得られた比較例２、実施例１および実施例２それぞれの場合の再構成画像のＳＳＩＭの時間変化を示すグラフである。ＳＳＩＭ（Structural Similarity Index）は、画像の輝度およびコントラストならびに構造の変化を定量化する指標であり、値が高いほど良好な画質であることを意味する。ＭＳＥ、ＰＳＮＲおよびＳＳＩＭの何れの指標も、比較例２と比べて実施例１，２の方が性能が優れていることを示している。

以上のとおり、本実施形態では、放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する。特に、順投影計算、逆投影計算および画像更新計算を繰り返し行って断層画像を作成する逐次近似的画像再構成処理において、画像更新計算の際に、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いてノイズ除去を図る。これにより、放射線断層撮影装置により収集されたリストデータに基づいて高性能にノイズ除去された断層画像を作成することができる。

なお、放射線断層撮影装置２は、上記の実施形態ではＰＥＴ装置であるとしたが、ＳＰＥＣＴ装置であってもよい。

１…放射線断層撮影システム、２…放射線断層撮影装置、１０…画像処理装置、１１…サイノグラム作成部、１２…ガイダンス画像作成部、１３…初期画像作成部、１４…順投影計算部、１５…逆投影計算部、１６…画像更新計算部、１７…記憶部。

Claims (7)

1. 放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する装置であって、
   前記リストデータを収集順に複数のフレームに分割し、前記リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成部と、
   初期画像を作成する初期画像作成部と、
   前記複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算部と、
   前記複数のフレームそれぞれについて前記順投影計算部により作成されたサイノグラムと前記サイノグラム作成部により作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算部と、
   前記複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、前記逆投影計算部による逆投影計算結果に基づいて画像更新計算を行う画像更新計算部と、
   を備え、
   前記順投影計算部による順投影計算、前記逆投影計算部による逆投影計算および前記画像更新計算部による画像更新計算を繰り返し行い、
   前記順投影計算部は、その繰り返し処理の初回においては前記初期画像を前記入力画像とし、以降においては前記画像更新計算部による画像更新計算の後の画像を前記入力画像とする、
   画像処理装置。
2. 前記画像更新計算部は、前記ガイディド・フィルタ処理の際のガイダンス画像として、前記リストデータのうち前記サイノグラム作成部において各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いる、
   請求項１に記載の画像処理装置。
3. 前記画像更新計算部は、前記ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いる、
   請求項１に記載の画像処理装置。
4. 被検体の断層画像を再構成するためのリストデータを収集する放射線断層撮影装置と、
   前記放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する請求項１〜３の何れか１項に記載の画像処理装置と、
   を備える放射線断層撮影システム。
5. 放射線断層撮影装置により収集されたリストデータに基づいてＭＡＰ-ＥＭ法により画像再構成を行って断層画像を作成する方法であって、
   前記リストデータを収集順に複数のフレームに分割し、前記リストデータのうち各フレームに含まれるデータ群を用いて、該フレームのサイノグラムを作成するサイノグラム作成ステップと、
   初期画像を作成する初期画像作成ステップと、
   前記複数のフレームそれぞれについて入力画像に基づいて順投影計算を行ってサイノグラムを作成する順投影計算ステップと、
   前記複数のフレームそれぞれについて前記順投影計算ステップにおいて作成されたサイノグラムと前記サイノグラム作成ステップにおいて作成されたサイノグラムとの比較に基づいて逆投影計算を行う逆投影計算ステップと、
   前記複数のフレームそれぞれについて、更新前の画像に対してガイディド・フィルタ処理を行った後の画像に基づく正則化項を用いて、前記逆投影計算ステップにおける逆投影計算結果に基づいて画像更新計算を行う画像更新計算ステップと、
   を備え、
   前記順投影計算ステップ、前記逆投影計算ステップおよび前記画像更新計算ステップを繰り返し行い、
   前記順投影計算ステップにおいて、その繰り返し処理の初回においては前記初期画像を前記入力画像とし、以降においては前記画像更新計算ステップにおける画像更新計算の後の画像を前記入力画像とする、
   画像処理方法。
6. 前記画像更新計算ステップは、前記ガイディド・フィルタ処理の際のガイダンス画像として、前記リストデータのうち前記サイノグラム作成ステップにおいて各フレームのサイノグラムを作成する際に用いたデータ群より多いデータ数のデータ群を用いて再構成した画像を用いる、
   請求項５に記載の画像処理方法。
7. 前記画像更新計算ステップは、前記ガイディド・フィルタ処理の際のガイダンス画像としてＭＲＩ画像を用いる、
   請求項５に記載の画像処理方法。

Images