JP2004309319A - Image reconstruction processing program and recording medium - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide an image reconstruction processing program capable of performing brain image reconstruction without performing subtraction, when using a twice (or more than twice) divided dose method. <P>SOLUTION: The S/N ratio of a load image in the Tc-99mECD twice divided dose method can be improved by integrating influence of intracerebral residual radioactivity into a calibration expression of a sequential image reconfiguration method. Actually, the influence of the intracerebral residual radioactivity is integrated by adding the total count number y<SB>i</SB><SP>[h-1]</SP>until the preceding time to the denominator of the expression 14. For example, in the case of the twice split dosage, the first-time count number y<SB>i</SB><SP>[1]</SP>is added to the denominator of the expression 14. Namely, when reconfiguring SPECT data in the second time, the brain image reconfiguration can be performed without performing subtraction, by supposing that the intracerebral residual radioactivity y<SB>i</SB><SP>[1]</SP>by first-time dosage exists as a base line in the denominator of the expression 14. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

ここで、λ_ｊ ^{０ ［１］}は前記第１回推定ＲＩ濃度初期値記録部に設定され、ｃ_ｉｊは前記検出確率記録部に記録され、ｙ_ｉ ^［１］は前記第１回測定データ記録部に記録されたものであり、
ｈ＝２からＬまで、第ｈ回目の分割投与後の再構成画像を、ｋ回目の逐次近似における再構成画像上の画素ｊの推定ＲＩ濃度λ_ｊ ^{ｋ ［ｈ］}を式（２）を用いてｋ＝０から所定の第２収束条件が成立するまで繰返し計算することにより求める第ｈ回再構成画像生成ステップ、
【数９】

ここで、λ_ｊ ^{０ ［ｈ］}は前記第ｈ回推定ＲＩ濃度初期値記録部に設定され、ｙ_ｉ ^{［ ｈ］}は前記第ｈ回測定データ記録部に記録されたものであるというステップを実行させるための画像再構成処理プログラムである。
【００１０】
ここで、この発明の画像再構成処理プログラムにおいて、ｈ＝１からＬについて、Ｓ_ｉ ^［ｈ］を第ｈ回目の分割投与後における投影データ上の１画素ｉにおける散乱補正項とし、ここでＳ_ｉ ^［ｈ］はｉ＝１から前記画素数記録部に記録された画素ｉの数ｎまで第ｈ回散乱補正項記録部に記録されており、
【数１０】

前記第１回再構成画像生成ステップにおいて式（１）の替わりに式（３）を用い、前記第ｈ回再構成画像生成ステップにおいて式（２）の替わりに式（４）を用いることができる。
【００１１】
ここで、この発明の画像再構成処理プログラムにおいて、前記所定の第１収束条件は、式（５）で求められる値が、
【数１１】

所定の値以下となる条件とすることができる。
【００１２】
ここで、この発明の画像再構成処理プログラムにおいて、前記所定の第２収束条件は、式（６）で求められる値が、
【数１２】

所定の値以下となる条件とすることができる。
【００１３】
ここで、この発明の画像再構成処理プログラムにおいて、前記所定の第１収束条件は、式（７）で求められる値が、
【数１３】

所定の値以下となる条件とすることができる。
【００１４】
ここで、この発明の画像再構成処理プログラムにおいて、前記所定の第２収束条件は、式（８）で求められる値が、
【数１４】

所定の値以下となる条件とすることができる。
【００１５】
ここで、この発明のいずれかの画像再構成処理プログラムにおいて、前記検出確率記録部に記録されたｃ_ｉｊは、再構成画像上の画素ｊの投影方向に投影データ上の画素ｉが存在する場合は１とし、他の場合は０とすることができる。
【００１６】
ここで、この発明のいずれかの画像再構成処理プログラムにおいて、前記第１回推定ＲＩ濃度初期値記録部に記録されたλ_ｊ ^{０ ［ｈ］}は、予めｊ＝１から前記画素数記録部に記録された画素ｊの数ｍまで１と記録しておくことができる。
【００１７】
この発明の記録媒体は、本発明のいずれかの画像再構成処理プログラムを記録したコンピュータ読取り可能な記録媒体である。
【００１８】
【発明の実施の形態】
以下、本発明の実施の形態について図面を参照して詳細に説明する。
【００１９】
実施の形態１．
本発明の画像再構成処理プログラムは、逐次近似型画像再構成法に前回（例えば１回）までのＥＣＤ投与によるベースラインの影響を予め組み込んだアルゴリズムを用いることが特徴である。そこで、まず逐次近似型画像再構成法について説明する。
【００２０】
未知数である体内のＲＩ分布を、収集された投影データから統計的手法により最も可能性の高いケースとして推定する方法として、最尤推定−期待値最大化（Ｍａｘｉｍｕｍ Ｌｉｋｅｌｉｈｏｏｄ − Ｅｘｐｅｃｔａｔｉｏｎ Ｍａｘｉｍｉｚａｔｉｏｎ ： ＭＬ−ＥＭ）再構成法がある。このＭＬ−ＥＭ法は最尤推定法と期待値最大化法とを組合せた方法である。ＭＬ−ＥＭ法は、推定した体内のＲＩ分布から計算した投影データを実際に収集した投影データにできるだけ近づけるように、逐次近似法を用いて体内のＲＩ分布を繰返し修正していく方法である。
【００２１】
図１は、最尤推定法の画像再構成への適用を説明するものである。図１において、符号１０は脳画像の断層画像、３は断層画像１０のマトリクスにおけるＲＩ濃度がλ（ＭＢｑ）の１画素、３０は検出器であり投影データ２０上の１画素を示す。ｃは画素３０で光子が検出される検出確率を示し、ｘは投影データ２０上の１画素３０（または画素を検出器と考えて、１検出器３０）における画素３からの光子の計数を示す。ＲＩ濃度がλ（ＭＢｑ）のＲＩ（画素３）から放出された光子を検出器３０（画素３０）でｎ回測定した場合、各測定毎の計数をｘｉとすると、計数ｘの最尤推定値は式９で与えられる。
【００２２】
【数１５】

【００２３】
ＲＩ濃度λの最尤推定値は式１０で与えられる。
【００２４】
【数１６】

【００２５】
次に、ＲＩをｎ方向の検出器で測定した場合について説明する。図２は、ＲＩをｎ方向の検出器で測定した場合における最尤推定法の画像再構成への適用を説明するものである。図２で図１と同じ符号を付した箇所は同じ要素を示すため説明は省略する。図２において、符号５は断層画像１０のマトリクスにおけるＲＩ濃度がλ_ｊ（ＭＢｑ）の１画素ｊ（ｊ＝１〜断層画像１０のマトリクスにおける画素の総数ｍ）を示す。符号３１ないし３３は検出器であり各々投影データ２０等上の１画素を示す。検出器はｎ個あるが図面の都合上３個のみ示されている。ｃ_ｉｊ（ｉ＝１〜検出器の数ｎまたは投影データ上の各画素の数ｎ、ｊ＝１〜ｍ）は画素ｊ（５）から放出された光子が画素ｉ（または検出器ｉ）で検出される検出確率を示し、ｘ_ｉｊは画素ｊ（５）から放出され画素ｉ（または検出器ｉ）で検出される光子数を示す。以上の条件の場合、画素ｊ（５）の濃度λ_ｊの最尤推定値は式１１で与えられる。
【００２６】
【数１７】

【００２７】
実際には画素ｊ（５）の周辺にもＲＩが存在するため、光子数ｘ_ｉｊは直接測定することはできない。図３は、画素ｊ（５）の周辺のＲＩを考慮した場合における最尤推定法の画像再構成への適用を説明するものである。図３で図２と同じ符号を付した箇所は同じ要素を示すため説明は省略する。図３において、ｙ_ｉは検出器ｉにおける各画素ｊ（５）等からの計測値の総和（投影カウント）を示す。総和ｙ_ｉから光子数ｘ_ｉｊを推定する必要がある。ここでは最大事後確率推定（ＭＡＰ推定）を用いると、光子数ｘ_ｉｊのＭＡＰ推定値は式１２で与えられる。
【００２８】
【数１８】

【００２９】
つまり、最尤推定法により求めたい未知数である画素ｊ（５）のＲＩ濃度λ_ｊを検出確率ｃ_ｉｊと光子数ｘ_ｉｊとから推定する場合、光子数ｘ_ｉｊは直接測定できないため、光子数ｘ_ｉｊのＭＡＰ推定値により求めることになる。そこで、式１１の右辺の光子数ｘ_ｉｊに式１２の光子数ｘ_ｉｊのＭＡＰ推定値を代入すると式１１の２を得ることができる。
【００３０】
【数１９】

【００３１】
式１１の２は右辺にもＲＩ濃度λ_ｊがあるためλ_ｊについて解けない。このためＥＭ法により逐次式を作成する。λ_ｊ ^ｋをｋ回目の逐次近似における推定値とすると、式１１の２より式１３を得ることができる。
【００３２】
【数２０】

【００３３】
図４は、式１３を用いてＲＩ濃度λ_ｊを推定する計算例を示す。図４に示されるように、まず投影方向数ｎ＝２とし、断層画像のマトリクスサイズを３×３（ｍ＝３×３＝９）とし、各検出確率ｃ_ｉｊは投影方向にのみ１とし他の方向へは０と仮定する（ステップＳ１０）。つまり投影方向以外の画素ｊ（５）からは計測されないものと仮定する。図４のステップＳ１２に示されるように、真の値をｍ＝１から９へかけて１ないし９とすると、各投影方向への投影カウント数ｙ_ｉは投影方向の画素ｊの値の総和となる。例えばｍ＝１、４および７の方向では画素ｊの値は各々１、４、７であるため投影カウント数ｙ_ｉ＝１＋４＋７＝１２となる。他の投影カウント数ｙ_ｉについても同様であるため、説明は省略する。推定値についてはステップＳ１２に示されるように、ｍ＝１から９へかけてλ_１ないしλ_９となる。次に、ｋ＝０のすべてのλ_ｊ ^０を１と仮定して式１３により推定計算を行なう（ステップＳ１４）。詳細はステップＳ１６に示されるようになる。例えばλ_１ ^１については、λ_１ ^０＝１．０であり、ｃ_１１λ_１ ^０からｃ_１９λ_９ ^０までの総和に関してはｃ_１１、ｃ_１４、ｃ_１７だけが１であって他は０であり、ｙ_１ｃ_１１＝１２×１、ｙ_２ｃ_２１＝６×１であり、ｃ_１１＋ｃ_２１＝２であるため、λ_１ ^１＝３．０となる。
【００３４】
１回目の推定値としてステップＳ１８に示されるように、ｍ＝１ないし９にかけて、λ_１ ^１＝３．０、λ_２ ^１＝３．５、λ_３ ^１＝４．０、λ_４ ^１＝４．５、λ_５ ^１＝５．０、λ_６ ^１＝５．５、λ_７ ^１＝６．０、λ_８ ^１＝６．５、λ_９ ^１＝７．０と求められる。
次に、今求められたλ_ｊ ^１を用いて式１３により同様にλ_ｊ ^２を求める（ステップＳ２０）。２回目の推定値としてステップＳ２２に示されるように、ｍ＝１ないし９にかけて、λ_１ ^２＝２．２、λ_２ ^２＝２．８、λ_３ ^２＝３．３、λ_４ ^２＝４．３、λ_５ ^２＝５．０、λ_６ ^２＝５．８、λ_７ ^２＝６．４、λ_８ ^２＝７．３、λ_９ ^２＝８．１と求められる。以上の計算を指定した回数だけ繰り返す（ステップＳ２４）。具体的には所望の収束条件により繰返しの終了判定を行なう。結果として真の値（または近似値）を求めることができる（ステップＳ２６）。
【００３５】
本願発明の画像再構成処理プログラムは、上述した逐次近似型画像再構成法を用いるものである。図５および図６は、本願発明の画像再構成処理プログラムの一連のフローチャートを示す。本フローチャートでは説明の便宜上、分割投与の回数ｈを２回とする。図５は第１回目（ｈ＝１）の分割投与後の再構成画像の生成のフローチャート（第１回再構成画像生成ステップ）であり、図６は図５中のステップＳ４０ａから先の処理、すなわち第２回目（ｈ＝２）の分割投与後の再構成画像の生成のフローチャート（第２回再構成画像生成ステップ）である。図５に示されるように、まず再構成画像上の画素の総数ｍと、投影データ上の画素の総数ｎとを画素数記録部４０（後述。図１８参照）から読み出す（ステップＳ３０ａ）。再構成画像上の画素ｊ＝１とし（ステップＳ３２ａ）、上述の繰返しの回数ｋ＝０とする（ステップＳ３４ａ）。再構成画像上の画素ｊの推定ＲＩ濃度λ_ｊ ^{ｋ ［ １ ］}＝を第１回（ｈ＝１）推定ＲＩ濃度初期値記録部５０−１（後述。図１８参照）から読み出す（ステップＳ３６ａ）。投影データ上の画素ｉ＝１からｎについて、検出確率ｃ_ｉｊを検出確率記録部４２（後述。図１８参照）から適宜読み出し、投影データ上の画素ｉで検出された光子のカウント数ｙ_ｉ ^［１］を第１（ｈ＝１）回測定データ記録部６０−１（後述。図１８参照）から適宜読み出して、上述の例のように式１３を計算する（ステップＳ３８ａ）。ただし、分割投与の回数ｈ＝１を考慮して式１３は式１３の２のようになっている。式１３の２で［１］はｈ＝１を意味する。
【００３６】
【数２１】

【００３７】
所望の収束条件（第１収束条件）が成立したかどうかを判断し（ステップＳ４０ａ）、成立していない場合は、繰返しの回数ｋを１増加し、式１３で求められたλ_ｊ ^{１ ［ １ ］}を用いてステップＳ３８ａと同様に新たなｋにつき式１３を計算する。ステップＳ４０ａで収束条件が成立した場合は、第１回目（ｈ＝１）の分割投与後の再構成画像の生成は終了し、次にＡへ進む。
【００３８】
図６に示されるように、まず再構成画像上の画素の総数ｍと、投影データ上の画素の総数ｎとを画素数記録部４０（後述。図１８参照）から読み出す（ステップＳ３０ｂ）。このステップＳ３０ｂはステップＳ３０ａと同様であるため省略してもよい。再構成画像上の画素ｊ＝１とし（ステップＳ３２ｂ）、上述の繰返しの回数ｋ＝０とする（ステップＳ３４ｂ）。再構成画像上の画素ｊの推定ＲＩ濃度λ_ｊ ^{ｋ ［ ２ ］}＝を第２回（ｈ＝２）推定ＲＩ濃度初期値記録部５０−２（後述。図１８参照）から読み出す（ステップＳ３６ｂ）。投影データ上の画素ｉ＝１からｎについて、検出確率ｃ_ｉｊを検出確率記録部４２（後述。図１８参照）から適宜読み出し、投影データ上の画素ｉで検出された光子のカウント数ｙ_ｉ ^{［ ２ ］}を第２（ｈ＝２）回測定データ記録部６０−２（後述。図１８参照）から適宜読み出して、以下の式１４を計算する（ステップＳ３８ｂ）。
【００３９】
【数２２】

【００４０】
本願発明の画像再構成処理プログラムの特徴は、式１４の分母中に前回までの総カウント数ｙ_ｉ ^{［ｈ−１］}を加えている点にある。したがって２回分割投与の場合、ｈ＝１〜２となるため、式１４の分母中に１回目のカウント数ｙ_ｉ ^［１］が加えられる。すなわち、２回目のＳＰＥＣＴデータを再構成する際に、式１４の分母に１回目の投与による脳内の残留放射能ｙ_ｉ ^［１］がベースラインとして存在しているものと考える。３回の分割投与を行なう場合はｈ＝１〜３となるため、式１４の分母中に２回目のカウント数ｙ_ｉ ^［２］が加えられる。ｈ回の分割投与を行なう場合は、分母中にｈ−１回目の総カウント数ｙ_ｉ ^{［ｈ−１］}が加えられる。
【００４１】
次に、所望の収束条件（第２収束条件）が成立したかどうかを判断し（ステップＳ４０ｂ）、成立していない場合は、繰返しの回数ｋを１増加し、式１４で求められたλ_ｊ ^ｋ［ｈ］を用いてステップＳ３８ｂと同様に新たなｋにつき式１４を計算する。ステップＳ４０ｂで収束条件が成立した場合は、第２回目（ｈ＝２）の分割投与後の再構成画像の生成は終了する。
【００４２】
上述の第１収束条件として、式１３の２の分子から分母を引いた差の２乗和（２乗誤差：χ^２）を用いることができる。例えば第１の収束条件として式１５を計算し、式１５のχ^２の値が所定の値以下になった時（理想的には０になった時）に、繰返し計算を終了させることができる。
【００４３】
【数２３】

【００４４】
上述の第２収束条件として、式１４の分子から分母を引いた差の２乗和（２乗誤差：χ^２）を用いることができる。例えば第２の収束条件として式１６を計算し、式１６のχ^２値が所定の値以下になった時（理想的には０になった時）に、繰返し計算を終了させることができる。
【００４５】
【数２４】

【００４６】
検証．
本願発明の画像再構成処理プログラムの有用性を検証するため、モンテカルロシミュレーションを用いた。モンテカルロシミュレーションでは、（１）「シミュレーション体系の定義」として、楕円球、円柱または楕円柱等の幾何形状を組み合わせることにより人体の体系の区画を定義する。人体形状は数値ファントム化されている場合が多い。（２）「物質データの定義」として所定の元素または水等の分子等を用い、物質に固有の線減弱係数μを仮定する。（３）「線源の定義」として光子、電子等を定義する。（４）「計測体系の定義」として所望の計測体系を定義する。ＳＰＥＣＴではγ線の入射方向を制限するために検出器の前面にコリメータ（ｃｏｌｌｉｍａｔｏｒ）を装着している。これらの装置もモンテカルロシミュレーションの中に組み込まれる。
【００４７】
モンテカルロシミュレーションでは、１４０ＫｅＶの光子発生に対して、吸収・散乱、コリメータによる分解能の劣化、および統計ノイズをシミュレートした。コリメータは低エネルギーコリメータ（例えばＬＥＨＲコリメータ）を想定した。マトリクスサイズは１２８×１２８（ｍ＝１６３８４）とし、投影方向数ｎ＝１２０とした。ピクセルサイズは３．１２５ｍｍである。光電ピークを収集するためのメインウィンドウは１４０ＫｅＶを中心とした２０％とし、メインウィンドウの下に設定するサブウィンドウは１２４ＫｅＶを中心とした５％とした。数値ファントムは単純な円柱ファントム（１８ｃｍφ）を使用し、一様吸収体（μ＝０．１５ｃｍ^−１）と仮定した。
【００４８】
投影データは、総カウント数が２倍毎に増加するように８種類（１．３６、２．７０、５．４０、１０．８、２１．６、４３．３、８６．５Ｍｃｏｕｎｔｓ／ｓｌｉｃｅ）作成し、半減期による減衰はないものとした。本シミュレーションでは１回目と２回目との投与量が１：１になるように、２倍毎に増加したデータをペアにしてできる７組の投影データを用いて、上述の式１３の２、式１４により画像再構成を行った。繰返し回数ｋは３回、ｓｕｂｓｅｔ＝２０とした。前処理フィルターはＢｕｔｔｅｒｗｏｒｔｈフィルターを用い、カットオフ周波数ｆｃを０．３６［ｃｙｃｌｅｓ／ｐｉｘｅｌ］、次数を８とした。吸収補正は仮定したμマップ（線減弱係数μを２次元の画像にしたもの）を用いて検出確率に組み込んだ。散乱成分はＴｒｉｐｌｅ
Ｅｎｅｒｇｙ Ｗｉｎｄｏｗ （ＴＥＷ）法を用いて推定した。
【００４９】
従来の方法（Ｓｕｂｔｒａｃｔｉｏｎ法という）と本願発明の画像再構成処理プログラム用いられている方法（Ａｄｄｉｔｉｏｎ法という）との比較は、視認検査（Ｖｉｓｕａｌ ｉｎｓｐｅｃｔｉｏｎ）による評価に加え、円柱ファントムでは中心に設定したＲＯＩ（５ｃｍφ）内の％ＣＯＶ（＝標準偏差／ＲＯＩ値 ×１００％）を計算することによって行った。
【００５０】
最適投与量比の評価．
上述した円柱ファントムシミュレーションによって、％ＣＯＶが１回目と２回目とで同一になる最適投与量比をＳｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法との両方で決定した。２回目と１回目との投与量比を０．５、１．０、１．５、２．０、２．５、３．０、３．５、４．０の８種類とし、投与量の総量が常に一定になるようにした。今回の場合、２回目の投影データはどの場合でも５．４０Ｍｃｏｕｎｔｓ／ｓｌｉｃｅとした。再構成条件、および、データ処理は上記と全く同じである。
【００５１】
求められた最適投与量比の条件下で、ディジタル化されたＨｏｆｆｍａｎファントムの１スライスを用いてモンテカルロシミュレーションを行い、Ｓｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法とで再構成した。灰白質と白質の放射能濃度比を４：１と仮定し、一様吸収体（μ＝０．１５ｃｍ^−１）を仮定した。
【００５２】
ノイズを含んだベースラインの影響の評価．
上述したように、Ａｄｄｉｔｉｏｎ法（式１３の２または１４）の分母において、ノイズを含んだベースラインが存在する場合は、ＭＬ法による収束性は保証されない。よって、どの程度の誤差を持つのかを検討するために、ディジタル化されたＨｏｆｆｍａｎファントムを用いて、モンテカルロシミュレーションを行った。１回目（ベースライン）および２回目の総光子数は、投与量比を１：１とし、それぞれの総カウントを６．８４Ｍｃｏｕｎｔｓ／ｓｌｉｃｅ、１３．７Ｍｃｏｕｎｔｓ／ｓｌｉｃｅとした。式１３の２または１４による再構成を行なう時のベースラインとして、ノイズフリーの理想的なベースラインを用いた場合と比較するために、３２倍の光子数の投影データを作成し規格化して用いた。
【００５３】
臨床データ．
３例の臨床データで評価を行った。ＳＰＥＣＴ測定は東芝ＧＣＡ−７２００で、コリメータはＬＥＧＰコリメータを用いた。マトリクスサイズは６４×６４（ｍ＝４０９６）とし、投影方向数ｎ＝６０とした。ピクセルサイズは４．３ｍｍである。メインウィンドウは１４０ＫｅＶを中心とした２０％とし、メインウィンドウの下に設定するサブウィンドウは１２２ＫｅＶを中心とした８％とした。１回目のＥＣＤ投与（３３３ＭＢｑ）から、５分後に７．５分間のＳＰＥＣＴ測定を行った。１回目のＳＰＥＣＴ測定直後に２回目のＥＣＤ投与（２２２ＭＢｑ）を行い、最初の投与から１２．５分後から２回目のＳＰＥＣＴ測定を同一条件で行った。ＡＣＺ
（１０００ｍｇ）の負荷は、１回目のＳＰＥＣＴ直前に行った。収集条件・プロトコルの詳細な情報は省略する。
【００５４】
データ処理はＳｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法とで画像再構成を行った。繰返し回数ｋは３回、ｓｕｂｓｅｔ＝１５とした。前処理フィルターはＢｕｔｔｅｒｗｏｒｔｈフィルター（カットオフ周波数ｆｃ＝０．２４ｃｙｃｌｅｓ／ｐｉｘｅｌ、次＝８）を用いた。吸収補正は、閾値による輪郭抽出を行って一様なμ値（０．１４ｃｍ^−１）を仮定し、検出確率に組み込んで行った。散乱成分の推定は、シミュレーションと同様にＴＥＷ法によって行った。ＲＯＩ解析によりＳｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法との比較を行った。
【００５５】
結果．
図７（Ａ）、（Ｂ）は、円柱ファントムシミュレーション（総カウント：１．３６Ｍ ｃｏｕｎｔｓ／ｓｌｉｃｅ）による収束特性を示すグラフである。図７（Ａ）、（Ｂ）において、横軸は繰返し回数ｋ、縦軸は上述の２乗誤差（χ^２：１回目は式１５、２回目は式１６による）である。Ｓｕｂｔｒａｃｔｉｏｎ法でもＡｄｄｉｔｉｏｎ法でも最終的な誤差の絶対値はほぼ同じであったが、収束特性には明らかな違いが見られた。すなわち、Ａｄｄｉｔｉｏｎ法では、２回目のＳＰＥＣＴデータの収束速度がやや遅くなることが明らかになった。しかし、ｓｕｂｓｅｔと繰返し回数ｋの積が２０以上では、両者はほぼ同じであり、今回用いた再構成条件では特性の違いはないものと考えられた。
【００５６】
図８は、投与量比が１：１の場合の円柱ファントムシミュレーション（２．７０Ｍｃｏｕｎｔｓ／ｓｌｉｃｅ）による再構成画像の比較である。Ｓｕｂｔｒａｃｔｉｏｎ法（図８で左側）では、２回目（２ｎｄ）の再構成画像のＳ／Ｎ比が悪いが、一方Ａｄｄｉｔｉｏｎ法（図８で右側）では著しく改善され、その結果として１回目（１ｓｔ）とほぼ同じ画質になっている。
【００５７】
図９（Ａ）、（Ｂ）は、投与量比が１：１の場合の総カウント数に対する％ＣＯＶの変化を示したグラフである。図９（Ａ）、（Ｂ）において、横軸はスライス（ｓｌｉｃｅ）当たりの総カウント数、縦軸は％ＣＯＶである。図９（Ａ）に示されるＳｕｂｔｒａｃｔｉｏｎ法の％ＣＯＶは、１回目に比べ２回目が明らかに悪く（総カウント２．７０Ｍｃｏｕｎｔｓ／ｓｌｉｃｅの時、１回目で５０．４％、２回目で８０．５％）、明らかにサブトラクションによる影響がみられる。一方、図９（Ｂ）に示されるＡｄｄｉｔｉｏｎ法では、１回目と２回目の＆￥％ＣＯＶはほぼ同じであり（総カウント２．７０Ｍｃｏｕｎｔｓ／ｓｌｉｃｅの時、１回目で４３．１％、２回目で４１．４％）、本発明のＡｄｄｉｔｉｏｎ法の有用性が裏付けられた。Ｓｕｂｔｒａｃｔｉｏｎ法の１回目とＡｄｄｉｔｉｏｎ法の１回目とは、散乱補正が施されていない場合は全く一致するが、本シミュレーションでは散乱補正をしているためにＡｄｄｉｔｉｏｎ法の方が若干良くなっている。
【００５８】
一方、図１０は横軸をスライス（ｓｌｉｃｅ）当たりの総カウント数、縦軸をＲＯＩ値の２回目／１回目の比に取った場合の総カウント数に対する変化を示したグラフである。本シミュレーションでは、この比は１．０になるはずであるが、Ｓｕｂｔｒａｃｔｉｏｎ法でもＡｄｄｉｔｉｏｎ法でも、ほぼ仮定した値になった。
【００５９】
図１１（Ａ）、（Ｂ）は横軸を投与量比（２回目／１回目）（実際には発生光子数比）、縦軸を％ＣＯＶとした場合のＳｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法との比較を示す。図１１（Ａ）に示されるようにＳｕｂｔｒａｃｔｉｏｎ法において、１回目および２回目の％ＣＯＶが同一になったのは投与量比が約２．０である。一方、図１１（Ｂ）に示されるようにＡｄｄｉｔｉｏｎ法において、１回目および２回目の％ＣＯＶが同一になったのは投与量比が約１．０である。最適投与量における％ＣＯＶの絶対値は、Ｓｕｂｔｒａｃｔｉｏｎ法で約５４％、Ａｄｄｉｔｉｏｎ法では約４３％であった。最適投与量比においても、Ａｄｄｉｔｉｏｎ法の方がＳ／Ｎ比が良いという結果が得られた。これは、全体の投与量（シミュレーションでは発生光子数）を同一に保っているために、Ｓｕｂｔｒａｃｔｉｏｎ法の１回目の投与量はＡｄｄｉｔｉｏｎ法の投与量の２／３になるためである。
【００６０】
図１２は最適投与量比におけるＳｕｂｔｒａｃｔｉｏｎ法、およびＡｄｄｉｔｉｏｎ法でのＨｏｆｆｍａｎファントムの再構成画像を示す。図１２に示されるように、円柱ファントムシミュレーションの結果と同様に、Ａｄｄｉｔｉｏｎ法（図１２の右側）の方がＳ／Ｎ比が良いことがわかる。
【００６１】
図１３は、図１２に示した画像の１回目と２回目との各画素値の散布図を示す。図１３（Ａ）、（Ｂ）で横軸は１回目ＳＰＥＣＴカウントであり、縦軸は２回目ＳＰＥＣＴカウントである。図１３（Ａ）に示されるように相関係数はＳｕｂｔｒａｃｔｉｏｎ法でｒ＝０．７６２となり、図１３（Ｂ）に示されるようにＡｄｄｉｔｉｏｎ法でｒ＝０．８１３となり、Ｓｕｂｔｒａｃｔｉｏｎ法のばらつきが大きかった。
【００６２】
図１４はＡｄｄｉｔｉｏｎ法の２回目再構成において、１回目の脳内残留放射能のベースラインに、１）ノイズを含んだベースラインの場合（図１４の左側）、２）ノイズフリーの理想的なベースラインを用いた場合（図１４の右側）の２種類の再構成画像を示す。視覚的評価ではほぼ両者は一致している。
【００６３】
図１５は、図１４の両者の画素値を比較した散布図を示す。図１５において横軸は図１４のノイズフリーの理想的なベースライン、縦軸は図１４のノイズを含んだベースラインである。両者に系統的な誤差は認められなかった。
【００６４】
図１６は、臨床データに適用した負荷（２回目）の再構成画像での比較を示す。図１６で左側はＳｕｂｔｒａｃｔｉｏｎ法の場合を示し、右側はＡｄｄｉｔｉｏｎ法の場合を示している。
【００６５】
図１７は３例の再構成画像にとったＲＯＩ（１５５４ＲＯＩｓ）における、２つの方法の比較を示す。図１７（Ａ）は１回目（安静時）、図１７（Ｂ）は２回目（ＡＣＺ負荷時）である。図１７（Ａ）、（Ｂ）において横軸はＲＯＩカウント（Ａｄｄｉｔｉｏｎ法）、縦軸はＲＯＩカウント（Ｓｕｂｔｒａｃｔｉｏｎ法）である。視覚的に評価するとＳｕｂｔｒａｃｔｉｏｎ法はノイジーな画像であったが、ＲＯＩ値の比較では、やや２回目のばらつきが存在するものの、大きな違いは見られなかった。
【００６６】
以上説明したように、本発明の実施の形態１によれば、逐次型画像再構成法の計算式の中に脳内残留放射能の影響を組み込むことにより、Ｔｃ−９９ｍＥＣＤ２回分割投与法における負荷画像のＳ／Ｎ比を向上させることができる。具体的には、式１４の分母中に前回までの総和カウント数ｙ_ｉ ^{［ｈ−１］}を加えることにより脳内残留放射能の影響を組み込んでいる。例えば２回分割投与の場合、式１４の分母中に１回目のカウント数ｙ_ｉ ^［１］が加えられる。すなわち、２回目のＳＰＥＣＴデータを再構成する際に、式１４の分母に１回目の投与による脳内の残留放射能ｙ_ｉ ^［１］がベースラインとして存在しているものと考える。本発明の画像再構成処理プログラムによる方法をシミュレーションデータと臨床データとに適用して評価した。その結果、投与量比が１：１の場合でも、１回目と２回目との再構成データのＳ／Ｎ比はほぼ同一となり、従来のＳｕｂｔｒａｃｔｉｏｎ法に比べＳ／Ｎ比の向上を確認することができた。上述のシミュレーション結果によれば、理想的な投影データをベースラインとした場合と、ノイズを含んだ投影データをベースラインとして用いた場合とで、定量的には差がなかった。以上より、２回（または２回以上の）分割投与法を用いる場合に、サブトラクションを行なわずに脳画像再構成を行うことができる画像再構成処理プログラムを提供することができる。
【００６７】
実施の形態２．
上述の実施の形態１における式１３の２および式１４において、散乱補正を行なうこともできる。分割投与の回数ｈ＝１からＬについて、Ｓ_ｉ ^［ｈ］を第ｈ回目の分割投与後における投影データ上の１画素ｉにおける散乱補正項とする。ここでＳ_ｉ ^［ｈ］はｉ＝１から画素数記録部１０に記録された画素ｉの数ｎまで第ｈ回散乱補正項記録部に記録されているものとする。この場合、式１３の２の替りに式１７を用い、式１４の替りに式１８を用いればよい。
【００６８】
【数２５】

【００６９】
式１７および１８を用いた場合、散乱補正を考慮した第１収束条件は式１９のようになり、第２収束条件は式２０のようになる。
【００７０】
【数２６】

【００７１】
実施の形態３．
図１８は、上述した実施の形態を実現するための本発明のコンピュータ・プログラムを実行するコンピュータの内部回路７０を示すブロック図である。図１８に示されるように、ＣＰＵ７２、ＲＯＭ７４、ＲＡＭ７６、画像制御部７８、コントローラ８６、入力制御部９０および外部インタフェース（Ｉｎｔｅｒｆａｃｅ ： Ｉ／Ｆ）部９２はバス８２に接続されている。図１８において、上述の本発明のコンピュータ・プログラムは、ＲＯＭ７４、ディスク８４ａまたはＣＤ−ＲＯＭ８４ｎ等の記録媒体（脱着可能な記録媒体を含む）に記録されている。ディスク８４ａ等には、画素数記録部４０、検出確率記録部４２、第ｈ回推定ＲＩ濃度記録部５０−ｈ、第ｈ回測定データ記録部６０−ｈ等が記録されている。本発明のコンピュータ・プログラムは、ＲＯＭ７４からバス８２を介し、またはディスク８４ａ若しくはＣＤ−ＲＯＭ８４ｎ等の記録媒体からコントローラ８６を経由してバス８２を介しＲＡＭ７６へロードされる。画像制御部７８は脳画像データをＶＲＡＭ８０へ送出し、ディスプレイ等の表示装置である表示部４５はＶＲＡＭ８０から送出された脳画像データに基づいて脳画像を表示する。ＶＲＡＭ８０は表示部４５の一画面分のデータ容量に相当する容量を有している画像メモリである。入力操作部８８は画像処理装置４４に入力を行うためのマウス、テンキー等の入力装置であり、入力制御部９０は入力操作部８８と接続され入力制御等を行う。外部Ｉ／Ｆ部９２は、例えばインターネット等の外部の通信網（不図示）と接続する際のインタフェース機能を有している。
【００７２】
上述のようにＣＰＵ７２が本発明のコンピュータ・プログラムを実行することにより、本発明の目的を達成することができる。当該コンピュータ・プログラムは上述のようにＣＤ−ＲＯＭ８４ｎ等の記録媒体の形態でコンピュータＣＰＵ７２に供給することができ、当該コンピュータ・プログラムを記録したＣＤ−ＲＯＭ８４ｎ等の記録媒体も同様に本発明を構成することになる。当該コンピュータ・プログラムを記録した記録媒体としては上述された記録媒体の他に、例えばメモリ・カード、メモリ・スティック、ＤＶＤ、光ディスク、ＦＤ等を用いることができる。
【００７３】
【発明の効果】
以上説明したように、本発明の画像再構成処理プログラムによれば、逐次型画像再構成法の計算式の中に脳内残留放射能の影響を組み込むことにより、Ｔｃ−９９ｍＥＣＤ２回分割投与法における負荷画像のＳ／Ｎ比を向上させることができる。具体的には、式１４の分母中に前回までの総カウント数ｙ_ｉ ^{［ｈ−１］}を加えることにより脳内残留放射能の影響を組み込んでいる。例えば２回分割投与の場合、式１４の分母中に１回目のカウント数ｙ_ｉ ^［１］が加えられる。すなわち、２回目のＳＰＥＣＴデータを再構成する際に、式１４の分母に１回目の投与による脳内の残留放射能ｙ_ｉ ^［１］がベースラインとして存在しているものと考える。本発明の画像再構成処理プログラムによる方法をシミュレーションデータと臨床データとに適用して評価した結果、投与量比が１：１の場合でも、１回目と２回目との再構成データのＳ／Ｎ比はほぼ同一となり、従来のＳｕｂｔｒａｃｔｉｏｎ法に比べＳ／Ｎ比の向上を確認することができた。上述のシミュレーション結果によれば、理想的な投影データをベースラインとした場合と、ノイズを含んだ投影データをベースラインとして用いた場合とで、定量的には差がなかった。以上より、２回（または２回以上の）分割投与法を用いる場合に、サブトラクションを行なわずに脳画像再構成を行うことができる画像再構成処理プログラムを提供することができる。
【図面の簡単な説明】
【図１】最尤推定法の画像再構成への適用を説明する図である。
【図２】ＲＩをｎ方向の検出器で測定した場合における最尤推定法の画像再構成への適用を説明する図である。
【図３】画素ｊ（５）の周辺のＲＩを考慮した場合における最尤推定法の画像再構成への適用を説明する図である。
【図４】式１３を用いてＲＩ濃度λ_ｊを推定する計算例を示す図である。
【図５】第１回目（ｈ＝１）の分割投与後の再構成画像の生成の流れを示すフローチャートである。
【図６】第２回目（ｈ＝２）の分割投与後の再構成画像の生成の流れを示すフローチャートである。
【図７】円柱ファントムシミュレーション（総カウント：１．３６Ｍ ｃｏｕｎｔｓ／ｓｌｉｃｅ）による収束特性を示すグラフである。
【図８】投与量比が１：１の場合の円柱ファントムシミュレーション（２．７０Ｍｃｏｕｎｔｓ／ｓｌｉｃｅ）による再構成画像の比較を示す図である。
【図９】投与量比が１：１の場合の総カウント数に対する％ＣＯＶの変化を示すグラフである。
【図１０】横軸をスライス（ｓｌｉｃｅ）当たりの総カウント数、縦軸をＲＯＩ値の２回目／１回目の比に取った場合の総カウント数に対する変化を示すグラフである。
【図１１】横軸を投与量比（２回目／１回目）（実際には発生光子数比）、縦軸を％ＣＯＶとした場合のＳｕｂｔｒａｃｔｉｏｎ法とＡｄｄｉｔｉｏｎ法との比較を示す図である。
【図１２】最適投与量比におけるＳｕｂｔｒａｃｔｉｏｎ法、およびＡｄｄｉｔｉｏｎ法でのＨｏｆｆｍａｎファントムの再構成画像を示す図である。
【図１３】図１２に示した画像の１回目と２回目との各画素値の散布図である。
【図１４】Ａｄｄｉｔｉｏｎ法の２回目再構成において、１回目の脳内残留放射能のベースラインに、ノイズを含んだベースラインの場合とノイズフリーの理想的なベースラインを用いた場合の２種類の再構成画像を示す図である。
【図１５】図１４の両者の画素値を比較した散布図である。
【図１６】臨床データに適用した負荷（２回目）の再構成画像での比較を示す図である。
【図１７】３例の再構成画像にとったＲＯＩ（１５５４ＲＯＩｓ）における、２つの方法の比較を示す図である。
【図１８】本発明のコンピュータ・プログラムを実行する画像処理装置４４等のコンピュータの内部回路７０を示すブロック図である。
【図１９】２回分割投与法の概念を示す図である。
【符号の説明】
３，５ 画素、 １０ 断層画像、 ２０ 投影データ、 ３０，３１，３２，３３ 検出器、 ４０ 画素数記録部、 ４２ 検出確率記録部、 ５０−ｈ第ｈ回推定ＲＩ濃度初期値記録部、 ５２，６２ 散乱成分、 ５４，６４ プライマリ成分、 ６０−ｈ 第ｈ回測定データ記録部、 ７０ 内部回路、 ７２ ＣＰＵ、 ７４ ＲＯＭ、 ７６ ＲＡＭ、 ７８ ＶＲＡＭ、 ８０ 画像制御部、 ８２ バス、 ８４ａ ディスク、 ８４ｎ ＣＤ−ＲＯＭ、 ８６ コントローラ、 ８８ 入力操作部、 ９０ 入力制御部、 ９２ 外部Ｉ／Ｆ部。[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention provides a single photon emission computed tomography (SPECT) or a positron emission tomography (Position emission tomography: PET) based on each SPECT measurement data obtained by performing PET measurement data or each PET measurement data based on the measurement. The present invention also relates to an image reconstruction processing program that executes a process of obtaining a reconstructed image after each divided administration using a computer.
[0002]
[Prior art]
[Non-patent Document 1] Tsunehiko Nishimura, Revised version of the latest brain SPECT / PET clinical-brain function test method, page 12-page 17, page 65-page 66, Medical View Inc.
[0003]
2. Description of the Related Art In recent years, an image diagnostic apparatus for capturing an image of the inside of a body as an image without imposing a large burden on a patient or the like and performing an accurate diagnosis has become essential in a medical field. Such image diagnostic apparatuses include, for example, an X-ray tomography apparatus (Computer Tomography: CT), a magnetic resonance imaging apparatus (Magnetic Resonance Imaging: MRI), an ultrasonic diagnostic apparatus, and a radiation diagnostic apparatus. Image diagnosis using these image diagnostic machines is widely used as providing functional information such as early diagnosis of a patient's disease, selection of a treatment method, prediction and determination of a treatment effect, and the like.
[0004]
In a clinical field of nuclear medicine, for example, as described in Non-Patent Document 1, a radioisotope (Radio Isotope: RI) is introduced into a patient's body, and SPECT that uses γ-rays emitted therefrom is used. And PET are used by using respective devices. In the medical image processing apparatus as described above, various computer programs are prepared so that various image processing such as image reconstruction and image analysis can be performed from collected data.
[0005]
Tc (technesium) -99m ECD (ethyl cysteinate dimer) is frequently used as a short half-life tracer which is a radiopharmaceutical for cerebral blood flow SPECT. A feature of this tracer is that the stability after labeling, in which a part of the constituent elements of the compound is replaced with RI and marked, is good. Since the distribution in the brain is fixed within a few minutes after this tracer is taken into the brain tissue, it is a characteristic feature that the distribution is maintained for a long time. Taking advantage of this characteristic property, a split-dose method and a protocol for evaluating cerebral circulatory reserve due to acetazolamide (acetazolamide, ACZ) load have been clinically used. As described in Non-Patent Document 1, a twice-dose administration method is to administer a tracer of a dose usually used at one time in two divided doses, and to continuously administer two times under different conditions. This is a method of taking a brain image of the subject. The image of the local cerebral blood flow distribution under the second condition can be obtained by subtracting the first SPECT data (at rest) from the second SPECT data (at rest).
[0006]
FIG. 19 shows the concept of the two-dose administration method. FIG. 19A shows the first SPECT data, and FIG. 19B shows the second SPECT data. In FIG. 19A, reference numeral 54 denotes a primary component P of the first SPECT data._i ^[1st]And 52 is the scatter component S of the first SPECT data_i ^[1st]And y is the sum of the two_i ^[1st]Is the first SPECT data. In FIG. 19B, reference numeral 64 denotes a true primary component P of the second SPECT data._i ^{[2nd true]}And 62 is the true scattering component S of the second SPECT data_i ^{[2nd true]}The combination of the two is the second true SPECT data. As shown in FIG. 19B, the second SPECT data y_i ^[2nd]Shows the first SPECT data y, which is the residual radioactivity in the brain due to the first ECD administration._i ^[1st]Is added as a baseline. Therefore, the second SPECT data y_i ^[2nd]Is the first SPECT data y_i ^[1st]And the second true SPECT data. Therefore, in the conventional brain image processing in the twice-dose administration method, the second SPECT data y_i ^[2nd]First SPECT data y_i ^[1st]By subtracting, the reconstructed image of local cerebral blood flow distribution was obtained. Alternatively, subtraction is performed after generating a reconstructed image.
[0007]
[Problems to be solved by the invention]
However, the two-dose administration method has a problem that the signal-to-noise ratio (S / N ratio) of the load image is deteriorated by subtracting the first SPECT data from the second SPECT data. When the dose ratio is 1: 1, it is said that the variance of the second reconstructed image is reduced by a factor of 3 (that is, the standard deviation is reduced by a factor of 3) due to the subtraction. It is theoretically said that when the dose ratio is 1: 2, the first and second reconstructed images have the same S / N ratio. In order to avoid such image deterioration, it is necessary to increase the total ECD dose. However, an increase in the ECD dose results in an increase in exposure to the patient and an increase in examination cost, which is preferable. Not something.
[0008]
Therefore, an object of the present invention is to solve the above-described problem, and to perform brain image reconstruction without performing subtraction when using twice (or two or more) divided administration methods. It is another object of the present invention to provide an image reconstruction processing program which can perform the following.
[0009]
[Means for Solving the Problems]
An image reconstruction processing program according to the present invention provides a radioisotope (Radio Isotope: RI) -labeled tracer L (L is 2 or more) divided doses to a predetermined tissue, and a single photon emission computed tomography after each divided dose. Based on each SPECT measurement data or each PET measurement data based on each SPECT measurement data or each PET measurement data obtained by performing single photon emission computed tomography (SPECT) or positron emission tomography (PET) measurement. An image reconstruction processing program that executes the process of obtaining by using a computer,
Let j (= 1 to m) be the pixel on the reconstructed image, where m is the total number of pixels on the reconstructed image and release i (= 1 to n) from the RI administered to the given tissue The photon is a pixel on the projection data when projected to a detector that detects the photon, where m and n are recorded in the pixel number recording unit in advance,
c_ijIs the probability that a photon emitted from pixel j on the reconstructed image is detected at one pixel i on the projection data, where c_ijAre previously recorded in the detection probability recording unit from j = 1 to m recorded in the pixel number recording unit and from i = 1 to n recorded in the pixel number recording unit,
For h = 1 to L, λ_j ^{0 [H]}Is the estimated RI density of pixel j on the reconstructed image in the 0th successive approximation used for the h-th divided administration, where λ_j ^{0 [H]}Are previously recorded in the h-th estimated RI density initial value recording unit from j = 1 to the number m of the pixels j recorded in the pixel number recording unit,
For h = 1 to L, y_i ^[H]Is the count number of photons detected at one pixel i on the projection data after the h-th divided administration, where y_i ^[H]Are recorded in the h-th measurement data recording unit from i = 1 to the number n of the pixels i recorded in the pixel number recording unit after the measurement.
The reconstructed image after the first divided administration is converted to the estimated RI density λ of pixel j on the reconstructed image in the k-th successive approximation._j ^{k [1]}A first reconstructed image generation step of repeatedly calculating from equation (1) from k = 0 until a predetermined first convergence condition is satisfied;
(Equation 8)

Where λ_j ^{0 [1]}Is set in the first estimated RI concentration initial value recording section, and c_ijIs recorded in the detection probability recording unit, and y_i ^[1]Is recorded in the first measurement data recording unit,
From h = 2 to L, the reconstructed image after the h-th divided administration is converted to the estimated RI density λ of the pixel j on the reconstructed image in the k-th successive approximation._j ^{k [H]}H-th reconstructed image generation step of repeatedly calculating from k = 0 until a predetermined second convergence condition is satisfied using Expression (2).
(Equation 9)

Where λ_j ^{0 [H]}Is set in the h-th estimated RI concentration initial value recording unit, and y_i ^{[ h]}Is an image reconstruction processing program for executing the step of being recorded in the h-th measurement data recording section.
[0010]
Here, in the image reconstruction processing program of the present invention, for h = 1 to L, S_i ^[H]Is the scattering correction term at one pixel i on the projection data after the h-th divided administration, where S_i ^[H]Are recorded in the h-th scattering correction term recording unit from i = 1 to the number n of the pixels i recorded in the pixel number recording unit;
(Equation 10)

Equation (3) can be used in place of Equation (1) in the first reconstructed image generation step, and Equation (4) can be used in place of Equation (2) in the h-th reconstructed image generation step. .
[0011]
Here, in the image reconstruction processing program according to the present invention, the predetermined first convergence condition is such that a value obtained by Expression (5) is:
(Equation 11)

The condition may be equal to or less than a predetermined value.
[0012]
Here, in the image reconstruction processing program according to the present invention, the predetermined second convergence condition is such that a value obtained by Expression (6) is:
(Equation 12)

The condition may be equal to or less than a predetermined value.
[0013]
Here, in the image reconstruction processing program according to the present invention, the predetermined first convergence condition is such that a value obtained by Expression (7) is:
(Equation 13)

The condition may be equal to or less than a predetermined value.
[0014]
Here, in the image reconstruction processing program according to the present invention, the predetermined second convergence condition is such that a value obtained by Expression (8) is:
[Equation 14]

The condition may be equal to or less than a predetermined value.
[0015]
Here, in any one of the image reconstruction processing programs of the present invention, c recorded in the detection probability recording unit_ijCan be set to 1 when the pixel i on the projection data exists in the projection direction of the pixel j on the reconstructed image, and can be set to 0 otherwise.
[0016]
Here, in any one of the image reconstruction processing programs according to the present invention, the λ recorded in the first estimated RI density initial value recording unit is stored._j ^{0 [H]}Can be previously recorded as 1 from j = 1 to the number m of pixels j recorded in the pixel number recording unit.
[0017]
A recording medium according to the present invention is a computer-readable recording medium recording any one of the image reconstruction processing programs according to the present invention.
[0018]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.
[0019]
Embodiment 1 FIG.
The image reconstruction processing program according to the present invention is characterized by using an algorithm in which the influence of the baseline by the previous ECD administration up to the previous time (for example, once) is incorporated in the successive approximation type image reconstruction method in advance. Therefore, the successive approximation type image reconstruction method will be described first.
[0020]
As a method of estimating the unknown RI distribution in the body as the most likely case from the collected projection data by a statistical method, a maximum likelihood estimation-expected value maximization (Maximum Likelihood-Expectation Maximization: ML-EM) is used. There is a reconstruction method. This ML-EM method is a method combining the maximum likelihood estimation method and the expectation value maximization method. The ML-EM method is a method of repeatedly correcting the RI distribution in the body by using the successive approximation method so that the projection data calculated from the estimated RI distribution in the body is made as close as possible to the actually collected projection data.
[0021]
FIG. 1 illustrates the application of the maximum likelihood estimation method to image reconstruction. In FIG. 1, reference numeral 10 denotes a tomographic image of a brain image, 3 denotes one pixel having a RI density of λ (MBq) in the matrix of the tomographic image 10, and 30 denotes a detector, which is one pixel on the projection data 20. c indicates a detection probability that a photon is detected at the pixel 30, and x indicates a count of photons from the pixel 3 in one pixel 30 (or one pixel 30 assuming that the pixel is a detector) on the projection data 20. . When the photon emitted from the RI (pixel 3) having the RI concentration of λ (MBq) is measured n times by the detector 30 (pixel 30), the maximum likelihood estimation value of the count x is defined as xi where xi is a count for each measurement. Is given by equation 9.
[0022]
[Equation 15]

[0023]
The maximum likelihood estimation value of the RI concentration λ is given by Expression 10.
[0024]
(Equation 16)

[0025]
Next, a case where RI is measured by an n-direction detector will be described. FIG. 2 illustrates application of the maximum likelihood estimation method to image reconstruction when RI is measured by an n-direction detector. In FIG. 2, the portions denoted by the same reference numerals as those in FIG. In FIG. 2, reference numeral 5 indicates that the RI density in the matrix of the tomographic image 10 is λ._j(MBq) indicates one pixel j (j = 1 to the total number m of pixels in the matrix of the tomographic image 10). Reference numerals 31 to 33 denote detectors, each of which indicates one pixel on the projection data 20 or the like. Although there are n detectors, only three are shown for the sake of illustration. c_ij(I = 1 to the number n of the detectors or the number n of each pixel on the projection data, j = 1 to m), the photon emitted from the pixel j (5) is detected by the pixel i (or the detector i). X indicates the detection probability_ijDenotes the number of photons emitted from pixel j (5) and detected at pixel i (or detector i). Under the above conditions, the density λ of the pixel j (5)_jIs given by equation 11.
[0026]
[Equation 17]

[0027]
Actually, since RI exists also around the pixel j (5), the number of photons x_ijCannot be measured directly. FIG. 3 illustrates the application of the maximum likelihood estimation method to image reconstruction when considering the RI around the pixel j (5). In FIG. 3, the portions denoted by the same reference numerals as those in FIG. In FIG. 3, y_iDenotes the total sum (projection count) of the measurement values from each pixel j (5) in the detector i. Sum y_iFrom the number of photons x_ijNeeds to be estimated. Here, when the maximum posterior probability estimation (MAP estimation) is used, the number of photons x_ijIs given by Equation 12.
[0028]
(Equation 18)

[0029]
That is, the RI density λ of the pixel j (5), which is an unknown to be obtained by the maximum likelihood estimation method,_jIs the detection probability c_ijAnd the number of photons x_ijAnd the number of photons x_ijCannot be measured directly, so the number of photons x_ijFrom the MAP estimated value. Therefore, the number of photons x on the right side of Equation 11_ijThe number of photons x in Equation 12 is_ijBy substituting the MAP estimated value of the above, 2 of Expression 11 can be obtained.
[0030]
[Equation 19]

[0031]
Equation 2 in Equation 11 indicates that the RI concentration λ_jBecause λ_jI can not solve Therefore, a sequential equation is created by the EM method. λ_j ^kIs the estimated value in the k-th successive approximation, Expression 13 can be obtained from Expression 11-2.
[0032]
(Equation 20)

[0033]
FIG. 4 shows the RI concentration λ_jA calculation example for estimating is shown. As shown in FIG. 4, first, the number of projection directions is set to n = 2, the matrix size of the tomographic image is set to 3 × 3 (m = 3 × 3 = 9), and each detection probability c_ijIs assumed to be 1 only in the projection direction and 0 in other directions (step S10). That is, it is assumed that measurement is not performed from the pixel j (5) other than the projection direction. As shown in step S12 in FIG. 4, if the true value is from 1 to 9 from m = 1 to 9, the projection count number y in each projection direction_iIs the sum of the values of pixel j in the projection direction. For example, in the direction of m = 1, 4, and 7, the values of the pixel j are 1, 4, and 7, respectively, so the projection count number y_i= 1 + 4 + 7 = 12. Other projection count number y_iIs also the same, and the description is omitted. As shown in step S12, the estimated value is λ from m = 1 to 9.₁Or λ₉Becomes Next, for all k = 0, λ_j ⁰Is assumed to be 1 and the estimation calculation is performed by Expression 13 (Step S14). Details are as shown in step S16. For example, λ₁ ¹For λ₁ ⁰= 1.0 and c₁₁λ₁ ⁰From c₁₉λ₉ ⁰For the sum up to₁₁, C₁₄, C₁₇Is only 1 and the others are 0, y₁c₁₁= 12 × 1, y₂c₂₁= 6 × 1 and c₁₁+ C₂₁= 2, λ₁ ¹= 3.0.
[0034]
As shown in step S18 as the first estimated value, from m = 1 to 9, λ₁ ¹= 3.0, λ₂ ¹= 3.5, λ₃ ¹= 4.0, λ₄ ¹= 4.5, λ₅ ¹= 5.0, λ₆ ¹= 5.5, λ₇ ¹= 6.0, λ₈ ¹= 6.5, λ₉ ¹= 7.0.
Next, λ_j ¹Similarly, λ_j ²Is obtained (step S20). As shown in step S22 as the second estimation value, from m = 1 to 9, λ₁ ²= 2.2, λ₂ ²= 2.8, λ₃ ²= 3.3, λ₄ ²= 4.3, λ₅ ²= 5.0, λ₆ ²= 5.8, λ₇ ²= 6.4, λ₈ ²= 7.3, λ₉ ²= 8.1. The above calculation is repeated a specified number of times (step S24). More specifically, it is determined whether or not the repetition is completed according to a desired convergence condition. As a result, a true value (or an approximate value) can be obtained (step S26).
[0035]
An image reconstruction processing program according to the present invention uses the above-described successive approximation type image reconstruction method. 5 and 6 show a series of flowcharts of the image reconstruction processing program of the present invention. In this flowchart, the number h of divided doses is set to two for convenience of explanation. FIG. 5 is a flowchart (first reconstructed image generation step) of generating a reconstructed image after the first (h = 1) divided administration, and FIG. 6 is a flowchart of the processing after step S40a in FIG. That is, this is a flowchart (second reconstructed image generation step) of generating a reconstructed image after the second (h = 2) divided administration. As shown in FIG. 5, first, the total number m of pixels on the reconstructed image and the total number n of pixels on the projection data are read from the pixel number recording unit 40 (described later; see FIG. 18) (step S30a). The pixel j on the reconstructed image is set to 1 (step S32a), and the number of repetitions k is set to k = 0 (step S34a). Estimated RI density λ of pixel j on reconstructed image_j ^{k [ 1 ]}Is read from the first (h = 1) estimated RI concentration initial value recording unit 50-1 (described later; see FIG. 18) (step S36a). For pixels i = 1 to n on the projection data, detection probability c_ijIs read out from the detection probability recording unit 42 (described later; see FIG. 18), and the count y of photons detected at the pixel i on the projection data_i ^[1]Is read out from the first (h = 1) -times measurement data recording unit 60-1 (described later; see FIG. 18), and Expression 13 is calculated as in the above-described example (step S38a). However, in consideration of the number of divided doses h = 1, Expression 13 is expressed as Expression 13-2. In [2] of Expression 13, [1] means h = 1.
[0036]
(Equation 21)

[0037]
It is determined whether or not a desired convergence condition (first convergence condition) is satisfied (step S40a). If not, the number of repetitions k is increased by one, and the λ obtained by Expression 13 is determined._j ^{1 [ 1 ]}Is used to calculate Equation 13 for a new k in the same manner as in Step S38a. When the convergence condition is satisfied in Step S40a, the generation of the first (h = 1) reconstructed image after the divided administration is completed, and the process proceeds to A.
[0038]
As shown in FIG. 6, first, the total number m of pixels on the reconstructed image and the total number n of pixels on the projection data are read from the pixel number recording unit 40 (described later; see FIG. 18) (step S30b). Step S30b is the same as step S30a and may be omitted. The pixel j on the reconstructed image is set to 1 (step S32b), and the number of repetitions k is set to k = 0 (step S34b). Estimated RI density λ of pixel j on reconstructed image_j ^{k [ 2 ]}Is read from the second (h = 2) estimated RI concentration initial value recording unit 50-2 (described later; see FIG. 18) (step S36b). For pixels i = 1 to n on the projection data, detection probability c_ijIs read out from the detection probability recording unit 42 (described later; see FIG. 18), and the count y of photons detected at the pixel i on the projection data_i ^{[ 2 ]}Is read from the second (h = 2) times measurement data recording unit 60-2 (described later; see FIG. 18), and the following Expression 14 is calculated (step S38b).
[0039]
(Equation 22)

[0040]
The feature of the image reconstruction processing program of the present invention is that the total count y_i ^[H-1]The point is to add. Therefore, in the case of two divided doses, h = 1 to 2, so that the first count number y in the denominator of Expression 14 is obtained._i ^[1]Is added. That is, when reconstructing the second SPECT data, the residual radioactivity y in the brain due to the first administration is added to the denominator of Equation 14_i ^[1]Is considered to exist as a baseline. In the case where three divided doses are performed, h = 1 to 3; therefore, the second count number y in the denominator of Expression 14 is obtained._i ^[2]Is added. If h divided doses are given, the total count y of the (h-1) th time in the denominator_i ^[H-1]Is added.
[0041]
Next, it is determined whether or not a desired convergence condition (second convergence condition) is satisfied (step S40b). If not, the number of repetitions k is increased by one, and the λ obtained by Expression 14 is determined._j ^{k [h]}Is used to calculate Equation 14 for a new k in the same manner as in step S38b. If the convergence condition is satisfied in step S40b, the generation of the reconstructed image after the second (h = 2) divided administration is completed.
[0042]
As the above-mentioned first convergence condition, the sum of squares of the difference obtained by subtracting the denominator from the numerator of Equation 13 (square error: χ²) Can be used. For example, Equation 15 is calculated as the first convergence condition, and²When the value of is less than or equal to a predetermined value (ideally, when it becomes 0), iterative calculation can be terminated.
[0043]
(Equation 23)

[0044]
As the above-described second convergence condition, the sum of squares of the difference obtained by subtracting the denominator from the numerator of Expression 14 (square error: χ²) Can be used. For example, Equation 16 is calculated as the second convergence condition, and χ²When the value becomes equal to or less than a predetermined value (ideally, when it becomes 0), the repetitive calculation can be terminated.
[0045]
[Equation 24]

[0046]
Verification.
In order to verify the usefulness of the image reconstruction processing program of the present invention, Monte Carlo simulation was used. In the Monte Carlo simulation, (1) As a “definition of a simulation system”, a section of a human body system is defined by combining geometric shapes such as an elliptical sphere, a cylinder, and an elliptic cylinder. The human body shape is often converted into a numerical phantom. (2) A predetermined element or a molecule such as water is used as the “definition of substance data”, and a linear attenuation coefficient μ specific to the substance is assumed. (3) Photons, electrons, and the like are defined as “definition of the radiation source”. (4) A desired measurement system is defined as “definition of measurement system”. In SPECT, a collimator is mounted on the front of the detector in order to limit the direction of incidence of gamma rays. These devices are also incorporated into the Monte Carlo simulation.
[0047]
In the Monte Carlo simulation, absorption / scattering, degradation of resolution by a collimator, and statistical noise were simulated for 140 KeV photon generation. The collimator was assumed to be a low energy collimator (for example, a LEHR collimator). The matrix size was 128 × 128 (m = 16384), and the number of projection directions was n = 120. The pixel size is 3.125 mm. The main window for collecting photoelectric peaks was 20% centered on 140 KeV, and the sub-window set below the main window was 5% centered on 124 KeV. The numerical phantom uses a simple cylindrical phantom (18 cmφ) and has a uniform absorber (μ = 0.15 cm).^-1).
[0048]
Eight types (1.36, 2.70, 5.40, 10.8, 21.6, 43.3, 86.5 Mcounts / slice) of projection data are created so that the total count number increases every two times. However, there was no decay due to half-life. In this simulation, using the seven sets of projection data that can be paired with the data increased every two times so that the doses of the first and second doses are 1: 1, the above-mentioned Expression 13-2 and Expression 2 are used. 14, image reconstruction was performed. The number of repetitions k was set to 3 and subset = 20. As the pre-processing filter, a Butterworth filter was used. The cutoff frequency fc was 0.36 [cycles / pixel], and the order was 8. The absorption correction was incorporated into the detection probability using the assumed μ map (two-dimensional image of the linear attenuation coefficient μ). The scattering component is Triple
Estimation was performed using the Energy Window (TEW) method.
[0049]
The comparison between the conventional method (referred to as Subtraction method) and the method using the image reconstruction processing program of the present invention (referred to as Addition method) is based on the evaluation by the visual inspection, and the center is set in the cylindrical phantom. The calculation was performed by calculating% COV (= standard deviation / ROI value × 100%) within the ROI (5 cmφ).
[0050]
Evaluation of optimal dose ratio.
By the above-described cylindrical phantom simulation, the optimal dose ratio at which the% COV becomes the same between the first time and the second time was determined by both the subtraction method and the addition method. The dose ratio between the second and first doses was set to eight types of 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0, and the The total amount was always constant. In this case, the second projection data was 5.40 Mcounts / slice in any case. The reconstruction conditions and data processing are exactly the same as above.
[0051]
Under the conditions of the determined optimal dose ratio, Monte Carlo simulation was performed using one slice of the digitized Hoffman phantom, and reconstructed by the Subtraction method and the Addition method. Assuming that the radioactivity concentration ratio between gray matter and white matter is 4: 1, a uniform absorber (μ = 0.15 cm^-1) Was assumed.
[0052]
Evaluate the effects of a noisy baseline.
As described above, in the denominator of the Addition method (2 or 14 of Expression 13), if there is a baseline including noise, convergence by the ML method is not guaranteed. Therefore, in order to examine the degree of error, a Monte Carlo simulation was performed using a digitized Hoffman phantom. The total number of photons for the first (baseline) and the second time was a dose ratio of 1: 1 and the total counts were 6.84 Mcounts / slice and 13.7 Mcounts / slice, respectively. In order to compare with the case where an ideal noise-free baseline is used as the baseline when performing the reconstruction according to 2 or 14 in Equation 13, projection data having 32 times the number of photons is created, standardized, and used. Was.
[0053]
Clinical data.
Evaluation was performed on clinical data of three cases. SPECT measurement was performed using Toshiba GCA-7200, and a LEGP collimator was used as a collimator. The matrix size was 64 × 64 (m = 4096), and the number of projection directions n = 60. The pixel size is 4.3 mm. The main window was set to 20% centering on 140 KeV, and the subwindow set below the main window was set to 8% centering on 122 KeV. After 5 minutes from the first ECD administration (333 MBq), a 7.5 minute SPECT measurement was performed. A second ECD administration (222 MBq) was performed immediately after the first SPECT measurement, and a second SPECT measurement was performed under the same conditions 12.5 minutes after the first administration. ACZ
A load of (1000 mg) was performed immediately before the first SPECT. Detailed information on collection conditions and protocols is omitted.
[0054]
In the data processing, the image was reconstructed by the subtraction method and the addition method. The number of repetitions k was set to 3 and subset = 15. As a pre-processing filter, a Butterworth filter (cutoff frequency fc = 0.24 cycles / pixel, next = 8) was used. The absorption correction is performed by performing contour extraction based on a threshold to obtain a uniform μ value (0.14 cm^-1) Was assumed and incorporated into the detection probability. The scattering component was estimated by the TEW method as in the simulation. The ROI analysis was used to compare the subtraction method and the addition method.
[0055]
result.
FIGS. 7A and 7B are graphs showing convergence characteristics obtained by a cylindrical phantom simulation (total count: 1.36 M counts / slice). 7A and 7B, the horizontal axis represents the number of repetitions k, and the vertical axis represents the square error (χ²The first time is given by Equation 15, and the second time is given by Equation 16.) Although the absolute value of the final error was substantially the same in both the subtraction method and the addition method, a clear difference was observed in the convergence characteristics. That is, it has been clarified that the convergence speed of the second SPECT data is slightly reduced in the Addition method. However, when the product of the subset and the number of repetitions k is 20 or more, the two are almost the same, and it is considered that there is no difference in the characteristics under the reconstruction conditions used this time.
[0056]
FIG. 8 is a comparison of reconstructed images by a cylindrical phantom simulation (2.70 Mcounts / slice) when the dose ratio is 1: 1. In the Subtraction method (left side in FIG. 8), the S / N ratio of the second (2nd) reconstructed image is poor. On the other hand, in the Addition method (right side in FIG. 8), the S / N ratio is significantly improved, and as a result, the first time (1st) The image quality is almost the same as.
[0057]
FIGS. 9A and 9B are graphs showing changes in% COV with respect to the total count number when the dose ratio is 1: 1. 9A and 9B, the horizontal axis represents the total number of counts per slice, and the vertical axis represents% COV. The% COV of the subtraction method shown in FIG. 9A is clearly worse in the second time than in the first time (50.4% at the first time and 80.5 at the second time when the total count is 2.70 Mcounts / slice). %), Which is clearly affected by subtraction. On the other hand, in the Addition method shown in FIG. 9 (B), the &% COV of the first time and the second time are almost the same (when the total count is 2.70 Mcounts / slice, the first time is 43.1%, the second time). 41.4%), confirming the usefulness of the Addition method of the present invention. The first time of the subtraction method and the first time of the addition method completely match when the scattering correction is not performed, but in the present simulation, the addition method is slightly better because the scattering correction is performed.
[0058]
On the other hand, FIG. 10 is a graph showing the change with respect to the total count when the horizontal axis represents the total count per slice, and the vertical axis represents the ratio of the second / first ROI value. In this simulation, this ratio should be 1.0, but the value was almost assumed in both the subtraction method and the addition method.
[0059]
FIGS. 11A and 11B show a comparison between the Subtraction method and the Addition method when the abscissa is the dose ratio (second time / first time) (actually, the ratio of the number of generated photons) and the ordinate is% COV. Is shown. As shown in FIG. 11 (A), in the Subtraction method, the first and second% COV became the same when the dose ratio was about 2.0. On the other hand, as shown in FIG. 11 (B), in the Addition method, the first and second% COV became the same when the dose ratio was about 1.0. The absolute value of% COV at the optimal dose was about 54% by the Subtraction method and about 43% by the Addition method. Also at the optimal dose ratio, the result that the Addition method had a better S / N ratio was obtained. This is because the total dose (the number of generated photons in the simulation) is kept the same, so that the first dose of the Subtraction method is 2/3 of the dose of the Addition method.
[0060]
FIG. 12 shows reconstructed images of the Hoffman phantom by the Subtraction method and the Addition method at the optimal dose ratio. As shown in FIG. 12, it is understood that the Addition method (right side in FIG. 12) has a better S / N ratio, similarly to the result of the cylindrical phantom simulation.
[0061]
FIG. 13 shows a scatter diagram of each pixel value of the first and second times of the image shown in FIG. 13A and 13B, the horizontal axis represents the first SPECT count, and the vertical axis represents the second SPECT count. As shown in FIG. 13A, the correlation coefficient becomes r = 0.762 by the subtraction method, and as shown in FIG. 13B, r becomes 0.813 by the addition method, and the variation of the subtraction method is large. Was.
[0062]
FIG. 14 shows that in the second reconstruction of the Addition method, the first baseline of residual radioactivity in the brain is 1) a baseline including noise (left side of FIG. 14), and 2) an ideal noise-free baseline. 14 shows two types of reconstructed images when a baseline is used (right side in FIG. 14). In visual evaluation, the two are almost the same.
[0063]
FIG. 15 is a scatter diagram comparing the pixel values of both pixels in FIG. In FIG. 15, the horizontal axis is the ideal noise-free baseline in FIG. 14, and the vertical axis is the noise-containing baseline in FIG. No systematic error was observed in either case.
[0064]
FIG. 16 shows a comparison of the reconstructed image of the load (second time) applied to the clinical data. In FIG. 16, the left side shows the case of the Subtraction method, and the right side shows the case of the Addition method.
[0065]
FIG. 17 shows a comparison of the two methods on ROIs (1554 ROIs) for three reconstructed images. FIG. 17A shows the first time (at rest), and FIG. 17B shows the second time (at the time of ACZ load). 17A and 17B, the horizontal axis represents the ROI count (Addition method), and the vertical axis represents the ROI count (Subtraction method). When visually evaluated, the Subtraction method was a noisy image, but in the comparison of ROI values, although there was a slight variation for the second time, no significant difference was found.
[0066]
As described above, according to the first embodiment of the present invention, by incorporating the effect of residual radioactivity in the brain into the calculation formula of the sequential image reconstruction method, the load in the Tc-99mECD two-dose administration method is reduced. The S / N ratio of an image can be improved. Specifically, in the denominator of Expression 14, the total count y_i ^[H-1]To incorporate the effects of residual radioactivity in the brain. For example, in the case of two divided doses, the first count number y in the denominator of Expression 14_i ^[1]Is added. That is, when reconstructing the second SPECT data, the residual radioactivity y in the brain due to the first administration is added to the denominator of Equation 14_i ^[1]Is considered to exist as a baseline. The method according to the image reconstruction processing program of the present invention was evaluated by applying it to simulation data and clinical data. As a result, even when the dose ratio is 1: 1, the S / N ratios of the first and second reconstruction data are almost the same, and it is confirmed that the S / N ratio is improved as compared with the conventional Subtraction method. Was completed. According to the above simulation results, there was no quantitative difference between the case where ideal projection data was used as the baseline and the case where projection data containing noise was used as the baseline. As described above, it is possible to provide an image reconstruction processing program capable of performing brain image reconstruction without performing subtraction when using two (or two or more) divided administration methods.
[0067]
Embodiment 2 FIG.
In Equations 13-2 and 14 in Embodiment 1 described above, scattering correction can also be performed. For the number of divided doses h = 1 to L, S_i ^[H]Is the scattering correction term for one pixel i on the projection data after the h-th divided administration. Where S_i ^[H]Is stored in the h-th scattering correction term recording unit from i = 1 to the number n of pixels i recorded in the pixel number recording unit 10. In this case, Equation 17 may be used instead of Equation 2 and Equation 18 may be used instead of Equation 14.
[0068]
(Equation 25)

[0069]
When Expressions 17 and 18 are used, the first convergence condition in consideration of the scattering correction is expressed by Expression 19, and the second convergence condition is expressed by Expression 20.
[0070]
(Equation 26)

[0071]
Embodiment 3 FIG.
FIG. 18 is a block diagram showing an internal circuit 70 of a computer that executes a computer program of the present invention for implementing the above-described embodiment. As shown in FIG. 18, the CPU 72, the ROM 74, the RAM 76, the image control unit 78, the controller 86, the input control unit 90, and the external interface (I / F) unit 92 are connected to the bus 82. In FIG. 18, the above-described computer program of the present invention is recorded on a recording medium (including a removable recording medium) such as the ROM 74, the disk 84a, or the CD-ROM 84n. On the disk 84a and the like, a pixel number recording unit 40, a detection probability recording unit 42, an h-th estimated RI density recording unit 50-h, an h-th measurement data recording unit 60-h, and the like are recorded. The computer program of the present invention is loaded into the RAM 76 from the ROM 74 via the bus 82 or from a recording medium such as the disk 84a or the CD-ROM 84n via the controller 86 via the bus 82. The image control unit 78 sends the brain image data to the VRAM 80, and the display unit 45, which is a display device such as a display, displays the brain image based on the brain image data sent from the VRAM 80. The VRAM 80 is an image memory having a capacity corresponding to the data capacity of one screen of the display unit 45. The input operation unit 88 is an input device such as a mouse and a numeric keypad for inputting to the image processing device 44. The input control unit 90 is connected to the input operation unit 88 and performs input control and the like. The external I / F unit 92 has an interface function for connecting to an external communication network (not shown) such as the Internet, for example.
[0072]
The object of the present invention can be achieved when the CPU 72 executes the computer program of the present invention as described above. The computer program can be supplied to the computer CPU 72 in the form of a recording medium such as a CD-ROM 84n as described above, and a recording medium such as a CD-ROM 84n that stores the computer program also constitutes the present invention. Will be. As the recording medium on which the computer program is recorded, for example, a memory card, a memory stick, a DVD, an optical disk, an FD, or the like can be used in addition to the recording medium described above.
[0073]
【The invention's effect】
As described above, according to the image reconstruction processing program of the present invention, by incorporating the effect of residual radioactivity in the brain into the calculation formula of the sequential image reconstruction method, the Tc-99mECD in the two-dose administration method can be used. The S / N ratio of the load image can be improved. Specifically, the total count number y up to the previous time in the denominator of Expression 14 is_i ^[H-1]To incorporate the effects of residual radioactivity in the brain. For example, in the case of two divided doses, the first count number y in the denominator of Expression 14_i ^[1]Is added. That is, when reconstructing the second SPECT data, the residual radioactivity y in the brain due to the first administration is added to the denominator of Equation 14_i ^[1]Is considered to exist as a baseline. As a result of applying the method according to the image reconstruction processing program of the present invention to simulation data and clinical data and evaluating the results, the S / N ratio of the first and second reconstruction data was obtained even when the dose ratio was 1: 1. The ratios were almost the same, and it was confirmed that the S / N ratio was improved compared to the conventional Subtraction method. According to the above simulation results, there was no quantitative difference between the case where ideal projection data was used as the baseline and the case where projection data containing noise was used as the baseline. As described above, it is possible to provide an image reconstruction processing program capable of performing brain image reconstruction without performing subtraction when using two (or two or more) divided administration methods.
[Brief description of the drawings]
FIG. 1 is a diagram for explaining application of a maximum likelihood estimation method to image reconstruction.
FIG. 2 is a diagram illustrating the application of maximum likelihood estimation to image reconstruction when RI is measured by an n-direction detector.
FIG. 3 is a diagram illustrating application of the maximum likelihood estimation method to image reconstruction when considering RI around a pixel j (5).
FIG. 4 is a graph showing RI concentration λ using equation 13._jFIG. 9 is a diagram illustrating a calculation example for estimating.
FIG. 5 is a flowchart showing a flow of generating a reconstructed image after the first (h = 1) divided administration.
FIG. 6 is a flowchart showing a flow of generation of a reconstructed image after the second (h = 2) divided administration.
FIG. 7 is a graph showing convergence characteristics obtained by a cylindrical phantom simulation (total count: 1.36 M counts / slice).
FIG. 8 is a diagram showing a comparison of reconstructed images by a cylindrical phantom simulation (2.70 Mcounts / slice) when the dose ratio is 1: 1.
FIG. 9 is a graph showing the change in% COV with respect to the total count number when the dose ratio is 1: 1.
FIG. 10 is a graph showing a change with respect to the total count when the horizontal axis is the total count per slice and the vertical axis is the ratio of the ROI value to the second / first time.
FIG. 11 is a graph showing a comparison between the subtraction method and the addition method when the dose ratio (second time / first time) (actually, the ratio of the number of generated photons) is plotted on the horizontal axis and% COV is plotted on the vertical axis.
FIG. 12 is a diagram showing reconstructed images of a Hoffman phantom by the Subtraction method and the Addition method at the optimal dose ratio.
13 is a scatter diagram of each pixel value of the first and second times of the image shown in FIG. 12;
FIG. 14 shows two types of cases in which a noise-containing baseline and a noise-free ideal baseline are used as the first baseline of residual radioactivity in the brain in the second reconstruction of the Addition method. FIG. 3 is a diagram showing a reconstructed image of FIG.
FIG. 15 is a scatter diagram comparing the two pixel values of FIG. 14;
FIG. 16 is a diagram showing a comparison of a load (second time) applied to clinical data in a reconstructed image.
FIG. 17 is a diagram showing a comparison between two methods in ROIs (1554 ROIs) obtained for three reconstructed images.
FIG. 18 is a block diagram showing an internal circuit 70 of a computer such as an image processing device 44 that executes a computer program of the present invention.
FIG. 19 is a diagram showing the concept of a two-dose administration method.
[Explanation of symbols]
3, 5 pixels, 10 tomographic images, 20 projection data, 30, 31, 32, 33 detector, 40 pixel number recording unit, 42 detection probability recording unit, 50-h h-th estimated RI density initial value recording unit, 52 , 62 scattering components, 54, 64 primary components, 60-h hth measurement data recording section, 70 internal circuit, 72 CPU, 74 ROM, 76 RAM, 78 VRAM, 80 image control section, 82 bus, 84a disk, 84n CD-ROM, 86 controller, 88 input operation unit, 90 input control unit, 92 external I / F unit.

Claims (9)

A tracer labeled with a radioisotope (Radio Isotope: RI) is divided into L (L is 2 or more) divided doses to a predetermined tissue, and after each divided dose, single photon emission computed tomography (SPECT) is performed. Alternatively, a process of obtaining a reconstructed image after each divided administration is performed using a computer based on each SPECT measurement data or each PET measurement data obtained by performing Positron Emission Tomography (PET) measurement. An image reconstruction processing program,
Let j (= 1 to m) be the pixel on the reconstructed image, where m is the total number of pixels on the reconstructed image and release i (= 1 to n) from the RI administered to the given tissue The photon is a pixel on the projection data when projected to a detector that detects the photon, where m and n are recorded in the pixel number recording unit in advance,
and the probability that the photons emitted from the pixel j in the reconstructed image c _ij are detected on a pixel i on the projection data, wherein c _ij is recorded from previously j = 1 to the number of pixels recorded portion m And from i = 1 to n recorded in the pixel number recording section are recorded in the detection probability recording section,
For h = 1 to L, let λ _j ^{0 [h]} be the estimated RI density of pixel j on the reconstructed image in the 0 th successive approximation used for the h th divided dose, where λ _j ^{0 [ h]} are previously recorded in the h-th estimated RI density initial value recording unit from j = 1 to the number m of the pixels j recorded in the pixel number recording unit in advance,
For h = 1 to L, y _i ^[h] is the count of photons detected at one pixel i on the projection data after the h-th divided administration, where y _i ^[h] is i = 1 to the number n of the pixels i recorded in the pixel number recording unit are recorded in the h-th measurement data recording unit.
The reconstructed image after the divided doses of the first round, the estimated RI concentration lambda _{j k} pixel j in the reconstructed image in k-th iterative ^[1] Equation (1) from k = 0 predetermined by using A first reconstructed image generation step obtained by repeatedly calculating until a first convergence condition is satisfied;

Here, λ _j ^{0 [1]} is set in the first estimated RI concentration initial value recording unit, c _ij is recorded in the detection probability recording unit, and y _i ^[1] is the first measurement data recording. Is recorded in the department,
from h = 2 to L, and the h-th reconstructed image after the divided doses of the estimated RI concentration lambda _{j k} pixel j in the reconstructed image ^[h] in the successive approximation of the k-th using the equation (2) An h-th reconstructed image generation step of repeatedly calculating from k = 0 until a predetermined second convergence condition is satisfied;

Here, lambda _{j 0} ^[h] is set to the first h times estimated RI concentration initial value recording unit, executes the steps of y i ^_[h] are those recorded in the first h measurements data recording unit Image reconstruction processing program for causing In claim 1, wherein the image reconstruction processing program, the L from h = 1, the scatter correction term in 1 pixel i on the projection data S i ^_[h] after divided doses of the h-th, where S _i ^[H] is recorded in the h-th scattering correction term recording unit from i = 1 to the number n of the pixels i recorded in the pixel number recording unit;

Equation (3) is used in place of Equation (1) in the first reconstructed image generation step, and Equation (4) is used in place of Equation (2) in the h-th reconstructed image generation step. Image reconstruction processing program. 2. The image reconstruction processing program according to claim 1, wherein the predetermined first convergence condition is a value obtained by Expression (5):

An image reconstruction processing program, characterized in that the condition is not more than a predetermined value. 2. The image reconstruction processing program according to claim 1, wherein the predetermined second convergence condition is a value obtained by Expression (6):

An image reconstruction processing program, characterized in that the condition is not more than a predetermined value. 3. The image reconstruction processing program according to claim 2, wherein the predetermined first convergence condition is:

An image reconstruction processing program, characterized in that the condition is not more than a predetermined value. 3. The image reconstruction processing program according to claim 2, wherein the predetermined second convergence condition is a value obtained by Expression (8):

An image reconstruction processing method, characterized in that the condition is not more than a predetermined value. In claims 1 to image reconstruction processing program according to any one of 6, c _ij recorded on the detection probability recording unit, there is a pixel i on the projection data in the projection direction of the pixel j in the reconstructed image An image reconstruction processing program characterized in that the value is set to 1 when the processing is performed and set to 0 otherwise. In the image reconstruction processing program according to any one of claims 1 to 7, wherein the first time estimation RI concentration initial value recorded in the recording unit the λ _{j 0} ^[h] is the number of pixels recorded in advance from j = 1 An image reconstruction processing program, wherein 1 is recorded up to several m of pixels j recorded in a section. A computer-readable recording medium on which the program according to claim 1 is recorded.

Images