知財求人 - 知財ポータルサイト「IP Force」
特許6761610吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6761610
(24)【登録日】2020年9月9日
(45)【発行日】2020年9月30日
(54)【発明の名称】吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置
(51)【国際特許分類】
G01T 1/161 20060101AFI20200917BHJP
【ＦＩ】
G01T1/161 A
G01T1/161 C
【請求項の数】11
【全頁数】28
(21)【出願番号】特願2019-521549(P2019-521549)
(86)(22)【出願日】2017年5月29日
(86)【国際出願番号】JP2017019944
(87)【国際公開番号】WO2018220686
(87)【国際公開日】20181206
【審査請求日】2019年6月18日
(73)【特許権者】
【識別番号】000001993
【氏名又は名称】株式会社島津製作所
(74)【代理人】
【識別番号】100093056
【弁理士】
【氏名又は名称】杉谷 勉
(74)【代理人】
【識別番号】100142930
【弁理士】
【氏名又は名称】戸高 弘幸
(74)【代理人】
【識別番号】100175020
【弁理士】
【氏名又は名称】杉谷 知彦
(74)【代理人】
【識別番号】100180596
【弁理士】
【氏名又は名称】栗原 要
(74)【代理人】
【識別番号】100195349
【弁理士】
【氏名又は名称】青野 信喜
(72)【発明者】
【氏名】小林 哲哉
【審査官】
右田 純生
(56)【参考文献】
【文献】
国際公開第２０１６／１９７１２７（ＷＯ，Ａ１）
【文献】
米国特許出願公開第２０１５／０１９３９２７（ＵＳ，Ａ１）
【文献】
特表２０１５−５１９５５５（ＪＰ，Ａ）
【文献】
特開２００７−１３２８００（ＪＰ，Ａ）
【文献】
国際公開第２０１５／０６８２０４（ＷＯ，Ａ１）
【文献】
特表２０１６−５２５６７２（ＪＰ，Ａ）
【文献】
特開２０１５−１４５８２８（ＪＰ，Ａ）
【文献】
米国特許出願公開第２０１５／００５７５３５（ＵＳ，Ａ１）
【文献】
国際公開第２０１８／１６３３６２（ＷＯ，Ａ１）
【文献】
LI, Yusheng et al.，"Transmission-less attenuation estimation from time-of-flight PET histo-images using consistensy equations"，Physics in Medicine and Biology，２０１５年 ８月２１日，Volume 60, Number 16，pp.6563-6583
【文献】
LI, Quanzheng et al.，"Joint estimation of activity image and attenuation sinogram using time-of-flight positron emission tomography data consistency condition filtering"，Journal of Medical Imaging，２０１７年 ４月２６日，Volume 4, Issue 2，pp.023502-1〜023502-7
(58)【調査した分野】（Int.Cl.，ＤＢ名）
Ｇ０１Ｔ １／１６１− １／１６６
Ａ６１Ｂ ６／０３
Ｇ０１Ｎ ２３／０４６
Ｇ０１Ｔ １／２９
(57)【特許請求の範囲】
【請求項1】
消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データから吸収係数画像を推定する方法であって、
μ’を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、前記計測データに関する評価関数の最適化に基づいて、前記画像μ’を計算する再構成計算工程と、
前記計測データに基づいて投影データ空間における被検体マスクデータである被検体マスク投影データを算出するマスク算出工程と、
μｏｆｆを不均一なオフセット画像としたときに、オフセット画像μｏｆｆの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、前記オフセット画像μｏｆｆを推定するオフセット推定工程と、
Ωを既知の吸収係数値で近似可能な領域としたときに、前記計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、前記領域Ωを少なくとも１つ以上抽出する参照領域抽出工程と、
αを係数としたときに、前記領域Ωにおける前記画像μ’の値と既知の吸収係数値との誤差を減少させる前記係数αを算出する係数算出工程と、
前記画像μ’の値に、前記オフセット画像μｏｆｆを前記係数α倍したα×μｏｆｆを加算して得られた値を吸収係数値として補正する吸収係数値補正工程と
を備える、
吸収係数画像推定方法。
μ’を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、前記計測データに関する評価関数の最適化に基づいて、前記画像μ’を計算する再構成計算工程と、
前記計測データに基づいて投影データ空間における被検体マスクデータである被検体マスク投影データを算出するマスク算出工程と、
μｏｆｆを不均一なオフセット画像としたときに、オフセット画像μｏｆｆの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、前記オフセット画像μｏｆｆを推定するオフセット推定工程と、
Ωを既知の吸収係数値で近似可能な領域としたときに、前記計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、前記領域Ωを少なくとも１つ以上抽出する参照領域抽出工程と、
αを係数としたときに、前記領域Ωにおける前記画像μ’の値と既知の吸収係数値との誤差を減少させる前記係数αを算出する係数算出工程と、
前記画像μ’の値に、前記オフセット画像μｏｆｆを前記係数α倍したα×μｏｆｆを加算して得られた値を吸収係数値として補正する吸収係数値補正工程と
を備える、
吸収係数画像推定方法。
【請求項2】
請求項１に記載の吸収係数画像推定方法において、
前記再構成計算工程を、(a)前記画像μ’を未知数に含む計算アルゴリズムで実施する、または(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび前記吸収係数投影データを再構成した画像を前記画像μ’とするアルゴリズムの組み合わせで実施する、
吸収係数画像推定方法。
前記再構成計算工程を、(a)前記画像μ’を未知数に含む計算アルゴリズムで実施する、または(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび前記吸収係数投影データを再構成した画像を前記画像μ’とするアルゴリズムの組み合わせで実施する、
吸収係数画像推定方法。
【請求項3】
請求項１または請求項２に記載の吸収係数画像推定方法において、
(A)前記マスク算出工程は、
前記画像μ’の二値化画像を被検体マスク画像として算出する工程と、
前記被検体マスク画像の投影データを算出する工程と、
前記被検体マスク画像の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
からなる、
または
(B)前記マスク算出工程は、
前記画像μ’の投影データを算出する工程と、
前記画像μ’の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
からなる、
吸収係数画像推定方法。
(A)前記マスク算出工程は、
前記画像μ’の二値化画像を被検体マスク画像として算出する工程と、
前記被検体マスク画像の投影データを算出する工程と、
前記被検体マスク画像の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
からなる、
または
(B)前記マスク算出工程は、
前記画像μ’の投影データを算出する工程と、
前記画像μ’の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
からなる、
吸収係数画像推定方法。
【請求項4】
請求項１または請求項２に記載の吸収係数画像推定方法において、
前記マスク算出工程は、
前記計測データを投影データ形式に変換したものを二値化したデータを前記被検体マスク投影データとして算出する工程からなる、
吸収係数画像推定方法。
前記マスク算出工程は、
前記計測データを投影データ形式に変換したものを二値化したデータを前記被検体マスク投影データとして算出する工程からなる、
吸収係数画像推定方法。
【請求項5】
請求項１または請求項２に記載の吸収係数画像推定方法において、
(C)前記マスク算出工程は、
前記計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、
前記放射能画像の投影データを算出する工程と、
前記投影データを二値化したデータを前記被検体マスク投影データとして算出する工程と
からなる、
または
(D)前記マスク算出工程は、
前記計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、
前記放射能画像の二値化画像を算出する工程と、
前記二値化画像の投影データを算出する工程と、
前記二値化画像の投影データを二値化したデータを前記被検体マスク投影データとして算出する工程と
からなる、
吸収係数画像推定方法。
(C)前記マスク算出工程は、
前記計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、
前記放射能画像の投影データを算出する工程と、
前記投影データを二値化したデータを前記被検体マスク投影データとして算出する工程と
からなる、
または
(D)前記マスク算出工程は、
前記計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、
前記放射能画像の二値化画像を算出する工程と、
前記二値化画像の投影データを算出する工程と、
前記二値化画像の投影データを二値化したデータを前記被検体マスク投影データとして算出する工程と
からなる、
吸収係数画像推定方法。
【請求項6】
請求項１から請求項５のいずれかに記載の吸収係数画像推定方法において、
前記オフセット推定工程で実施する再構成処理は、解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施する、
吸収係数画像推定方法。
前記オフセット推定工程で実施する再構成処理は、解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施する、
吸収係数画像推定方法。
【請求項7】
請求項１から請求項６のいずれかに記載の吸収係数画像推定方法において、
前記参照領域抽出工程において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である、
吸収係数画像推定方法。
前記参照領域抽出工程において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である、
吸収係数画像推定方法。
【請求項8】
請求項１から請求項７のいずれかに記載の吸収係数画像推定方法において、
K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の前記領域Ωとし、μｎｋｎｏｗｎ(n=1,…,K)を前記領域Ωｎの既知の吸収係数値とし、S(X;Ωｎ)を、前記領域Ωｎにおける画像Xの統計量または統計量から算出される値の代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値の代表値とし、αｎを前記領域Ωｎにおける前記係数αとしたときに、前記係数算出工程における前記係数αは、α= T(α1,α2,…,αK)，αｎ=(μｎｋｎｏｗｎ-S(μ’; Ωｎ))/S(μｏｆｆ; Ωｎ) (n=1,…,K)である、
吸収係数画像推定方法。
K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の前記領域Ωとし、μｎｋｎｏｗｎ(n=1,…,K)を前記領域Ωｎの既知の吸収係数値とし、S(X;Ωｎ)を、前記領域Ωｎにおける画像Xの統計量または統計量から算出される値の代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値の代表値とし、αｎを前記領域Ωｎにおける前記係数αとしたときに、前記係数算出工程における前記係数αは、α= T(α1,α2,…,αK)，αｎ=(μｎｋｎｏｗｎ-S(μ’; Ωｎ))/S(μｏｆｆ; Ωｎ) (n=1,…,K)である、
吸収係数画像推定方法。
【請求項9】
請求項１から請求項７のいずれかに記載の吸収係数画像推定方法において、
K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の前記領域Ωとし、μｎｋｎｏｗｎ(n=1,…,K)を前記領域Ωｎ内に既知の吸収係数を設定した画像とし、DΩｎ(X,Y)を前記領域Ωｎ内に関する画像Xおよび画像Yの誤差評価関数とし、ｗｎ(n=1,…,K)を０以上１以下の係数としたときに、前記係数算出工程における前記係数αは、関数f(α)= Σn=1,…,K[ｗｎ×DΩｎ(μｎｋｎｏｗｎ，μ’+α×μｏｆｆ)]を最小化するαである、
吸収係数画像推定方法。
K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の前記領域Ωとし、μｎｋｎｏｗｎ(n=1,…,K)を前記領域Ωｎ内に既知の吸収係数を設定した画像とし、DΩｎ(X,Y)を前記領域Ωｎ内に関する画像Xおよび画像Yの誤差評価関数とし、ｗｎ(n=1,…,K)を０以上１以下の係数としたときに、前記係数算出工程における前記係数αは、関数f(α)= Σn=1,…,K[ｗｎ×DΩｎ(μｎｋｎｏｗｎ，μ’+α×μｏｆｆ)]を最小化するαである、
吸収係数画像推定方法。
【請求項10】
請求項１から請求項９のいずれかに記載の吸収係数画像推定方法をコンピュータに実行させる、吸収係数画像推定プログラム。
【請求項11】
請求項１０に記載の吸収係数画像推定プログラムを搭載したポジトロンＣＴ装置において、
当該吸収係数画像推定プログラムを実行する演算手段を備える、ポジトロンＣＴ装置。
当該吸収係数画像推定プログラムを実行する演算手段を備える、ポジトロンＣＴ装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、ポジトロンＣＴ装置（陽電子放射断層撮影装置）の計測データから吸収係数画像を推定する、吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置に関する。【背景技術】
【0002】
ポジトロンＣＴ装置、すなわちＰＥＴ(Positron Emission Tomography)装置は、陽電子（Positron）の消滅によって発生する２本のγ線を、複数個の検出器で同時に検出したときのみ（つまり同時計数したときのみ）有効な信号と見なして計測し、計測データに基づいて被検体の断層画像を再構成するように構成されている。具体的には、陽電子放出核種を含んだ放射性薬剤を被検体に投与して、投与された被検体から放出される511keVの対消滅γ線を多数の検出器素子（例えばシンチレータ）群からなる検出器で検出する。そして、２つの検出器で一定時間内にγ線を検出した場合に“同時”に検出したとして、それを一対の対消滅γ線として計数し、さらに対消滅発生位置を検出した２つの検出器を結ぶ直線（LOR: Line Of Response）を特定する。このように検出された同時計数情報を蓄積して再構成処理を行って、陽電子放出核種画像（すなわち断層画像）を得る。【0003】
ポジトロンＣＴ（ＰＥＴ）において、被検体内の放射能濃度を定量計測するためには、様々なデータ補正処理が必要である。代表的な補正処理には、感度補正，散乱補正，ランダム補正，減衰補正，デッドタイム補正，吸収補正がある。本発明は、放射性薬剤（ラジオアイソトープ）から放出されるγ線の吸収による画像の定量性低下を防止するための吸収補正に関する。吸収補正を行うためには被検体内の吸収係数分布を画像化した吸収係数画像を推定する必要がある。【0004】
吸収補正では、推定した吸収係数画像からγ線の透過率を求め、ＰＥＴの計測データから透過率を除算することでγ線の吸収の影響が排除されたデータに変換する。もしくは、推定した吸収係数画像を画像再構成の計算式に組み込んでγ線の吸収の影響が排除された再構成画像を取得する。【0005】
吸収係数画像を推定するためには、陽電子放出核種の外部線源を照射して得られたトランスミッション(Transmission)データが必要である。もしくは、トランスミッションデータの替わりにＸ線ＣＴ(Computed Tomography)装置から得られたＣＴデータを用いて、吸収係数画像を推定することが可能である。【0006】
近年では、このようなトランスミッションデータが不要な再構成アルゴリズムがある（例えば、非特許文献１，２参照）。非特許文献１，２では、消滅放射線の検出時間差（「飛行時間差」とも呼ばれる）(TOF: Time Of Flight)情報を計測するＰＥＴ（以下、「TOF-PET」と表記する）の計測データから、吸収係数サイノグラムの分布形状を推定することができる。そして、放射能画像および吸収係数サイノグラムを同時に推定することができる。非特許文献１，２における再構成アルゴリズムは、放射能画像および吸収係数に関するデータ（例えば吸収係数サイノグラム）を同時に推定することから「同時再構成アルゴリズム」とも呼ばれている。特に、放射能画像および吸収係数サイノグラムを同時に推定する同時再構成アルゴリズムは、MLACF法とも呼ばれている。【0007】
MLACF法における原理について、図９の概念図を参照して定性的に説明する。被検体の対象領域において放射能濃度が、図９に示すように分布しているとする。TOF情報t，半径方向sからなる２次元分布に放射能分布を展開すると、投影方向θ=0°のときには図９の上図のようになり、投影方向θ=90°のときには図９の右図のようになる。投影方向θを0°-180°の範囲で、それぞれの投影角度θ毎の２次元分布が計測データとして得られる。【0008】
もしγ線の吸収がなければ、投影角度θ毎の２次元分布を復元すると元の放射能分布に変換することができる。しかし、実際にはγ線の吸収により、元の放射能分布と異なる放射能分布に復元されてしまう。また、投影方向に応じてγ線の吸収の度合いが異なる。図９に示すように、投影方向θ=0°のときの吸収の度合いをAとし、投影方向θ=90°のときの吸収の度合いをA'とすると、通常はA≠A'となる。そこで、TOF情報の計測データを利用して、その計測データと、未知の吸収係数値と投影データとの乗算で得られた推定データとの誤差を最小にすることで、放射能画像と同時に、吸収係数値をサイノグラム形式に変換した吸収係数サイノグラムを推定することができる。【0009】
ただし、推定できるのは分布形状（つまり、相対値）のみで、絶対値を推定できない。吸収係数サイノグラムの絶対値が判らないと吸収補正を正しく行えない。したがって、定量的な放射能画像を取得できない。これは臨床上問題であり、実用化の大きな障壁となる。【0010】
同時再構成アルゴリズムに基づくPETイメージング技術を実用化するためには、下記のような技術が必要不可欠である。すなわち、分布形状は正しいが、絶対値が誤っている吸収係数サイノグラムから、絶対値が正しい吸収係数画像を推定する技術が必要不可欠である。【0011】
そこで、非特許文献２の従来技術として、吸収係数サイノグラムの直接的な定量化を行わずに、放射能画像を直接的に定量化するという方法がある。処理ステップは以下の通りである。(1)吸収補正なしの放射能画像を作成し、当該放射能画像から被検体マスク画像を作成する。
(2)被検体マスク画像の各画素値に既知の吸収係数（例えば、水の吸収係数や骨の吸収係数）を代入し、疑似的な吸収係数画像（以下、「疑似吸収係数画像」と略記する）を作成する。
(3)疑似吸収係数画像を用いて、従来の再構成アルゴリズムで放射能画像を作成し、体軸方向スライス毎に合計画素値を計算する。図１０に示すようにzを体軸方向とし、xy平面を体軸方向zに直交する平面（アキシャル面）とし、P(x,y)を画素値とし、S(z)を体軸方向zスライス毎の合計画素値とすると、S(z)=ΣyΣxP(x,y)となる。横軸を体軸方向z，縦軸を合計画素値S(z)とすると、図１０の右図のようなプロファイルが作成される。
(4)同時再構成アルゴリズムによって、定量的でない放射能画像および定量的でない吸収係数サイノグラムを推定する。
(5)手順(4)で求めた定量的でない放射能画像に対しても、手順(3)と同様にして、体軸方向スライス毎に合計画素値を計算する。ここでは、S'(z)を体軸方向zスライス毎の合計画素値とする。
(6)手順(5)で求めた合計画素値S'(z)と手順(3)で求めた合計画素値S(z)とが体軸方向スライス毎に一致するように、手順(4)で求めた定量的でない放射能画像をスケーリングし、定量化する。具体的には、手順(4)で求めた定量的でない放射能画像の各画素値にS(z)/ S'(z)を体軸方向スライス毎にそれぞれ乗算することで、当該放射能画像をスケーリングする。
【先行技術文献】
【非特許文献】
【0012】
【非特許文献1】A. Rezaei, M. Defrise and J. Nuyts, “ML-reconstruction for TOF-PET with simultaneous estimation of the Attenuation Factors”, IEEE Trans. Med. Imag., 33 (7), 1563-1572, 2014
【0013】
【非特許文献2】V. Y. Panin, H. Bal, M. Defrise, C. Hayden and M. E. Casey, “Transmission-less Brain TOF PET Imaging using MLACF”, 2013 IEEE Nuclear Science Symposium and Medical Imaging Conference, Seoul, Republic of Korea, 2013.
【発明の概要】
【発明が解決しようとする課題】
【0014】
しかしながら、上述した非特許文献２の従来技術では、最終的に求まる放射能画像が定量的であることを原理的に保証できないという問題点がある。【0015】
本発明は、このような事情に鑑みてなされたものであって、定量的な吸収係数画像を作成することができる吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置を提供することを目的とする。【課題を解決するための手段】
【0016】
発明者は、上記の問題を解決するために鋭意研究した結果、次のような知見を得た。すなわち、上述した従来技術は、手順(3)で求めた体軸方向スライス毎の合計画素値が真である（すなわち正しい）ということを仮定している。
【0017】
しかし、手順(3)で求めた体軸方向スライス毎の合計画素値が正しいという保証はない。その正しい保証のない結果に基づいて定量化しようとしているので、手順(6)で求まる放射能画像の定量性も当然ながら保証できない。したがって、従来技術では定量的な放射能画像を作成することができないという問題点が生じる。つまり、従来技術の問題の根源は、経験的な方法であり、吸収係数の定量化という課題の背景にある数学的な理論に基づいていないことであるという知見を得た。【0018】
このような知見に基づく本発明は、次のような構成をとる。すなわち、本発明に係る吸収係数画像推定方法は、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データから吸収係数画像を推定する方法であって、μ'を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、前記計測データに関する評価関数の最適化に基づいて、前記画像μ'を計算する再構成計算工程と、前記計測データに基づいて投影データ空間における被検体マスクデータである被検体マスク投影データを算出するマスク算出工程と、μoffを不均一なオフセット画像としたときに、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、前記オフセット画像μoffを推定するオフセット推定工程と、Ωを既知の吸収係数値で近似可能な領域としたときに、前記計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、前記領域Ωを少なくとも１つ以上抽出する参照領域抽出工程と、αを係数としたときに、前記領域Ωにおける前記画像μ'の値と既知の吸収係数値との誤差を減少させる前記係数αを算出する係数算出工程と、前記画像μ'の値に、前記オフセット画像μoffを前記係数α倍したα×μoffを加算して得られた値を吸収係数値として補正する吸収係数値補正工程とを備えるものである。
【0019】
本発明に係る吸収係数画像推定方法によれば、再構成計算工程では、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データ（TOF-PETの計測データ）に関する評価関数の最適化に基づいて、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像を計算する。それと同時に、λ'を非定量的な放射能画像としたときに、再構成計算工程では放射能画像λ'を計算する。ここでの再構成計算工程における再構成アルゴリズムは、上述した同時再構成アルゴリズムである。ここで、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像をμ'とすると、画像μ'も吸収係数画像となるが、最終的に求まる定量的な吸収係数画像と区別して、単に「画像μ'」とする。【0020】
一方、マスク算出工程では、計測データに基づいて投影データ空間における被検体マスクデータ（以下、「被検体マスク投影データ」と呼ぶ）を算出する。μoffを不均一なオフセット画像としたときに、オフセット推定工程では、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。【0021】
一方、Ωを既知の吸収係数値で近似可能な領域としたときに、計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、参照領域抽出工程では、領域Ωを少なくとも１つ以上抽出する。【0022】
ここで、「計測データに基づいて計算された、被検体領域が認識可能な画像」とは、例えば、上述した画像μ'，上述した放射能画像λ'，上述した再構成計算工程における再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズムで推定した非定量的な画像，上述した被検体マスク投影データから推定した被検体マスク画像である。もちろん、これらの例示した画像に限定されず、計測データに基づいて計算された、被検体領域が認識可能な画像であればよい。また、「計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する」ことは、単一の画像情報を用いて領域Ωを抽出することのみならず、複数の画像情報の組み合わせ（例えば論理和や論理積）を用いて領域Ωを抽出することも含まれる。【0023】
上述した各工程で得られた画像μ'，オフセット画像μoffおよび領域Ωを用いて、下記の係数算出工程および吸収係数値補正工程を実施する。μを未知の真の吸収係数画像の値（吸収係数値）とし、αを係数としたときに、係数算出工程では、領域Ωにおける画像μ'の値と既知の吸収係数値との誤差を減少させる係数αを算出する。そして、吸収係数値補正工程では、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値として補正する。このように、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能である（つまり、μ≒μ'+α×μoff）という数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0024】
上述した再構成計算工程を、(a)画像μ'を未知数に含む計算アルゴリズムで実施してもよい。または、上述した再構成計算工程を、(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムの組み合わせで実施してもよい。【0025】
前者の(a)のアルゴリズムは、非定量的な放射能画像λ'および画像μ'（非定量的な吸収係数画像）を同時に再構成するMLAA法である。後者の(b)のアルゴリズムは、非定量的な放射能画像λ'および非定量的な吸収係数投影データ（例えば吸収係数サイノグラム）を同時に推定する、非特許文献１で述べたMLACF法と、吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムとの組み合わせである。後者の(b)における、吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムについては、再構成アルゴリズムであれば特に限定されない。【0026】
上述したマスク算出工程の一例は、画像μ'の二値化画像を被検体マスク画像として算出する工程と、被検体マスク画像の投影データを算出する工程と、被検体マスク画像の投影データの二値化データを被検体マスク投影データとして算出する工程とからなる（(A)の態様）。また、マスク算出工程は、画像μ'の投影データを算出する工程と、画像μ'の投影データの二値化データを被検体マスク投影データとして算出する工程とからなる(B)の態様でもよい。上記(A)の態様では、画像μ'を先に二値化した後に投影データを算出し、その投影データを二値化した。上記(B)の態様では、画像μ'の投影データを先に算出した後に二値化する。画像μ'は定量的でない吸収係数画像であるが、吸収係数画像以外の画像や、投影データを二値化しても、被検体マスク投影データを算出することができる。【0027】
マスク算出工程の他の一例は、上述したTOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出する工程からなる。この一例の場合には、TOF-PETの計測データを直接的に用いて、投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出することができる。【0028】
また、マスク算出工程のさらなる他の一例は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の投影データを算出する工程と、（放射能画像の）投影データを二値化したデータを被検体マスク投影データとして算出する工程とからなる（(C)の態様）。また、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の二値化画像を算出する工程と、二値化画像の投影データを算出する工程と、二値化画像の投影データを二値化したデータを被検体マスク投影データとして算出する工程とからなる(D)の態様でもよい。上記(C)の態様では、放射能画像の投影データを先に算出した後に二値化していた。上記(D)の態様では、放射能画像を先に二値化した後に投影データを算出し、その投影データを二値化する。上記(C)の態様または上記(D)の態様のいずれにおいても、この一例の場合には、放射能画像を算出する再構成アルゴリズムについては、特に限定されない。上述した同時再構成アルゴリズム（MLACF法やMLAA法）を用いた場合には、同時再構成アルゴリズムで推定した上述の放射能画像λ'を利用して、被検体マスク投影データを算出することができる。また、上述した同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した放射能画像を利用して、被検体マスク投影データを算出してもよい。【0029】
上述したオフセット推定工程で実施する再構成処理は、解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施すればよい。解析的再構成として例えばFBP(Filtered Back Projection)法がある。統計的再構成として例えば上述したML-EM(Maximum Likelihood - Expectation Maximization)法がある。代数的再構成として例えばART(Algebraic Reconstruction Technique)法がある。【0030】
上述した参照領域抽出工程において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である。ここで、「吸収係数を既知と見なせる組織の領域」とは、例えば、水で近似可能な領域，空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域，肺組織で近似可能な領域，軟部組織で近似可能な領域などである。もちろん、これらの例示した組織に限定されず、吸収係数の近似値が判る組織であればよい。なお、脳組織の場合には水で近似可能な領域と見なされるので、水の吸収係数値である0.0096/mmを用いることができる。【0031】
上述した係数算出工程における係数αについては、下記のように代表値に基づいて算出する手法（前者の手法）と、下記のように誤差評価関数に基づいて算出する手法（後者の手法）とがある。【0032】
前者の手法においては、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の領域Ωとし、μｎｋｎｏｗｎ(n=1,…,K)を領域Ωｎの既知の吸収係数値とし、S(X；Ωｎ)を、領域Ωｎにおける画像Xの統計量または統計量から算出される値の代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値の代表値とし、αｎを領域Ωｎにおける係数αとする。ここで、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。もちろん、これらの例示した値に限定されず、「統計量または統計量から算出される値」であればよい。【0033】
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnの既知の吸収係数値μnknownは、領域Ωnにおける画像μ'の代表値S(μ';Ωn)に、領域Ωnにおけるオフセット画像μoff'の代表値S(μoff;Ωn)に領域Ωnにおける係数αnを乗算したαn×S(μoff;Ωn)を加算して得られた値と等しい。つまり、μnknown=S(μ';Ωn)+αn×S(μoff;Ωn)が成立する。したがって、αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn)によって、領域Ωnにおける係数αnをそれぞれ算出する。K個の値x1,x2,…,xKをα1,α2,…,αKにそれぞれ置き換えれば、T(α1,α2,…,αK)は係数α1,α2,…,αKの代表値である。以上をまとめると、n=1,…,Kで各領域Ω1,…,ΩKにおける係数α1,…,αKをそれぞれ算出した後に、係数α1,…,αKの代表値T(α1,α2,…,αK)を係数αとすることで、係数αを算出することができる。【0034】
後者の手法においては、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとし、μnknown (n=1,…,K)を領域Ωn内に既知の吸収係数を設定した画像とし、DΩn(X,Y)を領域Ωn内に関する画像Xおよび画像Yの誤差評価関数とし、wn (n=1,…,K)を０以上１以下の係数とする。画像Xを領域Ωn内に既知の吸収係数を設定した画像μnknownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩn(μnknown,μ'+α×μoff)は、領域Ωn内に関する既知の吸収係数を設定した画像μnknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。n=1,…,Kで各領域Ω1,…,ΩK毎の係数w1,…,wKによるDΩ1(μ1known,μ'+α×μoff),…,DΩK(μKknown,μ'+α×μoff)の重み付き加算で得られた、係数αを変数とした関数f(α)(=Σn=1,…,K [wn×DΩn(μnknown,μ'+α×μoff)])を最小化するαを算出することで、係数αを算出することができる。【0035】
また、本発明に係る吸収係数画像推定プログラムは、本発明に係る吸収係数画像推定方法をコンピュータに実行させるものである。【0036】
本発明に係る吸収係数画像推定プログラムによれば、本発明に係る吸収係数画像推定方法をコンピュータに実行させることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0037】
また、本発明に係るポジトロンＣＴ装置は、本発明に係る吸収係数画像推定プログラムを搭載したポジトロンＣＴ装置において、当該吸収係数画像推定プログラムを実行する演算手段を備えるものである。【0038】
本発明に係るポジトロンＣＴ装置によれば、本発明に係る吸収係数画像推定プログラムを実行する演算手段を備えることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【発明の効果】
【0039】
本発明に係る吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置によれば、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【図面の簡単な説明】
【0040】
【図1】各実施例に係るＰＥＴ装置の概略斜視図およびブロック図である。
【図2】γ線検出器の概略斜視図である。
【図3】実施例１に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図4】実施例２に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図5】MLACF法で推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図6】MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図7】実施例４に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図8】実施例５に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
【図9】MLACF法における原理の説明に供する概念図である。
【図10】体軸方向スライス毎の合計画素値、および横軸を体軸方向，縦軸を合計画素値としたときのプロファイルの概略図である。
【実施例1】
【0041】
以下、図面を参照して本発明の実施例１を説明する。図１は、各実施例に係るＰＥＴ装置の概略斜視図およびブロック図であり、図２は、γ線検出器の概略斜視図である。また、図１および図２は各実施例とも共通の構成である。【0042】
図１に示すように、ＰＥＴ装置１は、被検体の周囲を取り囲む検出器リング２を被検体の体軸方向に積層配置して備えている。検出器リング２内には複数のγ線検出器３が埋設されている。ＰＥＴ装置１は、本発明におけるポジトロンＣＴ装置に相当する。また、γ線検出器３は、本発明における検出器に相当する。【0043】
その他にも、ＰＥＴ装置１は、同時計数回路４と演算回路５とを備えている。図１では、γ線検出器３から同時計数回路４への結線を２つのみ図示しているが、実際には、γ線検出器３の光電子増倍管(PMT: Photo Multiplier Tube)３３（図２を参照）の総チャンネル数分、同時計数回路４に接続されている。演算回路５は、吸収係数画像推定プログラム６による、後述する図３に示す吸収係数画像推定方法の処理を実行する。演算回路５は、本発明における演算手段に相当する。【0044】
放射性薬剤が投与された被検体（図示省略）から発生したγ線をγ線検出器３のシンチレータブロック３１（図２を参照）が光に変換して、変換されたその光をγ線検出器３の光電子増倍管（ＰＭＴ）３３（図２を参照）は増倍させて電気信号に変換する。その電気信号を同時計数回路４に送り込み、カウント値の検出信号データを生成する。【0045】
具体的には、被検体（図示省略）に放射性薬剤を投与すると、ポジトロン放出型のＲＩのポジトロンが消滅することにより、２本のγ線が発生する。同時計数回路４は、シンチレータブロック３１（図２を参照）の位置とγ線の入射タイミングとをチェックし、被検体の両側にある２つのシンチレータブロック３１でγ線が同時に入射したときのみ、送り込まれた電気信号を適正なデータと判定する。一方のシンチレータブロック３１のみにγ線が入射したときには、同時計数回路４は棄却する。つまり、同時計数回路４は、上述した電気信号に基づいて、２つのγ線検出器３においてγ線が同時観測（すなわち同時計数）されたことを検出する。【0046】
同時計数回路４により同時計数と判定された適正なデータにより構成された検出信号データ（カウント値）を、演算回路５に送り込む。演算回路５は、後述するステップＳ１〜Ｓ８（図３を参照）を行って、ＰＥＴ装置１によって得られた被検体（図示省略）の検出信号データから吸収係数画像を推定する。演算回路５の具体的な機能については後述する。【0047】
なお、ＲＯＭ（Read-only Memory）などに代表される記憶媒体（図示省略）に吸収係数画像推定プログラム６が記憶されており、当該吸収係数画像推定プログラム６を記憶媒体から演算回路５に読み出して、吸収係数画像推定プログラム６を演算回路５が実行することで、図３のフローチャートに示す吸収係数画像推定方法の処理を行う。演算回路５は、ＧＰＵ（Graphics Processing Unit），中央演算処理装置（ＣＰＵ）あるいはプログラムデータに応じて内部の使用するハードウェア回路（例えば論理回路）が変更可能なプログラマブルデバイス（例えばＦＰＧＡ(Field Programmable Gate Array)）などで構成されている。【0048】
γ線検出器３は、図２に示すようにシンチレータブロック３１と、そのシンチレータブロック３１に対して光学的に結合されたライトガイド３２と、そのライトガイド３２に対して光学的に結合された光電子増倍管（以下、単に「ＰＭＴ」と略記する）３３とを備えている。シンチレータブロック３１を構成する各シンチレータ素子は、γ線の入射に伴って発光することでγ線から光に変換する。この変換によってシンチレータ素子はγ線を検出する。シンチレータ素子において発光した光がシンチレータブロック３１で十分に拡散されて、ライトガイド３２を介してＰＭＴ３３に入力される。ＰＭＴ３３は、シンチレータブロック３１で変換された光を増倍させて電気信号に変換する。その電気信号は画素値として同時計数回路４（図１を参照）に送り込まれる。【0049】
また、γ線検出器３は、図２に示すように、３次元的に配置されたシンチレータ素子からなり、深さ方向に複数の層からなるＤＯＩ検出器である。図２では、４層のＤＯＩ検出器を図示しているが、層の数については、複数であれば特に限定されない。【0050】
ここで、ＤＯＩ検出器は、各々のシンチレータ素子を放射線の深さ方向に積層して構成されたものであり、相互作用を起こした深さ(DOI: Depth of Interaction)方向と横方向（入射面に平行な方向）との座標情報を重心演算により求める。ＤＯＩ検出器を用いることにより深さ方向の空間分解能をより一層向上させることができる。よって、ＤＯＩ検出器の層の数は、深さ方向に積層されたシンチレータ素子の層の数である。【0051】
次に、演算回路５の具体的な機能について、図３を参照して説明する。図３は、実施例１に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。【0052】
先ず、図１に示すＰＥＴ装置１によって被検体の撮影を行い、同時計数回路４（図１を参照）によってリストモードデータを取得する。リストモードデータには、検出した光子のエネルギー情報が記録されている。【0053】
通常のエネルギーウィンドウ（例えば400keV-600keV）、すなわち再構成データ用のエネルギーウィンドウ、およびTOF計測データのTOF方向の計測範囲とビン幅とをそれぞれ設定する。なお、ここで、ビン(bin)とは、離散化すること（区切ること）を意味する。画像の場合には画素がビンに対応する。TOFビンはTOF情報の時間的区切りを意味し、例えばTOFビンが100[ps]の場合には、検出時間差は100[ps]の精度で時間的に区切られる。【0054】
その設定にしたがって、リストモードデータから、TOF-PETの計測データを作成する。【0055】
（ステップＳ１）MLACFμ'を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、λ'を非定量的な放射能画像とする。計測データに関する評価関数の最適化に基づいて、画像μ'および放射能画像λ'を同時に計算する。「課題を解決するための手段」の欄でも述べたように、画像μ'も吸収係数画像となるが、最終的に求まる定量的な吸収係数画像と区別して、単に「画像μ'」とする。
【0056】
ステップＳ１での再構成アルゴリズムは、非特許文献１での同時再構成アルゴリズムである。後述する実施例２，３も含めて、本実施例１では、同時再構成アルゴリズムとして、非特許文献１で述べたMLACF法を適用する。MLACF法の具体的な手法については、上述した非特許文献１を参照されたい。【0057】
A'を吸収係数サイノグラムとする。MLACF法によって、放射能画像λ'および吸収係数サイノグラムA'を推定する。吸収係数サイノグラムA'は、本発明における吸収係数投影データに相当する。【0058】
（ステップＳ２）ML-TR or ML-EM吸収係数サイノグラムA'を再構成アルゴリズム（例えばML-TR法やML-EM法）で再構成し、定量的でない画像μ'を作成する。ML-EM法を利用する場合は、吸収係数サイノグラムA’を事前にログ変換する。ML-TR法の具体的な手法については、参考文献１を参照されたい（参考文献１：Erdo?an H, Fessler JA: Ordered subsets algorithms for transmission tomography. Phys Med Biol 44: 2835-2851, 1999）。また、ML-EM法の具体的な手法については、参考文献２を参照されたい（参考文献２：L.A. Shepp and Y. Vardi. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging, Vol. 1, pp. 113-122, 1982）。ステップＳ１，Ｓ２は、本発明における再構成計算工程に相当する。
【0059】
（ステップＳ３）Binarization Processing閾値処理によって画像μ'を二値化処理(Binarization Processing)する。そして、被検体領域を“１”，その他の領域を“０”となる二値化画像を被検体マスク画像として算出する。mimgを被検体マスク画像とする。
【0060】
（ステップＳ４）Projection + Binarization Processing被検体マスク画像mimgの線積分データ（投影データ）を算出する（Projection）。そして、閾値処理によって被検体マスク画像mimgの投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データとして算出する。mprojを被検体マスク投影データとする。ステップＳ３，Ｓ４は、本発明におけるマスク算出工程に相当する。
【0061】
（ステップＳ５）ML-EMμoffを不均一なオフセット画像とする。オフセット画像μoffの順投影データが被検体マスク投影データmprojを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。つまり、被検体マスク投影データmprojを再構成アルゴリズム（例えばML-EM法）で画像データに変換する。この画像データをオフセット画像μoffとする。ステップＳ５は、本発明におけるオフセット推定工程に相当する。
【0062】
（ステップＳ６）ExtractionΩを既知の吸収係数値で近似可能な領域とする。計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する（Extraction）。後述する実施例２〜４も含めて、本実施例１では、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとしたときに、n=1,…,Kで各領域Ω1,…,ΩKをそれぞれ抽出する。例えば、被検体の頭部を撮影する場合には、空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域で各領域Ω1,…,ΩKをそれぞれ抽出する。
【0063】
計測データに基づいて計算された、被検体領域が認識可能な画像としては、例えば、ステップＳ２で作成した画像μ'，ステップＳ１で求まった放射能画像λ'，ステップＳ１での再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した非定量的な画像，上述した被検体マスク投影データmprojから推定した被検体マスク画像mimgである。ステップＳ６は、本発明における参照領域抽出工程に相当する。【0064】
（ステップＳ７）αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)μnknown(n=1,…,K)を領域Ωnの既知の吸収係数値とし、S(X;Ωn)を、領域Ωnにおける画像Xの統計量または統計量から算出される値を代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値を代表値とし、αnを領域Ωnにおける係数αとする。なお、αは係数である。ここで、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。ここで、「トリムド平均値」とは、値が極端に大きい／小さいデータを取り除いた残りのデータの平均値を意味する。「トリムド中央値」とは、値が極端に大きい／小さいデータを取り除いた残りのデータの中央値を意味する。
【0065】
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnにおける係数αnは下記（１）式のように表される。αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn) …（１）
【0066】
K個の値x1,x2,…,xKをα1,α2,…,αKにそれぞれ置き換えれば、T(α1,α2,…,αK)は係数α1,α2,…,αKの代表値である。したがって、係数αは下記（２）式のように表される。ステップＳ７は、本発明における係数算出工程に相当する。α=T(α1,α2,…,αK) …（２）
【0067】
（ステップＳ８）μ=μ'+α×μoffμを未知の真の吸収係数画像の値（吸収係数値）とする。吸収係数値μは下記（３）式のように表される。
μ=μ'+α×μoff …（３）
【0068】
上記（３）式のように、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値μとして補正する。ステップＳ８は、本発明における吸収係数値補正工程に相当する。【0069】
ステップＳ８で求まった吸収係数値μを有する吸収係数画像を用いて吸収補正を行う。「背景技術」の欄でも述べたように、推定した吸収係数画像からγ線の透過率を求め、ＰＥＴの計測データから透過率を除算することでγ線の吸収の影響が排除されたデータに変換することによって、吸収補正を行う。もしくは、推定した吸収係数画像を画像再構成の計算式に組み込んでγ線の吸収の影響が排除された再構成画像を取得することによって、吸収補正を行う。【0070】
本実施例１に係る吸収係数画像推定方法によれば、ステップＳ１，Ｓ２では、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データ（TOF-PETの計測データ）に関する評価関数の最適化に基づいて、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像μ'を計算する。それと同時に、ステップＳ１では放射能画像λ'を計算する。ステップＳ１における再構成アルゴリズムは、上述した同時再構成アルゴリズムである。【0071】
一方、ステップＳ３，Ｓ４では、計測データに基づいて投影データ空間における被検体マスクデータ（すなわち、被検体マスク投影データmproj）を算出する。ステップＳ５では、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。【0072】
一方、計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、ステップＳ６では、領域Ωを少なくとも１つ以上抽出する。【0073】
ここで、「課題を解決するための手段」の欄でも述べたように、「計測データに基づいて計算された、被検体領域が認識可能な画像」とは、例えば、ステップＳ２で作成した画像μ'，ステップＳ１で求まった放射能画像λ'，ステップＳ１での再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した非定量的な画像，上述した被検体マスク投影データmprojから推定した被検体マスク画像mimgである。もちろん、これらの例示した画像に限定されず、計測データに基づいて計算された、被検体領域が認識可能な画像であればよい。また、「課題を解決するための手段」の欄でも述べたように、「計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する」ことは、単一の画像情報を用いて領域Ωを抽出することのみならず、複数の画像情報の組み合わせ（例えば論理和や論理積）を用いて領域Ωを抽出することも含まれる。【0074】
上述した各ステップＳ１〜Ｓ６で得られた画像μ'，オフセット画像μoffおよび領域Ωを用いて、ステップＳ７，Ｓ８を実施する。ステップＳ７では、領域Ωにおける画像μ'の値と既知の吸収係数値との誤差を減少させる係数αを算出する。そして、ステップＳ８では、上記（３）式のように、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値μとして補正する。このように、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能である（つまり、μ≒μ'+α×μoff）という数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0075】
後述する実施例２，３も含めて、本実施例１では、上述した再構成計算工程に相当するステップＳ１，Ｓ２を、(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムの組み合わせで実施している。【0076】
(b)のアルゴリズムは、非定量的な放射能画像λ'および非定量的な吸収係数投影データ（各実施例１〜３では吸収係数サイノグラムA'）を同時に推定する、非特許文献１で述べたMLACF法と、吸収係数投影データ（吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムとの組み合わせである。「課題を解決するための手段」の欄でも述べたように、(b)における、吸収係数投影データ（吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムについては、再構成アルゴリズムであれば特に限定されない。図３では、ステップＳ２で述べたように、ML-TR法またはML-EM法で吸収係数サイノグラムA'を再構成し、画像μ'を作成する。【0077】
本実施例１では、上述したマスク算出工程に相当するステップＳ３，Ｓ４を実施している。すなわち、マスク算出工程は、画像μ'の二値化画像を被検体マスク画像mimgとして算出する工程（ステップＳ３）と、被検体マスク画像mimgの投影データを算出する工程（ステップＳ４の前半の「Projection」）と、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出する工程（ステップＳ４の後半の「Binarization Processing」）とからなる。画像μ'は定量的でない吸収係数画像であるが、後述する実施例２，３のように、吸収係数画像以外の画像（例えば実施例３の放射能画像λ'，または放射能画像λ２'）や、投影データを二値化しても、被検体マスク投影データmprojを算出することができる。【0078】
上述したオフセット推定工程に相当するステップＳ５で実施する再構成処理は、本実施例１ではML-EM(Maximum Likelihood - Expectation Maximization)法等に代表される統計的再構成の計算方式で実施している。もちろん、ステップＳ５で実施する再構成処理は、本実施例１のような統計的再構成に限定されない。解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施すればよい。解析的再構成として例えばFBP(Filtered Back Projection)法がある。代数的再構成として例えばART(Algebraic Reconstruction Technique)法がある。【0079】
上述した参照領域抽出工程に相当するステップＳ６において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である。ここで、「課題を解決するための手段」の欄でも述べたように、「吸収係数を既知と見なせる組織の領域」とは、例えば、水で近似可能な領域，空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域，肺組織で近似可能な領域，軟部組織で近似可能な領域などである。もちろん、これらの例示した組織に限定されず、吸収係数の近似値が判る組織であればよい。なお、脳組織の場合には水で近似可能な領域と見なされるので、水の吸収係数値である0.0096/mmを用いることができる。【0080】
本実施例１では、上述した係数算出工程に相当するステップＳ７を実施している。すなわち、後述する実施例２〜４も含めて、本実施例１では、代表値に基づいて係数αを算出している。上述したように、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。「課題を解決するための手段」の欄でも述べたように、これらの例示した値に限定されず、「統計量または統計量から算出される値」であればよい。【0081】
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnの既知の吸収係数値μnknownは、領域Ωnにおける画像μ'の代表値S(μ';Ωn)に、領域Ωnにおけるオフセット画像μoff'の代表値S(μoff;Ωn)に領域Ωnにおける係数αnを乗算したαn×S(μoff;Ωn)を加算して得られた値と等しい。つまり、μnknown= S(μ';Ωn)+αn×S(μoff;Ωn)が成立する。したがって、上記（１）式のように、αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn)によって、領域Ωnにおける係数αnをそれぞれ算出する。K個の値x1,x2,…,xKをα1,α2,…,αKにそれぞれ置き換えれば、T(α1,α2,…,αK)は係数α1,α2,…,αKの代表値である。以上をまとめると、n=1,…,Kで各領域Ω1,…,ΩKにおける係数α1,…,αKをそれぞれ算出した後に、上記（２）式のように、係数α1,…,αKの代表値T(α1,α2,…,αK)を係数αとすることで、係数αを算出することができる。【0082】
本実施例１に係る吸収係数画像推定プログラム６（図１を参照）によれば、本実施例１に係る吸収係数画像推定プログラムをコンピュータ（各実施例では図１に示す演算回路５を構成するＧＰＵ，ＣＰＵまたはプログラマブルデバイス）に実行させることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0083】
本実施例１に係るＰＥＴ装置１（図１を参照）によれば、本実施例１に係る吸収係数画像推定プログラム６を実行する演算手段（各実施例では図１に示す演算回路５を構成するＧＰＵ，ＣＰＵまたはプログラマブルデバイス）を備えることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【実施例2】
【0084】
次に、図面を参照して本発明の実施例２を説明する。図４は、実施例２に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。【0085】
上述した実施例１では、画像μ'の二値化画像を被検体マスク画像mimgとして算出し、被検体マスク画像mimgの投影データを算出し、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出していた。これに対して、本実施例２では、TOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出する。【0086】
（ステップＳ１１）MLACF図４のステップＳ１１は、上述した実施例１のステップＳ１と同じであるので、その説明については省略する。
【0087】
（ステップＳ１２）ML-TR or ML-EM図４のステップＳ１２は、上述した実施例１のステップＳ２と同じであるので、その説明については省略する。ステップＳ１１，Ｓ１２は、本発明における再構成計算工程に相当する。
【0088】
（ステップＳ１４）Projection + Binarization Processing上述した実施例１では、画像μ'から被検体マスク画像mimgを算出するのに、図３のステップＳ３，Ｓ４を実施していた。本実施例２では、画像μ'の替わりに、計測データから被検体マスク画像mimgを算出するために、図４のステップＳ１４を実施する。先ず、計測データを投影データ形式に変換する（Conversion）。そして、閾値処理によって計測データの投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ１４は、本発明におけるマスク算出工程に相当する。
【0089】
（ステップＳ１５）ML-EM図４のステップＳ１５は、上述した実施例１のステップＳ５と同じであるので、その説明については省略する。ステップＳ１５は、本発明におけるオフセット推定工程に相当する。
【0090】
（ステップＳ１６）Extraction図４のステップＳ１６は、上述した実施例１のステップＳ６と同じであるので、その説明については省略する。ステップＳ１６は、本発明における参照領域抽出工程に相当する。
【0091】
（ステップＳ１７）αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)図４のステップＳ１７は、上述した実施例１のステップＳ７と同じであるので、その説明については省略する。ステップＳ１７は、本発明における係数算出工程に相当する。
【0092】
（ステップＳ１８）μ=μ'+α×μoff図４のステップＳ１８は、上述した実施例１のステップＳ８と同じであるので、その説明については省略する。ステップＳ１８は、本発明における吸収係数値補正工程に相当する。
【0093】
本実施例２に係る吸収係数画像推定方法によれば、上述した実施例１と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ１８で吸収係数値を補正する。したがって、ステップＳ１８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0094】
本実施例２では、上述したマスク算出工程に相当するステップＳ１４を実施している。すなわち、マスク算出工程は、上述したTOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出する工程（ステップＳ１４）からなる。本実施例２の場合には、TOF-PETの計測データを直接的に用いて、投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出することができる。【0095】
本実施例２に係る吸収係数画像推定プログラム６（図１を参照）や本実施例２に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１と同じであるので、その説明については省略する。【実施例3】
【0096】
次に、図面を参照して本発明の実施例３を説明する。図５は、MLACF法で推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートであり、図６は、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。【0097】
上述した実施例１では、画像μ'の二値化画像を被検体マスク画像mimgとして算出し、被検体マスク画像mimgの投影データを算出し、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出していた。また、上述した実施例２では、TOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出していた。これに対して、本実施例３では、TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出し、放射能画像の投影データを算出し、（放射能画像の）投影データを二値化したデータを被検体マスク投影データmprojとして算出する。【0098】
先ず、MLACF法で推定した放射能画像を利用して被検体マスク投影データmprojを算出する場合について、図５を参照して説明する。【0099】
（ステップＳ２１）MLACF図５のステップＳ２１は、上述した実施例１のステップＳ１，上述した実施例２のステップＳ１１と同じであるので、その説明については省略する。
【0100】
（ステップＳ２２）ML-TR or ML-EM図５のステップＳ２２は、上述した実施例１のステップＳ２，上述した実施例２のステップＳ１２と同じであるので、その説明については省略する。ステップＳ２１，Ｓ２２は、本発明における再構成計算工程に相当する。
【0101】
（ステップＳ２４）Projection + Binarization Processing上述した実施例１では、画像μ'から被検体マスク画像mimgを算出するのに、図３のステップＳ３，Ｓ４を実施していた。また、上述した実施例２では、計測データから被検体マスク画像mimgを算出するのに、図４のステップＳ１４を実施していた。本実施例３では、実施例１の画像μ'，実施例２の計測データの替わりに、計測データに関する評価関数の最適化に基づいて推定した放射能画像から被検体マスク画像mimgを算出するために、図５のステップＳ２４を実施する。
【0102】
先ず、図５では、MLACF法における計測データに関する評価関数の最適化に基づいて推定した放射能画像λ'の線積分データ（投影データ）を算出する（Projection）。つまり、図５では、ステップＳ２１のMLACF法で推定した放射能画像λ'を利用して放射能画像λ'の投影データを算出すればよい。そして、閾値処理によって放射能画像λ'の投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ２４は、本発明におけるマスク算出工程に相当する。【0103】
（ステップＳ２５）ML-EM図５のステップＳ２５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５と同じであるので、その説明については省略する。ステップＳ２５は、本発明におけるオフセット推定工程に相当する。
【0104】
（ステップＳ２６）Extraction図５のステップＳ２６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６と同じであるので、その説明については省略する。ステップＳ２６は、本発明における参照領域抽出工程に相当する。
【0105】
（ステップＳ２７）αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)図５のステップＳ２７は、上述した実施例１のステップＳ７，上述した実施例２のステップＳ１７と同じであるので、その説明については省略する。ステップＳ２７は、本発明における係数算出工程に相当する。
【0106】
（ステップＳ２８）μ=μ'+α×μoff図５のステップＳ２８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８と同じであるので、その説明については省略する。ステップＳ２８は、本発明における吸収係数値補正工程に相当する。
【0107】
次に、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データmprojを算出する場合について、図６を参照して説明する。【0108】
（ステップＳ３１）MLACF図６のステップＳ３１は、図５のステップＳ２１と同じであるので、その説明については省略する。
【0109】
（ステップＳ３２）ML-TR or ML-EM図６のステップＳ３２は、図５のステップＳ２２と同じであるので、その説明については省略する。ステップＳ３１，Ｓ３２は、本発明における再構成計算工程に相当する。
【0110】
（ステップＳ３３）ML-EM図５ではMLACF法で推定した放射能画像λ'を利用して被検体マスク投影データmprojを算出するのに、図５のステップＳ２４を実施していた。図６では、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データmprojを算出するために、図６のステップＳ３３，Ｓ３４を実施する。
【0111】
先ず、図６では、MLACF法とは異なる再構成アルゴリズム（例えばML-EM法）における計測データに関する評価関数の最適化に基づいて、放射能画像を推定する。図６ではMLACF法で推定した放射能画像λ'と区別して、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像をλ2'とする。【0112】
（ステップＳ３４）Projection + Binarization ProcessingMLACF法とは異なる再構成アルゴリズムにおける計測データに関する評価関数の最適化に基づいて推定した放射能画像λ2'の線積分データ（投影データ）を算出する（Projection）。そして、閾値処理によって放射能画像λ2'の投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ３３，Ｓ３４は、本発明におけるマスク算出工程に相当する。
【0113】
（ステップＳ３５）ML-EM図６のステップＳ３５は、図５のステップＳ２５と同じであるので、その説明については省略する。ステップＳ３５は、本発明におけるオフセット推定工程に相当する。
【0114】
（ステップＳ３６）Extraction図６のステップＳ３６は、図５のステップＳ２６と同じであるので、その説明については省略する。ステップＳ３６は、本発明における参照領域抽出工程に相当する。
【0115】
（ステップＳ３７）αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)図６のステップＳ３７は、図５のステップＳ２７と同じであるので、その説明については省略する。ステップＳ３７は、本発明における係数算出工程に相当する。
【0116】
（ステップＳ３８）μ=μ'+α×μoff図６のステップＳ３８は、図５のステップＳ２８と同じであるので、その説明については省略する。ステップＳ３８は、本発明における吸収係数値補正工程に相当する。
【0117】
本実施例３に係る吸収係数画像推定方法によれば、上述した実施例１，２と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、図５のステップＳ２８，または図６のステップＳ３８で吸収係数値を補正する。したがって、図５のステップＳ２８，または図６のステップＳ３８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0118】
本実施例３では、上述したマスク算出工程に相当する図５のステップＳ２４，または図６のステップＳ３３，Ｓ３４を実施している。すなわち、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像（図５ではλ'，図６ではλ2'）を算出する工程（図５ではステップＳ２１，図６ではステップＳ３３）と、放射能画像（図５ではλ'，図６ではλ2'）の投影データを算出する工程（図５ではステップＳ２４の前半の「Projection」，図６ではステップＳ３４の前半の「Projection」）と、（放射能画像の）投影データを二値化したデータを被検体マスク投影データmprojとして算出する工程（図５ではステップＳ２４の後半の「Binarization Processing」，図６ではステップＳ３４の後半の「Binarization Processing」）とからなる。本実施例３の場合には、放射能画像を算出する再構成アルゴリズムについては、特に限定されない。図５のように、上述した同時再構成アルゴリズム（MLACF法やMLAA法）を用いた場合には、同時再構成アルゴリズムで推定した上述の放射能画像λ'を利用して、被検体マスク投影データmprojを算出することができる。また、図６のように、上述した同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した放射能画像λ2'を利用して、被検体マスク投影データmprojを算出してもよい。【0119】
本実施例３に係る吸収係数画像推定プログラム６（図１を参照）や本実施例３に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１，２と同じであるので、その説明については省略する。【実施例4】
【0120】
次に、図面を参照して本発明の実施例４を説明する。図７は、実施例４に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。【0121】
上述した実施例１〜３では、MLACF法と、吸収係数投影データ（各実施例１〜３では吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムとの組み合わせで、上述した再構成計算工程を実施していた。これに対して、本実施例４では、(a)画像μ'を未知数に含む計算アルゴリズム（MLAA法）で実施する。【0122】
（ステップＳ４１）MLAA上述した実施例１〜３では、画像μ'および放射能画像λ'を同時に計算するために、MLACF法によって、放射能画像λ'および吸収係数サイノグラムA'を推定したステップＳ１（実施例２ではステップＳ１１，実施例３ではステップＳ２１またはＳ３１）の後に、吸収係数サイノグラムA'を再構成した画像を画像μ'とするステップＳ２（実施例２ではステップＳ１２，実施例３ではステップＳ２２またはＳ３２）を実施していた。本実施例４では、MLAA法によって、放射能画像λ'および画像μ'を同時に計算する。
【0123】
つまり、MLACF法を用いた２つのステップＳ１・Ｓ２（実施例２ではステップＳ１１・Ｓ１２，実施例３ではステップＳ２１・Ｓ２２またはＳ３１・Ｓ３２）が、本実施例４では、MLAA法を用いた場合には１つのステップＳ４１に統合される。MLAA法の具体的な手法については、参考文献３を参照されたい（参考文献３：A. Rezaei (K.U. Leuven), M. Defrise, G. Bal et. al., “Simultaneous reconstruction of activity and attenuation in time-of-flight PET”, IEEE Trans. Med. Imag., 31 (12), 2224-2233, 2012）。ステップＳ４１は、本発明における再構成計算工程に相当する。【0124】
（ステップＳ４３）Binarization Processing図７のステップＳ４３は、上述した実施例１のステップＳ３と同じであるので、その説明については省略する。
【0125】
（ステップＳ４４）Projection + Binarization Processing図７のステップＳ４４は、上述した実施例１のステップＳ４と同じであるので、その説明については省略する。ステップＳ４３，Ｓ４４は、本発明におけるマスク算出工程に相当する。
【0126】
（ステップＳ４５）ML-EM図７のステップＳ４５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５，上述した実施例３のステップ（図５ではステップＳ２５，図６ではステップＳ３５）と同じであるので、その説明については省略する。ステップＳ４５は、本発明におけるオフセット推定工程に相当する。
【0127】
（ステップＳ４６）Extraction図７のステップＳ４６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６，上述した実施例３のステップ（図５ではステップＳ２６，図６ではステップＳ３６）と同じであるので、その説明については省略する。ステップＳ４６は、本発明における参照領域抽出工程に相当する。
【0128】
（ステップＳ４７）αn=(μnknown-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)図７のステップＳ４７は、上述した実施例１のステップＳ７，上述した実施例２のステップＳ１７，上述した実施例３のステップ（図５ではステップＳ２７，図６ではステップＳ３７）と同じであるので、その説明については省略する。ステップＳ４７は、本発明における係数算出工程に相当する。
【0129】
（ステップＳ４８）μ=μ'+α×μoff図７のステップＳ４８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８，上述した実施例３のステップ（図５ではステップＳ２８，図６ではステップＳ３８）と同じであるので、その説明については省略する。ステップＳ４８は、本発明における吸収係数値補正工程に相当する。
【0130】
なお、図７は、上述した実施例１の図３におけるステップＳ１・Ｓ２をステップＳ４３に統合した場合に適用したときの実施例４のフローチャートである。もちろん、上述した実施例２の図４におけるステップＳ１１・Ｓ１２をステップＳ４３に統合した場合に実施例４を適用してもよいし、上述した実施例３の図５におけるステップＳ２１・Ｓ２２または図６におけるステップＳ３１・Ｓ３２に統合した場合に実施例４を適用してもよい。【0131】
本実施例４に係る吸収係数画像推定方法によれば、上述した実施例１〜３と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ４８で吸収係数値を補正する。したがって、ステップＳ４８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0132】
本実施例４では、上述した再構成計算工程に相当するステップＳ４１を、(a)画像μ'を未知数に含む計算アルゴリズムで実施している。(a)のアルゴリズムは、非定量的な放射能画像λ'および画像μ'（非定量的な吸収係数画像）を同時に再構成するMLAA法である。【0133】
本実施例４に係る吸収係数画像推定プログラム６（図１を参照）や本実施例４に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１〜３と同じであるので、その説明については省略する。【実施例5】
【0134】
次に、図面を参照して本発明の実施例５を説明する。図８は、実施例５に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。【0135】
上述した実施例１〜４では、代表値に基づいて係数αを算出していた。これに対して、本実施例５では、誤差評価関数に基づいて係数αを算出する。【0136】
（ステップＳ５１）MLACF図８のステップＳ５１は、上述した実施例１のステップＳ１，上述した実施例２のステップＳ１１，上述した実施例３のステップ（図５ではステップＳ２１，図６ではステップＳ３１）と同じであるので、その説明については省略する。
【0137】
（ステップＳ５２）ML-TR or ML-EM図８のステップＳ５２は、上述した実施例１のステップＳ２，上述した実施例２のステップＳ１２，上述した実施例３のステップ（図５ではステップＳ２２，図６ではステップＳ３２）と同じであるので、その説明については省略する。ステップＳ５１，Ｓ５２は、本発明における再構成計算工程に相当する。
【0138】
（ステップＳ５３）Binarization Processing図８のステップＳ５３は、上述した実施例１のステップＳ３，上述した実施例４のステップＳ４３と同じであるので、その説明については省略する。
【0139】
（ステップＳ５４）Projection + Binarization Processing図８のステップＳ５４は、上述した実施例１のステップＳ４，上述した実施例４のステップＳ４４と同じであるので、その説明については省略する。ステップＳ５３，Ｓ５４は、本発明におけるマスク算出工程に相当する。
【0140】
（ステップＳ５５）ML-EM図８のステップＳ５５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５，上述した実施例３のステップ（図５ではステップＳ２５，図６ではステップＳ３５），上述した実施例４のステップＳ４５と同じであるので、その説明については省略する。ステップＳ５５は、本発明におけるオフセット推定工程に相当する。
【0141】
（ステップＳ５６）Extraction図８のステップＳ５６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６，上述した実施例３のステップ（図５ではステップＳ２６，図６ではステップＳ３６），上述した実施例４のステップＳ４６と同じであるので、その説明については省略する。ステップＳ５６は、本発明における参照領域抽出工程に相当する。
【0142】
（ステップＳ５７）f(α)=DΩ(μknown,μ'+α×μoff)上述した実施例１〜４では、代表値に基づいて係数αを算出していた。本実施例５では、誤差評価関数に基づいて係数αを算出する。なお、K(≧1)を既知の吸収係数値で近似可能な領域数とし、wn (n=1,…,K)を０以上１以下の係数としたときに、図８では、K=1（つまり、n=1のみ）として領域数を１つのみとし、係数w1=1として説明する。一般化した式については、後述する。
【0143】
具体的には、μknownを領域Ω内に既知の吸収係数を設定した画像とし、DΩ(X,Y)を領域Ω内に関する画像Xおよび画像Yの誤差評価関数とする。画像Xを領域Ω内に既知の吸収係数を設定した画像μknownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩ(μknown,μ'+α×μoff)は、領域Ω内に関する既知の吸収係数を設定した画像μknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。f(α)を、係数αを変数とした関数とすると、関数f(α)は下記（４）式のように表される。f(α)=DΩ(μknown,μ'+α×μoff) …（４）
【0144】
上記（４）式で表された関数f(α)を最小化するαを算出する。誤差評価関数DΩ(X,Y)は、例えば、二次関数DΩ(X,Y)=Σn∈Ω(X-Y)2である。もちろん、画像Xおよび画像Yの誤差を評価する関数であれば、絶対値の差分和DΩ(X,Y)=Σn∈Ω|X-Y|に例示されるように、誤差評価関数DΩ(X,Y)については特に限定されない。ステップＳ５７は、本発明における係数算出工程に相当する。【0145】
（ステップＳ５８）μ=μ'+α×μoff図８のステップＳ４８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８，上述した実施例３のステップ（図５ではステップＳ２８，図６ではステップＳ３８），上述した実施例４のステップＳ４８と同じであるので、その説明については省略する。ステップＳ５８は、本発明における吸収係数値補正工程に相当する。
【0146】
なお、図８は、上述した実施例１の図３におけるステップＳ７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に適用したときの実施例５のフローチャートである。もちろん、上述した実施例２の図４におけるステップＳ１７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよいし、上述した実施例３の図５におけるステップＳ２７または図６におけるステップＳ３７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよいし、上述した実施例４の図７におけるステップＳ４７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよい。【0147】
本実施例５に係る吸収係数画像推定方法によれば、上述した実施例１〜４と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ５８で吸収係数値を補正する。したがって、ステップＳ５８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。【0148】
本実施例５では、誤差評価関数に基づいて係数αを算出している。上記（４）式を一般化した場合について説明する。K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとし、μnknown(n=1,…,K)を領域Ωn内に既知の吸収係数を設定した画像とし、DΩn(X,Y)を領域Ωn内に関する画像Xおよび画像Yの誤差評価関数とし、wn(n=1,…,K)を０以上１以下の係数とする。画像Xを領域Ωn内に既知の吸収係数を設定した画像μnknownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩn(μnknown,μ'+α×μoff)は、領域Ωn内に関する既知の吸収係数を設定した画像μnknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。このとき、係数αを変数とした関数f(α)は、n=1,…,Kで各領域Ω1,…,ΩK毎の係数w1,…,wKによるDΩ1(μ1known,μ'+α×μoff),…,DΩK(μKknown,μ'+α×μoff)の重み付き加算で表される。すなわち、一般化した場合には、関数f(α)は下記（５）式のように表される。f(α)=Σn=1,…,K[wn×DΩn(μnknown,μ'+α×μoff)] …（５）
【0149】
上記（５）式で表された関数f(α)を最小化するαを算出する。【0150】
本実施例５に係る吸収係数画像推定プログラム６（図１を参照）や本実施例５に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１〜４と同じであるので、その説明については省略する。【0151】
本発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。【0152】
（１）上述した各実施例において撮影対象については、特に限定されない。被検体の全身を撮影する装置や、被検体の頭部を撮影する装置や、被検体の乳房を撮影する装置に適用してもよい。【0153】
（２）上述した各実施例では、ＤＯＩ検出器であったが、深さ方向を弁別しない検出器に適用してもよい。【0154】
（３）上述した各実施例では、検出器リングを被検体の体軸方向に積層配置する構成であったが、検出器リングが１つのみの構成であってもよい。【0155】
（４）上述した各実施例では、データ形式がサイノグラムのときであったが、サイノグラムに限定されない。投影データであれば、例えばデータ形式がヒストグラムであってもよい。したがって、上述した実施例１〜３では、MLACF法で推定した吸収係数投影データは吸収係数サイノグラムA'であったが、MLACF法で推定した吸収係数投影データは吸収係数ヒストグラムであってもよい。【0156】
（５）上述した実施例１では、画像μ'を先に二値化した後に投影データを算出し、その投影データを二値化したが、下記のように被検体マスク投影データmprojを算出してもよい。すなわち、マスク算出工程は、画像μ'の投影データを算出する工程と、画像μ'の投影データの二値化データを被検体マスク投影データmprojとして算出する工程とからなる態様でもよい。この態様の場合には、画像μ'の投影データを先に算出した後に二値化している。【0157】
（６）上述した実施例３では、放射能画像（図５ではλ'，図６ではλ2'）の投影データを先に算出した後に二値化していたが、下記のように被検体マスク投影データmprojを算出してもよい。すなわち、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の二値化画像を算出する工程と、二値化画像の投影データを算出する工程と、二値化画像の投影データを二値化したデータを被検体マスク投影データmprojとして算出する工程とからなる態様でもよい。この態様の場合には、放射能画像を先に二値化した後に投影データを算出し、その投影データを二値化している。上述した実施例３と同様に、この態様の場合においても、放射能画像を算出する再構成アルゴリズムについては、特に限定されず、図５のような同時再構成アルゴリズム（MLACF法やMLAA法）あるいは、図６のような同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）を用いてもよい。【産業上の利用可能性】
【0158】
以上のように、本発明は、TOF計測型ＰＥＴ装置で計測された計測データを用いた吸収係数画像の推定に適している。【符号の説明】
【0159】
１ … ＰＥＴ装置３ … γ線検出器
５ … 演算回路
６ … 吸収係数画像推定プログラム
λ'，λ2' … 放射能画像
A' … 吸収係数サイノグラム
μ' … （定量的な吸収係数画像に対して不均一なオフセット値が加算された）画像
mimg … 被検体マスク画像
mproj … 被検体マスク投影データ
μoff … （不均一な）オフセット画像
Ω … （既知の吸収係数値で近似可能な）領域
K … 既知の吸収係数値で近似可能な領域数
Ωn … n番目の領域Ω
μnknown … 領域Ωnの既知の吸収係数値
S(μ';Ωn) … 領域Ωnにおける画像μ'の代表値
S(μoff;Ωn) … 領域Ωnにおけるオフセット画像μoff'の代表値
α … 係数
αn … 領域Ωnにおける係数α
T(α1,α2,…,αK) … 係数α1,α2,…,αKの代表値
μ … 真の吸収係数画像の値（吸収係数値）
DΩ(μknown,μ'+α×μoff) … 領域Ω内に既知の吸収係数を設定した画像μknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数
DΩn(μnknown,μ'+α×μoff) … 領域Ωn内に既知の吸収係数を設定した画像μnknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数
wn … （０以上１以下の）係数
f(α) … （係数αを変数とした）関数