JP6761610B2 - 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置 - Google Patents

吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置 Download PDF

Info

Publication number
JP6761610B2
JP6761610B2 JP2019521549A JP2019521549A JP6761610B2 JP 6761610 B2 JP6761610 B2 JP 6761610B2 JP 2019521549 A JP2019521549 A JP 2019521549A JP 2019521549 A JP2019521549 A JP 2019521549A JP 6761610 B2 JP6761610 B2 JP 6761610B2
Authority
JP
Japan
Prior art keywords
image
absorption coefficient
value
projection data
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2019521549A
Other languages
English (en)
Other versions
JPWO2018220686A1 (ja
Inventor
哲哉 小林
哲哉 小林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shimadzu Corp
Original Assignee
Shimadzu Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shimadzu Corp filed Critical Shimadzu Corp
Publication of JPWO2018220686A1 publication Critical patent/JPWO2018220686A1/ja
Application granted granted Critical
Publication of JP6761610B2 publication Critical patent/JP6761610B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G06T7/0014—Biomedical image inspection using an image reference approach
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10104—Positron emission tomography [PET]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Medical Informatics (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Optics & Photonics (AREA)
  • Quality & Reliability (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine (AREA)

Description

本発明は、ポジトロンＣＴ装置（陽電子放射断層撮影装置）の計測データから吸収係数画像を推定する、吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置に関する。
ポジトロンＣＴ装置、すなわちＰＥＴ(Positron Emission Tomography)装置は、陽電子（Positron）の消滅によって発生する２本のγ線を、複数個の検出器で同時に検出したときのみ（つまり同時計数したときのみ）有効な信号と見なして計測し、計測データに基づいて被検体の断層画像を再構成するように構成されている。具体的には、陽電子放出核種を含んだ放射性薬剤を被検体に投与して、投与された被検体から放出される511keVの対消滅γ線を多数の検出器素子（例えばシンチレータ）群からなる検出器で検出する。そして、２つの検出器で一定時間内にγ線を検出した場合に“同時”に検出したとして、それを一対の対消滅γ線として計数し、さらに対消滅発生位置を検出した２つの検出器を結ぶ直線（LOR: Line Of Response）を特定する。このように検出された同時計数情報を蓄積して再構成処理を行って、陽電子放出核種画像（すなわち断層画像）を得る。
ポジトロンＣＴ（ＰＥＴ）において、被検体内の放射能濃度を定量計測するためには、様々なデータ補正処理が必要である。代表的な補正処理には、感度補正，散乱補正，ランダム補正，減衰補正，デッドタイム補正，吸収補正がある。本発明は、放射性薬剤（ラジオアイソトープ）から放出されるγ線の吸収による画像の定量性低下を防止するための吸収補正に関する。吸収補正を行うためには被検体内の吸収係数分布を画像化した吸収係数画像を推定する必要がある。
吸収補正では、推定した吸収係数画像からγ線の透過率を求め、ＰＥＴの計測データから透過率を除算することでγ線の吸収の影響が排除されたデータに変換する。もしくは、推定した吸収係数画像を画像再構成の計算式に組み込んでγ線の吸収の影響が排除された再構成画像を取得する。
吸収係数画像を推定するためには、陽電子放出核種の外部線源を照射して得られたトランスミッション(Transmission)データが必要である。もしくは、トランスミッションデータの替わりにＸ線ＣＴ(Computed Tomography)装置から得られたＣＴデータを用いて、吸収係数画像を推定することが可能である。
近年では、このようなトランスミッションデータが不要な再構成アルゴリズムがある（例えば、非特許文献１，２参照）。非特許文献１，２では、消滅放射線の検出時間差（「飛行時間差」とも呼ばれる）(TOF: Time Of Flight)情報を計測するＰＥＴ（以下、「TOF-PET」と表記する）の計測データから、吸収係数サイノグラムの分布形状を推定することができる。そして、放射能画像および吸収係数サイノグラムを同時に推定することができる。非特許文献１，２における再構成アルゴリズムは、放射能画像および吸収係数に関するデータ（例えば吸収係数サイノグラム）を同時に推定することから「同時再構成アルゴリズム」とも呼ばれている。特に、放射能画像および吸収係数サイノグラムを同時に推定する同時再構成アルゴリズムは、MLACF法とも呼ばれている。
MLACF法における原理について、図９の概念図を参照して定性的に説明する。被検体の対象領域において放射能濃度が、図９に示すように分布しているとする。TOF情報t，半径方向sからなる２次元分布に放射能分布を展開すると、投影方向θ=0°のときには図９の上図のようになり、投影方向θ=90°のときには図９の右図のようになる。投影方向θを0°-180°の範囲で、それぞれの投影角度θ毎の２次元分布が計測データとして得られる。
もしγ線の吸収がなければ、投影角度θ毎の２次元分布を復元すると元の放射能分布に変換することができる。しかし、実際にはγ線の吸収により、元の放射能分布と異なる放射能分布に復元されてしまう。また、投影方向に応じてγ線の吸収の度合いが異なる。図９に示すように、投影方向θ=0°のときの吸収の度合いをAとし、投影方向θ=90°のときの吸収の度合いをA'とすると、通常はA≠A'となる。そこで、TOF情報の計測データを利用して、その計測データと、未知の吸収係数値と投影データとの乗算で得られた推定データとの誤差を最小にすることで、放射能画像と同時に、吸収係数値をサイノグラム形式に変換した吸収係数サイノグラムを推定することができる。
ただし、推定できるのは分布形状（つまり、相対値）のみで、絶対値を推定できない。吸収係数サイノグラムの絶対値が判らないと吸収補正を正しく行えない。したがって、定量的な放射能画像を取得できない。これは臨床上問題であり、実用化の大きな障壁となる。
同時再構成アルゴリズムに基づくPETイメージング技術を実用化するためには、下記のような技術が必要不可欠である。すなわち、分布形状は正しいが、絶対値が誤っている吸収係数サイノグラムから、絶対値が正しい吸収係数画像を推定する技術が必要不可欠である。
そこで、非特許文献２の従来技術として、吸収係数サイノグラムの直接的な定量化を行わずに、放射能画像を直接的に定量化するという方法がある。処理ステップは以下の通りである。
(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)を体軸方向スライス毎にそれぞれ乗算することで、当該放射能画像をスケーリングする。
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
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.
しかしながら、上述した非特許文献２の従来技術では、最終的に求まる放射能画像が定量的であることを原理的に保証できないという問題点がある。
本発明は、このような事情に鑑みてなされたものであって、定量的な吸収係数画像を作成することができる吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置を提供することを目的とする。
発明者は、上記の問題を解決するために鋭意研究した結果、次のような知見を得た。
すなわち、上述した従来技術は、手順(3)で求めた体軸方向スライス毎の合計画素値が真である（すなわち正しい）ということを仮定している。
しかし、手順(3)で求めた体軸方向スライス毎の合計画素値が正しいという保証はない。その正しい保証のない結果に基づいて定量化しようとしているので、手順(6)で求まる放射能画像の定量性も当然ながら保証できない。したがって、従来技術では定量的な放射能画像を作成することができないという問題点が生じる。つまり、従来技術の問題の根源は、経験的な方法であり、吸収係数の定量化という課題の背景にある数学的な理論に基づいていないことであるという知見を得た。
このような知見に基づく本発明は、次のような構成をとる。
すなわち、本発明に係る吸収係数画像推定方法は、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データから吸収係数画像を推定する方法であって、μ'を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、前記計測データに関する評価関数の最適化に基づいて、前記画像μ'を計算する再構成計算工程と、前記計測データに基づいて投影データ空間における被検体マスクデータである被検体マスク投影データを算出するマスク算出工程と、μoffを不均一なオフセット画像としたときに、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、前記オフセット画像μoffを推定するオフセット推定工程と、Ωを既知の吸収係数値で近似可能な領域としたときに、前記計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、前記領域Ωを少なくとも１つ以上抽出する参照領域抽出工程と、αを係数としたときに、前記領域Ωにおける前記画像μ'の値と既知の吸収係数値との誤差を減少させる前記係数αを算出する係数算出工程と、前記画像μ'の値に、前記オフセット画像μoffを前記係数α倍したα×μoffを加算して得られた値を吸収係数値として補正する吸収係数値補正工程とを備えるものである。
本発明に係る吸収係数画像推定方法によれば、再構成計算工程では、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データ（TOF-PETの計測データ）に関する評価関数の最適化に基づいて、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像を計算する。それと同時に、λ'を非定量的な放射能画像としたときに、再構成計算工程では放射能画像λ'を計算する。ここでの再構成計算工程における再構成アルゴリズムは、上述した同時再構成アルゴリズムである。ここで、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像をμ'とすると、画像μ'も吸収係数画像となるが、最終的に求まる定量的な吸収係数画像と区別して、単に「画像μ'」とする。
一方、マスク算出工程では、計測データに基づいて投影データ空間における被検体マスクデータ（以下、「被検体マスク投影データ」と呼ぶ）を算出する。μoffを不均一なオフセット画像としたときに、オフセット推定工程では、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。
一方、Ωを既知の吸収係数値で近似可能な領域としたときに、計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、参照領域抽出工程では、領域Ωを少なくとも１つ以上抽出する。
ここで、「計測データに基づいて計算された、被検体領域が認識可能な画像」とは、例えば、上述した画像μ'，上述した放射能画像λ'，上述した再構成計算工程における再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズムで推定した非定量的な画像，上述した被検体マスク投影データから推定した被検体マスク画像である。もちろん、これらの例示した画像に限定されず、計測データに基づいて計算された、被検体領域が認識可能な画像であればよい。また、「計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する」ことは、単一の画像情報を用いて領域Ωを抽出することのみならず、複数の画像情報の組み合わせ（例えば論理和や論理積）を用いて領域Ωを抽出することも含まれる。
上述した各工程で得られた画像μ'，オフセット画像μoffおよび領域Ωを用いて、下記の係数算出工程および吸収係数値補正工程を実施する。μを未知の真の吸収係数画像の値（吸収係数値）とし、αを係数としたときに、係数算出工程では、領域Ωにおける画像μ'の値と既知の吸収係数値との誤差を減少させる係数αを算出する。そして、吸収係数値補正工程では、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値として補正する。このように、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能である（つまり、μ≒μ'+α×μoff）という数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
上述した再構成計算工程を、(a)画像μ'を未知数に含む計算アルゴリズムで実施してもよい。または、上述した再構成計算工程を、(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムの組み合わせで実施してもよい。
前者の(a)のアルゴリズムは、非定量的な放射能画像λ'および画像μ'（非定量的な吸収係数画像）を同時に再構成するMLAA法である。後者の(b)のアルゴリズムは、非定量的な放射能画像λ'および非定量的な吸収係数投影データ（例えば吸収係数サイノグラム）を同時に推定する、非特許文献１で述べたMLACF法と、吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムとの組み合わせである。後者の(b)における、吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムについては、再構成アルゴリズムであれば特に限定されない。
上述したマスク算出工程の一例は、画像μ'の二値化画像を被検体マスク画像として算出する工程と、被検体マスク画像の投影データを算出する工程と、被検体マスク画像の投影データの二値化データを被検体マスク投影データとして算出する工程とからなる（(A)の態様）。また、マスク算出工程は、画像μ'の投影データを算出する工程と、画像μ'の投影データの二値化データを被検体マスク投影データとして算出する工程とからなる(B)の態様でもよい。上記(A)の態様では、画像μ'を先に二値化した後に投影データを算出し、その投影データを二値化した。上記(B)の態様では、画像μ'の投影データを先に算出した後に二値化する。画像μ'は定量的でない吸収係数画像であるが、吸収係数画像以外の画像や、投影データを二値化しても、被検体マスク投影データを算出することができる。
マスク算出工程の他の一例は、上述したTOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出する工程からなる。この一例の場合には、TOF-PETの計測データを直接的に用いて、投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出することができる。
また、マスク算出工程のさらなる他の一例は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の投影データを算出する工程と、（放射能画像の）投影データを二値化したデータを被検体マスク投影データとして算出する工程とからなる（(C)の態様）。また、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の二値化画像を算出する工程と、二値化画像の投影データを算出する工程と、二値化画像の投影データを二値化したデータを被検体マスク投影データとして算出する工程とからなる(D)の態様でもよい。上記(C)の態様では、放射能画像の投影データを先に算出した後に二値化していた。上記(D)の態様では、放射能画像を先に二値化した後に投影データを算出し、その投影データを二値化する。上記(C)の態様または上記(D)の態様のいずれにおいても、この一例の場合には、放射能画像を算出する再構成アルゴリズムについては、特に限定されない。上述した同時再構成アルゴリズム（MLACF法やMLAA法）を用いた場合には、同時再構成アルゴリズムで推定した上述の放射能画像λ'を利用して、被検体マスク投影データを算出することができる。また、上述した同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した放射能画像を利用して、被検体マスク投影データを算出してもよい。
上述したオフセット推定工程で実施する再構成処理は、解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施すればよい。解析的再構成として例えばFBP(Filtered Back Projection)法がある。統計的再構成として例えば上述したML-EM(Maximum Likelihood - Expectation Maximization)法がある。代数的再構成として例えばART(Algebraic Reconstruction Technique)法がある。
上述した参照領域抽出工程において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である。ここで、「吸収係数を既知と見なせる組織の領域」とは、例えば、水で近似可能な領域，空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域，肺組織で近似可能な領域，軟部組織で近似可能な領域などである。もちろん、これらの例示した組織に限定されず、吸収係数の近似値が判る組織であればよい。なお、脳組織の場合には水で近似可能な領域と見なされるので、水の吸収係数値である0.0096/mmを用いることができる。
上述した係数算出工程における係数αについては、下記のように代表値に基づいて算出する手法（前者の手法）と、下記のように誤差評価関数に基づいて算出する手法（後者の手法）とがある。
前者の手法においては、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の領域Ωとし、μｎ ｋｎｏｗｎ(n=1,…,K)を領域Ωｎの既知の吸収係数値とし、S(X；Ωｎ)を、領域Ωｎにおける画像Xの統計量または統計量から算出される値の代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値の代表値とし、αｎを領域Ωｎにおける係数αとする。ここで、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。もちろん、これらの例示した値に限定されず、「統計量または統計量から算出される値」であればよい。
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnの既知の吸収係数値μn knownは、領域Ωnにおける画像μ'の代表値S(μ';Ωn)に、領域Ωnにおけるオフセット画像μoff'の代表値S(μoff;Ωn)に領域Ωnにおける係数αnを乗算したαn×S(μoff;Ωn)を加算して得られた値と等しい。つまり、μn known=S(μ';Ωn)+αn×S(μoff;Ωn)が成立する。したがって、αn=(μn known-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)を係数αとすることで、係数αを算出することができる。
後者の手法においては、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとし、μn known (n=1,…,K)を領域Ωn内に既知の吸収係数を設定した画像とし、DΩn(X,Y)を領域Ωn内に関する画像Xおよび画像Yの誤差評価関数とし、wn (n=1,…,K)を０以上１以下の係数とする。画像Xを領域Ωn内に既知の吸収係数を設定した画像μn knownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩn(μn known,μ'+α×μoff)は、領域Ωn内に関する既知の吸収係数を設定した画像μn knownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。n=1,…,Kで各領域Ω1,…,ΩK毎の係数w1,…,wKによるDΩ1(μ1 known,μ'+α×μoff),…,DΩK(μK known,μ'+α×μoff)の重み付き加算で得られた、係数αを変数とした関数f(α)(=Σn=1,…,K [wn×DΩn(μn known,μ'+α×μoff)])を最小化するαを算出することで、係数αを算出することができる。
また、本発明に係る吸収係数画像推定プログラムは、本発明に係る吸収係数画像推定方法をコンピュータに実行させるものである。
本発明に係る吸収係数画像推定プログラムによれば、本発明に係る吸収係数画像推定方法をコンピュータに実行させることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
また、本発明に係るポジトロンＣＴ装置は、本発明に係る吸収係数画像推定プログラムを搭載したポジトロンＣＴ装置において、当該吸収係数画像推定プログラムを実行する演算手段を備えるものである。
本発明に係るポジトロンＣＴ装置によれば、本発明に係る吸収係数画像推定プログラムを実行する演算手段を備えることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本発明に係る吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンＣＴ装置によれば、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、吸収係数値補正工程で吸収係数値を補正する。したがって、吸収係数値補正工程で補正された吸収係数値を有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
各実施例に係るＰＥＴ装置の概略斜視図およびブロック図である。 γ線検出器の概略斜視図である。 実施例１に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 実施例２に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 MLACF法で推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 実施例４に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 実施例５に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。 MLACF法における原理の説明に供する概念図である。 体軸方向スライス毎の合計画素値、および横軸を体軸方向，縦軸を合計画素値としたときのプロファイルの概略図である。
以下、図面を参照して本発明の実施例１を説明する。図１は、各実施例に係るＰＥＴ装置の概略斜視図およびブロック図であり、図２は、γ線検出器の概略斜視図である。また、図１および図２は各実施例とも共通の構成である。
図１に示すように、ＰＥＴ装置１は、被検体の周囲を取り囲む検出器リング２を被検体の体軸方向に積層配置して備えている。検出器リング２内には複数のγ線検出器３が埋設されている。ＰＥＴ装置１は、本発明におけるポジトロンＣＴ装置に相当する。また、γ線検出器３は、本発明における検出器に相当する。
その他にも、ＰＥＴ装置１は、同時計数回路４と演算回路５とを備えている。図１では、γ線検出器３から同時計数回路４への結線を２つのみ図示しているが、実際には、γ線検出器３の光電子増倍管(PMT: Photo Multiplier Tube)３３（図２を参照）の総チャンネル数分、同時計数回路４に接続されている。演算回路５は、吸収係数画像推定プログラム６による、後述する図３に示す吸収係数画像推定方法の処理を実行する。演算回路５は、本発明における演算手段に相当する。
放射性薬剤が投与された被検体（図示省略）から発生したγ線をγ線検出器３のシンチレータブロック３１（図２を参照）が光に変換して、変換されたその光をγ線検出器３の光電子増倍管（ＰＭＴ）３３（図２を参照）は増倍させて電気信号に変換する。その電気信号を同時計数回路４に送り込み、カウント値の検出信号データを生成する。
具体的には、被検体（図示省略）に放射性薬剤を投与すると、ポジトロン放出型のＲＩのポジトロンが消滅することにより、２本のγ線が発生する。同時計数回路４は、シンチレータブロック３１（図２を参照）の位置とγ線の入射タイミングとをチェックし、被検体の両側にある２つのシンチレータブロック３１でγ線が同時に入射したときのみ、送り込まれた電気信号を適正なデータと判定する。一方のシンチレータブロック３１のみにγ線が入射したときには、同時計数回路４は棄却する。つまり、同時計数回路４は、上述した電気信号に基づいて、２つのγ線検出器３においてγ線が同時観測（すなわち同時計数）されたことを検出する。
同時計数回路４により同時計数と判定された適正なデータにより構成された検出信号データ（カウント値）を、演算回路５に送り込む。演算回路５は、後述するステップＳ１〜Ｓ８（図３を参照）を行って、ＰＥＴ装置１によって得られた被検体（図示省略）の検出信号データから吸収係数画像を推定する。演算回路５の具体的な機能については後述する。
なお、ＲＯＭ（Read-only Memory）などに代表される記憶媒体（図示省略）に吸収係数画像推定プログラム６が記憶されており、当該吸収係数画像推定プログラム６を記憶媒体から演算回路５に読み出して、吸収係数画像推定プログラム６を演算回路５が実行することで、図３のフローチャートに示す吸収係数画像推定方法の処理を行う。演算回路５は、ＧＰＵ（Graphics Processing Unit），中央演算処理装置（ＣＰＵ）あるいはプログラムデータに応じて内部の使用するハードウェア回路（例えば論理回路）が変更可能なプログラマブルデバイス（例えばＦＰＧＡ(Field Programmable Gate Array)）などで構成されている。
γ線検出器３は、図２に示すようにシンチレータブロック３１と、そのシンチレータブロック３１に対して光学的に結合されたライトガイド３２と、そのライトガイド３２に対して光学的に結合された光電子増倍管（以下、単に「ＰＭＴ」と略記する）３３とを備えている。シンチレータブロック３１を構成する各シンチレータ素子は、γ線の入射に伴って発光することでγ線から光に変換する。この変換によってシンチレータ素子はγ線を検出する。シンチレータ素子において発光した光がシンチレータブロック３１で十分に拡散されて、ライトガイド３２を介してＰＭＴ３３に入力される。ＰＭＴ３３は、シンチレータブロック３１で変換された光を増倍させて電気信号に変換する。その電気信号は画素値として同時計数回路４（図１を参照）に送り込まれる。
また、γ線検出器３は、図２に示すように、３次元的に配置されたシンチレータ素子からなり、深さ方向に複数の層からなるＤＯＩ検出器である。図２では、４層のＤＯＩ検出器を図示しているが、層の数については、複数であれば特に限定されない。
ここで、ＤＯＩ検出器は、各々のシンチレータ素子を放射線の深さ方向に積層して構成されたものであり、相互作用を起こした深さ(DOI: Depth of Interaction)方向と横方向（入射面に平行な方向）との座標情報を重心演算により求める。ＤＯＩ検出器を用いることにより深さ方向の空間分解能をより一層向上させることができる。よって、ＤＯＩ検出器の層の数は、深さ方向に積層されたシンチレータ素子の層の数である。
次に、演算回路５の具体的な機能について、図３を参照して説明する。図３は、実施例１に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
先ず、図１に示すＰＥＴ装置１によって被検体の撮影を行い、同時計数回路４（図１を参照）によってリストモードデータを取得する。リストモードデータには、検出した光子のエネルギー情報が記録されている。
通常のエネルギーウィンドウ（例えば400keV-600keV）、すなわち再構成データ用のエネルギーウィンドウ、およびTOF計測データのTOF方向の計測範囲とビン幅とをそれぞれ設定する。なお、ここで、ビン(bin)とは、離散化すること（区切ること）を意味する。画像の場合には画素がビンに対応する。TOFビンはTOF情報の時間的区切りを意味し、例えばTOFビンが100[ps]の場合には、検出時間差は100[ps]の精度で時間的に区切られる。
その設定にしたがって、リストモードデータから、TOF-PETの計測データを作成する。
（ステップＳ１）MLACF
μ'を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、λ'を非定量的な放射能画像とする。計測データに関する評価関数の最適化に基づいて、画像μ'および放射能画像λ'を同時に計算する。「課題を解決するための手段」の欄でも述べたように、画像μ'も吸収係数画像となるが、最終的に求まる定量的な吸収係数画像と区別して、単に「画像μ'」とする。
ステップＳ１での再構成アルゴリズムは、非特許文献１での同時再構成アルゴリズムである。後述する実施例２，３も含めて、本実施例１では、同時再構成アルゴリズムとして、非特許文献１で述べたMLACF法を適用する。MLACF法の具体的な手法については、上述した非特許文献１を参照されたい。
A'を吸収係数サイノグラムとする。MLACF法によって、放射能画像λ'および吸収係数サイノグラムA'を推定する。吸収係数サイノグラムA'は、本発明における吸収係数投影データに相当する。
（ステップＳ２）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）。ステップＳ１，Ｓ２は、本発明における再構成計算工程に相当する。
（ステップＳ３）Binarization Processing
閾値処理によって画像μ'を二値化処理(Binarization Processing)する。そして、被検体領域を“１”，その他の領域を“０”となる二値化画像を被検体マスク画像として算出する。mimgを被検体マスク画像とする。
（ステップＳ４）Projection + Binarization Processing
被検体マスク画像mimgの線積分データ（投影データ）を算出する（Projection）。そして、閾値処理によって被検体マスク画像mimgの投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データとして算出する。mprojを被検体マスク投影データとする。ステップＳ３，Ｓ４は、本発明におけるマスク算出工程に相当する。
（ステップＳ５）ML-EM
μoffを不均一なオフセット画像とする。オフセット画像μoffの順投影データが被検体マスク投影データmprojを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。つまり、被検体マスク投影データmprojを再構成アルゴリズム（例えばML-EM法）で画像データに変換する。この画像データをオフセット画像μoffとする。ステップＳ５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ６）Extraction
Ωを既知の吸収係数値で近似可能な領域とする。計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する（Extraction）。後述する実施例２〜４も含めて、本実施例１では、K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとしたときに、n=1,…,Kで各領域Ω1,…,ΩKをそれぞれ抽出する。例えば、被検体の頭部を撮影する場合には、空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域で各領域Ω1,…,ΩKをそれぞれ抽出する。
計測データに基づいて計算された、被検体領域が認識可能な画像としては、例えば、ステップＳ２で作成した画像μ'，ステップＳ１で求まった放射能画像λ'，ステップＳ１での再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した非定量的な画像，上述した被検体マスク投影データmprojから推定した被検体マスク画像mimgである。ステップＳ６は、本発明における参照領域抽出工程に相当する。
（ステップＳ７）αn=(μn known-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)
μn known(n=1,…,K)を領域Ωnの既知の吸収係数値とし、S(X;Ωn)を、領域Ωnにおける画像Xの統計量または統計量から算出される値を代表値とし、T(x1,x2,…,xK)を任意のK個の値x1,x2,…,xKの統計量または統計量から算出される値を代表値とし、αnを領域Ωnにおける係数αとする。なお、αは係数である。ここで、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。ここで、「トリムド平均値」とは、値が極端に大きい／小さいデータを取り除いた残りのデータの平均値を意味する。「トリムド中央値」とは、値が極端に大きい／小さいデータを取り除いた残りのデータの中央値を意味する。
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnにおける係数αnは下記（１）式のように表される。
αn=(μn known-S(μ';Ωn))/S(μoff;Ωn) …（１）
K個の値x1,x2,…,xKをα1,α2,…,αKにそれぞれ置き換えれば、T(α1,α2,…,αK)は係数α1,α2,…,αKの代表値である。したがって、係数αは下記（２）式のように表される。ステップＳ７は、本発明における係数算出工程に相当する。
α=T(α1,α2,…,αK) …（２）
（ステップＳ８）μ=μ'+α×μoff
μを未知の真の吸収係数画像の値（吸収係数値）とする。吸収係数値μは下記（３）式のように表される。
μ=μ'+α×μoff …（３）
上記（３）式のように、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値μとして補正する。ステップＳ８は、本発明における吸収係数値補正工程に相当する。
ステップＳ８で求まった吸収係数値μを有する吸収係数画像を用いて吸収補正を行う。「背景技術」の欄でも述べたように、推定した吸収係数画像からγ線の透過率を求め、ＰＥＴの計測データから透過率を除算することでγ線の吸収の影響が排除されたデータに変換することによって、吸収補正を行う。もしくは、推定した吸収係数画像を画像再構成の計算式に組み込んでγ線の吸収の影響が排除された再構成画像を取得することによって、吸収補正を行う。
本実施例１に係る吸収係数画像推定方法によれば、ステップＳ１，Ｓ２では、消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データ（TOF-PETの計測データ）に関する評価関数の最適化に基づいて、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像μ'を計算する。それと同時に、ステップＳ１では放射能画像λ'を計算する。ステップＳ１における再構成アルゴリズムは、上述した同時再構成アルゴリズムである。
一方、ステップＳ３，Ｓ４では、計測データに基づいて投影データ空間における被検体マスクデータ（すなわち、被検体マスク投影データmproj）を算出する。ステップＳ５では、オフセット画像μoffの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、オフセット画像μoffを推定する。
一方、計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、ステップＳ６では、領域Ωを少なくとも１つ以上抽出する。
ここで、「課題を解決するための手段」の欄でも述べたように、「計測データに基づいて計算された、被検体領域が認識可能な画像」とは、例えば、ステップＳ２で作成した画像μ'，ステップＳ１で求まった放射能画像λ'，ステップＳ１での再構成アルゴリズム（同時再構成アルゴリズム）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した非定量的な画像，上述した被検体マスク投影データmprojから推定した被検体マスク画像mimgである。もちろん、これらの例示した画像に限定されず、計測データに基づいて計算された、被検体領域が認識可能な画像であればよい。また、「課題を解決するための手段」の欄でも述べたように、「計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、領域Ωを少なくとも１つ以上抽出する」ことは、単一の画像情報を用いて領域Ωを抽出することのみならず、複数の画像情報の組み合わせ（例えば論理和や論理積）を用いて領域Ωを抽出することも含まれる。
上述した各ステップＳ１〜Ｓ６で得られた画像μ'，オフセット画像μoffおよび領域Ωを用いて、ステップＳ７，Ｓ８を実施する。ステップＳ７では、領域Ωにおける画像μ'の値と既知の吸収係数値との誤差を減少させる係数αを算出する。そして、ステップＳ８では、上記（３）式のように、画像μ'の値に、オフセット画像μoffを係数α倍したα×μoffを加算して得られた値を吸収係数値μとして補正する。このように、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能である（つまり、μ≒μ'+α×μoff）という数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
後述する実施例２，３も含めて、本実施例１では、上述した再構成計算工程に相当するステップＳ１，Ｓ２を、(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび吸収係数投影データを再構成した画像を画像μ'とするアルゴリズムの組み合わせで実施している。
(b)のアルゴリズムは、非定量的な放射能画像λ'および非定量的な吸収係数投影データ（各実施例１〜３では吸収係数サイノグラムA'）を同時に推定する、非特許文献１で述べたMLACF法と、吸収係数投影データ（吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムとの組み合わせである。「課題を解決するための手段」の欄でも述べたように、(b)における、吸収係数投影データ（吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムについては、再構成アルゴリズムであれば特に限定されない。図３では、ステップＳ２で述べたように、ML-TR法またはML-EM法で吸収係数サイノグラムA'を再構成し、画像μ'を作成する。
本実施例１では、上述したマスク算出工程に相当するステップＳ３，Ｓ４を実施している。すなわち、マスク算出工程は、画像μ'の二値化画像を被検体マスク画像mimgとして算出する工程（ステップＳ３）と、被検体マスク画像mimgの投影データを算出する工程（ステップＳ４の前半の「Projection」）と、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出する工程（ステップＳ４の後半の「Binarization Processing」）とからなる。画像μ'は定量的でない吸収係数画像であるが、後述する実施例２，３のように、吸収係数画像以外の画像（例えば実施例３の放射能画像λ'，または放射能画像λ２'）や、投影データを二値化しても、被検体マスク投影データmprojを算出することができる。
上述したオフセット推定工程に相当するステップＳ５で実施する再構成処理は、本実施例１ではML-EM(Maximum Likelihood - Expectation Maximization)法等に代表される統計的再構成の計算方式で実施している。もちろん、ステップＳ５で実施する再構成処理は、本実施例１のような統計的再構成に限定されない。解析的再構成，統計的再構成，代数的再構成のいずれかの計算方式で実施すればよい。解析的再構成として例えばFBP(Filtered Back Projection)法がある。代数的再構成として例えばART(Algebraic Reconstruction Technique)法がある。
上述した参照領域抽出工程に相当するステップＳ６において抽出される少なくとも１つ以上の領域Ωは、吸収係数を既知と見なせる組織の領域である。ここで、「課題を解決するための手段」の欄でも述べたように、「吸収係数を既知と見なせる組織の領域」とは、例えば、水で近似可能な領域，空気で近似可能な領域，脳組織で近似可能な領域，骨で近似可能な領域，肺組織で近似可能な領域，軟部組織で近似可能な領域などである。もちろん、これらの例示した組織に限定されず、吸収係数の近似値が判る組織であればよい。なお、脳組織の場合には水で近似可能な領域と見なされるので、水の吸収係数値である0.0096/mmを用いることができる。
本実施例１では、上述した係数算出工程に相当するステップＳ７を実施している。すなわち、後述する実施例２〜４も含めて、本実施例１では、代表値に基づいて係数αを算出している。上述したように、代表値としては、例えば、平均値，中央値，トリムド平均値，トリムド中央値またはそれらの２つ以上の値の重み付き平均値である。「課題を解決するための手段」の欄でも述べたように、これらの例示した値に限定されず、「統計量または統計量から算出される値」であればよい。
画像Xを画像μ'に置き換えれば、S(μ';Ωn)は領域Ωnにおける画像μ'の代表値であり、画像Xをオフセット画像μoff'に置き換えれば、S(μoff;Ωn)は領域Ωnにおけるオフセット画像μoff'の代表値である。領域Ωnの既知の吸収係数値μn knownは、領域Ωnにおける画像μ'の代表値S(μ';Ωn)に、領域Ωnにおけるオフセット画像μoff'の代表値S(μoff;Ωn)に領域Ωnにおける係数αnを乗算したαn×S(μoff;Ωn)を加算して得られた値と等しい。つまり、μn known= S(μ';Ωn)+αn×S(μoff;Ωn)が成立する。したがって、上記（１）式のように、αn=(μn known-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)を係数αとすることで、係数αを算出することができる。
本実施例１に係る吸収係数画像推定プログラム６（図１を参照）によれば、本実施例１に係る吸収係数画像推定プログラムをコンピュータ（各実施例では図１に示す演算回路５を構成するＧＰＵ，ＣＰＵまたはプログラマブルデバイス）に実行させることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本実施例１に係るＰＥＴ装置１（図１を参照）によれば、本実施例１に係る吸収係数画像推定プログラム６を実行する演算手段（各実施例では図１に示す演算回路５を構成するＧＰＵ，ＣＰＵまたはプログラマブルデバイス）を備えることによって、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ８で吸収係数値を補正する。したがって、ステップＳ８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
次に、図面を参照して本発明の実施例２を説明する。図４は、実施例２に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
上述した実施例１では、画像μ'の二値化画像を被検体マスク画像mimgとして算出し、被検体マスク画像mimgの投影データを算出し、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出していた。これに対して、本実施例２では、TOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出する。
（ステップＳ１１）MLACF
図４のステップＳ１１は、上述した実施例１のステップＳ１と同じであるので、その説明については省略する。
（ステップＳ１２）ML-TR or ML-EM
図４のステップＳ１２は、上述した実施例１のステップＳ２と同じであるので、その説明については省略する。ステップＳ１１，Ｓ１２は、本発明における再構成計算工程に相当する。
（ステップＳ１４）Projection + Binarization Processing
上述した実施例１では、画像μ'から被検体マスク画像mimgを算出するのに、図３のステップＳ３，Ｓ４を実施していた。本実施例２では、画像μ'の替わりに、計測データから被検体マスク画像mimgを算出するために、図４のステップＳ１４を実施する。先ず、計測データを投影データ形式に変換する（Conversion）。そして、閾値処理によって計測データの投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ１４は、本発明におけるマスク算出工程に相当する。
（ステップＳ１５）ML-EM
図４のステップＳ１５は、上述した実施例１のステップＳ５と同じであるので、その説明については省略する。ステップＳ１５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ１６）Extraction
図４のステップＳ１６は、上述した実施例１のステップＳ６と同じであるので、その説明については省略する。ステップＳ１６は、本発明における参照領域抽出工程に相当する。
（ステップＳ１７）αn=(μn known-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)
図４のステップＳ１７は、上述した実施例１のステップＳ７と同じであるので、その説明については省略する。ステップＳ１７は、本発明における係数算出工程に相当する。
（ステップＳ１８）μ=μ'+α×μoff
図４のステップＳ１８は、上述した実施例１のステップＳ８と同じであるので、その説明については省略する。ステップＳ１８は、本発明における吸収係数値補正工程に相当する。
本実施例２に係る吸収係数画像推定方法によれば、上述した実施例１と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ１８で吸収係数値を補正する。したがって、ステップＳ１８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本実施例２では、上述したマスク算出工程に相当するステップＳ１４を実施している。すなわち、マスク算出工程は、上述したTOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データとして算出する工程（ステップＳ１４）からなる。本実施例２の場合には、TOF-PETの計測データを直接的に用いて、投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出することができる。
本実施例２に係る吸収係数画像推定プログラム６（図１を参照）や本実施例２に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１と同じであるので、その説明については省略する。
次に、図面を参照して本発明の実施例３を説明する。図５は、MLACF法で推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートであり、図６は、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データを算出する場合における、実施例３に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
上述した実施例１では、画像μ'の二値化画像を被検体マスク画像mimgとして算出し、被検体マスク画像mimgの投影データを算出し、被検体マスク画像mimgの投影データの二値化データを被検体マスク投影データmprojとして算出していた。また、上述した実施例２では、TOF-PETの計測データを投影データ形式に変換したものを二値化したデータを被検体マスク投影データmprojとして算出していた。これに対して、本実施例３では、TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出し、放射能画像の投影データを算出し、（放射能画像の）投影データを二値化したデータを被検体マスク投影データmprojとして算出する。
先ず、MLACF法で推定した放射能画像を利用して被検体マスク投影データmprojを算出する場合について、図５を参照して説明する。
（ステップＳ２１）MLACF
図５のステップＳ２１は、上述した実施例１のステップＳ１，上述した実施例２のステップＳ１１と同じであるので、その説明については省略する。
（ステップＳ２２）ML-TR or ML-EM
図５のステップＳ２２は、上述した実施例１のステップＳ２，上述した実施例２のステップＳ１２と同じであるので、その説明については省略する。ステップＳ２１，Ｓ２２は、本発明における再構成計算工程に相当する。
（ステップＳ２４）Projection + Binarization Processing
上述した実施例１では、画像μ'から被検体マスク画像mimgを算出するのに、図３のステップＳ３，Ｓ４を実施していた。また、上述した実施例２では、計測データから被検体マスク画像mimgを算出するのに、図４のステップＳ１４を実施していた。本実施例３では、実施例１の画像μ'，実施例２の計測データの替わりに、計測データに関する評価関数の最適化に基づいて推定した放射能画像から被検体マスク画像mimgを算出するために、図５のステップＳ２４を実施する。
先ず、図５では、MLACF法における計測データに関する評価関数の最適化に基づいて推定した放射能画像λ'の線積分データ（投影データ）を算出する（Projection）。つまり、図５では、ステップＳ２１のMLACF法で推定した放射能画像λ'を利用して放射能画像λ'の投影データを算出すればよい。そして、閾値処理によって放射能画像λ'の投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ２４は、本発明におけるマスク算出工程に相当する。
（ステップＳ２５）ML-EM
図５のステップＳ２５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５と同じであるので、その説明については省略する。ステップＳ２５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ２６）Extraction
図５のステップＳ２６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６と同じであるので、その説明については省略する。ステップＳ２６は、本発明における参照領域抽出工程に相当する。
（ステップＳ２７）αn=(μn known-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)
図５のステップＳ２７は、上述した実施例１のステップＳ７，上述した実施例２のステップＳ１７と同じであるので、その説明については省略する。ステップＳ２７は、本発明における係数算出工程に相当する。
（ステップＳ２８）μ=μ'+α×μoff
図５のステップＳ２８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８と同じであるので、その説明については省略する。ステップＳ２８は、本発明における吸収係数値補正工程に相当する。
次に、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データmprojを算出する場合について、図６を参照して説明する。
（ステップＳ３１）MLACF
図６のステップＳ３１は、図５のステップＳ２１と同じであるので、その説明については省略する。
（ステップＳ３２）ML-TR or ML-EM
図６のステップＳ３２は、図５のステップＳ２２と同じであるので、その説明については省略する。ステップＳ３１，Ｓ３２は、本発明における再構成計算工程に相当する。
（ステップＳ３３）ML-EM
図５ではMLACF法で推定した放射能画像λ'を利用して被検体マスク投影データmprojを算出するのに、図５のステップＳ２４を実施していた。図６では、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像を利用して被検体マスク投影データmprojを算出するために、図６のステップＳ３３，Ｓ３４を実施する。
先ず、図６では、MLACF法とは異なる再構成アルゴリズム（例えばML-EM法）における計測データに関する評価関数の最適化に基づいて、放射能画像を推定する。図６ではMLACF法で推定した放射能画像λ'と区別して、MLACF法とは異なる再構成アルゴリズムで推定した放射能画像をλ2'とする。
（ステップＳ３４）Projection + Binarization Processing
MLACF法とは異なる再構成アルゴリズムにおける計測データに関する評価関数の最適化に基づいて推定した放射能画像λ2'の線積分データ（投影データ）を算出する（Projection）。そして、閾値処理によって放射能画像λ2'の投影データを二値化処理(Binarization Processing)することで、被検体を通過する投影線を“１”，その他の投影線を“０”となる二値化データを被検体マスク投影データmprojとして算出する。ステップＳ３３，Ｓ３４は、本発明におけるマスク算出工程に相当する。
（ステップＳ３５）ML-EM
図６のステップＳ３５は、図５のステップＳ２５と同じであるので、その説明については省略する。ステップＳ３５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ３６）Extraction
図６のステップＳ３６は、図５のステップＳ２６と同じであるので、その説明については省略する。ステップＳ３６は、本発明における参照領域抽出工程に相当する。
（ステップＳ３７）αn=(μn known-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)
図６のステップＳ３７は、図５のステップＳ２７と同じであるので、その説明については省略する。ステップＳ３７は、本発明における係数算出工程に相当する。
（ステップＳ３８）μ=μ'+α×μoff
図６のステップＳ３８は、図５のステップＳ２８と同じであるので、その説明については省略する。ステップＳ３８は、本発明における吸収係数値補正工程に相当する。
本実施例３に係る吸収係数画像推定方法によれば、上述した実施例１，２と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、図５のステップＳ２８，または図６のステップＳ３８で吸収係数値を補正する。したがって、図５のステップＳ２８，または図６のステップＳ３８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本実施例３では、上述したマスク算出工程に相当する図５のステップＳ２４，または図６のステップＳ３３，Ｓ３４を実施している。すなわち、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像（図５ではλ'，図６ではλ2'）を算出する工程（図５ではステップＳ２１，図６ではステップＳ３３）と、放射能画像（図５ではλ'，図６ではλ2'）の投影データを算出する工程（図５ではステップＳ２４の前半の「Projection」，図６ではステップＳ３４の前半の「Projection」）と、（放射能画像の）投影データを二値化したデータを被検体マスク投影データmprojとして算出する工程（図５ではステップＳ２４の後半の「Binarization Processing」，図６ではステップＳ３４の後半の「Binarization Processing」）とからなる。本実施例３の場合には、放射能画像を算出する再構成アルゴリズムについては、特に限定されない。図５のように、上述した同時再構成アルゴリズム（MLACF法やMLAA法）を用いた場合には、同時再構成アルゴリズムで推定した上述の放射能画像λ'を利用して、被検体マスク投影データmprojを算出することができる。また、図６のように、上述した同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）で推定した放射能画像λ2'を利用して、被検体マスク投影データmprojを算出してもよい。
本実施例３に係る吸収係数画像推定プログラム６（図１を参照）や本実施例３に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１，２と同じであるので、その説明については省略する。
次に、図面を参照して本発明の実施例４を説明する。図７は、実施例４に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
上述した実施例１〜３では、MLACF法と、吸収係数投影データ（各実施例１〜３では吸収係数サイノグラムA'）を再構成した画像を画像μ'とするアルゴリズムとの組み合わせで、上述した再構成計算工程を実施していた。これに対して、本実施例４では、(a)画像μ'を未知数に含む計算アルゴリズム（MLAA法）で実施する。
（ステップＳ４１）MLAA
上述した実施例１〜３では、画像μ'および放射能画像λ'を同時に計算するために、MLACF法によって、放射能画像λ'および吸収係数サイノグラムA'を推定したステップＳ１（実施例２ではステップＳ１１，実施例３ではステップＳ２１またはＳ３１）の後に、吸収係数サイノグラムA'を再構成した画像を画像μ'とするステップＳ２（実施例２ではステップＳ１２，実施例３ではステップＳ２２またはＳ３２）を実施していた。本実施例４では、MLAA法によって、放射能画像λ'および画像μ'を同時に計算する。
つまり、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）。ステップＳ４１は、本発明における再構成計算工程に相当する。
（ステップＳ４３）Binarization Processing
図７のステップＳ４３は、上述した実施例１のステップＳ３と同じであるので、その説明については省略する。
（ステップＳ４４）Projection + Binarization Processing
図７のステップＳ４４は、上述した実施例１のステップＳ４と同じであるので、その説明については省略する。ステップＳ４３，Ｓ４４は、本発明におけるマスク算出工程に相当する。
（ステップＳ４５）ML-EM
図７のステップＳ４５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５，上述した実施例３のステップ（図５ではステップＳ２５，図６ではステップＳ３５）と同じであるので、その説明については省略する。ステップＳ４５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ４６）Extraction
図７のステップＳ４６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６，上述した実施例３のステップ（図５ではステップＳ２６，図６ではステップＳ３６）と同じであるので、その説明については省略する。ステップＳ４６は、本発明における参照領域抽出工程に相当する。
（ステップＳ４７）αn=(μn known-S(μ';Ωn))/S(μoff;Ωn), α=T(α1,α2,…,αK)
図７のステップＳ４７は、上述した実施例１のステップＳ７，上述した実施例２のステップＳ１７，上述した実施例３のステップ（図５ではステップＳ２７，図６ではステップＳ３７）と同じであるので、その説明については省略する。ステップＳ４７は、本発明における係数算出工程に相当する。
（ステップＳ４８）μ=μ'+α×μoff
図７のステップＳ４８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８，上述した実施例３のステップ（図５ではステップＳ２８，図６ではステップＳ３８）と同じであるので、その説明については省略する。ステップＳ４８は、本発明における吸収係数値補正工程に相当する。
なお、図７は、上述した実施例１の図３におけるステップＳ１・Ｓ２をステップＳ４３に統合した場合に適用したときの実施例４のフローチャートである。もちろん、上述した実施例２の図４におけるステップＳ１１・Ｓ１２をステップＳ４３に統合した場合に実施例４を適用してもよいし、上述した実施例３の図５におけるステップＳ２１・Ｓ２２または図６におけるステップＳ３１・Ｓ３２に統合した場合に実施例４を適用してもよい。
本実施例４に係る吸収係数画像推定方法によれば、上述した実施例１〜３と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ４８で吸収係数値を補正する。したがって、ステップＳ４８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本実施例４では、上述した再構成計算工程に相当するステップＳ４１を、(a)画像μ'を未知数に含む計算アルゴリズムで実施している。(a)のアルゴリズムは、非定量的な放射能画像λ'および画像μ'（非定量的な吸収係数画像）を同時に再構成するMLAA法である。
本実施例４に係る吸収係数画像推定プログラム６（図１を参照）や本実施例４に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１〜３と同じであるので、その説明については省略する。
次に、図面を参照して本発明の実施例５を説明する。図８は、実施例５に係る吸収係数画像推定方法の処理手順およびデータの流れを示したフローチャートである。
上述した実施例１〜４では、代表値に基づいて係数αを算出していた。これに対して、本実施例５では、誤差評価関数に基づいて係数αを算出する。
（ステップＳ５１）MLACF
図８のステップＳ５１は、上述した実施例１のステップＳ１，上述した実施例２のステップＳ１１，上述した実施例３のステップ（図５ではステップＳ２１，図６ではステップＳ３１）と同じであるので、その説明については省略する。
（ステップＳ５２）ML-TR or ML-EM
図８のステップＳ５２は、上述した実施例１のステップＳ２，上述した実施例２のステップＳ１２，上述した実施例３のステップ（図５ではステップＳ２２，図６ではステップＳ３２）と同じであるので、その説明については省略する。ステップＳ５１，Ｓ５２は、本発明における再構成計算工程に相当する。
（ステップＳ５３）Binarization Processing
図８のステップＳ５３は、上述した実施例１のステップＳ３，上述した実施例４のステップＳ４３と同じであるので、その説明については省略する。
（ステップＳ５４）Projection + Binarization Processing
図８のステップＳ５４は、上述した実施例１のステップＳ４，上述した実施例４のステップＳ４４と同じであるので、その説明については省略する。ステップＳ５３，Ｓ５４は、本発明におけるマスク算出工程に相当する。
（ステップＳ５５）ML-EM
図８のステップＳ５５は、上述した実施例１のステップＳ５，上述した実施例２のステップＳ１５，上述した実施例３のステップ（図５ではステップＳ２５，図６ではステップＳ３５），上述した実施例４のステップＳ４５と同じであるので、その説明については省略する。ステップＳ５５は、本発明におけるオフセット推定工程に相当する。
（ステップＳ５６）Extraction
図８のステップＳ５６は、上述した実施例１のステップＳ６，上述した実施例２のステップＳ１６，上述した実施例３のステップ（図５ではステップＳ２６，図６ではステップＳ３６），上述した実施例４のステップＳ４６と同じであるので、その説明については省略する。ステップＳ５６は、本発明における参照領域抽出工程に相当する。
（ステップＳ５７）f(α)=DΩ(μknown,μ'+α×μoff)
上述した実施例１〜４では、代表値に基づいて係数αを算出していた。本実施例５では、誤差評価関数に基づいて係数αを算出する。なお、K(≧1)を既知の吸収係数値で近似可能な領域数とし、wn (n=1,…,K)を０以上１以下の係数としたときに、図８では、K=1（つまり、n=1のみ）として領域数を１つのみとし、係数w1=1として説明する。一般化した式については、後述する。
具体的には、μknownを領域Ω内に既知の吸収係数を設定した画像とし、DΩ(X,Y)を領域Ω内に関する画像Xおよび画像Yの誤差評価関数とする。画像Xを領域Ω内に既知の吸収係数を設定した画像μknownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩ(μknown,μ'+α×μoff)は、領域Ω内に関する既知の吸収係数を設定した画像μknownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。f(α)を、係数αを変数とした関数とすると、関数f(α)は下記（４）式のように表される。
f(α)=DΩ(μknown,μ'+α×μoff) …（４）
上記（４）式で表された関数f(α)を最小化するαを算出する。誤差評価関数DΩ(X,Y)は、例えば、二次関数DΩ(X,Y)=Σn∈Ω(X-Y)2である。もちろん、画像Xおよび画像Yの誤差を評価する関数であれば、絶対値の差分和DΩ(X,Y)=Σn∈Ω|X-Y|に例示されるように、誤差評価関数DΩ(X,Y)については特に限定されない。ステップＳ５７は、本発明における係数算出工程に相当する。
（ステップＳ５８）μ=μ'+α×μoff
図８のステップＳ４８は、上述した実施例１のステップＳ８，上述した実施例２のステップＳ１８，上述した実施例３のステップ（図５ではステップＳ２８，図６ではステップＳ３８），上述した実施例４のステップＳ４８と同じであるので、その説明については省略する。ステップＳ５８は、本発明における吸収係数値補正工程に相当する。
なお、図８は、上述した実施例１の図３におけるステップＳ７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に適用したときの実施例５のフローチャートである。もちろん、上述した実施例２の図４におけるステップＳ１７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよいし、上述した実施例３の図５におけるステップＳ２７または図６におけるステップＳ３７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよいし、上述した実施例４の図７におけるステップＳ４７による代表値に基づく係数αの算出を、ステップＳ５７による誤差評価関数に基づく係数αの算出に置き換えた場合に実施例５を適用してもよい。
本実施例５に係る吸収係数画像推定方法によれば、上述した実施例１〜４と同様に、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ５８で吸収係数値を補正する。したがって、ステップＳ５８で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。
本実施例５では、誤差評価関数に基づいて係数αを算出している。上記（４）式を一般化した場合について説明する。K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωnをn番目の領域Ωとし、μn known(n=1,…,K)を領域Ωn内に既知の吸収係数を設定した画像とし、DΩn(X,Y)を領域Ωn内に関する画像Xおよび画像Yの誤差評価関数とし、wn(n=1,…,K)を０以上１以下の係数とする。画像Xを領域Ωn内に既知の吸収係数を設定した画像μn knownに置き換え、画像Yを、オフセット画像μoff'に係数α倍したα×μoffを画像μ'に加算して得られた値(μ'+α×μoff)に置き換えれば、DΩn(μn known,μ'+α×μoff)は、領域Ωn内に関する既知の吸収係数を設定した画像μn knownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数となる。このとき、係数αを変数とした関数f(α)は、n=1,…,Kで各領域Ω1,…,ΩK毎の係数w1,…,wKによるDΩ1(μ1 known,μ'+α×μoff),…,DΩK(μK known,μ'+α×μoff)の重み付き加算で表される。すなわち、一般化した場合には、関数f(α)は下記（５）式のように表される。
f(α)=Σn=1,…,K[wn×DΩn(μn known,μ'+α×μoff)] …（５）
上記（５）式で表された関数f(α)を最小化するαを算出する。
本実施例５に係る吸収係数画像推定プログラム６（図１を参照）や本実施例５に係るＰＥＴ装置１（図１を参照）の作用・効果については、上述した実施例１〜４と同じであるので、その説明については省略する。
本発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
（１）上述した各実施例において撮影対象については、特に限定されない。被検体の全身を撮影する装置や、被検体の頭部を撮影する装置や、被検体の乳房を撮影する装置に適用してもよい。
（２）上述した各実施例では、ＤＯＩ検出器であったが、深さ方向を弁別しない検出器に適用してもよい。
（３）上述した各実施例では、検出器リングを被検体の体軸方向に積層配置する構成であったが、検出器リングが１つのみの構成であってもよい。
（４）上述した各実施例では、データ形式がサイノグラムのときであったが、サイノグラムに限定されない。投影データであれば、例えばデータ形式がヒストグラムであってもよい。したがって、上述した実施例１〜３では、MLACF法で推定した吸収係数投影データは吸収係数サイノグラムA'であったが、MLACF法で推定した吸収係数投影データは吸収係数ヒストグラムであってもよい。
（５）上述した実施例１では、画像μ'を先に二値化した後に投影データを算出し、その投影データを二値化したが、下記のように被検体マスク投影データmprojを算出してもよい。すなわち、マスク算出工程は、画像μ'の投影データを算出する工程と、画像μ'の投影データの二値化データを被検体マスク投影データmprojとして算出する工程とからなる態様でもよい。この態様の場合には、画像μ'の投影データを先に算出した後に二値化している。
（６）上述した実施例３では、放射能画像（図５ではλ'，図６ではλ2'）の投影データを先に算出した後に二値化していたが、下記のように被検体マスク投影データmprojを算出してもよい。すなわち、マスク算出工程は、（上述した）TOF-PETの計測データに関する評価関数の最適化に基づいて、放射能画像を算出する工程と、放射能画像の二値化画像を算出する工程と、二値化画像の投影データを算出する工程と、二値化画像の投影データを二値化したデータを被検体マスク投影データmprojとして算出する工程とからなる態様でもよい。この態様の場合には、放射能画像を先に二値化した後に投影データを算出し、その投影データを二値化している。上述した実施例３と同様に、この態様の場合においても、放射能画像を算出する再構成アルゴリズムについては、特に限定されず、図５のような同時再構成アルゴリズム（MLACF法やMLAA法）あるいは、図６のような同時再構成アルゴリズム（MLACF法やMLAA法）とは異なる再構成アルゴリズム（例えばML-EM法）を用いてもよい。
以上のように、本発明は、TOF計測型ＰＥＴ装置で計測された計測データを用いた吸収係数画像の推定に適している。
１ … ＰＥＴ装置
３ … γ線検出器
５ … 演算回路
６ … 吸収係数画像推定プログラム
λ'，λ2' … 放射能画像
A' … 吸収係数サイノグラム
μ' … （定量的な吸収係数画像に対して不均一なオフセット値が加算された）画像
mimg … 被検体マスク画像
mproj … 被検体マスク投影データ
μoff … （不均一な）オフセット画像
Ω … （既知の吸収係数値で近似可能な）領域
K … 既知の吸収係数値で近似可能な領域数
Ωn … n番目の領域Ω
μn known … 領域Ω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(μn known,μ'+α×μoff) … 領域Ωn内に既知の吸収係数を設定した画像μn knownと、定量的な吸収係数画像に対して不均一なオフセット値が加算された画像(μ'+α×μoff)との誤差評価関数
wn … （０以上１以下の）係数
f(α) … （係数αを変数とした）関数

Claims (11)

  1. 消滅放射線の飛行時間差(Time Of Flight)情報を含んだポジトロンＣＴの計測データから吸収係数画像を推定する方法であって、
    μ’を定量的な吸収係数画像に対して不均一なオフセット値が加算された画像とし、前記計測データに関する評価関数の最適化に基づいて、前記画像μ’を計算する再構成計算工程と、
    前記計測データに基づいて投影データ空間における被検体マスクデータである被検体マスク投影データを算出するマスク算出工程と、
    μｏｆｆを不均一なオフセット画像としたときに、オフセット画像μｏｆｆの順投影データが被検体マスク投影データを近似するように構成された再構成アルゴリズムによって、前記オフセット画像μｏｆｆを推定するオフセット推定工程と、
    Ωを既知の吸収係数値で近似可能な領域としたときに、前記計測データに基づいて計算された、被検体領域が認識可能な画像を用いて、前記領域Ωを少なくとも１つ以上抽出する参照領域抽出工程と、
    αを係数としたときに、前記領域Ωにおける前記画像μ’の値と既知の吸収係数値との誤差を減少させる前記係数αを算出する係数算出工程と、
    前記画像μ’の値に、前記オフセット画像μｏｆｆを前記係数α倍したα×μｏｆｆを加算して得られた値を吸収係数値として補正する吸収係数値補正工程と
    を備える、
    吸収係数画像推定方法。
  2. 請求項１に記載の吸収係数画像推定方法において、
    前記再構成計算工程を、(a)前記画像μ’を未知数に含む計算アルゴリズムで実施する、または(b)吸収係数投影データを未知数に含む計算アルゴリズムおよび前記吸収係数投影データを再構成した画像を前記画像μ’とするアルゴリズムの組み合わせで実施する、
    吸収係数画像推定方法。
  3. 請求項１または請求項２に記載の吸収係数画像推定方法において、
    (A)前記マスク算出工程は、
    前記画像μ’の二値化画像を被検体マスク画像として算出する工程と、
    前記被検体マスク画像の投影データを算出する工程と、
    前記被検体マスク画像の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
    からなる、
    または
    (B)前記マスク算出工程は、
    前記画像μ’の投影データを算出する工程と、
    前記画像μ’の投影データの二値化データを前記被検体マスク投影データとして算出する工程と
    からなる、
    吸収係数画像推定方法。
  4. 請求項１または請求項２に記載の吸収係数画像推定方法において、
    前記マスク算出工程は、
    前記計測データを投影データ形式に変換したものを二値化したデータを前記被検体マスク投影データとして算出する工程からなる、
    吸収係数画像推定方法。
  5. 請求項１または請求項２に記載の吸収係数画像推定方法において、
    (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)である、
    吸収係数画像推定方法。
  9. 請求項１から請求項７のいずれかに記載の吸収係数画像推定方法において、
    K(≧1)を既知の吸収係数値で近似可能な領域数とし、Ωｎをｎ番目の前記領域Ωとし、μｎ ｋｎｏｗｎ(n=1,…,K)を前記領域Ωｎ内に既知の吸収係数を設定した画像とし、DΩｎ(X,Y)を前記領域Ωｎ内に関する画像Xおよび画像Yの誤差評価関数とし、ｗｎ(n=1,…,K)を０以上１以下の係数としたときに、前記係数算出工程における前記係数αは、関数f(α)= Σn=1,…,K[ｗｎ×DΩｎ(μｎ ｋｎｏｗｎ，μ’+α×μｏｆｆ)]を最小化するαである、
    吸収係数画像推定方法。
  10. 請求項１から請求項９のいずれかに記載の吸収係数画像推定方法をコンピュータに実行させる、吸収係数画像推定プログラム。
  11. 請求項１０に記載の吸収係数画像推定プログラムを搭載したポジトロンＣＴ装置において、
    当該吸収係数画像推定プログラムを実行する演算手段を備える、ポジトロンＣＴ装置。
JP2019521549A 2017-05-29 2017-05-29 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置 Active JP6761610B2 (ja)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/019944 WO2018220686A1 (ja) 2017-05-29 2017-05-29 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置

Publications (2)

Publication Number Publication Date
JPWO2018220686A1 JPWO2018220686A1 (ja) 2019-11-07
JP6761610B2 true JP6761610B2 (ja) 2020-09-30

Family

ID=64455240

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2019521549A Active JP6761610B2 (ja) 2017-05-29 2017-05-29 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置

Country Status (5)

Country Link
US (1) US20200085397A1 (ja)
EP (1) EP3633411A4 (ja)
JP (1) JP6761610B2 (ja)
CN (1) CN110691992A (ja)
WO (1) WO2018220686A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108986182B (zh) * 2018-07-10 2022-11-25 上海联影医疗科技股份有限公司 一种重建ct图像的方法、系统及存储介质
JP7247782B2 (ja) * 2019-06-24 2023-03-29 株式会社島津製作所 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置

Family Cites Families (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4604970B2 (ja) * 2005-11-10 2011-01-05 株式会社島津製作所 核医学診断装置
JP2009047602A (ja) * 2007-08-21 2009-03-05 Toshiba Corp 陽電子放出コンピュータ断層撮影装置、減弱マップ作成装置および減弱マップ作成プログラム
US8299438B2 (en) * 2009-07-16 2012-10-30 Siemens Medical Solutions Usa, Inc. Model based estimation of a complete or partial positron emission tomography attenuation map using maximum likelihood expectation maximization
US9342903B2 (en) * 2012-03-28 2016-05-17 National Institute Of Radiological Services Method for generating image for PET attenuation correction from MR image and computer program
JP6235564B2 (ja) * 2012-05-04 2017-11-22 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 陽電子放出断層撮影における散乱された同時発生を用いた減衰マップ
EP2662022A1 (en) * 2012-05-07 2013-11-13 National Cerebral and Cardiovascular Center Method of extracting contour of tomogram
WO2014066629A1 (en) * 2012-10-26 2014-05-01 The Regents Of The University Of California Image-based point-spread-function modelling in time-of-flight positron-emission-tomography iterative list-mode reconstruction
US9928617B2 (en) * 2012-11-08 2018-03-27 The General Hospital Corporation System and method for multi-modality time-of-flight attenuation correction
US9507033B2 (en) * 2013-02-05 2016-11-29 Siemens Medical Solutions Usa, Inc. Method and apparatus for compensating for scattering of emission gamma photons for PET imaging
US20150057535A1 (en) * 2013-05-08 2015-02-26 Arkadiusz Sitek Systems and methods for attenuation correction in time-of-flight positron emission tomography
US9155514B2 (en) * 2013-08-01 2015-10-13 Siemens Medical Solutions Usa, Inc. Reconstruction with partially known attenuation information in time of flight positron emission tomography
US9750475B2 (en) * 2013-11-05 2017-09-05 Shimadzu Corporation Contour image generating device and nuclear medicine diagnosis apparatus
US9600910B2 (en) * 2014-01-08 2017-03-21 Rensselaer Polytechnic Institute Attenuation map reconstruction from TOF PET data
JP2015145828A (ja) * 2014-02-03 2015-08-13 株式会社東芝 画像処理装置及び画像処理プログラム
US10537299B2 (en) * 2015-06-04 2020-01-21 Rensselaer Polytechnic Institute Attenuation map reconstruction from TOF PET data
US10311604B2 (en) * 2016-07-08 2019-06-04 Shanghai United Imaging Healthcare Co., Ltd. System and method for generating attenuation map
CN110446946A (zh) * 2017-03-09 2019-11-12 株式会社岛津制作所 散射估计方法、散射估计程序以及搭载有该散射估计程序的正电子ct装置

Also Published As

Publication number Publication date
US20200085397A1 (en) 2020-03-19
EP3633411A4 (en) 2021-03-03
CN110691992A (zh) 2020-01-14
JPWO2018220686A1 (ja) 2019-11-07
EP3633411A1 (en) 2020-04-08
WO2018220686A1 (ja) 2018-12-06

Similar Documents

Publication Publication Date Title
US7718954B2 (en) Component method and system for PET detector efficiency normalization
EP1934942B1 (en) Iterative reconstruction with enhanced noise control filtering
Kacperski et al. Iterative deconvolution of simultaneous 99mTc and 201Tl projection data measured on a CdZnTe-based cardiac SPECT scanner
JP4208284B2 (ja) 核像形成方法及び装置
JP7118048B2 (ja) 局所的に修正された飛行時間（ｔｏｆ）カーネルを使用するｔｏｆ ｐｅｔ画像再構成
US10215864B2 (en) System and method to improve image quality of emission tomography when using advanced radionuclides
JP6711450B2 (ja) 散乱推定方法、散乱推定プログラム並びにそれを搭載したポジトロンｃｔ装置
US10036817B2 (en) Solving outside-field of view scatter correction problem in positron emission tomography via digital experimentation
US8712124B2 (en) Artifact removal in nuclear images
JP7317586B2 (ja) 医用画像処理装置、方法及びプログラム
JP6761610B2 (ja) 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置
JP2021173755A (ja) 医用画像処理装置、医用画像処理方法及びプログラム
Bouwens et al. Image-correction techniques in SPECT
JP7302737B2 (ja) リストモード画像再構成方法および核医学診断装置
JP7247782B2 (ja) 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置
JPWO2020059135A1 (ja) データ処理方法、プログラム、データ処理装置および陽電子放出断層撮像装置
US20240037814A1 (en) Convolutional neural network for dynamic pet frame clustering
US20230206516A1 (en) Scatter estimation for pet from image-based convolutional neural network
EP4318400A1 (en) Image processing apparatus, image processing method and program
US20230260171A1 (en) Event property-dependent point spread function modeling and image reconstruction for pet
Ljungberg Instrumentation, Calibration, Quantitative Imaging, and Quality Control
CN113253330A (zh) 伽玛射线放射成像装置及能量校准方法
Wallstén Correction for partial volume effects in PET imaging
Lam et al. Correction of inter-crystal scatter effect in iterative image reconstruction of the jPET-D4
Ben-Haim Iterative deconvolution of simultaneous Tc-99m and Tl-201 projection data measured on a CdZnTe-based cardiac SPECT...

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20190618

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20200630

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20200727

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20200807

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20200820

R151 Written notification of patent or utility model registration

Ref document number: 6761610

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151