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

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

Info

Publication number
WO2018220686A1
WO2018220686A1 PCT/JP2017/019944 JP2017019944W WO2018220686A1 WO 2018220686 A1 WO2018220686 A1 WO 2018220686A1 JP 2017019944 W JP2017019944 W JP 2017019944W WO 2018220686 A1 WO2018220686 A1 WO 2018220686A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
absorption coefficient
value
region
projection data
Prior art date
Application number
PCT/JP2017/019944
Other languages
English (en)
French (fr)
Inventor
哲哉 小林
Original Assignee
株式会社島津製作所
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 株式会社島津製作所 filed Critical 株式会社島津製作所
Priority to CN201780091422.6A priority Critical patent/CN110691992A/zh
Priority to EP17912063.9A priority patent/EP3633411A4/en
Priority to JP2019521549A priority patent/JP6761610B2/ja
Priority to US16/499,466 priority patent/US20200085397A1/en
Priority to PCT/JP2017/019944 priority patent/WO2018220686A1/ja
Publication of WO2018220686A1 publication Critical patent/WO2018220686A1/ja

Links

Images

Classifications

    • 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
    • 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
    • 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]

Definitions

  • the present invention relates to an absorption coefficient image estimation method, an absorption coefficient image estimation program for estimating an absorption coefficient image from measurement data of a positron CT apparatus (positron emission tomography apparatus), and a positron CT apparatus equipped with the absorption coefficient image estimation program.
  • a positron CT device or PET (Positron Emission Tomography) device, is effective only when two ⁇ -rays generated by the annihilation of positron (Positron) are detected simultaneously by multiple detectors (that is, only when they are counted simultaneously). It is configured to measure the signal as a simple signal and reconstruct a tomographic image of the subject based on the measurement data.
  • a radiopharmaceutical containing a positron emitting nuclide is administered to a subject, and 511 keV paired annihilation gamma rays released from the administered subject are detected by a group of detector elements (for example, scintillators). Detect with instrument.
  • positron CT In positron CT (PET), various data correction processes are required to quantitatively measure the radioactivity concentration in the subject. Typical correction processes include sensitivity correction, scattering correction, random correction, attenuation correction, dead time correction, and absorption correction.
  • the present invention relates to an absorption correction for preventing deterioration in image quantitativeness due to absorption of ⁇ rays emitted from a radiopharmaceutical (radioisotope). In order to perform the absorption correction, it is necessary to estimate an absorption coefficient image obtained by imaging the absorption coefficient distribution in the subject.
  • the transmittance of ⁇ rays is obtained from the estimated absorption coefficient image, and the transmittance is divided from the measurement data of PET, and converted to data from which the influence of absorption of ⁇ rays has been eliminated.
  • the estimated absorption coefficient image is incorporated into an image reconstruction calculation formula to obtain a reconstructed image from which the influence of ⁇ -ray absorption is eliminated.
  • transmission data obtained by irradiating an external source of positron emitting nuclides is necessary.
  • CT data obtained from an X-ray CT (Computed Tomography) apparatus instead of transmission data.
  • Non-Patent Documents 1 and 2 there are reconstruction algorithms that do not require such transmission data (for example, see Non-Patent Documents 1 and 2).
  • TOF-PET the measurement data of PET
  • TOF time difference
  • the distribution shape of the absorption coefficient sinogram can be estimated.
  • the radioactivity image and the absorption coefficient sinogram can be estimated simultaneously.
  • the reconstruction algorithms in Non-Patent Documents 1 and 2 are also called “simultaneous reconstruction algorithms” because they simultaneously estimate the radioactivity image and the data related to the absorption coefficient (for example, absorption coefficient sinogram).
  • the simultaneous reconstruction algorithm that simultaneously estimates the radioactivity image and the absorption coefficient sinogram is also called the MLACF method.
  • the radioactivity concentration is distributed as shown in FIG. 9 in the target region of the subject.
  • the radioactivity distribution is developed into a two-dimensional distribution consisting of TOF information t and radial direction s
  • the projection direction ⁇ is in the range of 0 ° to 180 °, a two-dimensional distribution for each projection angle ⁇ is obtained as measurement data.
  • the following technologies are indispensable. That is, a technique for estimating an absorption coefficient image with a correct absolute value from an absorption coefficient sinogram with a correct distribution shape but an incorrect absolute value is indispensable.
  • Non-Patent Document 2 there is a method of directly quantifying a radioactivity image without directly quantifying an absorption coefficient sinogram.
  • the processing steps are as follows. (1) A radioactivity image without absorption correction is created, and an object mask image is created from the radioactivity image. (2) Substituting a known absorption coefficient (for example, water absorption coefficient or bone absorption coefficient) into each pixel value of the subject mask image, and abbreviated as a pseudo absorption coefficient image (hereinafter referred to as “pseudo absorption coefficient image”). Create). (3) Using the pseudo absorption coefficient image, a radioactivity image is created by a conventional reconstruction algorithm, and the total pixel value is calculated for each slice in the body axis direction. As shown in FIG.
  • z is the body axis direction
  • the xy plane is a plane (axial plane) orthogonal to the body axis direction z
  • P (x, y) is the pixel value
  • S (z) is the body axis direction z.
  • the total pixel value is calculated for each slice in the body axis direction in the same manner as in step (3).
  • S ′ (z) is the total pixel value for each z-slice in the body axis direction.
  • Step (4) so that the total pixel value S ′ (z) obtained in step (5) matches the total pixel value S (z) obtained in step (3) for each slice in the body axis direction.
  • Non-Patent Document 2 has a problem that it cannot be guaranteed in principle that the finally obtained radioactivity image is quantitative.
  • the present invention has been made in view of such circumstances, and an absorption coefficient image estimation method, an absorption coefficient image estimation program capable of creating a quantitative absorption coefficient image, and a positron CT apparatus equipped with the absorption coefficient image estimation program are provided.
  • the purpose is to provide.
  • the inventors have obtained the following knowledge. That is, the above-described conventional technique assumes that the total pixel value for each slice in the body axis direction obtained in the procedure (3) is true (that is, correct).
  • the absorption coefficient image estimation method is a method for estimating an absorption coefficient image from positron CT measurement data including time-of-flight information of annihilation radiation, and quantitatively determines ⁇ ′.
  • a reconstruction calculation step for calculating the image ⁇ ′ based on optimization of an evaluation function related to the measurement data, and an image obtained by adding an uneven offset value to the absorption coefficient image
  • a mask calculation step for calculating subject mask projection data, which is subject mask data in the projection data space, and forward projection data of the offset image ⁇ off is subject mask projection data when ⁇ off is a non-uniform offset image.
  • the configured reconstruction algorithm to approximate the offset estimation step of estimating the offset image mu off, the ⁇ with known absorption coefficients
  • an absorption coefficient value correcting step of correcting, as an absorption coefficient value, a value obtained by adding ⁇ ⁇ ⁇ off obtained by multiplying off by the coefficient ⁇ .
  • the evaluation function of the positron CT measurement data (TOF-PET measurement data) including the time-of-flight difference (Time Of Flight) information of annihilation radiation is calculated. Based on the optimization, an image in which a non-uniform offset value is added to the quantitative absorption coefficient image is calculated.
  • ⁇ ′ is a non-quantitative radioactivity image
  • the radioactivity image ⁇ ′ is calculated in the reconstruction calculation step.
  • the reconstruction algorithm in the reconstruction calculation step here is the simultaneous reconstruction algorithm described above.
  • subject mask projection data in the projection data space (hereinafter referred to as “subject mask projection data”) is calculated based on the measurement data.
  • subject mask projection data in the projection data space
  • the offset estimation step the reconstruction algorithm forward projection data of the offset image mu off is configured to approximate the object mask projection data, an offset image mu off presume.
  • the reference region extraction step uses an image that is calculated based on the measurement data and can recognize the subject region. Extract one or more.
  • the image calculated based on the measurement data and capable of recognizing the subject region is, for example, the above-described image ⁇ ′, the above-described radioactivity image ⁇ ′, and the reconstruction in the above-described reconstruction calculation step.
  • These are non-quantitative images estimated by a reconstruction algorithm different from the algorithm (simultaneous reconstruction algorithm), and object mask images estimated from the above-described object mask projection data.
  • the image is not limited to these exemplified images, and may be any image that can be recognized based on the measurement data and that can recognize the subject region.
  • extracting at least one region ⁇ using an image that can be recognized by the subject region calculated based on the measurement data means that the region ⁇ is extracted using single image information. In addition to this, extracting the region ⁇ using a combination of a plurality of pieces of image information (for example, logical sum or logical product) is also included.
  • the following coefficient calculation step and absorption coefficient value correction step are performed using the image ⁇ ′, the offset image ⁇ off and the region ⁇ obtained in the above-described steps.
  • ⁇ is an unknown true absorption coefficient image value (absorption coefficient value)
  • ⁇ is a coefficient
  • the coefficient calculation process reduces the error between the image ⁇ ′ value and the known absorption coefficient value in the region ⁇ .
  • the coefficient ⁇ to be calculated is calculated.
  • the absorption coefficient value correcting step a value obtained by adding ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by a coefficient ⁇ to the value of the image ⁇ ′ is corrected as an absorption coefficient value.
  • the difference between the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (the value of the true absorption coefficient image) is approximated by ⁇ ⁇ ⁇ off, which is obtained by multiplying the offset image ⁇ off by a coefficient ⁇ .
  • the absorption coefficient value is corrected in the absorption coefficient value correction step based on a mathematical relationship that is possible (that is, ⁇ ′ + ⁇ ⁇ ⁇ off ). Therefore, the absorption coefficient image having the absorption coefficient value corrected in the absorption coefficient value correction step has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • the above-described reconstruction calculation process may be performed with a calculation algorithm including (a) the image ⁇ ′ as an unknown.
  • the above-described reconstruction calculation step may be performed by a combination of (b) a calculation algorithm including absorption coefficient projection data as an unknown and an algorithm in which an image obtained by reconstructing absorption coefficient projection data is an image ⁇ ′.
  • the former algorithm (a) is an MLAA method in which a non-quantitative radioactivity image ⁇ ′ and an image ⁇ ′ (non-quantitative absorption coefficient image) are simultaneously reconstructed.
  • the latter algorithm (b) uses the MLACF method described in Non-Patent Document 1 for simultaneous estimation of non-quantitative radioactivity image ⁇ ′ and non-quantitative absorption coefficient projection data (for example, absorption coefficient sinogram), and absorption.
  • This is a combination with an algorithm in which an image obtained by reconstructing coefficient projection data is an image ⁇ ′.
  • the algorithm in which the image ⁇ ′ is an image obtained by reconstructing the absorption coefficient projection data in the latter (b) is not particularly limited as long as it is a reconstruction algorithm.
  • An example of the mask calculating step described above includes a step of calculating a binarized image of the image ⁇ ′ as a subject mask image, a step of calculating projection data of the subject mask image, and a projection data of the subject mask image. And a step of calculating the value data as subject mask projection data (aspect (A)). Further, the mask calculation step may be a mode of (B) including a step of calculating projection data of the image ⁇ ′ and a step of calculating binarized data of the projection data of the image ⁇ ′ as subject mask projection data. .
  • the projection data is calculated after binarizing the image ⁇ ′ first, and the projection data is binarized.
  • the projection data of the image ⁇ ′ is binarized after being calculated first.
  • the image ⁇ ′ is a non-quantitative absorption coefficient image, but the object mask projection data can be calculated by binarizing the image other than the absorption coefficient image and the projection data.
  • Another example of the mask calculation step includes a step of calculating data obtained by binarizing the above-described TOF-PET measurement data converted into a projection data format as subject mask projection data.
  • the TOF-PET measurement data is directly used, and binarized data converted into the projection data format can be calculated as the subject mask projection data.
  • the mask calculation step is a step of calculating a radioactivity image based on optimization of an evaluation function related to the measurement data of TOF-PET (described above) and a projection data of the radioactivity image. And a step of calculating data obtained by binarizing projection data (of a radioactivity image) as subject mask projection data (aspect (C)).
  • the mask calculation step is a step of calculating a radioactivity image based on optimization of the evaluation function related to the measurement data of TOF-PET (described above), a step of calculating a binarized image of the radioactivity image,
  • the aspect of (D) may include a step of calculating projection data of a binarized image and a step of calculating data obtained by binarizing the projection data of the binarized image as subject mask projection data.
  • the projection data of the radioactivity image is binarized after previously calculated.
  • the projection data is calculated after binarizing the radioactivity image first, and the projection data is binarized.
  • the reconstruction algorithm for calculating the radioactivity image is not particularly limited.
  • the subject mask projection data can be calculated using the above-described radioactivity image ⁇ ′ estimated by the simultaneous reconstruction algorithm.
  • the subject mask projection data may be calculated using a radioactivity image estimated by a reconstruction algorithm (for example, ML-EM method) different from the above-described simultaneous reconstruction algorithm (MLACF method or MLAA method). .
  • the reconstruction processing performed in the offset estimation step described above may be performed by any one of analytical reconstruction, statistical reconstruction, and algebraic reconstruction.
  • analytical reconstruction for example, there is an FBP (Filtered Back Projection) method.
  • statistical reconstruction for example, the above-described ML-EM (Maximum Likelihood- ⁇ Expectation Maximization) method is available.
  • algebraic reconstruction for example, there is an ART (Algebraic Reconstruction Technique) method.
  • the at least one region ⁇ extracted in the above-described reference region extraction step is a tissue region in which the absorption coefficient can be regarded as known.
  • the “region of tissue whose absorption coefficient can be regarded as known” means, for example, a region that can be approximated by water, a region that can be approximated by air, a region that can be approximated by brain tissue, a region that can be approximated by bone, and lung tissue The region that can be approximated by, and the region that can be approximated by soft tissue.
  • any structure can be used as long as an approximate value of the absorption coefficient can be understood.
  • the water absorption coefficient value of 0.0096 / mm can be used.
  • K ( ⁇ 1) is the number of regions that can be approximated by a known absorption coefficient value
  • ⁇ n is the n-th region ⁇
  • n is a known absorption coefficient value
  • S (X; ⁇ n ) is a statistic of image X in region ⁇ n or a value calculated from the statistic is a representative value
  • the representative value is, for example, an average value, a median value, a trimmed average value, a trimmed median value, or a weighted average value of two or more values thereof.
  • the value is not limited to these exemplified values, and may be “a statistic or a value calculated from the statistic”.
  • S ( ⁇ ′; ⁇ n ) is a representative value of image ⁇ ′ in region ⁇ n
  • S ( ⁇ off ; ⁇ n ) is a representative value of the offset image ⁇ off ′ in the region ⁇ n .
  • n 1, ..., each region K Omega 1, ..., coefficient alpha 1 in Omega K, ..., an alpha K after calculating each coefficient ⁇ 1, ..., ⁇ K representative value T
  • the coefficient ⁇ can be calculated by setting ⁇ 1 , ⁇ 2 ,..., ⁇ K ) as the coefficient ⁇ .
  • K ( ⁇ 1) is the number of regions that can be approximated with a known absorption coefficient value
  • ⁇ n is the nth region ⁇
  • image set the known absorption coefficient in n, the error evaluation function D ⁇ n (X, Y) image X and image Y about the region ⁇ n, w n (n 1, ..., K) 0
  • the coefficient is 1 or less.
  • Image X is obtained by replacing image X n known by setting a known absorption coefficient in region ⁇ n and adding image Y to image ⁇ ′ by multiplying offset image ⁇ off ′ by ⁇ ⁇ ⁇ off.
  • D ⁇ n ( ⁇ n known, ⁇ value ⁇ )' + ⁇ ⁇ ⁇ off ) includes an image mu n known set of known absorption coefficient for region Omega n, This is an error evaluation function with an image ( ⁇ ′ + ⁇ ⁇ ⁇ off ) obtained by adding a non-uniform offset value to a quantitative absorption coefficient image.
  • the coefficient ⁇ can be calculated by calculating ⁇ that minimizes known , ⁇ ′ + ⁇ ⁇ ⁇ off )]).
  • the absorption coefficient image estimation program according to the present invention causes a computer to execute the absorption coefficient image estimation method according to the present invention.
  • the absorption coefficient image estimation program by causing a computer to execute the absorption coefficient image estimation method according to the present invention, a non-quantitative image ⁇ ′ value in a region ⁇ and a known absorption coefficient value (true The absorption coefficient value is corrected in the absorption coefficient value correcting step based on a mathematical relationship that the difference from the absorption coefficient image value) can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by a coefficient ⁇ . . Therefore, the absorption coefficient image having the absorption coefficient value corrected in the absorption coefficient value correction step has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • the positron CT apparatus is a positron CT apparatus equipped with the absorption coefficient image estimation program according to the present invention, and includes a calculation means for executing the absorption coefficient image estimation program.
  • the positron CT apparatus by including a calculation unit that executes the absorption coefficient image estimation program according to the present invention, the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (true The absorption coefficient value is corrected in the absorption coefficient value correcting step based on a mathematical relationship that the difference from the absorption coefficient image value) can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by a coefficient ⁇ . . Therefore, the absorption coefficient image having the absorption coefficient value corrected in the absorption coefficient value correction step has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (true
  • the absorption coefficient value is corrected in the absorption coefficient value correction step based on the mathematical relationship that the offset image ⁇ off can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by the coefficient ⁇ . Therefore, the absorption coefficient image having the absorption coefficient value corrected in the absorption coefficient value correction step has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • FIG. 3 is a flowchart illustrating a processing procedure and a data flow of an absorption coefficient image estimation method according to the first embodiment.
  • 10 is a flowchart illustrating a processing procedure and a data flow of an absorption coefficient image estimation method according to the second embodiment.
  • 10 is a flowchart showing a processing procedure and a data flow of an absorption coefficient image estimation method according to Example 3 when subject mask projection data is calculated using a radioactivity image estimated by the MLACF method.
  • FIG. 14 is a flowchart illustrating a processing procedure and a data flow of an absorption coefficient image estimation method according to a fourth embodiment.
  • 10 is a flowchart showing a processing procedure and data flow of an absorption coefficient image estimation method according to Embodiment 5. It is a conceptual diagram with which the principle in the MLACF method is provided.
  • FIG. 6 is a schematic diagram of a profile when a total pixel value for each slice in the body axis direction, a horizontal axis is a body axis direction, and a vertical axis is a total pixel value.
  • FIG. 1 is a schematic perspective view and a block diagram of a PET apparatus according to each embodiment
  • FIG. 2 is a schematic perspective view of a ⁇ -ray detector. 1 and 2 have the same configuration in each embodiment.
  • the PET apparatus 1 includes a detector ring 2 that surrounds the periphery of the subject in a stacked arrangement in the body axis direction of the subject.
  • a plurality of ⁇ -ray detectors 3 are embedded in the detector ring 2.
  • the PET apparatus 1 corresponds to the positron CT apparatus in the present invention.
  • the ⁇ -ray detector 3 corresponds to the detector in the present invention.
  • the PET apparatus 1 includes a coincidence counting circuit 4 and an arithmetic circuit 5.
  • a coincidence counting circuit 4 includes a coincidence counting circuit 4 and an arithmetic circuit 5.
  • FIG. 1 only two connections from the ⁇ -ray detector 3 to the coincidence counting circuit 4 are shown, but actually, a photomultiplier tube (PMT: Photo Multiplier Tube) 33 of the ⁇ -ray detector 3 ( Are connected to the coincidence counting circuit 4 by the total number of channels (see FIG. 2).
  • the arithmetic circuit 5 executes processing of an absorption coefficient image estimation method shown in FIG. 3 described later by the absorption coefficient image estimation program 6.
  • the arithmetic circuit 5 corresponds to the arithmetic means in the present invention.
  • a ⁇ -ray generated from a subject (not shown) to which a radiopharmaceutical has been administered is converted into light by a scintillator block 31 (see FIG. 2) of the ⁇ -ray detector 3, and the converted light is converted into a ⁇ -ray detector.
  • 3 photomultiplier tube (PMT) 33 (see FIG. 2) is multiplied and converted into an electrical signal. The electrical signal is sent to the coincidence counting circuit 4 to generate count value detection signal data.
  • the coincidence circuit 4 checks the position of the scintillator block 31 (see FIG. 2) and the incident timing of the ⁇ -ray, and sends it only when ⁇ -rays are simultaneously incident on the two scintillator blocks 31 on both sides of the subject. The determined electrical signal is determined as appropriate data.
  • the coincidence counting circuit 4 rejects. That is, the coincidence counting circuit 4 detects that ⁇ rays are simultaneously observed (that is, coincidence counting) by the two ⁇ ray detectors 3 based on the above-described electrical signal.
  • the detection signal data (count value) composed of appropriate data determined to be coincidence by the coincidence circuit 4 is sent to the arithmetic circuit 5.
  • the arithmetic circuit 5 performs steps S1 to S8 (see FIG. 3) to be described later, and estimates an absorption coefficient image from detection signal data of a subject (not shown) obtained by the PET apparatus 1. Specific functions of the arithmetic circuit 5 will be described later.
  • An absorption coefficient image estimation program 6 is stored in a storage medium (not shown) represented by ROM (Read-only Memory), and the absorption coefficient image estimation program 6 is read from the storage medium to the arithmetic circuit 5. Then, the calculation circuit 5 executes the absorption coefficient image estimation program 6 to perform the process of the absorption coefficient image estimation method shown in the flowchart of FIG.
  • the arithmetic circuit 5 is a GPU (Graphics Processing Unit), a central processing unit (CPU), or a programmable device (for example, an FPGA (Field Programmable Gate) that can change a hardware circuit (for example, a logic circuit) used in accordance with program data. Array)).
  • the ⁇ -ray detector 3 includes a scintillator block 31, a light guide 32 optically coupled to the scintillator block 31, and photoelectrons optically coupled to the light guide 32.
  • a multiplier (hereinafter simply abbreviated as “PMT”) 33 is provided.
  • Each scintillator element constituting the scintillator block 31 converts ⁇ rays into light by emitting light with the incidence of ⁇ rays. By this conversion, the scintillator element detects ⁇ rays.
  • Light emitted from the scintillator element is sufficiently diffused by the scintillator block 31 and input to the PMT 33 via the light guide 32.
  • the PMT 33 multiplies the light converted by the scintillator block 31 and converts it into an electric signal.
  • the electric signal is sent to the coincidence counting circuit 4 (see FIG. 1) as a pixel value.
  • the ⁇ -ray detector 3 is a DOI detector composed of scintillator elements arranged three-dimensionally and composed of a plurality of layers in the depth direction.
  • a four-layer DOI detector is illustrated, but the number of layers is not particularly limited as long as it is plural.
  • the DOI detector is constructed by laminating the respective scintillator elements in the radiation depth direction, and the interaction depth (DOI: Depth of Interaction) direction and lateral direction (incident surface). Coordinate information with a direction parallel to the center of gravity).
  • DOI Depth of Interaction
  • lateral direction incident surface
  • FIG. 3 is a flowchart illustrating the processing procedure and the data flow of the absorption coefficient image estimation method according to the first embodiment.
  • the subject is imaged by the PET apparatus 1 shown in FIG. 1, and list mode data is acquired by the coincidence counting circuit 4 (see FIG. 1).
  • the list mode data energy information of detected photons is recorded.
  • ⁇ a normal energy window (for example, 400 keV-600 keV), that is, an energy window for reconstruction data, and a measurement range and bin width in the TOF direction of TOF measurement data.
  • bin means binarization (separation).
  • a pixel corresponds to a bin.
  • the TOF bin means a temporal separation of TOF information. For example, when the TOF bin is 100 [ps], the detection time difference is temporally separated with an accuracy of 100 [ps].
  • MLACF ⁇ ′ is an image obtained by adding a non-uniform offset value to a quantitative absorption coefficient image
  • ⁇ ′ is a non-quantitative radioactivity image. Based on the optimization of the evaluation function regarding the measurement data, the image ⁇ ′ and the radioactivity image ⁇ ′ are calculated simultaneously. As described in the section “Means for Solving the Problems”, the image ⁇ ′ is also an absorption coefficient image. However, it is simply referred to as “image ⁇ ′” to distinguish it from the quantitative absorption coefficient image finally obtained. .
  • the reconstruction algorithm in Step S1 is a simultaneous reconstruction algorithm in Non-Patent Document 1.
  • the MLACF method described in Non-Patent Document 1 is applied as the simultaneous reconstruction algorithm.
  • the specific method of the MLACF method see Non-Patent Document 1 described above.
  • a ′ is the absorption coefficient sinogram. Radioactivity image ⁇ ′ and absorption coefficient sinogram A ′ are estimated by MLACF method. The absorption coefficient sinogram A ′ corresponds to the absorption coefficient projection data in the present invention.
  • Step S2 ML-TR or ML-EM
  • the absorption coefficient sinogram A ′ is reconstructed by a reconstruction algorithm (for example, ML-TR method or ML-EM method) to create a non-quantitative image ⁇ ′.
  • a reconstruction algorithm for example, ML-TR method or ML-EM method
  • log conversion of the absorption coefficient sinogram A ′ is performed in advance.
  • Reference 1 Reference 1
  • Reference 2 Reference 2
  • Reference 2 Reference 2
  • Step 2 Reference 2: LA Shepp and Y. Vardi. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging, Vol. 1 , pp. 113-122, 1982.
  • Steps S1 and S2 correspond to a reconstruction calculation step in the present invention.
  • Step S3 Binarization Processing
  • the image ⁇ ′ is binarized by threshold processing. Then, a binarized image in which the subject region is “1” and the other regions are “0” is calculated as a subject mask image. Let m img be the subject mask image.
  • Step S4 Projection + Binarization Processing Line integral data (projection data) of the subject mask image m img is calculated (Projection). Then, the projection data of the subject mask image m img is binarized by threshold processing, so that the projection line passing through the subject becomes “1” and the other projection lines become “0”. The digitized data is calculated as subject mask projection data. Let m proj be subject mask projection data. Steps S3 and S4 correspond to a mask calculation step in the present invention.
  • Step S5 ML-EM Let ⁇ off be a non-uniform offset image.
  • the configured reconstruction algorithm as the forward projection data of the offset image mu off approximates the subject mask projection data m proj, estimates the offset image mu off. That is, the subject mask projection data m proj is converted into image data by a reconstruction algorithm (for example, ML-EM method). This image data is set as an offset image ⁇ off .
  • Step S5 corresponds to the offset estimation step in the present invention.
  • Step S6 Extraction Let ⁇ be a region that can be approximated by a known absorption coefficient value. At least one region ⁇ is extracted (Extraction) using an image recognizable from the subject region calculated based on the measurement data.
  • K ( ⁇ 1) is the number of regions that can be approximated by a known absorption coefficient value
  • ⁇ n is the n-th region ⁇
  • each region ⁇ 1 ,..., ⁇ K is extracted from a region that can be approximated by air, a region that can be approximated by brain tissue, and a region that can be approximated by bone.
  • Examples of the image that can be recognized based on the measurement data that can recognize the subject region include the image ⁇ ′ created in step S2, the radioactivity image ⁇ ′ obtained in step S1, and the reconstruction algorithm in step S1 ( A non-quantitative image estimated by a different reconstruction algorithm (for example, ML-EM method), and a subject mask image m img estimated from the above-described subject mask projection data m proj .
  • Step S6 corresponds to the reference region extraction step in the present invention.
  • ⁇ is a coefficient.
  • the representative value is, for example, an average value, a median value, a trimmed average value, a trimmed median value, or a weighted average value of two or more values thereof.
  • the “trimmed average value” means an average value of the remaining data from which data having extremely large / small values is removed.
  • the “trimmed median” means the median of the remaining data obtained by removing data with extremely large / small values.
  • Step S7 corresponds to a coefficient calculation step in the present invention.
  • Step S8 corresponds to an absorption coefficient value correcting step in the present invention.
  • the absorption correction is performed using the absorption coefficient image having the absorption coefficient value ⁇ obtained in step S8.
  • the transmittance of ⁇ rays is obtained from the estimated absorption coefficient image, and the transmittance is divided from the PET measurement data, thereby eliminating the influence of the absorption of ⁇ rays.
  • Absorption correction is performed by conversion.
  • the absorption correction is performed by incorporating the estimated absorption coefficient image into a calculation formula for image reconstruction to obtain a reconstructed image from which the influence of ⁇ -ray absorption is eliminated.
  • step S1 the evaluation of the positron CT measurement data (TOF-PET measurement data) including the time-of-flight difference (Time ⁇ Of ⁇ Flight) information of the annihilation radiation is performed.
  • TOF-PET measurement data the time-of-flight difference
  • step S1 a radioactivity image ⁇ ′ is calculated.
  • the reconstruction algorithm in step S1 is the above-described simultaneous reconstruction algorithm.
  • step S3 subject mask data in the projection data space (that is, subject mask projection data m proj ) is calculated based on the measurement data.
  • step S5 the offset image ⁇ off is estimated by a reconstruction algorithm configured such that the forward projection data of the offset image ⁇ off approximates the subject mask projection data.
  • At least one region ⁇ is extracted in step S6 using an image that can be recognized based on the measurement data and that can recognize the subject region.
  • the image that can be recognized based on the measurement data and in which the subject region can be recognized is, for example, the image created in step S2.
  • ⁇ ′, radioactivity image ⁇ ′ obtained in step S1 non-quantitative image estimated by a reconstruction algorithm (for example, ML-EM method) different from the reconstruction algorithm (simultaneous reconstruction algorithm) in step S1,
  • a subject mask image m img estimated from the subject mask projection data m proj is not limited to these exemplified images, and may be any image that can be recognized based on the measurement data and that can recognize the subject region.
  • At least one region ⁇ is extracted using an image recognizable from the subject region calculated based on the measurement data”. “Not only extracting the region ⁇ using a single image information but also extracting the region ⁇ using a combination of a plurality of pieces of image information (for example, logical sum or logical product).
  • Steps S7 and S8 are performed using the image ⁇ ′, the offset image ⁇ off and the region ⁇ obtained in the above steps S1 to S6.
  • step S7 a coefficient ⁇ that reduces the error between the value of the image ⁇ ′ in the region ⁇ and the known absorption coefficient value is calculated.
  • step S8 as described above (3), the value of the image mu ', corrects the value obtained by adding the alpha ⁇ mu off that coefficient alpha multiplying the offset image mu off as the absorption coefficient mu To do.
  • the difference between the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (the value of the true absorption coefficient image) is approximated by ⁇ ⁇ ⁇ off, which is obtained by multiplying the offset image ⁇ off by a coefficient ⁇ .
  • the absorption coefficient value is corrected in step S8. Therefore, the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S8 has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • steps S1 and S2 corresponding to the above-described reconstruction calculation step are performed as follows: (b) a calculation algorithm including absorption coefficient projection data in an unknown and absorption coefficient projection data This is implemented by a combination of algorithms that make the image reconstructed as an image ⁇ ′.
  • Non-Patent Document 1 which simultaneously estimates a non-quantitative radioactivity image ⁇ ′ and non-quantitative absorption coefficient projection data (absorption coefficient sinogram A ′ in each of Examples 1 to 3).
  • This is a combination of the MLACF method and an algorithm in which an image obtained by reconstructing the absorption coefficient projection data (absorption coefficient sinogram A ′) is an image ⁇ ′.
  • the algorithm for reconstructing the image obtained by reconstructing the absorption coefficient projection data (absorption coefficient sinogram A ′) in (b) is the image ⁇ ′. If it is an algorithm, it will not specifically limit.
  • step S2 the absorption coefficient sinogram A ′ is reconstructed by the ML-TR method or the ML-EM method, and an image ⁇ ′ is created.
  • the mask calculation step includes a step of calculating a binarized image of the image ⁇ ′ as the subject mask image m img (step S3) and a step of calculating projection data of the subject mask image m img (the first half of step S4).
  • the process of calculating the binarized data of the projection data of the subject mask image m img as the subject mask projection data m proj (“Binarization Processing” in the latter half of step S4).
  • the image ⁇ ′ is a non-quantitative absorption coefficient image.
  • an image other than the absorption coefficient image for example, the radioactivity image ⁇ ′ of Example 3 or the radioactivity image ⁇ 2 ′.
  • the object mask projection data m proj can also be calculated by binarizing the projection data.
  • the reconstruction process performed in step S5 corresponding to the offset estimation process described above is performed by a statistical reconstruction calculation method typified by the ML-EM (Maximum-Likelihood--Expectation-Maximization) method in the first embodiment. Yes.
  • the reconstruction process performed in step S5 is not limited to the statistical reconstruction as in the first embodiment. What is necessary is just to implement by any one of analytical reconstruction, statistical reconstruction, and algebraic reconstruction.
  • analytical reconstruction for example, there is an FBP (Filtered Back Projection) method.
  • an algebraic reconstruction for example, there is an ART (Algebraic Reconstruction Technique) method.
  • the at least one region ⁇ extracted in step S6 corresponding to the above-described reference region extraction step is a tissue region in which the absorption coefficient can be regarded as known.
  • the region of the tissue whose absorption coefficient can be regarded as known means, for example, a region that can be approximated by water, a region that can be approximated by air, These include areas that can be approximated by brain tissue, areas that can be approximated by bone, areas that can be approximated by lung tissue, and areas that can be approximated by soft tissue.
  • any structure can be used as long as an approximate value of the absorption coefficient can be understood.
  • the water absorption coefficient value of 0.0096 / mm can be used.
  • step S7 corresponding to the coefficient calculation step described above is performed. That is, the coefficient ⁇ is calculated based on the representative value in the first embodiment, including later-described second to fourth embodiments.
  • the representative value is, for example, an average value, a median value, a trimmed average value, a trimmed median value, or a weighted average value of two or more values thereof.
  • the value is not limited to these exemplified values, and may be “a statistic or a value calculated from the statistic”.
  • S ( ⁇ ′; ⁇ n ) is a representative value of image ⁇ ′ in region ⁇ n
  • S ( ⁇ off ; ⁇ n ) is a representative value of the offset image ⁇ off ′ in the region ⁇ n .
  • the coefficients ⁇ 1 ,. ..., ⁇ K can be calculated by setting the representative value T ( ⁇ 1 , ⁇ 2 ,..., ⁇ K ) as the coefficient ⁇ .
  • the absorption coefficient image estimation program 6 is configured by a computer (in each embodiment, the arithmetic circuit 5 shown in FIG. 1 is configured).
  • the absorption coefficient value is corrected in step S8 based on the mathematical relationship that it is possible to approximate with ⁇ ⁇ ⁇ off multiplied by ⁇ . Therefore, the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S8 has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • arithmetic means for executing the absorption coefficient image estimation program 6 according to the first embodiment (in each embodiment, the arithmetic circuit 5 shown in FIG. 1 is configured).
  • the difference between the non-quantitative image ⁇ ′ value in the region ⁇ and the known absorption coefficient value (the true absorption coefficient image value) is the coefficient of the offset image ⁇ off .
  • the absorption coefficient value is corrected in step S8 based on the mathematical relationship that it is possible to approximate with ⁇ ⁇ ⁇ off multiplied by ⁇ . Therefore, the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S8 has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • FIG. 4 is a flowchart illustrating a processing procedure and a data flow of the absorption coefficient image estimation method according to the second embodiment.
  • the binarized image of the image mu ' is calculated as a subject mask image m img, calculate the projection data of the object mask image m img, the projection data of the object mask image m img two The value data was calculated as the subject mask projection data m proj .
  • data obtained by binarizing the TOF-PET measurement data converted into the projection data format is calculated as the subject mask projection data m proj .
  • Step S11 MLACF Since step S11 in FIG. 4 is the same as step S1 in the first embodiment described above, description thereof is omitted.
  • Step S12 ML-TR or ML-EM Since step S12 in FIG. 4 is the same as step S2 in the first embodiment described above, description thereof is omitted. Steps S11 and S12 correspond to a reconstruction calculation step in the present invention.
  • Step S14 Projection + Binarization Processing
  • steps S3 and S4 in FIG. 3 are performed to calculate the subject mask image m img from the image ⁇ ′.
  • step S14 in FIG. 4 is performed in order to calculate the subject mask image m img from the measurement data instead of the image ⁇ ′.
  • the measurement data is converted into a projection data format (Conversion).
  • binarization processing of the projection data of the measurement data by threshold processing binarization data in which the projection line passing through the subject is “1” and the other projection lines are “0” is obtained.
  • subject mask projection data m proj Step S14 corresponds to a mask calculation step in the present invention.
  • Step S15 ML-EM Since step S15 in FIG. 4 is the same as step S5 in the first embodiment described above, description thereof is omitted. Step S15 corresponds to the offset estimation step in the present invention.
  • Step S16 Extraction Since step S16 in FIG. 4 is the same as step S6 in the first embodiment described above, description thereof is omitted. Step S16 corresponds to the reference region extraction step in the present invention.
  • the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (value of the true absorption coefficient image).
  • the absorption coefficient value is corrected in step S18 based on a mathematical relationship that the offset image ⁇ off can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by a coefficient ⁇ . Therefore, the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S18 has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • step S14 corresponding to the mask calculation process described above is performed. That is, the mask calculation step includes a step (step S14) of calculating, as object mask projection data, binary data obtained by converting the above-described TOF-PET measurement data into a projection data format.
  • step S14 data obtained by binarizing the data converted into the projection data format using the TOF-PET measurement data directly can be calculated as the subject mask projection data m proj .
  • the binarized image of the image mu ' is calculated as a subject mask image m img, calculate the projection data of the object mask image m img, the projection data of the object mask image m img two The value data was calculated as the subject mask projection data m proj .
  • data obtained by binarizing the TOF-PET measurement data converted into the projection data format is calculated as the subject mask projection data m proj .
  • the radioactivity image is calculated, the projection data of the radioactivity image is calculated, (of the radioactivity image)
  • Data obtained by binarizing the projection data is calculated as subject mask projection data m proj .
  • Step S21 MLACF Step S21 of FIG. 5 is the same as step S1 of the first embodiment described above and step S11 of the second embodiment described above, and thus description thereof is omitted.
  • Step S22 ML-TR or ML-EM Step S22 in FIG. 5 is the same as step S2 in the first embodiment described above and step S12 in the second embodiment described above, and therefore description thereof is omitted. Steps S21 and S22 correspond to a reconstruction calculation step in the present invention.
  • the line integral data (projection data) of the radioactivity image ⁇ ′ estimated based on the optimization of the evaluation function related to the measurement data in the MLACF method is calculated (Projection). That is, in FIG. 5, the projection data of the radioactivity image ⁇ ′ may be calculated using the radioactivity image ⁇ ′ estimated by the MLACF method in step S21. Then, by performing binarization processing on the projection data of the radioactivity image ⁇ ′ by threshold processing, the projection line passing through the subject is “1”, and the other projection lines are “0”. Data is calculated as subject mask projection data m proj . Step S24 corresponds to a mask calculation step in the present invention.
  • Step S25 ML-EM Step S25 in FIG. 5 is the same as step S5 in the first embodiment described above and step S15 in the second embodiment described above, and thus description thereof is omitted. Step S25 corresponds to an offset estimation step in the present invention.
  • Step S26 Extraction Step S26 in FIG. 5 is the same as step S6 in the first embodiment described above and step S16 in the second embodiment described above, and therefore description thereof is omitted. Step S26 corresponds to the reference region extraction step in the present invention.
  • Step S27 in FIG. 5 is the same as step S7 in the first embodiment described above and step S17 in the second embodiment described above, and thus description thereof is omitted.
  • Step S27 corresponds to a coefficient calculation step in the present invention.
  • Step S31 MLACF Since step S31 in FIG. 6 is the same as step S21 in FIG. 5, the description thereof is omitted.
  • Step S32 ML-TR or ML-EM Since step S32 in FIG. 6 is the same as step S22 in FIG. 5, the description thereof is omitted. Steps S31 and S32 correspond to a reconstruction calculation step in the present invention.
  • Step S33 ML-EM
  • step S24 in FIG. 5 is performed to calculate the subject mask projection data m proj using the radioactivity image ⁇ ′ estimated by the MLACF method.
  • steps S33 and S34 of FIG. 6 are performed in order to calculate the subject mask projection data m proj using a radioactivity image estimated by a reconstruction algorithm different from the MLACF method.
  • a radioactivity image is estimated based on optimization of an evaluation function related to measurement data in a reconstruction algorithm (for example, ML-EM method) different from the MLACF method.
  • a reconstruction algorithm for example, ML-EM method
  • the radioactivity image estimated by the reconstruction algorithm different from that of the MLACF method is defined as ⁇ 2 ′ in distinction from the radioactivity image ⁇ ′ estimated by the MLACF method.
  • Step S34 Projection + Binarization Processing
  • the line integral data (projection data) of the radioactivity image ⁇ 2 ′ estimated based on the optimization of the evaluation function for the measurement data in the reconstruction algorithm different from the MLACF method is calculated (Projection).
  • the projection data of the radioactivity image ⁇ 2 ′ is binarized by threshold processing, so that the projection line passing through the subject becomes “1” and the other projection lines become “0”.
  • the digitized data is calculated as subject mask projection data m proj .
  • Steps S33 and S34 correspond to a mask calculation step in the present invention.
  • Step S35 ML-EM Since step S35 in FIG. 6 is the same as step S25 in FIG. 5, the description thereof is omitted. Step S35 corresponds to an offset estimation step in the present invention.
  • Step S36 Extraction Since step S36 in FIG. 6 is the same as step S26 in FIG. 5, the description thereof is omitted. Step S36 corresponds to the reference region extraction step in the present invention.
  • Step S37 in FIG. 6 is the same as step S27 in FIG. Step S37 corresponds to a coefficient calculation step in the present invention.
  • the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value is based on the mathematical relationship that the offset image ⁇ off can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by the coefficient ⁇ , and the absorption coefficient value in step S28 of FIG. 5 or step S38 of FIG. Correct. Therefore, the systematic error of the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S28 of FIG. 5 or step S38 of FIG. 6 is reduced. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • step S24 in FIG. 5 or steps S33 and S34 in FIG. 6 corresponding to the above-described mask calculation step is performed. That is, the mask calculation step is a step of calculating a radioactivity image ( ⁇ ′ in FIG. 5 and ⁇ 2 ′ in FIG. 6) based on optimization of the evaluation function related to the measurement data of TOF-PET (described above) (FIG. 6).
  • 5 is step S21
  • FIG. 6 is step S33
  • a step of calculating projection data of the radioactivity image ⁇ ′ in FIG. 5 and ⁇ 2 ′ in FIG. 6) (“Projection” in the first half of step S24 in FIG. In FIG.
  • the reconstruction algorithm for calculating the radioactivity image is not particularly limited. As shown in FIG. 5, when the above-described simultaneous reconstruction algorithm (MLACF method or MLAA method) is used, subject mask projection data is obtained using the above-described radioactivity image ⁇ ′ estimated by the simultaneous reconstruction algorithm. m proj can be calculated. Further, as shown in FIG.
  • Radioactivity image lambda 2 'estimated in the, subject Mask projection data m proj may be calculated.
  • FIG. 7 is a flowchart illustrating the processing procedure and the data flow of the absorption coefficient image estimation method according to the fourth embodiment.
  • Step S41 MLAA
  • the radioactivity image ⁇ ′ and the absorption coefficient sinogram A ′ are estimated by the MLACF method in step S1 (in Example 2).
  • Step S11, Step S21 or S31 in Example 3 After Step S21 or S31), an image obtained by reconstructing the absorption coefficient sinogram A ′ is set as an image ⁇ ′ (Step S12 in Example 2, Step S22 or S32 in Example 3) Had been implemented.
  • the radioactivity image ⁇ ′ and the image ⁇ ′ are simultaneously calculated by the MLAA method.
  • steps S1 and S2 using the MLACF method steps S11 and S12 in the second embodiment, steps S21 and S22 or S31 and S32 in the third embodiment
  • the MLAA method is used. Is integrated into one step S41.
  • Reference 3 Reference 3: A. Rezaei (KU 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).
  • Step S41 corresponds to a reconstruction calculation step in the present invention.
  • Step S43 Binarization Processing Since step S43 in FIG. 7 is the same as step S3 in the first embodiment described above, description thereof is omitted.
  • Step S44 Projection + Binarization Processing Since step S44 of FIG. 7 is the same as step S4 of the first embodiment described above, description thereof is omitted. Steps S43 and S44 correspond to a mask calculation step in the present invention.
  • Step S45 ML-EM Step S45 in FIG. 7 is the same as Step S5 in Embodiment 1 described above, Step S15 in Embodiment 2 described above, and Step 3 in Embodiment 3 described above (Step S25 in FIG. 5 and Step S35 in FIG. 6). The description is omitted. Step S45 corresponds to an offset estimation step in the present invention.
  • Step S46 Extraction Step S46 in FIG. 7 is the same as step S6 in the first embodiment described above, step S16 in the second embodiment described above, and step S3 in the third embodiment described above (step S26 in FIG. 5 and step S36 in FIG. 6). The description is omitted. Step S46 corresponds to the reference region extraction step in the present invention.
  • Step S47 in FIG. 7 is the same as Step S7 in Embodiment 1 described above, Step S17 in Embodiment 2 described above, and Step 3 in Embodiment 3 described above (Step S27 in FIG. 5 and Step S37 in FIG. 6). The description is omitted.
  • Step S47 corresponds to a coefficient calculation step in the present invention.
  • FIG. 7 is a flowchart of the fourth embodiment when applied to the case where steps S1 and S2 in FIG. 3 of the first embodiment described above are integrated into step S43.
  • the fourth embodiment may be applied when the steps S11 and S12 in FIG. 4 of the second embodiment described above are integrated into the step S43, or the steps S21 and S22 in FIG.
  • the fourth embodiment may be applied when the steps are integrated into steps S31 and S32.
  • step S41 corresponding to the above-described reconstruction calculation step is performed by a calculation algorithm including (a) the image ⁇ ′ as an unknown.
  • the algorithm (a) is an MLAA method that simultaneously reconstructs a non-quantitative radioactivity image ⁇ ′ and an image ⁇ ′ (non-quantitative absorption coefficient image).
  • FIG. 8 is a flowchart illustrating the processing procedure and data flow of the absorption coefficient image estimation method according to the fifth embodiment.
  • the coefficient ⁇ is calculated based on the representative value.
  • the coefficient ⁇ is calculated based on the error evaluation function.
  • Step S51 MLACF Step S51 in FIG. 8 is the same as Step S1 in Embodiment 1 described above, Step S11 in Embodiment 2 described above, and Step 3 in Embodiment 3 described above (Step S21 in FIG. 5 and Step S31 in FIG. 6). The description is omitted.
  • Step S52 ML-TR or ML-EM Step S52 in FIG. 8 is the same as Step S2 in Embodiment 1 described above, Step S12 in Embodiment 2 described above, and Step 3 in Embodiment 3 described above (Step S22 in FIG. 5 and Step S32 in FIG. 6). The description is omitted. Steps S51 and S52 correspond to a reconstruction calculation step in the present invention.
  • Step S53 Binarization Processing Step S53 in FIG. 8 is the same as step S3 in the first embodiment described above and step S43 in the fourth embodiment described above, and thus description thereof is omitted.
  • Step S54 Projection + Binarization Processing Step S54 in FIG. 8 is the same as step S4 in the first embodiment described above and step S44 in the fourth embodiment described above, and therefore description thereof is omitted. Steps S53 and S54 correspond to the mask calculation step in the present invention.
  • Step S55 ML-EM Step S55 in FIG. 8 includes step S5 in the first embodiment described above, step S15 in the second embodiment described above, step S3 in the third embodiment described above (step S25 in FIG. 5, step S35 in FIG. 6), and the embodiment described above. Since this is the same as step S45 in FIG. Step S55 corresponds to the offset estimation step in the present invention.
  • Step S56 of Extraction Step S56 of FIG. 8 includes step S6 of the first embodiment described above, step S16 of the second embodiment described above, step S3 of the third embodiment described above (step S26 in FIG. 5 and step S36 in FIG. 6), and the embodiment described above. Since this is the same as step S46 in FIG. Step S56 corresponds to the reference region extraction step in the present invention.
  • f ( ⁇ ) D ⁇ ( ⁇ known , ⁇ ′ + ⁇ ⁇ ⁇ off )
  • the coefficient ⁇ is calculated based on the representative value.
  • the coefficient ⁇ is calculated based on the error evaluation function.
  • ⁇ known is an image in which a known absorption coefficient is set in the region ⁇
  • D ⁇ (X, Y) is an error evaluation function of the image X and the image Y related to the region ⁇ .
  • D ⁇ ( ⁇ known , ⁇ ' + ⁇ ⁇ ⁇ off ) is an image ⁇ known with a known absorption coefficient in the region ⁇ and a quantitative absorption coefficient
  • This is an error evaluation function for an image ( ⁇ ′ + ⁇ ⁇ ⁇ off ) obtained by adding a non-uniform offset value to the image.
  • f ( ⁇ ) D ⁇ ( ⁇ known , ⁇ '+ ⁇ ⁇ ⁇ off ) (4)
  • Step S57 corresponds to a coefficient calculation step in the present invention.
  • Step S48 in FIG. 8 includes step S8 in the above-described first embodiment, step S18 in the above-described second embodiment, step in the above-described third embodiment (step S28 in FIG. 5, step S38 in FIG. 6), and the above-described embodiment. Since this is the same as step S48 in FIG. Step S58 corresponds to the absorption coefficient value correcting step in the present invention.
  • FIG. 8 shows an implementation when the calculation of the coefficient ⁇ based on the representative value in step S7 in FIG. 3 of the first embodiment described above is replaced with the calculation of the coefficient ⁇ based on the error evaluation function in step S57.
  • 10 is a flowchart of Example 5.
  • the fifth embodiment may be applied when the calculation of the coefficient ⁇ based on the representative value in step S17 in FIG. 4 of the second embodiment described above is replaced with the calculation of the coefficient ⁇ based on the error evaluation function in step S57.
  • the calculation of the coefficient ⁇ based on the representative value in step S27 in FIG. 5 or the step S37 in FIG. 6 of the third embodiment described above is replaced with the calculation of the coefficient ⁇ based on the error evaluation function in step S57.
  • the calculation of the coefficient ⁇ based on the representative value in step S47 in FIG. 7 in the fourth embodiment described above is replaced with the calculation of the coefficient ⁇ based on the error evaluation function in step S57. May be applied.
  • the value of the non-quantitative image ⁇ ′ in the region ⁇ and the known absorption coefficient value (the true absorption coefficient image
  • the absorption coefficient value is corrected in step S58 on the basis of a mathematical relationship that the offset image ⁇ off can be approximated by ⁇ ⁇ ⁇ off obtained by multiplying the offset image ⁇ off by a coefficient ⁇ . Therefore, the absorption coefficient image having the absorption coefficient value ⁇ corrected in step S58 has a small systematic error. As a result, since a quantitative absorption coefficient image can be created, accurate absorption correction of the radioactivity image becomes possible.
  • the coefficient ⁇ is calculated based on the error evaluation function.
  • K ( ⁇ 1) be the number of regions that can be approximated by a known absorption coefficient value
  • ⁇ n is the nth region ⁇
  • w n (n 1, ..., K) and the 0 to 1 inclusive coefficients To do.
  • Image X is obtained by replacing image X n known by setting a known absorption coefficient in region ⁇ n and adding image Y to image ⁇ ′ by multiplying offset image ⁇ off ′ by ⁇ ⁇ ⁇ off.
  • 'be replaced with (+ ⁇ ⁇ ⁇ off, D ⁇ n ( ⁇ n known, ⁇ value ⁇ )' + ⁇ ⁇ ⁇ off ) includes an image mu n known set of known absorption coefficient for region Omega n, This is an error evaluation function with an image ( ⁇ ′ + ⁇ ⁇ ⁇ off ) obtained by adding a non-uniform offset value to a quantitative absorption coefficient image.
  • the present invention is not limited to the above embodiment, and can be modified as follows.
  • the subject to be photographed in each of the embodiments described above is not particularly limited. You may apply to the apparatus which image
  • the DOI detector is used, but it may be applied to a detector that does not distinguish the depth direction.
  • the detector rings are stacked in the body axis direction of the subject.
  • the detector ring may have only one configuration.
  • the data format is sinogram, but is not limited to sinogram. If it is projection data, for example, the data format may be a histogram. Therefore, in Examples 1 to 3 described above, the absorption coefficient projection data estimated by the MLACF method is the absorption coefficient sinogram A ′, but the absorption coefficient projection data estimated by the MLACF method may be an absorption coefficient histogram.
  • the projection data is calculated after binarizing the image ⁇ ′ first, and the projection data is binarized, but the object mask projection data m proj is calculated as follows. May be. That is, the mask calculation step may be an aspect including a step of calculating projection data of the image ⁇ ′ and a step of calculating binarized data of the projection data of the image ⁇ ′ as the subject mask projection data m proj . In this embodiment, the projection data of the image ⁇ ′ is binarized after being calculated first.
  • the radioactivity image in FIG. 5 lambda ', lambda 2 in Fig. 6'
  • the projection data m proj may be calculated. That is, the mask calculation step is a step of calculating a radioactivity image based on optimization of the evaluation function related to the measurement data of TOF-PET (described above), a step of calculating a binarized image of the radioactivity image,
  • the aspect which consists of the process of calculating the projection data of a binarized image and the process of calculating the data which binarized the projection data of the binarized image as object mask projection data m proj may be sufficient.
  • the projection data is calculated after binarizing the radioactivity image first, and the projection data is binarized.
  • the reconstruction algorithm for calculating the radioactivity image is not particularly limited, and a simultaneous reconstruction algorithm (MLACF method or MLAA method) as shown in FIG.
  • a reconstruction algorithm for example, ML-EM method
  • MCACF method or MLAA method different from the simultaneous reconstruction algorithm (MLACF method or MLAA method) as shown in FIG. 6 may be used.
  • the present invention is suitable for estimation of an absorption coefficient image using measurement data measured by a TOF measurement type PET apparatus.

Landscapes

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

Abstract

本発明の吸収係数画像推定方法では、領域Ωにおける定量的でない画像μ'の値と既知の吸収係数値（真の吸収係数画像の値）との差分は、オフセット画像μoffを係数α倍したα×μoffで近似可能であるという数学的な関係に基づいて、ステップＳ８(μ=μ'+α×μoff)で吸収係数値を補正する。したがって、ステップＳ８(μ=μ'+α×μoff)で補正された吸収係数値μを有する吸収係数画像は、系統誤差が小さくなる。その結果、定量的な吸収係数画像を作成することができるので、放射能画像の正確な吸収補正が可能となる。

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

Priority Applications (5)

Application Number Priority Date Filing Date Title
CN201780091422.6A CN110691992A (zh) 2017-05-29 2017-05-29 吸收系数图像估计方法、吸收系数图像估计程序以及安装有该吸收系数图像估计程序的正电子ct装置
EP17912063.9A EP3633411A4 (en) 2017-05-29 2017-05-29 CALCULATION PROCEDURE FOR ABSORPTION COEFFICIENT IMAGE, CALCULATION PROGRAM FOR ABSORPTION COEFFICIENT IMAGE AND THE POSITRON CT DEVICE EQUIPPED WITH IT
JP2019521549A JP6761610B2 (ja) 2017-05-29 2017-05-29 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置
US16/499,466 US20200085397A1 (en) 2017-05-29 2017-05-29 Attenuation coefficient image estimation method, attenuation coefficient image estimation program, and positron ct apparatus equipped with the same
PCT/JP2017/019944 WO2018220686A1 (ja) 2017-05-29 2017-05-29 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置

Applications Claiming Priority (1)

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

Publications (1)

Publication Number Publication Date
WO2018220686A1 true WO2018220686A1 (ja) 2018-12-06

Family

ID=64455240

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/019944 WO2018220686A1 (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)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021001841A (ja) * 2019-06-24 2021-01-07 株式会社島津製作所 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108986182B (zh) * 2018-07-10 2022-11-25 上海联影医疗科技股份有限公司 一种重建ct图像的方法、系统及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015519555A (ja) * 2012-05-04 2015-07-09 コーニンクレッカ フィリップス エヌ ヴェ 陽電子放出断層撮影における散乱された同時発生を用いた減衰マップ
US20150193927A1 (en) * 2014-01-08 2015-07-09 Rensselaer Polytechnic Institute Attenuation map reconstruction from tof pet data

Family Cites Families (15)

* 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
WO2013147013A1 (ja) * 2012-03-28 2013-10-03 独立行政法人放射線医学総合研究所 Ｍｒ画像からのｐｅｔ吸収補正画像生成方法及びコンピュータプログラム
EP2662022A1 (en) * 2012-05-07 2013-11-13 National Cerebral and Cardiovascular Center Method of extracting contour of tomogram
JP6335181B2 (ja) * 2012-10-26 2018-05-30 カリフォルニア大学Ｔｈｅ Ｒｅｇｅｎｔｓ ｏｆ ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｃａｌｉｆｏｒｎｉａ ポジトロン放射断層撮影（Ｐｏｓｉｔｒｏｎ−ＥｍｉｓｓｉｏｎＴｏｍｏｇｒａｐｈｙ：ＰＥＴ）画像のタイムオブフライト（Ｔｉｍｅ−Ｏｆ−Ｆｌｉｇｈｔ：ＴＯＦ）リストモード再構成のためのシステム行列を計算する方法及び装置
WO2014074148A1 (en) * 2012-11-08 2014-05-15 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
EP3067716A4 (en) * 2013-11-05 2016-11-30 Shimadzu Corp CONTOUR IMAGE GENERATING DEVICE AND DIAGNOSTIC DEVICE IN NUCLEAR MEDICINE
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
WO2018006419A1 (en) * 2016-07-08 2018-01-11 Shanghai United Imaging Healthcare Co., Ltd. System and method for generating attenuation map
WO2018163362A1 (ja) * 2017-03-09 2018-09-13 株式会社島津製作所 散乱推定方法、散乱推定プログラム並びにそれを搭載したポジトロンｃｔ装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2015519555A (ja) * 2012-05-04 2015-07-09 コーニンクレッカ フィリップス エヌ ヴェ 陽電子放出断層撮影における散乱された同時発生を用いた減衰マップ
US20150193927A1 (en) * 2014-01-08 2015-07-09 Rensselaer Polytechnic Institute Attenuation map reconstruction from tof pet data

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
LI, QUANZHENG ET AL.: "Joint estimation of activity image and attenuation sinogram using time-of-flight positron emission tomography data consistency condition filtering", JOURNAL MEDICAL IMAGING, vol. 4, no. 2, 26 April 2017 (2017-04-26), pages 023502, XP055634118, DOI: 10.1117/1.JMI.4.2.023502 *
LI, YUSHENG ET AL.: "Transmission-less attenuation estimation from time-off-flight PET histo-images using consitsensy equations", PHYSICS IN MEDICINE AND BIOLOGY, vol. 60, no. 16, 21 August 2015 (2015-08-21), pages 6563, XP020287693, DOI: 10.1088/0031-9155/60/16/6563 *
See also references of EP3633411A4 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021001841A (ja) * 2019-06-24 2021-01-07 株式会社島津製作所 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置
JP7247782B2 (ja) 2019-06-24 2023-03-29 株式会社島津製作所 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置

Also Published As

Publication number Publication date
JPWO2018220686A1 (ja) 2019-11-07
JP6761610B2 (ja) 2020-09-30
CN110691992A (zh) 2020-01-14
EP3633411A1 (en) 2020-04-08
US20200085397A1 (en) 2020-03-19
EP3633411A4 (en) 2021-03-03

Similar Documents

Publication Publication Date Title
EP1934942B1 (en) Iterative reconstruction with enhanced noise control filtering
US7718954B2 (en) Component method and system for PET detector efficiency normalization
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
JP7317586B2 (ja) 医用画像処理装置、方法及びプログラム
JP6256608B2 (ja) 画像再構成処理方法
JP6761610B2 (ja) 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンｃｔ装置
JP2021173755A (ja) 医用画像処理装置、医用画像処理方法及びプログラム
Bouwens et al. Image-correction techniques in SPECT
JP7302737B2 (ja) リストモード画像再構成方法および核医学診断装置
JP7247782B2 (ja) 吸収係数画像推定方法、吸収係数画像推定プログラム、および、ポジトロンｃｔ装置
Bentourkia et al. Simultaneous attenuation and scatter corrections in small animal PET imaging
JP2011002306A (ja) Ｐｅｔ装置の逐次近似画像再構成法
JPWO2020059135A1 (ja) データ処理方法、プログラム、データ処理装置および陽電子放出断層撮像装置
US20230260171A1 (en) Event property-dependent point spread function modeling and image reconstruction for pet
EP4318400A1 (en) Image processing apparatus, image processing method and program
US20240037814A1 (en) Convolutional neural network for dynamic pet frame clustering
Brard et al. Axially oriented crystal geometry applied to small-animal PET system: A proof of concept
Ljungberg Instrumentation, Calibration, Quantitative Imaging, and Quality Control
JP2023094592A (ja) 核医学診断装置、画像処理方法及びプログラム
Wallstén Correction for partial volume effects in PET imaging
CN113253330A (zh) 伽玛射线放射成像装置及能量校准方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17912063

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2019521549

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2017912063

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 2017912063

Country of ref document: EP

Effective date: 20200102