WO2016002084A1 - 画像再構成処理方法 - Google Patents

画像再構成処理方法 Download PDF

Info

Publication number
WO2016002084A1
WO2016002084A1 PCT/JP2014/067976 JP2014067976W WO2016002084A1 WO 2016002084 A1 WO2016002084 A1 WO 2016002084A1 JP 2014067976 W JP2014067976 W JP 2014067976W WO 2016002084 A1 WO2016002084 A1 WO 2016002084A1
Authority
WO
WIPO (PCT)
Prior art keywords
function
reconstruction processing
processing method
image reconstruction
image
Prior art date
Application number
PCT/JP2014/067976
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 JP2016530794A priority Critical patent/JP6256608B2/ja
Priority to CN201480080398.2A priority patent/CN106471392B/zh
Priority to PCT/JP2014/067976 priority patent/WO2016002084A1/ja
Priority to US15/320,244 priority patent/US10304218B2/en
Publication of WO2016002084A1 publication Critical patent/WO2016002084A1/ja

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/02Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computerised tomographs
    • A61B6/032Transmission computed tomography [CT]
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/02Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03Computerised tomographs
    • A61B6/037Emission tomography
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/42Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for detecting radiation specially adapted for radiation diagnosis
    • A61B6/4266Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment with arrangements for detecting radiation specially adapted for radiation diagnosis characterised by using a plurality of detector units
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50Clinical applications
    • A61B6/502Clinical applications involving diagnosis of breast, i.e. mammography
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/52Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/16Measuring radiation intensity
    • G01T1/17Circuit arrangements not adapted to a particular type of detector
    • G01T1/172Circuit arrangements not adapted to a particular type of detector with coincidence circuit arrangements
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914Measurement of spatial distribution of radiation
    • G01T1/2985In depth localisation, e.g. using positron emitters; Tomographic imaging (longitudinal and transverse section imaging; apparatus for radiation diagnosis sequentially in different planes, steroscopic radiation diagnosis)
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10104Positron emission tomography [PET]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10108Single photon emission computed tomography [SPECT]

Definitions

  • the present invention performs image reconstruction processing for reconstructing the physical quantity distribution of the subject related to the generation factor of the measurement data set as a multi-dimensional digital image from the measurement data set of the subject obtained by the radiation detection apparatus.
  • the present invention relates to a configuration processing method.
  • This image reconstruction processing method is used for an image reconstruction technique of a general tomographic imaging apparatus (CT (Computed Tomography) apparatus) having a radiation detection apparatus.
  • CT Computer Tomography
  • Examples of the tomographic imaging apparatus having a radiation detection apparatus include a nuclear medicine diagnostic apparatus and an X-ray computed tomography apparatus (X-ray CT apparatus).
  • X-ray CT apparatus X-ray computed tomography apparatus
  • the physical quantity distribution of the subject related to the generation factor of the measurement data set is reconstructed as a multi-dimensional digital image (such as a tomographic image or a three-dimensional reconstructed image) from the measurement data set of the subject obtained by the radiation detection apparatus. Perform reconfiguration processing.
  • a nuclear medicine diagnostic apparatus there are a positron emission tomography apparatus (PET (Positron Emission Tomography) apparatus) and a single photon emission tomography apparatus (SPECT (Single Photon Emission CT) apparatus).
  • PET Positron Emission Tomography
  • SPECT Single Photon Emission CT
  • the PET device detects only multiple rays ( ⁇ rays) generated by the annihilation of positrons (positrons) and detects the radiation ( ⁇ rays) with multiple detectors simultaneously (that is, only when they are counted simultaneously).
  • a detection signal is recorded, and a tomographic image of the subject is created by performing reconstruction processing on the detection signal (detection signals for a number of ⁇ rays).
  • the SPECT apparatus detects a single radiation ( ⁇ -ray) and performs a reconstruction process to create a tomographic image of the subject.
  • ML reconstruction method An emission CT image based on maximum likelihood estimation (ML: Maximum Likelihood) of Poisson distribution is used.
  • a reconstruction processing method (ML reconstruction method) has been proposed (see, for example, Non-Patent Document 1).
  • image reconstruction methods used in today's PET devices and SPECT devices although there are differences depending on the device manufacturer, the mathematical framework (the basis theory) of almost all methods is described in Non-Patent Document 1.
  • ML reconstruction method is described in Non-Patent Document 1.
  • Non-Patent Document 1 is a very famous academic paper on the ML reconstruction method, and almost all of today's image reconstruction methods are described in Non-Patent Document 1. It can be said that it is a derivation method of the method.
  • the measurement data set (that is, measurement data) of the subject obtained by the radiation detection apparatus includes a statistical error, and the statistical error distribution (error distribution) follows a Poisson distribution.
  • the ML reconstruction method described in Non-Patent Document 1 is a method for obtaining a solution (image) that maximizes the likelihood function derived from the Poisson property of measurement data as a likely radioactivity distribution image (physical quantity distribution).
  • the likelihood function is generally called a “data function”.
  • the likelihood function is maximized by using an iterative calculation algorithm (sequential approximation method).
  • the likelihood function of the Poisson distribution is L (x)
  • the likelihood function L (x) of the Poisson distribution is expressed by the following equation (1).
  • x is a reconstructed image vector (where the pixel value is non-negative)
  • I is the number of measurement data points
  • a i is the sensitivity distribution function of the i-th measurement data point (i-th row of system matrix A)
  • Y i is the prompt simultaneous clock value (count value) at the i-th measurement data point
  • r i is a count other than the prompt simultaneous clock value (count value) at the i-th measurement data point ( It is an estimated value of a count value (count value) of accidental coincidence count and scattering coincidence count).
  • Non-Patent Document 1 An image reconstruction method (ML reconstruction method) described in Non-Patent Document 1 and an image reconstruction method derived from the method have been established as mathematical theories. However, when these methods are applied directly to actual measurement data, (I) streak artifacts (streak artifacts) that are linear false images (linear noise) occur in the image; (Ii) the spatial resolution of the image is lower than the value predicted from the device parameters (mainly the size of the radiation detection element); Such a problem may occur. This is a problem (problem) that arises in a nuclear medicine diagnostic apparatus (emission CT apparatus), but even when it is extended to an X-ray computed tomography apparatus (X-ray CT apparatus), artifacts and degradation of the spatial resolution of the image may occur. It is thought to occur.
  • X-ray computed tomography apparatus X-ray computed tomography apparatus
  • the present invention has been made in view of such circumstances, and an object thereof is to provide an image reconstruction processing method capable of suppressing artifacts generated in an image or improving the spatial resolution of the image. .
  • the present invention has the following configuration. That is, according to the image reconstruction processing method of the present invention, the physical quantity distribution of the subject related to the generation factor of the measurement data set is obtained as a multi-dimensional digital image from the measurement data set of the subject obtained by the radiation detection apparatus.
  • An image reconstruction processing method for performing reconstruction processing for reconstruction wherein the first multivariable function having the digital image as an unknown is configured based on (1) an error distribution of element data constituting the measurement data set A data function represented by the sum of the partial functions generated, or (2) a data function represented by the sum of the partial functions configured based on the error distribution of the element data constituting the measurement data set, and the physical quantity distribution And a second multivariable function configured based on the prior information, and the image reconstruction processing method applies to the reconstructed image of the element data corresponding to the partial function.
  • Weighting the partial function with a projection weighting factor, the weighted data function, or the first function comprising the sum of the weighted data function or the second multivariable function configured based on prior information of the physical quantity distribution A reconstruction processing step for performing the reconstruction processing based on optimization calculation of a variable function is provided.
  • weighting is applied in the conventional reconstruction processing step. That is, weighting is performed when reconstruction processing is performed based on optimization calculation of a multivariable function having a digital image as an unknown, including a data function that generalizes the likelihood function of Poisson distribution of Non-Patent Document 1.
  • a multivariable function having an unknown digital image is a “first multivariable function”
  • the first multivariable function is expressed by the following (1) or (2). That is, the first multivariable function is (1) data represented by the sum of partial functions constructed based on the error distribution of element data constituting the measurement data set (of the subject obtained by the radiation detection apparatus). It is a function.
  • the first multivariable function is a multivariable function configured based on (2) the data function described in (1) and the prior information on the physical quantity distribution (of the subject related to the generation factor of the measurement data set).
  • the multivariable function configured based on the prior information of the physical quantity distribution is referred to as a “second multivariable function” in order to distinguish it from the first multivariable function.
  • the partial function is weighted with a back projection weight coefficient for the reconstructed image of the element data corresponding to the partial function described above.
  • This weighting coefficient is a non-negative coefficient (also called “influence adjustment coefficient”) that adjusts the degree of influence of the element data on the reconstructed image.
  • the radiation detection apparatus includes a positron emission tomography apparatus (PET apparatus), a single photon emission tomography apparatus (SPECT apparatus), and an X-ray computed tomography apparatus (X-ray CT apparatus).
  • PET apparatus positron emission tomography apparatus
  • SPECT apparatus single photon emission tomography apparatus
  • X-ray CT apparatus X-ray computed tomography apparatus
  • the radiation detection apparatus is any one of a PET apparatus, a SPECT apparatus, and an X-ray CT apparatus
  • the weight coefficient is not set as in the prior art. This can be attributed to the fact that reconstruction processing is performed using a function. Therefore, by setting the weighting factor based on the directivity of the linear noise generated in the reconstructed image when the weighting factor is a constant that does not depend on element data, linear noise (streak artifact) can be suppressed. .
  • a weighting factor for element data along the running direction of linear noise (however, a weighting factor greater than 0) is used as a weighting factor for element data not along the running direction of linear noise (greater than 0).
  • a weighting factor for element data not along the running direction of linear noise (greater than 0).
  • linear noise is more likely to occur when a detector unit that is partially open is used rather than a full ring type detector unit.
  • a detector unit that is partially opened radiation that passes through the opening (missing portion) is not detected. Therefore, noise having a strong spatial correlation is generated in the reconstructed image due to the loss of a part of the projection data.
  • linear noise when a plurality of detector units separated from each other are used, it is considered that linear noise (streak artifact) occurs along a linear direction connecting detection elements in the same detector unit.
  • the radiation detectors constituting the radiation detection apparatus are composed of a plurality of detector units separated from each other, in a linear direction connecting the detection elements in the same detector unit.
  • the weighting factor for the element data along is smaller than the weighting factor (the weighting factor greater than 0) for the element data along the linear direction connecting the detection elements in different detector units.
  • the radiation detection apparatus is either a positron emission tomography apparatus (PET apparatus) or a single photon emission tomography apparatus (SPECT apparatus).
  • PET apparatus positron emission tomography apparatus
  • SPECT apparatus single photon emission tomography apparatus
  • the measurement data set (of the subject obtained by the radiation detection apparatus) described above is one of sinogram data, histogram data, and list mode data. is there.
  • the radiation detection device is either a PET device or a SPECT device and the measurement data set is any of sinogram data, histogram data, and list mode data.
  • the radiation detector constituting the radiation detection apparatus is configured to measure radiation detection depth position information, that is, each detection element is configured to be stacked in the radiation depth direction.
  • a DOI detector When a DOI detector is used, the following phenomenon occurs. That is, the width of the sensitivity distribution function of the pair of detection elements (shallow DOI layer) closer to the measurement object and the pair of detection elements (deep DOI layer) on the far side becomes wider. The reliability is lower than the former. Therefore, the spatial resolution of the image can be improved by setting the weighting coefficient depending on the detected depth position information corresponding to the measurement data set.
  • the radiation detector is configured to measure N-stage detection depth position information (that is, a DOI detector),
  • the detection elements constituting the radiation detector have g, h (number of detection depth stage numbers of the two detection elements for which the coincidence count is measured so that the numbers increase from the shallow stage to the deep stage, respectively. 1 ⁇ g, h ⁇ N)
  • the weighting coefficient is a two-dimensional function with the stage numbers g and h as discrete variables, and the two-dimensional function is for the other variable obtained when one variable is fixed.
  • the one-dimensional function is a non-increasing function.
  • the error distribution described above is either a Poisson distribution or a Gaussian distribution.
  • the error distribution is a Poisson distribution, it is used for a nuclear medicine diagnostic apparatus (emission CT apparatus), and when the error distribution is a Gaussian distribution, it is used for an X-ray computed tomography apparatus (X-ray CT apparatus).
  • weighting is performed when performing reconstruction processing based on optimization calculation of a multivariable function including a data function or the like with a digital image as an unknown number. Then, by weighting the partial function with a weighting factor of back projection for the reconstructed image of the element data corresponding to the partial function, artifacts generated in the image can be suppressed or the spatial resolution of the image can be improved.
  • It is a schematic perspective view of a gamma ray detector. 3 is a flowchart of image reconstruction processing according to the first embodiment.
  • (A) is a table showing an example of a likelihood weighting factor for a pair of DOI layers, and (b) is a weight for the other variable when one variable is fixed when the number of detection depth steps is a discrete variable. It is a graph of the non-increasing function of a coefficient. It is the side view and block diagram of the mammography apparatus which concern on a modification.
  • (A), (b) is a schematic front view of the detector unit which concerns on the further modification.
  • FIG. 1 is a schematic perspective view and a block diagram showing an embodiment of a gamma ray detector arrangement of a partial ring type PET apparatus according to each example
  • FIG. 2 is a schematic perspective view of the gamma ray detector.
  • PET apparatus positron emission tomography apparatus
  • the partial ring type PET apparatus 1 includes detector units 2A and 2B.
  • a plurality of ⁇ -ray detectors 3 are embedded in the detector units 2A and 2B.
  • the partial ring type PET apparatus 1 corresponds to the radiation detection apparatus in the present invention, and also corresponds to the positron emission tomography apparatus in the present invention.
  • the detector units 2A and 2B correspond to detector units, and the ⁇ -ray detector 3 corresponds to a radiation detector in the present invention.
  • the detector units 2A and 2B are configured to be partially opened. That is, an open region (opening region) where the ⁇ -ray detector 3 does not exist exists between the detector units 2A and 2B.
  • the detector units 2A and 2B have a vertically arranged geometry.
  • the open region (opening region) is not limited to the YZ plane direction, and the ⁇ -ray detector 3 may be arranged so that the open region (opening region) exists along the XZ plane direction (in this case, detection is performed).
  • the unit 2A, 2B has a left-right geometry. Moreover, you may arrange
  • the partial ring type PET apparatus 1 includes a coincidence counting circuit 4 and an arithmetic circuit 5.
  • a coincidence counting circuit 4 and an arithmetic circuit 5.
  • PMT Photo Multiplier Tube
  • 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 electric signal is sent to the coincidence counting circuit 4.
  • 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 electrical signal sent to the coincidence counting circuit 4 is sent to the arithmetic circuit 5.
  • the arithmetic circuit 5 performs steps S1 to S6 (see FIG. 3) to be described later, and a measurement data set (here, measurement data of count values) of the subject (not shown) obtained by the partial ring type PET apparatus 1. That is, the physical quantity distribution (here, radioactivity distribution image) of the subject related to the generation factor of the measurement data set (measurement data) (here, generation of ⁇ -rays by administration of the radiopharmaceutical) from the count data in each embodiment.
  • a reconstruction process is performed to reconstruct a multi-dimensional digital image (here, a reconstructed image). Specific functions of the arithmetic circuit 5 will be described later.
  • 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
  • the spatial resolution in the depth direction can be further improved. Therefore, the number of DOI detector layers is the number of scintillator element layers stacked in the depth direction.
  • a weighting factor described later is set based on detection depth position information of a detection element (scintillator element) of the DOI detector.
  • FIG. 3 is a flowchart of image reconstruction processing according to the first embodiment
  • FIG. 4 is a schematic front view of a partial ring type detector unit
  • FIG. 5 is a partial ring type detector unit and a weighting factor.
  • FIG. 6 is a schematic diagram showing the relationship between streak artifacts and the setting of the weighting factor.
  • the detector units 2A and 2B have a vertically arranged geometry, and as shown in FIG. 4, the relationship between the arrangement of the top plate and bed (both not shown) on which the subject M is placed.
  • a subject M is arranged in proximity to the upper and upper detector units 2A.
  • the ⁇ -ray detector 3 (not shown in FIG. 4) includes a plurality of detector units separated from each other (two detector units in FIGS. 1 and 4). 2A, 2B).
  • the count data measured in the pair of detection elements (scintillator elements) in the same detector unit 2A the count data of the pairs of detection elements (scintillator elements) in the detector units 2A and 2B different from each other. Also, by applying a small weighting factor and performing steps S1 to S6 described later, streak artifacts can be suppressed.
  • the element data constituting the measurement data set is the i-th measurement data point in the first embodiment including the second and third embodiments described later.
  • the linear direction connecting the detection elements (scintillator elements) in the same detector unit 2A is LOR AA
  • the element data (i-th measurement data point) along the linear direction LOR AA The weighting factor is w AA
  • the linear direction connecting the detector elements (scintillator elements) in different detector units 2A and 2B is LOR AB
  • element data (i-th measurement data point) along the linear direction LOR AB ) Is assumed to be w AB .
  • the same detector weighting factors w AA for detecting element element data along a linear direction LOR AA connecting the (scintillator elements) (i th measurement data points) in unit 2A, different detector units 2A the weighting factor for sensing element element data along a linear direction LOR AB connecting the (scintillator elements) (i th measurement data points) in 2B is smaller than w AB.
  • the linear direction LOR AB is not limited to one direction as long as the detection elements (scintillator elements) in different detector units 2A and 2B are connected to each other, and the detection elements in different detector units 2A and 2B ( Various linear directions corresponding to the linear direction connecting the scintillator elements become the linear direction LOR AB .
  • the weighting factor w AA for the element data (i-th measurement data point) along the linear direction LOR AA connecting the detection elements (scintillator elements) in the same detector unit 2A is larger than 0,
  • the value is set to a value smaller than 1 means that the weighting factor w AA on the linear direction LOR AA is set to a value larger than 0 and smaller than 1. . That is, even if it is parallel to the linear direction LOR AA (see circle in FIG. 5) as shown in FIG. 5, the detection elements in the detector units 2A and 2B are parallel to the linear direction LOR AA and are different from each other.
  • the weighting coefficient w AB riding on the linear direction LOR AB is set to 1.
  • the detector units are not limited to the detector units 2A and 2B that are separated from each other. Even if one detector unit is used, if it is a detector unit that is partially opened, it passes through the opening (missing portion). Gamma rays to be detected are not detected. Therefore, noise having a strong spatial correlation (for example, streak artifact) occurs in the reconstructed image due to the loss of a part of the projection data. Therefore, assuming that the streak artifact is SA as shown in FIG. 6, the streak artifact SA, which is a linear noise, is reconstructed using a partial function without setting a weighting coefficient as in the prior art. It is thought to be caused. Therefore, the weighting factor may be set based on the directivity of the streak artifact SA generated in the reconstructed image when the weighting factor is a constant that does not depend on the element data (i-th measurement data point).
  • the solid line shown in FIG. 6 is the streak artifact SA
  • the broken line shown in FIG. 6 is a straight line parallel to the streak artifact SA (see ⁇ in FIG. 6), and a straight line that does not correspond to the streak artifact SA ( ⁇ in FIG. 6). Reference).
  • the weighting factor for the element data (i-th measurement data point) along the travel direction of the streak artifact SA is w SA
  • w EX be the weighting factor for the measurement data point).
  • the weight coefficient w SA for the element data (i-th measurement data point) along the travel direction of the streak artifact SA is used as the weight coefficient w SA for the element data (i-th measurement data point) not along the streak artifact SA.
  • the weight coefficient w SA is set to a value larger than 0 and smaller than 1 (0 ⁇ w SA ⁇ 1).
  • the weighting factor riding on the straight line is the streak artifact SA.
  • a weighting factor w EX respect have not element data (i-th measured data points) that along, set to 1 for weighting coefficient w EX.
  • the weighting coefficient set in this way is a non-negative coefficient (influence adjustment coefficient) that adjusts the degree of influence of element data (i-th measurement data point) on the reconstructed image, as is apparent from the above reason. Therefore, the weighting factor set in this way is also a weighting factor for back projection of the element data (i-th measurement data point) on the reconstructed image. Using this weighting coefficient, a partial function described later is weighted.
  • measurement data When applied to a nuclear medicine diagnostic apparatus (emission CT apparatus) typified by a positron emission tomography apparatus (PET apparatus) as in the first embodiment, including later-described second and third embodiments, measurement data
  • the error distribution of element data (i-th measurement data point) constituting the set (measurement data) is a Poisson distribution
  • the data function is a likelihood function. Therefore, weighting is performed on the likelihood function of the Poisson distribution.
  • the weighted likelihood function is referred to as a “weighted likelihood function”.
  • the measurement data set (measurement data) is sinogram data or histogram data
  • the weighted likelihood function of the Poisson distribution is L (x) as in the conventional case
  • the weighted likelihood function L of the Poisson distribution is L (X) is represented by the following formula (2).
  • x is a reconstructed image vector (where the pixel value is non-negative)
  • I is the number of measurement data points
  • a i is the sensitivity of the i-th measurement data point.
  • Distribution function (i-th row vector of the system matrix A)
  • y i is an immediate simultaneous clock value (count value) at the i-th measurement data point
  • r i is an immediate coincidence at the i-th measurement data point This is an estimated value of a count value (count value) of a count (coincident coincidence count and scattered coincidence count) other than the clock value (count value).
  • w i in the above equation (2) is a weighting coefficient (likelihood weighting coefficient) of the i-th measurement data point.
  • each measurement data point is weighted by the weighting factor w i. It is a point. That is, the role of the likelihood weighting factor is to adjust the degree of influence of each measurement data point on the reconstructed image.
  • a content of ⁇ in the right hand side of equation (2) "w i a i ⁇ x-w i y i log (w i a i + r i ) ”can be defined as a partial function constructed based on the error distribution of the element data (i-th measurement data point).
  • the data function (likelihood function in each embodiment) is represented by the sum of partial functions configured based on the error distribution of the element data (i-th measurement data point). Further, when a multivariable function having an unknown number of digital images after reconstruction processing is referred to as a “first multivariable function”, the first multivariable function is a data function (likelihood function).
  • the reconstruction process is performed based on the optimization calculation of the weighted likelihood function of the above-described weighted expression (2).
  • the weighted likelihood function of the above equation (2) can be maximized by the following iterative calculation algorithm (sequential approximation method), and a digital image after reconstruction processing can be acquired. .
  • the specific reconstruction process is as shown in FIG.
  • Step S1 Setting of Weighting Factor
  • Step S2 Setting of initial image A non-negative image is set as an initial image x (0) .
  • a nuclear medicine diagnostic apparatus emission CT apparatus
  • PET apparatus positron emission tomography apparatus
  • 0 is set. Except for x (0) > 0.
  • the initial image x (0) may be a reconstructed image having a uniform pixel value, for example.
  • Step S3 Initialization of the counter variable of the number of iterations
  • Step S4 Update Image Calculation
  • the (k + 1) -th update image x (k + 1) is calculated using the following equation (3).
  • J is the number of pixels of the reconstructed image.
  • the likelihood weight coefficient w i is also applied to the following equation (3) in order to maximize the weighted likelihood function of the above equation (2).
  • Step S5 Increment of counter variable of iteration count
  • the counter variable k is incremented (k ⁇ k + 1).
  • K ⁇ k + 1 means that (k + 1) on the right side is substituted for k on the left side.
  • Step S6 k ⁇ N iter ? The number of iterations to terminate the iterative calculation algorithm and N iter, whether the counter variable k has reached the number of iterations N iter? Judge whether or not. Note that the user may set the number of iterations N iter in advance. If k ⁇ N iter , the process returns to step S4 to continue the iterative calculation algorithm. If k ⁇ N iter , the iterative calculation algorithm is completed and a series of calculations is terminated.
  • the update image x (k + 1) obtained in this way is acquired as a reconstructed image.
  • the user observes the updated image x (k + 1) obtained at each update, and the user interrupts the iterative calculation algorithm based on the observation result.
  • the obtained update image x (k + 1) may be acquired as a reconstructed image.
  • steps S1 to S6 correspond to the reconstruction processing step in the present invention.
  • the image reconstruction processing method is characterized in that weighting is applied in the conventional reconstruction processing step. That is, weighting is performed when reconstruction processing is performed based on optimization calculation of a multivariable function having a digital image as an unknown, including a data function that generalizes the likelihood function of Poisson distribution of Non-Patent Document 1. .
  • the first multivariable function is a radiation detector (partial ring type PET apparatus in each embodiment).
  • Error distribution of element data (i-th measurement data point in each embodiment) constituting the measurement data set (here, measurement data of count values, that is, count data in each embodiment) obtained by 1) data functions (each embodiment is expressed by the sum of ( "w i a i ⁇ x-w i y i log (w i a i + r i) " in the above (2)) configured partial function on the basis of (Likelihood function). Then, the partial function is weighted with a weighting factor of back projection for the reconstructed image of the element data (i-th measurement data point) corresponding to the partial function described above.
  • This weight coefficient is a non-negative coefficient (also referred to as “influence adjustment coefficient”) for adjusting the influence of the element data (i-th measurement data point) on the reconstructed image as described above. By weighting, artifacts generated in the image can be suppressed and the spatial resolution of the image can be improved.
  • the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus).
  • PET apparatus positron emission tomography apparatus
  • SPECT apparatus single photon emission tomography apparatus
  • X-ray CT apparatus X-ray computed tomography apparatus
  • linear noise is more likely to occur when a detector unit that is partially open is used rather than a full ring type detector unit.
  • a detector unit that is partially opened radiation that passes through the opening (missing portion) is not detected. Therefore, noise having a strong spatial correlation is generated in the reconstructed image due to the loss of a part of the projection data.
  • linear noise is generated along the linear direction connecting the detection elements in the same detector unit.
  • the radiation detector ( ⁇ -ray detector 3 in each embodiment) constituting the radiation detection apparatus is changed into a plurality of detector units (two detectors) separated from each other.
  • the linear noise smooth artifact
  • the measurement data set (measurement data) of the subject M obtained by the radiation detection apparatus Is one of sinogram data, histogram data, and list mode data.
  • the measurement data set (measurement data) is sinogram data or histogram data
  • the weighted likelihood function L (x) of the Poisson distribution is expressed by the above equation (2) in the first embodiment.
  • the error distribution of the element data (i-th measurement data point) constituting the measurement data set (measurement data) is a Poisson distribution.
  • the error distribution is a Poisson distribution, it is used in a nuclear medicine diagnostic apparatus (emission CT apparatus) represented by a PET apparatus or the like.
  • FIG. 7 is a schematic diagram showing a relationship between a four-layer (four-stage) DOI detector and a detection depth stage number
  • FIG. 8A shows an example of a likelihood weighting coefficient for a pair of DOI layers
  • FIG. 8B is a graph of a non-increasing function of the weighting coefficient for the other variable when one variable is fixed when the number of detection depth stage numbers is a discrete variable.
  • FIG. 7 shows only the scintillator block 31 (see FIG. 2) for the ⁇ -ray detector 3 composed of the DOI detector for simplification of illustration, and the light guide 32 and the other configurations.
  • the illustration of the PMT 33 (see FIG. 2 for both) is omitted, and only four scintillator blocks 31 in the horizontal direction are shown.
  • the weighting factor is set based on the directivity of the linear noise generated in the reconstructed image when the weighting factor is a constant that does not depend on the element data (i-th measurement data point).
  • the weighting coefficient is set based on the detection depth position information of the detection element (scintillator element) of the DOI detector.
  • a PET apparatus DOI-PET apparatus having a DOI detector in which detection elements (scintillator elements) are stacked in multiple stages in the depth direction as shown in FIG.
  • detection elements sintillator elements
  • the latter sensitivity distribution function has a wider width, and the latter is less reliable than the former. . Therefore, the spatial resolution of the image is improved by setting the weighting coefficient depending on the detection depth position information corresponding to the measurement data set (measurement data).
  • the ⁇ -ray detector 3 is configured to measure detection depth position information of four stages (ie, four layers). . That is, the ⁇ -ray detector 3 is composed of a four-layer DOI detector.
  • the number of steps of the detection depth of the two detection elements (scintillator elements) for which the coincidence count was measured is g and h (1 ⁇ g, h ⁇ N), respectively, so that the number increases from the shallow level toward the deep level.
  • the number of DOI detector stages that is, the number of layers
  • the number of DOI detector stages is not particularly limited as long as it is a natural number (that is, a plurality) of 2 or more.
  • the weighting coefficient w is expressed by a two-dimensional function having the stage number numbers g and h as discrete variables. As described above, the reliability of the count data measured for the deep DOI layer pair is lower than the reliability of the count data measured for the shallow DOI layer pair.
  • the weighting coefficient w is set so that the coefficient w is smaller than the weighting coefficient w in the shallow DOI layer pair. As shown in FIG.
  • FIG. 8A shows an example of the likelihood weighting factor for the pair of DOI layers.
  • the two-dimensional function related to the weighting coefficient w is represented by a one-dimensional function for the other variable h as shown in FIG.
  • the non-increasing function is 25.
  • a one-dimensional function (non-increasing function) for the other variable obtained when one variable is fixed, and a one-dimensional function (non-increasing) for one variable obtained when the other variable is fixed. Function but not necessarily the same.
  • a one-dimensional function (non-increasing function) with respect to may be set with different functions.
  • the weighting factor set in this way is applied to the above equation (2) in the same manner as in the first embodiment.
  • the weighted likelihood function of the above equation (2) can be maximized by the following iterative calculation algorithm (sequential approximation method) as in the first embodiment, and a digital image after reconstruction processing is acquired. can do.
  • the specific reconstruction process is as shown in FIG. 3 as in the first embodiment. Note that the reference numerals of the steps in the second embodiment are the same as those in the first embodiment (S1 to S6).
  • Step S1 Setting of Weighting Factor
  • likelihood weighting factors w i for all data points are set.
  • the weight coefficient is set based on the detected depth position information of the detection element (scintillator element) of the DOI detector as described above in the second embodiment. For example, as shown in FIGS. 2 and 7, in the case of a 4-layer (4-stage) DOI detector, the likelihood weighting coefficient shown in FIG. 8 is set for each pair of DOI layers.
  • Step S2 to (Step S6) Steps S2 to S6 are the same as steps S2 to S6 of the first embodiment described above, and a description thereof will be omitted. As described above, steps S1 to S6 correspond to the reconstruction processing step in the present invention.
  • the optimal multivariable function including a data function (likelihood function in each embodiment) and the like with a digital image as an unknown number is used. Weighting is performed when the reconstruction process is performed based on the calculation.
  • the subfunction above (2) "w i a i ⁇ x-w i y i log (w i a i + r i) " in the equation) corresponding element data in the (i-th measured data points in each of the embodiments
  • the spatial resolution of the image can be improved by weighting the partial function with the weighting factor of the back projection for the reconstructed image.
  • the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus) as in the first embodiment described above and the third embodiment described later.
  • Measurement data of the subject M obtained by the radiation detection apparatus (partial ring type PET apparatus 1 in each embodiment) when the radiation detection apparatus is a PET apparatus as in each embodiment or a SPECT apparatus.
  • the set here, count value measurement data, that is, count data in each embodiment
  • the set is any of sinogram data, histogram data, and list mode data.
  • the weighted likelihood function L (x) of the Poisson distribution in the second embodiment is It is represented by the above formula (2).
  • the radiation detection apparatus is a PET apparatus as in each embodiment or a SPECT apparatus
  • the partial function is not set without setting the weighting coefficient as in the conventional case.
  • reconstruction processing is performed using
  • the radiation detector ( ⁇ -ray detector 3 in each embodiment) constituting the radiation detection apparatus is configured to measure radiation detection depth position information, that is, each detection element is irradiated with radiation.
  • a DOI detector constructed by laminating in the depth direction is used, the following phenomenon occurs. That is, the width of the sensitivity distribution function of the pair of detection elements (shallow DOI layer) closer to the measurement object and the pair of detection elements (deep DOI layer) on the far side becomes wider. The reliability is lower than the former. Therefore, the spatial resolution of the image can be improved by setting the weighting coefficient depending on the detection depth position information corresponding to the measurement data set (measurement data).
  • the radiation detector ( ⁇ -ray detector 3) has N detection depth positions (four steps). It is configured to measure information (ie, configured with a DOI detector).
  • the numbers are g and h (1 ⁇ g, h ⁇ N), respectively.
  • the weighting factor is a two-dimensional function with the stage numbers g and h as discrete variables, and the two-dimensional function is a non-increasing function with respect to the other variable obtained when one variable is fixed. .
  • the count data measured by the pair of deep DOI layers with low reliability is weighted by multiplying the count data measured by the pair of shallow DOI layers with high reliability by a weighting factor.
  • the spatial resolution of the image can be improved.
  • the error distribution of element data (i-th measurement data point) constituting the measurement data set (measurement data) is a Poisson distribution.
  • the error distribution is a Poisson distribution, it is used in a nuclear medicine diagnostic apparatus (emission CT apparatus) represented by a PET apparatus or the like.
  • the weighted likelihood function of the above equation (2) is the case where the measurement data is sinogram data or histogram data.
  • a weighted likelihood function is set when the measurement data is list mode data.
  • the list mode data is data in which detection event information (detector number, detection time, ⁇ -ray energy, etc.) obtained by the radiation detector of the PET apparatus is stored in time series.
  • the weighted likelihood function L (x) of Poisson distribution is expressed by the following equation (4).
  • N is the number of events (number of lists), and i (n) is the number of the measurement data point at which the nth event is measured (1 ⁇ i (n) ⁇ N).
  • the weighting coefficient set in this way is applied to the above equation (4).
  • the weighted likelihood function of the above equation (4) can be maximized by the iterative calculation algorithm (sequential approximation method) shown in FIG. 3 as in the first and second embodiments, and the digital image after the reconstruction process Can be obtained.
  • steps S1 to S6 shown in FIG. 3 are the same as steps S1 to S6 of the first and second embodiments described above, description thereof will be omitted.
  • step S4 updated image calculation
  • the data format is list mode data
  • steps S1 to S6 correspond to the reconstruction processing step in the present invention.
  • the weighting coefficient in the above equation (4) is based on the directivity of linear noise generated in the reconstructed image when the weighting coefficient is a constant that does not depend on element data as in the first embodiment.
  • a weighting factor may be set, or the weighting factor may be set based on the detection depth position information of the detection element of the DOI detector as in the second embodiment.
  • the present invention is not limited to the above embodiment, and can be modified as follows.
  • the positron emission tomography apparatus has been described as an example of the radiation detection apparatus.
  • the apparatus is a device that acquires a measurement data set of a subject based on detection of radiation. If there is, it will not be specifically limited. You may apply to a single photon emission tomography apparatus (SPECT apparatus), an X-ray computed tomography apparatus (X-ray CT apparatus), etc.
  • SPECT apparatus single photon emission tomography apparatus
  • X-ray CT apparatus X-ray computed tomography apparatus
  • the present invention may be applied to an apparatus that captures the entire body of the subject, an apparatus that captures the head of the subject, and an apparatus that captures the breast of the subject.
  • the partial ring type PET apparatus 1 as shown in FIG. 1 is used. May be.
  • the detector units 2A and 2B in FIG. 1 have the same configuration as that in FIG. 1 except that the detector units 2A and 2B are replaced with a breast examination unit 2C.
  • the breast inspection unit 2C has a notch, and the breast is inspected by being sandwiched by this notch.
  • a plurality of ⁇ -ray detectors 3 (not shown in FIG. 9) are arranged in parallel in the breast examination unit 2C in accordance with the notches.
  • the DOI detector is used.
  • the DOI detector may be applied to a radiation detector that does not distinguish the depth direction.
  • the weighting factor is set based on the directivity of the linear noise generated in the reconstructed image when the weighting factor is a constant that does not depend on element data as in the first embodiment described above, the above-described implementation is performed.
  • the weighting coefficient can be set without using the detection depth position information of the detection element of the DOI detector as in Example 2.
  • the PET apparatus when the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus), the PET apparatus (part in FIG. 1) includes a detector unit that is partially opened.
  • the ring-type PET apparatus 1) is not necessarily a detector unit that is partially opened, but may be applied to a normal full-ring type PET apparatus.
  • the weighting coefficient when the weighting coefficient is set based on the detection depth position information of the detection element of the DOI detector as in the second embodiment, the linear noise directivity as in the first embodiment described above is set. A weighting factor can be set without using it.
  • the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus)
  • PET apparatus positron emission tomography apparatus
  • two detector units 2A and 2B separated from each other see FIGS. 1, 4 and 5).
  • one detector unit 2D having a part opened as shown in FIG. 10A may be used.
  • the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus)
  • PET apparatus positron emission tomography apparatus
  • two detector units 2A and 2B separated from each other see FIGS. 1, 4 and 5).
  • the radiation detection apparatus is a positron emission tomography apparatus (PET apparatus)
  • PET apparatus positron emission tomography apparatus
  • detector units 2E, 2F, 2G, 2H separated from each other as shown in FIG. It may be a container unit.
  • the first multivariable function constitutes a measurement data set (measurement data) of the subject M obtained by the radiation detection apparatus (partial ring type PET apparatus 1 in each embodiment).
  • the data function (likelihood function in each embodiment) represented by the sum of partial functions configured based on the error distribution of element data (i-th measurement data point in each embodiment)
  • the variable function may be the sum of the above-described data function (likelihood function) and the second multivariable function configured based on the prior information on the physical quantity distribution.
  • the weighted likelihood function L (x) of the above equation (2) or (4) and the prior information on the physical quantity distribution An image reconstruction method based on maximization of L (x) + U (x) consisting of the sum of other functions U (x) consisting of the second multivariable function configured based on the above is also conceivable.
  • L (x) is a function derived on the basis of statistical properties, whereas U (x) is based on prior information of a subject to be photographed (a deterministic property that the reconstructed image x may have). Is a function defined.
  • U (x) is generally called “PenaltyalFunction”. An example of the successive approximation calculation formula when the penalty function is used is shown below (see the following formula (6)).
  • c (k) j is the j-th pixel value of the approximate curvature image of the regularization function U (x) in the vicinity of the k-th estimated solution x (k) .
  • the error distribution of the element data (i-th measurement data point) constituting the measurement data set (measurement data) is the Poisson distribution, but is not limited to the Poisson distribution.
  • a Gaussian distribution may be used.
  • the error distribution is a Gaussian distribution, it is used in an X-ray computed tomography apparatus (X-ray CT apparatus).
  • the present invention has a radiation detection apparatus such as a positron emission tomography apparatus (PET apparatus), a single photon emission tomography apparatus (SPECT apparatus), and an X-ray computed tomography apparatus (X-ray CT apparatus). It is suitable for the image reconstruction technology of the tomographic imaging apparatus in general.
  • PET apparatus positron emission tomography apparatus
  • SPECT apparatus single photon emission tomography apparatus
  • X-ray CT apparatus X-ray computed tomography apparatus

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • Public Health (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biomedical Technology (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Dentistry (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Pulmonology (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Quality & Reliability (AREA)
  • Nuclear Medicine (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Measurement Of Radiation (AREA)

Abstract

 この発明の画像再構成処理方法では、従来の再構成処理工程において、重み付けを加えることを特徴とするものである。すなわち、ポアソン分布の尤度関数を一般化したデータ関数などからなる、デジタル画像を未知数とする多変数関数の最適化計算に基づいて再構成処理を行う際に、重み付けを行う。例えば、重み係数を要素データに依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定する、あるいはDOI検出器の検出素子の検出深さ位置情報に基づいて重み係数を設定する(ステップS1)。そして、部分関数に対応する要素データの再構成画像に対する逆投影の重み係数で部分関数を重み付けることにより、画像に生じるアーティファクトを抑制し、または画像の空間分解能を向上させることができる。

Description

画像再構成処理方法
 この発明は、放射線検出装置によって得られた被検体の計測データ集合から、計測データ集合の発生要因に関連する被検体の物理量分布を複数次元のデジタル画像として再構成する再構成処理を行う画像再構成処理方法に関する。
 この画像再構成処理方法は、放射線検出装置を有した断層画像撮影装置(CT(Computed Tomography)装置)全般の画像再構成技術に用いられる。放射線検出装置を有した断層画像撮影装置としては、例えば、核医学診断装置やX線コンピュータ断層撮影装置(X線CT装置)がある。放射線検出装置によって得られた被検体の計測データ集合から、計測データ集合の発生要因に関連する被検体の物理量分布を複数次元のデジタル画像(断層画像や3次元再構成画像など)として再構成する再構成処理を行う。
 また、核医学診断装置としては、陽電子放射断層撮影装置(PET(Positron Emission Tomography)装置)や単一光子放射断層撮影装置(SPECT(Single Photon Emission CT)装置)がある。PET装置は、陽電子(ポジトロン)の消滅によって発生する複数本の放射線(γ線)を検出して複数個の検出器で放射線(γ線)を同時に検出したときのみ(つまり同時計数したときのみ)検出信号を記録し、その検出信号(多数のγ線の検出信号)に対して再構成処理を行って被検体の断層画像を作成する。SPECT装置は、単一の放射線(γ線)を検出して再構成処理を行って被検体の断層画像を作成する。
 PET装置やSPECT装置などの核医学診断装置(エミッションCT装置)を例に採って説明すると、エミッションCT装置の分野においては、ポアソン分布の最尤推定(ML: Maximum Likelihood)に基づくエミッションCT画像の再構成処理手法(ML再構成方法)が提案されている(例えば、非特許文献1参照)。今日のPET装置やSPECT装置で利用されている画像再構成手法において、装置メーカにより差異はあるものの、ほぼ全ての手法の数学的な枠組み(土台となる理論)は非特許文献1に記載されているML再構成方法である。その意味において、エミッションCT装置の分野においては、非特許文献1はML再構成方法に関する非常に有名な学術論文であり、今日の画像再構成方法は、ほぼ全て非特許文献1に記載されている手法の派生手法と言える。
 放射線検出装置によって得られた被検体の計測データ集合(すなわち測定データ)は統計誤差を含み、その統計誤差の分布(誤差分布)はポアソン分布にしたがう。非特許文献1に記載のML再構成方法は、測定データのポアソン性から導かれる尤度関数を最大化する解(画像)を尤もらしい放射能分布画像(物理量分布)として求める方法である。なお、X線CT装置の分野などのように、ポアソン分布以外の誤差分布(例えばガウス分布)にまで拡げた場合には、一般的には尤度関数は「データ関数」とも呼ばれる。なお、尤度関数の最大化は反復計算アルゴリズム(逐次近似法)を利用して実施される。
 また、ポアソン分布の尤度関数をL(x)とすると、ポアソン分布の尤度関数L(x)は下記(1)式で表される。
Figure JPOXMLDOC01-appb-M000001
 ここで、xは再構成画像ベクトル(ただし画素値は非負)であり、Iは測定データ点の数であり、aはi番目の測定データ点の感度分布関数(システム行列Aの第i行ベクトル)であり、yはi番目の測定データ点での即発同時計数値(カウント値)であり、rはi番目の測定データ点での即発同時計数値(カウント値)以外の計数(偶発同時計数および散乱同時計数)の計数値(カウント値)の推定値である。
L. A. Shepp and Y. Vardi. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging, Vol. 1, pp.113-122, 1982.
 非特許文献1に記載の画像再構成方法(ML再構成方法)や当該手法から派生した画像再構成方法は、数学的な理論としては確立されている。しかしながら、それらの方法を実際の測定データに直接に適用すると、
 (i)画像に直線状の偽像(直線状ノイズ)であるスジ状アーティファクト(ストリークアーティファクト)が生じる、
 (ii)画像の空間分解能が装置パラメータ(主に放射線検出素子のサイズ)から予測される値よりも低い、
 といった問題が生じる場合がある。これは、核医学診断装置(エミッションCT装置)において生じる課題(問題)であるが、X線コンピュータ断層撮影装置(X線CT装置)にまで拡げた場合においてもアーティファクトや画像の空間分解能の劣化が生じると考えられる。
 この発明は、このような事情に鑑みてなされたものであって、画像に生じるアーティファクトを抑制し、または画像の空間分解能を向上させることができる画像再構成処理方法を提供することを目的とする。
 この発明は、このような目的を達成するために、次のような構成をとる。
 すなわち、この発明の画像再構成処理方法は、放射線検出装置によって得られた被検体の計測データ集合から、前記計測データ集合の発生要因に関連する前記被検体の物理量分布を複数次元のデジタル画像として再構成する再構成処理を行う画像再構成処理方法であって、前記デジタル画像を未知数とする第1多変数関数は、(1)前記計測データ集合を構成する要素データの誤差分布に基づいて構成された部分関数の和で表されるデータ関数,または(2)前記計測データ集合を構成する要素データの誤差分布に基づいて構成された部分関数の和で表されるデータ関数と、前記物理量分布の事前情報に基づいて構成された第2多変数関数との和であり、前記画像再構成処理方法は、前記部分関数に対応する要素データの再構成画像に対する逆投影の重み係数で部分関数を重み付けて、重み付けされた前記データ関数,または前記データ関数と、前記物理量分布の事前情報に基づいて構成された第2多変数関数との和からなる前記第1多変数関数の最適化計算に基づいて前記再構成処理を行う再構成処理工程を備えることを特徴とするものである。
 この発明の画像再構成処理方法によれば、従来の再構成処理工程において、重み付けを加えることを特徴とするものである。すなわち、非特許文献1のポアソン分布の尤度関数を一般化したデータ関数などからなる、デジタル画像を未知数とする多変数関数の最適化計算に基づいて再構成処理を行う際に、重み付けを行う。ここで、デジタル画像を未知数とする多変数関数を「第1多変数関数」とすると、第1多変数関数は下記(1)または(2)で表される。すなわち、第1多変数関数は、(1)(放射線検出装置によって得られた被検体の)計測データ集合を構成する要素データの誤差分布に基づいて構成された部分関数の和で表されるデータ関数である。または、第1多変数関数は、(2)(1)で述べたデータ関数と、(計測データ集合の発生要因に関連する被検体の)物理量分布の事前情報に基づいて構成された多変数関数(以下、物理量分布の事前情報に基づいて構成された多変数関数を第1多変数関数と区別するために「第2多変数関数」とする)との和である。そして、上述した部分関数に対応する要素データの再構成画像に対する逆投影の重み係数で部分関数を重み付ける。この重み係数は、要素データの再構成画像への影響度を調節する非負係数(「影響調節係数」とも呼ばれる)であり、この重み係数で部分関数を重み付けることにより、画像に生じるアーティファクトを抑制し、または画像の空間分解能を向上させることができる。
 この発明に係る画像再構成処理方法において、放射線検出装置は、陽電子放射断層撮影装置(PET装置),単一光子放射断層撮影装置(SPECT装置),X線コンピュータ断層撮影装置(X線CT装置)のいずれかである。
 放射線検出装置が、PET装置,SPECT装置,X線CT装置のいずれかである場合において、上述した直線状ノイズ(ストリークアーティファクト)が生じる場合には、従来のように重み係数を設定せずに部分関数を用いて再構成処理を行っていることに起因すると考えられる。そこで、重み係数を要素データに依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定することで、直線状ノイズ(ストリークアーティファクト)を抑制することができる。
 より具体的には、直線状ノイズの走行方向に沿った要素データに対する重み係数(ただし、0より大きい重み係数)を、直線状ノイズの走行方向に沿っていない要素データに対する重み係数(0より大きい重み係数)よりも小さな値に設定する。これにより、直線状ノイズ(ストリークアーティファクト)に沿った計測データ集合の影響を相対的に弱める。その結果、直線状ノイズ(ストリークアーティファクト)を抑制することができる。
 上述した直線状ノイズ(ストリークアーティファクト)は、フルリング型の検出器ユニットよりも、一部が開口して構成された検出器ユニットを用いた場合に生じやすい。一部が開口して構成された検出器ユニットの場合には、開口部分(抜け部分)を通過する放射線は検出されない。よって、投影データの一部分が欠損することに起因して、再構成画像に空間的相関性の強いノイズが生じる。特に、互いに分離した複数の検出器ユニットを用いた場合には、同一の検出器ユニット内の検出素子を結んだ直線方向に沿って直線状ノイズ(ストリークアーティファクト)が生じると考えられる。
 そこで、見方を変えて、放射線検出装置を構成する放射線検出器が、互いに分離した複数の検出器ユニットで構成されている場合には、同一の検出器ユニット内の検出素子を結んだ直線方向に沿った要素データに対する重み係数(ただし、0より大きい重み係数)を、互いに異なる検出器ユニット内の検出素子を結んだ直線方向に沿った要素データに対する重み係数(0より大きい重み係数)よりも小さな値に設定する。これにより、同一の検出器ユニット内の検出素子を結んだ直線方向に沿って現れた直線状ノイズ(ストリークアーティファクト)を抑制することができる。
 X線CT装置を除いた核医学診断装置(エミッションCT装置)の場合、放射線検出装置は、陽電子放射断層撮影装置(PET装置),単一光子放射断層撮影装置(SPECT装置)のいずれかである。放射線検出装置が、PET装置,SPECT装置のいずれかである場合に、上述した(放射線検出装置によって得られた被検体の)計測データ集合は、サイノグラムデータ,ヒストグラムデータ,リストモードデータのいずれかである。
 放射線検出装置が、PET装置,SPECT装置のいずれかであり、計測データ集合が、サイノグラムデータ,ヒストグラムデータ,リストモードデータのいずれかである場合において、上述した画像の空間分解能の劣化が生じる場合には、従来のように重み係数を設定せずに部分関数を用いて再構成処理を行っていることに起因すると考えられる。特に、放射線検出装置を構成する放射線検出器が、放射線の検出深さ位置情報を計測するように構成されている場合、すなわち、各々の検出素子を放射線の深さ方向に積層して構成されたDOI検出器を用いた場合、下記のような現象が起こる。つまり、測定対象に近い側の検出素子(浅いDOI層)のペアと、遠い側の検出素子(深いDOI層)のペアとでは、後者の感度分布関数の幅の方が広くなって、後者の方が前者よりも信頼度が低くなる。そこで、計測データ集合に対応した検出深さ位置情報に依存して重み係数を設定することで、画像の空間分解能を向上させることができる。
 より具体的には、Nを2以上の自然数としたときに、放射線検出器はN段の検出深さ位置情報を計測するように構成されており(すなわちDOI検出器で構成されており)、放射線検出器を構成する検出素子であって、浅い段から深い段に向かうにしたがって番号が大きくなるように、同時計数を計測した2つの検出素子の検出深さの段数番号をそれぞれg,h(1≦g,h≦N)としたときに、重み係数は段数番号g,hを離散変数とする2次元関数であり、2次元関数は一方の変数を固定したときに得られる他方の変数に対する1次元関数が非増加関数である。これにより、信頼度が低い深いDOI層のペアで計測された計数データに対して、信頼度が高い浅いDOI層のペアで計測された計数データよりも小さな重み係数を掛けて重み付けを行うことで、画像の空間分解能を向上させることができる。
 これらの発明に係る画像再構成処理方法において、上述した誤差分布は、ポアソン分布,ガウス分布のいずれかである。誤差分布がポアソン分布の場合には核医学診断装置(エミッションCT装置)に用いられ、誤差分布がガウス分布の場合にはX線コンピュータ断層撮影装置(X線CT装置)に用いられる。
 この発明に係る画像再構成処理方法によれば、データ関数などからなる、デジタル画像を未知数とする多変数関数の最適化計算に基づいて再構成処理を行う際に、重み付けを行う。そして、部分関数に対応する要素データの再構成画像に対する逆投影の重み係数で部分関数を重み付けることにより、画像に生じるアーティファクトを抑制し、または画像の空間分解能を向上させることができる。
各実施例に係る部分リング型PET装置のγ線検出器配置の一実施態様を示す概略斜視図およびブロック図である。 γ線検出器の概略斜視図である。 実施例1に係る画像再構成処理のフローチャートである。 部分リング型の検出器ユニットの概略正面図である。 部分リング型の検出器ユニットと重み係数の設定との関係を示した模式図である。 ストリークアーティファクトと重み係数の設定との関係を示した模式図である。 4層(4段)のDOI検出器と検出深さの段数番号との関係を示した模式図である。 (a)はDOI層のペアに対する尤度重み係数の一例を示すテーブル、(b)は検出深さの段数番号を離散変数とした場合において、一方の変数を固定したときの他方の変数に対する重み係数の非増加関数のグラフである。 変形例に係るマンモグラフィ装置の側面図およびブロック図である。 (a)、(b)は、さらなる変形例に係る検出器ユニットの概略正面図である。
 以下、図面を参照してこの発明の実施例1を説明する。図1は、各実施例に係る部分リング型PET装置のγ線検出器配置の一実施態様を示す概略斜視図およびブロックであり、図2は、γ線検出器の概略斜視図である。後述する実施例2、3も含めて本実施例1では、放射線検出装置として、陽電子放射断層撮影装置(PET装置)を例に採って説明する。また、図1および図2は各実施例とも共通の構成である。
 図1に示すように、部分リング型PET装置1は、検出器ユニット2A,2Bを備えている。検出器ユニット2A,2B内には複数のγ線検出器3が埋設されている。部分リング型PET装置1は、この発明における放射線検出装置に相当し、この発明における陽電子放射断層撮影装置にも相当する。また、検出器ユニット2A,2Bは、検出器ユニットに相当し、γ線検出器3は、この発明における放射線検出器に相当する。
 また、検出器ユニット2A,2Bは、一部が開口して構成されている。すなわち、検出器ユニット2A,2B間には、γ線検出器3が存在しない開放領域(開口領域)が存在する。図1の場合には、開放領域(開口領域)はYZ平面(図4および図5を参照)に沿って存在しているので、検出器ユニット2A,2Bは上下配置型のジオメトリとなる。もちろん、開放領域(開口領域)はYZ平面方向に限定されず、XZ平面方向に沿って開放領域(開口領域)が存在するようにγ線検出器3を配置してもよい(このときは検出器ユニット2A,2Bは左右配置型のジオメトリとなる)。また、YZ平面,XZ平面以外の平面に沿って開放領域(開口領域)が存在するようにγ線検出器3を配置してもよい。
 その他にも、部分リング型PET装置1は、同時計数回路4と演算回路5とを備えている。図1では、γ線検出器3から同時計数回路4への結線を2つのみ図示しているが、実際には、γ線検出器3の光電子増倍管(PMT: Photo Multiplier Tube)33(図2を参照)の総チャンネル数分、同時計数回路4に接続されている。
 放射性薬剤が投与された被検体(図示省略)から発生したγ線をγ線検出器3のシンチレータブロック31(図2を参照)が光に変換して、変換されたその光をγ線検出器3の光電子増倍管(PMT)33(図2を参照)は増倍させて電気信号に変換する。その電気信号を同時計数回路4に送り込む。
 具体的には、被検体(図示省略)に放射性薬剤を投与すると、ポジトロン放出型のRIのポジトロンが消滅することにより、2本のγ線が発生する。同時計数回路4は、シンチレータブロック31(図2を参照)の位置とγ線の入射タイミングとをチェックし、被検体の両側にある2つのシンチレータブロック31でγ線が同時に入射したときのみ、送り込まれた電気信号を適正なデータと判定する。一方のシンチレータブロック31のみにγ線が入射したときには、同時計数回路4は棄却する。つまり、同時計数回路4は、上述した電気信号に基づいて、2つのγ線検出器3においてγ線が同時観測(すなわち同時計数)されたことを検出する。
 同時計数回路4に送り込まれた電気信号を、演算回路5に送り込む。演算回路5は、後述するステップS1~S6(図3を参照)を行って、部分リング型PET装置1によって得られた被検体(図示省略)の計測データ集合(ここではカウント値の測定データ、すなわち各実施例では計数データ)から、計測データ集合(測定データ)の発生要因(ここでは放射性薬剤の投与によるγ線の発生)に関連する被検体の物理量分布(ここでは放射能分布画像)を複数次元のデジタル画像(ここでは再構成画像)として再構成する再構成処理を行う。演算回路5の具体的な機能については後述する。
 γ線検出器3は、図2に示すようにシンチレータブロック31と、そのシンチレータブロック31に対して光学的に結合されたライトガイド32と、そのライトガイド32に対して光学的に結合された光電子増倍管(以下、単に「PMT」と略記する)33とを備えている。シンチレータブロック31を構成する各シンチレータ素子は、γ線の入射に伴って発光することでγ線から光に変換する。この変換によってシンチレータ素子はγ線を検出する。シンチレータ素子において発光した光がシンチレータブロック31で十分に拡散されて、ライトガイド32を介してPMT33に入力される。PMT33は、シンチレータブロック31で変換された光を増倍させて電気信号に変換する。その電気信号は画素値として同時計数回路4(図1を参照)に送り込まれる。
 また、γ線検出器3は、図2に示すように、3次元的に配置されたシンチレータ素子からなり、深さ方向に複数の層からなるDOI検出器である。図2では、4層のDOI検出器を図示しているが、層の数については、複数であれば特に限定されない。
 ここで、DOI検出器は、各々のシンチレータ素子を放射線の深さ方向に積層して構成されたものであり、相互作用を起こした深さ(DOI: Depth of Interaction)方向と横方向(入射面に平行な方向)との座標情報を重心演算により求める。DOI検出器を用いることにより深さ方向の空間分解能をより一層向上させることができる。よって、DOI検出器の層の数は、深さ方向に積層されたシンチレータ素子の層の数である。後述する実施例2では、DOI検出器の検出素子(シンチレータ素子)の検出深さ位置情報に基づいて後述する重み係数を設定する。
 次に、演算回路5の具体的な機能について、図3~図6を参照して説明する。図3は、実施例1に係る画像再構成処理のフローチャートであり、図4は、部分リング型の検出器ユニットの概略正面図であり、図5は、部分リング型の検出器ユニットと重み係数の設定との関係を示した模式図であり、図6は、ストリークアーティファクトと重み係数の設定との関係を示した模式図である。
 図3のフローチャートの説明よりも前に、先ず重み係数の設定について説明する。図1でも述べたように、検出器ユニット2A,2Bは上下配置型のジオメトリであり、図4に示すように被検体Mを載置する天板やベッド(いずれも図示省略)の配置の関係上、上側の検出器ユニット2Aに近接して被検体Mが配置される。なお、本実施例では、図1および図4に示すようにγ線検出器3(図4では図示省略)は、互いに分離した複数の検出器ユニット(図1および図4では2つの検出器ユニット2A,2B)で構成されている。
 かかる部分リング型の検出器ユニット2A,2Bを用いた場合には、従来のML再構成方法により画像再構成処理を実施すると、同一の検出器ユニット内の検出素子(シンチレータ素子)を結んだ直線方向に沿って直線状ノイズ(ストリークアーティファクト)が生じる。特に、図4に示すように上側の検出器ユニット2Aに近接して被検体Mを配置した場合には、同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向に沿ってストリークアーティファクトが生じる。
 そこで、同一の検出器ユニット2A内の検出素子(シンチレータ素子)のペアにおいて計測された計数データに対して、互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)のペアの計数データよりも小さな重み係数を掛けて、後述するステップS1~S6を行うことで、ストリークアーティファクトを抑制することができる。なお、計測データ集合(測定データ、すなわち計数データ)を構成する要素データは、後述する実施例2、3も含めて本実施例1ではi番目の測定データ点である。
 図5に示すように、同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向をLORAAとし、当該直線方向LORAAに沿った要素データ(i番目の測定データ点)に対する重み係数をwAAとし、互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)を結んだ直線方向をLORABとし、当該直線方向LORABに沿った要素データ(i番目の測定データ点)に対する重み係数をwABとする。このとき、同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向LORAAに沿った要素データ(i番目の測定データ点)に対する重み係数wAAを、互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)を結んだ直線方向LORABに沿った要素データ(i番目の測定データ点)に対する重み係数をwABよりも小さくする。
 一般的に重み係数は、0よりも大きく、1以下であるので、重み係数wAAを、0よりも大きく、かつ1よりも小さい値に設定し(0<wAA<1)、重み係数wABを1に設定する(wAB=1)。なお、直線方向LORAAは同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向であれば一方向だけでなく、同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向に該当する様々な直線方向が直線方向LORAAとなる。同じく、直線方向LORABは互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)を結んだ直線方向であれば一方向だけでなく、互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)を結んだ直線方向に該当する様々な直線方向が直線方向LORABとなる。
 ここで、「同一の検出器ユニット2A内の検出素子(シンチレータ素子)を結んだ直線方向LORAAに沿った要素データ(i番目の測定データ点)に対する重み係数wAAを、0よりも大きく、かつ1よりも小さい値に設定する」とは、当該直線方向LORAAに乗った重み係数wAAを、0よりも大きく、かつ1よりも小さい値に設定するという意味であることに留意されたい。つまり、図5に示すような直線方向LORAA(図5中の○を参照)に平行であっても、当該直線方向LORAAに平行で、かつ互いに異なる検出器ユニット2A,2B内の検出素子(シンチレータ素子)を結んだ直線方向LORAB(図5中の○を参照)であれば、当該直線方向LORABに乗った重み係数wABについては1に設定する。
 また、互いに分離した検出器ユニット2A,2Bに限定されず、1つの検出器ユニットであっても、一部が開口して構成された検出器ユニットであれば、開口部分(抜け部分)を通過するγ線は検出されない。よって、投影データの一部分が欠損することに起因して、再構成画像に空間的相関性の強いノイズ(例えばストリークアーティファクト)が生じる。よって、図6に示すようにストリークアーティファクトをSAとすると、直線状ノイズであるストリークアーティファクトSAは、従来のように重み係数を設定せずに部分関数を用いて再構成処理を行っていることに起因するとも考えられる。そこで、重み係数を要素データ(i番目の測定データ点)に依存しない定数とした場合に再構成画像に生じるストリークアーティファクトSAの指向性に基づいて重み係数を設定してもよい。
 図6に示す実線をストリークアーティファクトSAとし、図6に示す破線を、ストリークアーティファクトSA(図6中の○を参照)に平行な直線で、ストリークアーティファクトSAに該当しない直線(図6中の○を参照)とする。なお、ストリークアーティファクトSAの走行方向に沿った要素データ(i番目の測定データ点)に対する重み係数をwSAとし、図6に示す破線も含めてストリークアーティファクトSAに沿っていない要素データ(i番目の測定データ点)に対する重み係数をwEXとする。このとき、ストリークアーティファクトSAの走行方向に沿った要素データ(i番目の測定データ点)に対する重み係数wSAを、ストリークアーティファクトSAに沿っていない要素データ(i番目の測定データ点)に対する重み係数wEXよりも小さくすることで、ストリークアーティファクトSAに沿った計測データ集合(測定データ)の影響を相対的に弱める。
 上述したように一般的に重み係数は、0よりも大きく、1以下であるので、重み係数wSAを、0よりも大きく、かつ1よりも小さい値に設定し(0<wSA<1)、重み係数wEXを1に設定する(wEX=1)。ここで、「ストリークアーティファクトSAの走行方向に沿った要素データ(i番目の測定データ点)に対する重み係数wSAを、0よりも大きく、かつ1よりも小さい値に設定する」とは、ストリークアーティファクトSAに乗った重み係数wSAを、0よりも大きく、かつ1よりも小さい値に設定するという意味であることに留意されたい。つまり、図6に示すようなストリークアーティファクトSAに平行であっても、ストリークアーティファクトSAに該当しない直線(図6の破線を参照)であれば、当該直線に乗った重み係数は、ストリークアーティファクトSAに沿っていない要素データ(i番目の測定データ点)に対する重み係数wEXであるとして、重み係数wEXについては1に設定する。
 このように設定された重み係数は、以上の理由から明らかなように、要素データ(i番目の測定データ点)の再構成画像への影響度を調節する非負係数(影響調節係数)である。よって、このように設定された重み係数は、要素データ(i番目の測定データ点)の再構成画像に対する逆投影の重み係数でもある。この重み係数を用いて、後述する部分関数を重み付ける。後述する実施例2、3も含めて、本実施例1のような陽電子放射断層撮影装置(PET装置)などに代表される核医学診断装置(エミッションCT装置)に適用する場合には、計測データ集合(測定データ)を構成する要素データ(i番目の測定データ点)の誤差分布はポアソン分布となり、データ関数は尤度関数となる。よって、ポアソン分布の尤度関数に対して重み付けを行う。
 以下、重み付けが行われた尤度関数を「重み付き尤度関数」と呼ぶ。計測データ集合(測定データ)が、サイノグラムデータやヒストグラムデータである場合には、従来と同様に、ポアソン分布の重み付き尤度関数をL(x)とすると、ポアソン分布の重み付き尤度関数L(x)は下記(2)式で表される。
Figure JPOXMLDOC01-appb-M000002
 ここで、従来の上記(1)と同様に、xは再構成画像ベクトル(ただし画素値は非負)であり、Iは測定データ点の数であり、aはi番目の測定データ点の感度分布関数(システム行列Aの第i行ベクトル)であり、yはi番目の測定データ点での即発同時計数値(カウント値)であり、rはi番目の測定データ点での即発同時計数値(カウント値)以外の計数(偶発同時計数および散乱同時計数)の計数値(カウント値)の推定値である。さらに、上記(2)式中のwは、i番目の測定データ点の重み係数(尤度重み係数)である。
 従来の尤度関数(上記(1)式)と重み付き尤度関数(上記(2)式)との相違点は、各測定データ点に対して重み係数wが掛かって重み付けが行われている点である。つまり、尤度重み係数の役割は、各測定データ点の再構成画像への影響度を調節することである。なお、上記(2)式の右辺のΣを全体で括ると、上記(2)式の右辺におけるΣの中身である「w・x-wlog(w+r)」を、要素データ(i番目の測定データ点)の誤差分布に基づいて構成された部分関数と定義づけることができる。すなわち、上記(2)式から、データ関数(各実施例では尤度関数)は、要素データ(i番目の測定データ点)の誤差分布に基づいて構成された部分関数の和で表される。また、再構成処理後のデジタル画像を未知数とする多変数関数を「第1多変数関数」とすると、第1多変数関数はデータ関数(尤度関数)である。
 このように重み付けされた上記(2)式の重み付き尤度関数の最適化計算に基づいて再構成処理を行う。具体的には、上記(2)式の重み付き尤度関数は、以下に示す反復計算アルゴリズム(逐次近似法)により最大化することができ、再構成処理後のデジタル画像を取得することができる。具体的な再構成処理については、図3に示す通りである。
 (ステップS1)重み係数の設定
 上記(2)式において、全てのデータ点に対する尤度重み係数wを設定する。具体的には、本実施例1では同一の検出器ユニット2A内の検出素子(シンチレータ素子)のペアによるデータ点に対してはw=α(0<α<1)、それ以外のデータ点に対してはw=1とする。
 (ステップS2)初期画像の設定
 非負画像を初期画像x(0)とする。後述する実施例2、3も含めて、本実施例1のような陽電子放射断層撮影装置(PET装置)などに代表される核医学診断装置(エミッションCT装置)に適用する場合には、0を除いてx(0)>0とする。初期画像x(0)については、例えば一様な画素値を有する再構成画像であればよい。
 (ステップS3)反復回数のカウンタ変数の初期化
 反復計算アルゴリズム(逐次近似法)における反復回数のカウンタ変数をkとして、反復回数のカウンタ変数kを初期化する(k=0)。
 (ステップS4)更新画像の計算
 下記(3)式を用いて、(k+1)回目の更新画像x(k+1)を計算する。ただし、Jは再構成画像の画素数である。下記(3)式から明らかなように、上記(2)式の重み付き尤度関数を最大化するために、尤度重み係数wが下記(3)式にも掛かる。
Figure JPOXMLDOC01-appb-M000003
 (ステップS5)反復回数のカウンタ変数のインクリメント
 カウンタ変数kをインクリメントする(k←k+1)。なお、「k←k+1」とは右辺の(k+1)を左辺のkに代入するという意味である。
 (ステップS6)k<Niter
 反復計算アルゴリズムを終了する反復回数をNiterとし、カウンタ変数kが反復回数Niterに達したか?否かを判断する。なお、反復回数Niterについてはユーザが予め設定すればよい。k<Niterの場合には反復計算アルゴリズムを継続するために、ステップS4に戻る。k<Niterの場合には反復計算アルゴリズムが終了したとして一連の計算を終了する。
 このようにして得られた更新画像x(k+1)を再構成画像として取得する。また、反復回数Niterを設定せずに、更新の度に得られた更新画像x(k+1)をユーザが観察して、観察結果に基づいて反復計算アルゴリズムをユーザが中断して、そのときに得られた更新画像x(k+1)を再構成画像として取得してもよい。以上のように、ステップS1~S6は、この発明における再構成処理工程に相当する。
 本実施例1に係る画像再構成処理方法によれば、従来の再構成処理工程において、重み付けを加えることを特徴とするものである。すなわち、非特許文献1のポアソン分布の尤度関数を一般化したデータ関数などからなる、デジタル画像を未知数とする多変数関数の最適化計算に基づいて再構成処理を行う際に、重み付けを行う。上述したように、デジタル画像を未知数とする多変数関数を「第1多変数関数」とすると、本実施例1では第1多変数関数は、放射線検出装置(各実施例では部分リング型PET装置1)によって得られた被検体Mの計測データ集合(ここではカウント値の測定データ、すなわち各実施例では計数データ)を構成する要素データ(各実施例ではi番目の測定データ点)の誤差分布に基づいて構成された部分関数(上記(2)式では「w・x-wlog(w+r)」)の和で表されるデータ関数(各実施例では尤度関数)である。そして、上述した部分関数に対応する要素データ(i番目の測定データ点)の再構成画像に対する逆投影の重み係数で部分関数を重み付ける。この重み係数は、上述したように要素データ(i番目の測定データ点)の再構成画像への影響度を調節する非負係数(「影響調節係数」とも呼ばれる)であり、この重み係数で部分関数を重み付けることにより、画像に生じるアーティファクトを抑制し、画像の空間分解能を向上させることができる。
 後述する実施例2、3も含めて、本実施例1では、放射線検出装置は、陽電子放射断層撮影装置(PET装置)である。放射線検出装置が、各実施例のようなPET装置である場合、あるいは単一光子放射断層撮影装置(SPECT装置),X線コンピュータ断層撮影装置(X線CT装置)のいずれかである場合において、直線状ノイズ(ストリークアーティファクト)が生じる場合には、従来のように重み係数を設定せずに部分関数を用いて再構成処理を行っていることに起因すると考えられる。そこで、重み係数を要素データ(i番目の測定データ点)に依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定することで、直線状ノイズ(ストリークアーティファクト)を抑制することができる。
 より具体的には、直線状ノイズ(図6ではストリークアーティファクトSA)の走行方向に沿った要素データ(i番目の測定データ点)に対する重み係数wSA(ただし、0より大きい重み係数)を、直線状ノイズ(ストリークアーティファクトSA)の走行方向に沿っていない要素データ(i番目の測定データ点)に対する重み係数wEX(0より大きい重み係数)よりも小さな値に設定する(例えば0<wSA<1,wEX=1)。これにより、直線状ノイズ(ストリークアーティファクトSA)に沿った計測データ集合(測定データ)の影響を相対的に弱める。その結果、直線状ノイズ(ストリークアーティファクトSA)を抑制することができる。
 上述した直線状ノイズ(ストリークアーティファクト)は、フルリング型の検出器ユニットよりも、一部が開口して構成された検出器ユニットを用いた場合に生じやすい。一部が開口して構成された検出器ユニットの場合には、開口部分(抜け部分)を通過する放射線は検出されない。よって、投影データの一部分が欠損することに起因して、再構成画像に空間的相関性の強いノイズが生じる。特に、後述する実施例2、3も含めて、本実施例1のように互いに分離した複数の検出器ユニット(図1および図4では2つの検出器ユニット2A,2B)を用いた場合には、同一の検出器ユニット内の検出素子を結んだ直線方向に沿って直線状ノイズ(ストリークアーティファクト)が生じると考えられる。
 そこで、見方を変えて、放射線検出装置(部分リング型PET装置1)を構成する放射線検出器(各実施例ではγ線検出器3)が、互いに分離した複数の検出器ユニット(2つの検出器ユニット2A,2B)で構成されている場合には、同一の検出器ユニット内の検出素子を結んだ直線方向に沿った要素データ(i番目の測定データ点)に対する重み係数(ただし、0より大きい重み係数)を、互いに異なる検出器ユニット内の検出素子を結んだ直線方向に沿った要素データ(i番目の測定データ点)に対する重み係数(0より大きい重み係数)よりも小さな値に設定する(例えば0<wAA<1,wAB=1)。これにより、同一の検出器ユニット内の検出素子を結んだ直線方向に沿って現れた直線状ノイズ(ストリークアーティファクト)を抑制することができる。
 放射線検出装置が、各実施例のようなPET装置である場合、あるいはSPECT装置である場合に、放射線検出装置(部分リング型PET装置1)によって得られた被検体Mの計測データ集合(測定データ)は、サイノグラムデータ,ヒストグラムデータ,リストモードデータのいずれかである。特に、計測データ集合(測定データ)が、サイノグラムデータやヒストグラムデータである場合には、本実施例1ではポアソン分布の重み付き尤度関数L(x)は上記(2)式で表される。
 後述する実施例2、3も含めて、本実施例1では、計測データ集合(測定データ)を構成する要素データ(i番目の測定データ点)の誤差分布は、ポアソン分布である。誤差分布がポアソン分布の場合には、PET装置などに代表される核医学診断装置(エミッションCT装置)に用いられる。
 次に、図面を参照してこの発明の実施例2を説明する。図7は、4層(4段)のDOI検出器と検出深さの段数番号との関係を示した模式図であり、図8(a)は、DOI層のペアに対する尤度重み係数の一例を示すテーブルであり、図8(b)は、検出深さの段数番号を離散変数とした場合において、一方の変数を固定したときの他方の変数に対する重み係数の非増加関数のグラフである。
 なお、図7では、図示の簡略化のために、DOI検出器からなるγ線検出器3については、シンチレータブロック31(図2を参照)のみを図示し、その他の構成であるライトガイド32やPMT33(いずれも図2を参照)については図示を省略し、横方向のシンチレータブロック31についても4つのみ図示する。
 上述した実施例1では、重み係数を要素データ(i番目の測定データ点)に依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定した。これに対して、本実施例2では、DOI検出器の検出素子(シンチレータ素子)の検出深さ位置情報に基づいて重み係数を設定する。
 具体的には、図2に示すような深さ方向に多段に検出素子(シンチレータ素子)を積層したDOI検出器を備えたPET装置(DOI-PET装置)では、測定対象に近い側の検出素子(浅いDOI層)のペアと、遠い側の検出素子(深いDOI層)のペアとでは、後者の感度分布関数の幅の方が広くなって、後者の方が前者よりも信頼度が低くなる。そこで、計測データ集合(測定データ)に対応した検出深さ位置情報に依存して重み係数を設定することで、画像の空間分解能を向上させる。
 Nを2以上の自然数としたときに、図7および図8ではN=4として、γ線検出器3は4段(すなわち4層)の検出深さ位置情報を計測するように構成されている。すなわち、γ線検出器3は4層のDOI検出器で構成されている。浅い段から深い段に向かうにしたがって番号が大きくなるように、同時計数を計測した2つの検出素子(シンチレータ素子)の検出深さの段数番号をそれぞれg,h(1≦g,h≦N)とする。図7および図8ではN=4であるので、浅い段から深い段に向かうにしたがって、図7に示すように、g=1,2,3,4となり、h=1,2,3,4となる。実施例1の図2でも述べたように、DOI検出器の段数(すなわち層の数)については、2以上の自然数(すなわち複数)であれば特に限定されない。
 段数番号g,hは自然数であるので、段数番号g,hを離散変数とすると、重み係数wは段数番号g,hを離散変数とする2次元関数で表される。上述したように、深いDOI層のペアで計測された計数データでの信頼度が、浅いDOI層のペアで計測された計数データでの信頼度よりも低いので、深いDOI層のペアでの重み係数wが、浅いDOI層のペアでの重み係数wよりも小さくなるように、重み係数wを設定する。図7に示すように、浅い段から深い段に向かうにしたがって、g=1,2,3,4として、h=1,2,3,4とした場合には、一方の変数gを固定すると重み係数wに関する2次元関数は、他方の変数hに対する1次元関数で表され、当該1次元関数は非増加関数となる。逆に、一方の変数hを固定した場合においても、同様に重み係数wに関する2次元関数は、他方の変数gに対する1次元関数で表され、当該1次元関数は非増加関数となる。
 図8(a)にDOI層のペアに対する尤度重み係数の一例を示す。一方の変数gを1に固定すると、図8(b)に示すように重み係数wに関する2次元関数は、他方の変数hに対する1次元関数で表される。そして、当該1次元関数は、h=1,2のときには重み係数wが1(図8(a)では「1.00」で表記)となり、h=3,4のときには重み係数wが0.25となる非増加関数である。
 重み係数については図8に示すような値に限定されない。また、図8では、gを固定したときのh=1,2での重み係数を同じ値で設定し、gを固定したときのh=3,4での重み係数を同じ値で設定し、hを固定したときのg=1,2での重み係数を同じ値で設定し、hを固定したときのg=3,4での重み係数を同じ値で設定したが、非増加関数であれば、g,hの値が増える度に段階的に減少するように重み係数を設定してもよい。
 また、図8では、一方の変数を固定したときに得られる他方の変数に対する1次元関数(非増加関数)と、他方の変数を固定したときに得られる一方の変数に対する1次元関数(非増加関数)とは同一であったが、必ずしも同一である必要はない。個々のγ線検出器3の特性に応じて、一方の変数を固定したときに得られる他方の変数に対する1次元関数(非増加関数)と、他方の変数を固定したときに得られる一方の変数に対する1次元関数(非増加関数)とを互いに異なる関数でそれぞれ設定してもよい。
 このように設定された重み係数を、上述した実施例1と同様に上記(2)式に適用する。上記(2)式の重み付き尤度関数は、上述した実施例1と同様に、以下に示す反復計算アルゴリズム(逐次近似法)により最大化することができ、再構成処理後のデジタル画像を取得することができる。具体的な再構成処理については、上述した実施例1と同様に図3に示す通りである。なお、本実施例2におけるステップの符号については、上述した実施例1と同じステップの符号(S1~S6)を付する。
 (ステップS1)重み係数の設定
 上記(2)式において、全てのデータ点に対する尤度重み係数wを設定する。上述した実施例1と相違するのは、本実施例2では上述したようにDOI検出器の検出素子(シンチレータ素子)の検出深さ位置情報に基づいて重み係数を設定する点である。例えば、図2および図7に示すように、4層(4段)のDOI検出器の場合には、各DOI層のペアに対して図8に示す尤度重み係数を設定する。
 (ステップS2)~(ステップS6)
 ステップS2~S6は、上述した実施例1のステップS2~S6と同じであるので、その説明については省略する。以上のように、ステップS1~S6は、この発明における再構成処理工程に相当する。
 本実施例2に係る画像再構成処理方法によれば、上述した実施例1と同様に、データ関数(各実施例では尤度関数)などからなる、デジタル画像を未知数とする多変数関数の最適化計算に基づいて再構成処理を行う際に、重み付けを行う。そして、部分関数(上記(2)式では「w・x-wlog(w+r)」)に対応する要素データ(各実施例ではi番目の測定データ点)の再構成画像に対する逆投影の重み係数で部分関数を重み付けることにより、画像の空間分解能を向上させることができる。
 上述した実施例1や後述する実施例3と同様に、本実施例2では、放射線検出装置は、陽電子放射断層撮影装置(PET装置)である。放射線検出装置が、各実施例のようなPET装置である場合、あるいはSPECT装置である場合に、放射線検出装置(各実施例では部分リング型PET装置1)によって得られた被検体Mの計測データ集合(ここではカウント値の測定データ、すなわち各実施例では計数データ)は、サイノグラムデータ,ヒストグラムデータ,リストモードデータのいずれかである。特に、計測データ集合(測定データ)が、サイノグラムデータやヒストグラムデータである場合には、上述した実施例1でも述べたように本実施例2ではポアソン分布の重み付き尤度関数L(x)は上記(2)式で表される。
 放射線検出装置が、各実施例のようなPET装置である場合、あるいはSPECT装置である場合において、画像の空間分解能の劣化が生じる場合には、従来のように重み係数を設定せずに部分関数を用いて再構成処理を行っていることに起因すると考えられる。特に、放射線検出装置を構成する放射線検出器(各実施例ではγ線検出器3)が、放射線の検出深さ位置情報を計測するように構成されている場合、すなわち、各々の検出素子を放射線の深さ方向に積層して構成されたDOI検出器を用いた場合、下記のような現象が起こる。つまり、測定対象に近い側の検出素子(浅いDOI層)のペアと、遠い側の検出素子(深いDOI層)のペアとでは、後者の感度分布関数の幅の方が広くなって、後者の方が前者よりも信頼度が低くなる。そこで、計測データ集合(測定データ)に対応した検出深さ位置情報に依存して重み係数を設定することで、画像の空間分解能を向上させることができる。
 より具体的には、Nを2以上の自然数(図7および図8ではN=4)としたときに、放射線検出器(γ線検出器3)はN段(4段)の検出深さ位置情報を計測するように構成されている(すなわちDOI検出器で構成されている)。放射線検出器(γ線検出器3)を構成する検出素子であって、浅い段から深い段に向かうにしたがって番号が大きくなるように、同時計数を計測した2つの検出素子の検出深さの段数番号をそれぞれg,h(1≦g,h≦N)とする。その際に、重み係数は段数番号g,hを離散変数とする2次元関数であり、2次元関数は一方の変数を固定したときに得られる他方の変数に対する1次元関数が非増加関数である。これにより、信頼度が低い深いDOI層のペアで計測された計数データに対して、信頼度が高い浅いDOI層のペアで計測された計数データよりも小さな重み係数を掛けて重み付けを行うことで、画像の空間分解能を向上させることができる。
 上述した実施例1や後述する実施例3と同様に、本実施例2では、計測データ集合(測定データ)を構成する要素データ(i番目の測定データ点)の誤差分布は、ポアソン分布である。誤差分布がポアソン分布の場合には、PET装置などに代表される核医学診断装置(エミッションCT装置)に用いられる。
 次に、この発明の実施例3を説明する。
 上述した実施例1、2では、上記(2)式の重み付き尤度関数は、測定データがサイノグラムデータまたはヒストグラムデータの場合であった。これに対して、本実施例3では、測定データがリストモードデータの場合において重み付き尤度関数を設定する。ここで、リストモードデータとは、PET装置の放射線検出器で得られた検出イベント情報(検出器番号、検出時間、γ線のエネルギ等)を時系列で保存したデータである。測定データがリストモードデータ(時系列データ)の場合には、ポアソン分布の重み付き尤度関数L(x)は下記(4)式で表される。
Figure JPOXMLDOC01-appb-M000004
 ここで、Nはイベント数(リスト数)であり、i(n)は、n番目のイベントを計測した測定データ点の番号(1≦i(n)≦N)である。このように設定された重み係数を、上記(4)式に適用する。上記(4)式の重み付き尤度関数は、上述した実施例1、2と同様に図3に示す反復計算アルゴリズム(逐次近似法)により最大化することができ、再構成処理後のデジタル画像を取得することができる。
 図3に示すステップS1~S6は、上述した実施例1、2のステップS1~S6と同じであるので、その説明については省略する。ただし、ステップS4(更新画像の計算)に関しては、データ形式がリストモードデータの場合には、下記(5)式を利用する。
Figure JPOXMLDOC01-appb-M000005
 以上のように、ステップS1~S6は、この発明における再構成処理工程に相当する。
 本実施例3に係る画像再構成処理方法の作用・効果についても、上述した実施例1、2に係る画像再構成処理方法の作用・効果と同じであるので、その説明については省略する。なお、上記(4)式中の重み係数については、上述した実施例1のように、重み係数を要素データに依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定してもよいし、上述した実施例2のように、DOI検出器の検出素子の検出深さ位置情報に基づいて重み係数を設定してもよい。
 この発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
 (1)上述した各実施例では、放射線検出装置として、陽電子放射断層撮影装置(PET装置)を例に採って説明したが、放射線の検出に基づいて被検体の計測データ集合を取得する装置であれば、特に限定されない。単一光子放射断層撮影装置(SPECT装置)やX線コンピュータ断層撮影装置(X線CT装置)などに適用してもよい。
 (2)上述した各実施例において撮影対象については、特に限定されない。特に、上述した各実施例の場合には、被検体の全身を撮影する装置や、被検体の頭部を撮影する装置や、被検体の乳房を撮影する装置に適用してもよい。
 (3)上述した各実施例では、図1のような部分リング型PET装置1であったが、人体Mの乳房を撮影するマンモグラフィ装置のように、放射線検出器を対向配置した装置にも適用してもよい。図1の検出器ユニット2A,2Bが、図9に示すように、乳房検査部2Cに置き換わるのを除けば、図1と同じ構成を有する。なお、図9の場合には、乳房検査部2Cは切り欠きとなっており、この切り欠きに脇で挟むことで乳房を検査する。また、γ線検出器3(図9では図示省略)は、この切り欠きに合わせて乳房検査部2C内に複数に並設されている。
 (4)上述した各実施例では、DOI検出器であったが、深さ方向を弁別しない放射線検出器に適用してもよい。特に、上述した実施例1のように、重み係数を要素データに依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて重み係数を設定する場合には、上述した実施例2のようなDOI検出器の検出素子の検出深さ位置情報を用いずに重み係数を設定することができる。
 (5)上述した各実施例では、放射線検出装置が陽電子放射断層撮影装置(PET装置)である場合において、一部が開口して構成された検出器ユニットを備えたPET装置(図1では部分リング型PET装置1)であったが、必ずしも一部が開口して構成された検出器ユニットでなくとも、通常のフルリング型PET装置に適用してもよい。特に、上述した実施例2のように、DOI検出器の検出素子の検出深さ位置情報に基づいて重み係数を設定する場合には、上述した実施例1のような直線状ノイズの指向性を用いずに重み係数を設定することができる。
 (6)上述した各実施例では、放射線検出装置が陽電子放射断層撮影装置(PET装置)である場合において、互いに分離した2つの検出器ユニット2A,2B(図1,図4および図5を参照)を備えたPET装置であったが、必ずしも互いに分離した複数の検出器ユニットでなくてもよい。実施例1でも述べたように、図10(a)に示すように一部が開口して構成された1つの検出器ユニット2Dでもよい。
 (7)上述した各実施例では、放射線検出装置が陽電子放射断層撮影装置(PET装置)である場合において、互いに分離した2つの検出器ユニット2A,2B(図1,図4および図5を参照)を備えたPET装置であったが、2つに限定されない。互いに分離した複数の検出器ユニットであれば、例えば、図10(b)に示すように互いに分離した4つの検出器ユニット2E,2F,2G,2Hなどのように互いに分離した3つ以上の検出器ユニットであってもよい。
 (8)上述した各実施例では、第1多変数関数は、放射線検出装置(各実施例では部分リング型PET装置1)によって得られた被検体Mの計測データ集合(測定データ)を構成する要素データ(各実施例ではi番目の測定データ点)の誤差分布に基づいて構成された部分関数の和で表されるデータ関数(各実施例では尤度関数)であったが、第1多変数関数は、上述したデータ関数(尤度関数)と、物理量分布の事前情報に基づいて構成された第2多変数関数との和であってもよい。各実施例のように陽電子放射断層撮影装置(PET装置)を用いた場合には、上記(2)式または上記(4)式の重み付き尤度関数L(x)と、物理量分布の事前情報に基づいて構成された第2多変数関数からなる他の関数U(x)との和からなるL(x)+U(x)の最大化による画像再構成法も考えられる。L(x)が統計的性質に基づいて導出される関数であるのに対して、U(x)は撮影対象の事前情報(再構成画像xが有するであろう決定的論的性質)に基づいて定義される関数である。U(x)は一般的に「ペナルティ関数(Penalty Function)」と呼ばれる。ペナルティ関数を利用した場合の逐次近似計算式の例を以下に示す(下記(6)式を参照)。
Figure JPOXMLDOC01-appb-M000006
 c(k) はk回目の推定解x(k)の近傍における、正則化関数U(x)の近似曲率画像のj番目の画素値である。
 (9)上述した各実施例では、計測データ集合(測定データ)を構成する要素データ(i番目の測定データ点)の誤差分布は、ポアソン分布であったが、ポアソン分布に限定されない。ガウス分布を用いてもよい。誤差分布がガウス分布の場合にはX線コンピュータ断層撮影装置(X線CT装置)に用いられる。
 以上のように、この発明は、陽電子放射断層撮影装置(PET装置)や単一光子放射断層撮影装置(SPECT装置)やX線コンピュータ断層撮影装置(X線CT装置)などの放射線検出装置を有した断層画像撮影装置全般の画像再構成技術に適している。
 1 … 部分リング型PET装置
 2A,2B … 検出器ユニット
 3 … γ線検出器
 L(x) … 重み付き尤度関数
 w … 重み係数(尤度重み係数)
 SA … ストリークアーティファクト
 h,g … 離散変数
 U(x) … ペナルティ関数
 M … 被検体

Claims (9)

  1.  放射線検出装置によって得られた被検体の計測データ集合から、前記計測データ集合の発生要因に関連する前記被検体の物理量分布を複数次元のデジタル画像として再構成する再構成処理を行う画像再構成処理方法であって、
     前記デジタル画像を未知数とする第1多変数関数は、
     (1)前記計測データ集合を構成する要素データの誤差分布に基づいて構成された部分関数の和で表されるデータ関数,
     または(2)前記計測データ集合を構成する要素データの誤差分布に基づいて構成された部分関数の和で表されるデータ関数と、前記物理量分布の事前情報に基づいて構成された第2多変数関数との和
     であり、
     前記画像再構成処理方法は、
     前記部分関数に対応する要素データの再構成画像に対する逆投影の重み係数で部分関数を重み付けて、重み付けされた前記データ関数,または前記データ関数と、前記物理量分布の事前情報に基づいて構成された第2多変数関数との和からなる前記第1多変数関数の最適化計算に基づいて前記再構成処理を行う再構成処理工程
     を備えることを特徴とする画像再構成処理方法。
  2.  請求項1に記載の画像再構成処理方法において、
     前記放射線検出装置は、陽電子放射断層撮影装置,単一光子放射断層撮影装置,X線コンピュータ断層撮影装置のいずれかであることを特徴とする画像再構成処理方法。
  3.  請求項2に記載の画像再構成処理方法において、
     前記重み係数は、前記重み係数を前記要素データに依存しない定数とした場合に再構成画像に生じる直線状ノイズの指向性に基づいて設定されていることを特徴とする画像再構成処理方法。
  4.  請求項3に記載の画像再構成処理方法において、
     前記直線状ノイズの走行方向に沿った要素データに対する重み係数は、前記直線状ノイズの走行方向に沿っていない要素データに対する重み係数よりも小さいことを特徴とする画像再構成処理方法。
  5.  請求項2に記載の画像再構成処理方法において、
     前記放射線検出装置を構成する放射線検出器は、互いに分離した複数の検出器ユニットで構成されており、
     同一の検出器ユニット内の検出素子を結んだ直線方向に沿った要素データに対する重み係数は、互いに異なる検出器ユニット内の検出素子を結んだ直線方向に沿った要素データに対する重み係数よりも小さいことを特徴とする画像再構成処理方法。
  6.  請求項1に記載の画像再構成処理方法において、
     前記放射線検出装置は、陽電子放射断層撮影装置,単一光子放射断層撮影装置のいずれかであり、
     前記計測データ集合は、サイノグラムデータ,ヒストグラムデータ,リストモードデータのいずれかであることを特徴とする画像再構成処理方法。
  7.  請求項6に記載の画像再構成処理方法において、
     前記放射線検出装置を構成する放射線検出器は、放射線の検出深さ位置情報を計測するように構成されており、
     前記重み係数は、前記計測データ集合に対応した前記検出深さ位置情報に依存することを特徴とする画像再構成処理方法。
  8.  請求項7に記載の画像再構成処理方法において、
     Nを2以上の自然数としたときに、前記放射線検出器はN段の検出深さ位置情報を計測するように構成されており、
     前記放射線検出器を構成する検出素子であって、浅い段から深い段に向かうにしたがって番号が大きくなるように、同時計数を計測した2つの検出素子の検出深さの段数番号をそれぞれg,h(1≦g,h≦N)としたときに、前記重み係数は前記段数番号g,hを離散変数とする2次元関数であり、前記2次元関数は一方の変数を固定したときに得られる他方の変数に対する1次元関数が非増加関数であることを特徴とする画像再構成処理方法。
  9.  請求項1から請求項8のいずれかに記載の画像再構成処理方法において、
     前記誤差分布は、ポアソン分布,ガウス分布のいずれかであることを特徴とする画像再構成処理方法。
PCT/JP2014/067976 2014-07-04 2014-07-04 画像再構成処理方法 WO2016002084A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2016530794A JP6256608B2 (ja) 2014-07-04 2014-07-04 画像再構成処理方法
CN201480080398.2A CN106471392B (zh) 2014-07-04 2014-07-04 图像重构处理方法
PCT/JP2014/067976 WO2016002084A1 (ja) 2014-07-04 2014-07-04 画像再構成処理方法
US15/320,244 US10304218B2 (en) 2014-07-04 2014-07-04 Image reconstruction processing method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2014/067976 WO2016002084A1 (ja) 2014-07-04 2014-07-04 画像再構成処理方法

Publications (1)

Publication Number Publication Date
WO2016002084A1 true WO2016002084A1 (ja) 2016-01-07

Family

ID=55018680

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/067976 WO2016002084A1 (ja) 2014-07-04 2014-07-04 画像再構成処理方法

Country Status (4)

Country Link
US (1) US10304218B2 (ja)
JP (1) JP6256608B2 (ja)
CN (1) CN106471392B (ja)
WO (1) WO2016002084A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108550172A (zh) * 2018-03-07 2018-09-18 浙江大学 一种基于非局部特性和全变分联合约束的pet图像重建方法

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107970037B (zh) * 2016-10-24 2021-01-05 上海东软医疗科技有限公司 成像方法和成像系统
JP6842694B2 (ja) * 2017-02-20 2021-03-17 国立研究開発法人量子科学技術研究開発機構 部分リングpet装置及びpet装置
WO2019033390A1 (en) * 2017-08-18 2019-02-21 Shenzhen United Imaging Healthcare Co., Ltd. SYSTEM AND METHOD FOR IMAGE RECONSTRUCTION
CN107909647B (zh) * 2017-11-22 2020-09-15 长春理工大学 基于空间复用的真实感虚拟3d场景光场投影图像绘制方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61225640A (ja) * 1985-03-30 1986-10-07 Toshiba Corp 断層撮影装置
JP2008245695A (ja) * 2007-03-29 2008-10-16 Natl Inst Of Radiological Sciences 断層撮影装置の画像再構成方法、故障診断方法、断層撮影装置、及び、システムマトリクスの管理プログラム
JP4660706B2 (ja) * 2007-04-23 2011-03-30 独立行政法人放射線医学総合研究所 エネルギーと位置情報を利用した放射線検出方法及び装置

Family Cites Families (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7035467B2 (en) * 2002-01-09 2006-04-25 Eastman Kodak Company Method and system for processing images for themed imaging services
DE10251448A1 (de) * 2002-11-05 2004-05-19 Siemens Ag Verfahren für die Computertomographie eines periodisch sich bewegenden Untersuchungsobjektes, sowie ein CT-Gerät zur Durchführung dieses Verfahrens
JP3930493B2 (ja) * 2004-05-17 2007-06-13 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像処理方法、画像処理装置およびx線ct装置
CN100403987C (zh) * 2005-08-19 2008-07-23 清华大学 一种基于遗传算法的x射线ct多相流检测方法
JP4611168B2 (ja) * 2005-10-07 2011-01-12 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像再構成方法、およびx線ct装置
CN101152104A (zh) * 2006-04-21 2008-04-02 美国西门子医疗解决公司 半自动主动脉瘤分析用的系统和方法
US8090182B2 (en) * 2006-11-13 2012-01-03 National University Corporation Kyoto Institute Of Technology Image reconstruction device, image reconstruction method, image reconstruction program, and CT apparatus
EP2149284B1 (en) * 2007-05-31 2014-08-13 General Electric Company Methods and systems to facilitate correcting gain fluctuations in image
CN100571637C (zh) * 2008-05-28 2009-12-23 华中科技大学 动态模型指导下的血管造影三维重建方法
US20120001077A1 (en) * 2009-03-25 2012-01-05 Yoshihiro Inoue Radiation tomography apparatus
JP5622487B2 (ja) * 2009-09-14 2014-11-12 株式会社東芝 放射線診断装置および画像再構成方法
US8731266B2 (en) * 2009-12-17 2014-05-20 General Electric Company Method and system for correcting artifacts in image reconstruction
JP5707745B2 (ja) * 2010-06-08 2015-04-30 ソニー株式会社 画像安定化装置、画像安定化方法、及びプログラム
JP2012061196A (ja) * 2010-09-17 2012-03-29 Fujifilm Corp 断層画像表示方法および装置
US9123156B2 (en) * 2010-10-14 2015-09-01 Hitachi Medical Corporation X-ray CT apparatus and image reconstruction method
US8233586B1 (en) * 2011-02-17 2012-07-31 Franz Edward Boas Iterative reduction of artifacts in computed tomography images using forward projection and an edge-preserving blur filter
DE102011006579A1 (de) * 2011-03-31 2012-10-04 Siemens Aktiengesellschaft Verfahren zur Erzeugung von Bilddaten eines Untersuchungsobjekts, Projektionsdatenverarbeitungseinrichtung, Röntgensystem und Computerprogramm
JP6125148B2 (ja) * 2012-03-14 2017-05-10 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像生成方法、画像生成装置および放射線断層撮影装置並びにプログラム
KR101932720B1 (ko) * 2012-05-07 2018-12-26 삼성전자주식회사 이미지를 복원하는 방법, 이를 수행하는 장치 및 토모그래피 장치
WO2014034618A1 (ja) * 2012-08-30 2014-03-06 株式会社東芝 医用画像処理装置及びx線コンピュータ断層撮影装置
US9486178B2 (en) * 2012-08-31 2016-11-08 Shimadzu Corporation Radiation tomographic image generating apparatus, and radiation tomographic image generating method

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS61225640A (ja) * 1985-03-30 1986-10-07 Toshiba Corp 断層撮影装置
JP2008245695A (ja) * 2007-03-29 2008-10-16 Natl Inst Of Radiological Sciences 断層撮影装置の画像再構成方法、故障診断方法、断層撮影装置、及び、システムマトリクスの管理プログラム
JP4660706B2 (ja) * 2007-04-23 2011-03-30 独立行政法人放射線医学総合研究所 エネルギーと位置情報を利用した放射線検出方法及び装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108550172A (zh) * 2018-03-07 2018-09-18 浙江大学 一种基于非局部特性和全变分联合约束的pet图像重建方法
CN108550172B (zh) * 2018-03-07 2020-05-19 浙江大学 一种基于非局部特性和全变分联合约束的pet图像重建方法

Also Published As

Publication number Publication date
CN106471392A (zh) 2017-03-01
JP6256608B2 (ja) 2018-01-10
US20170154444A1 (en) 2017-06-01
CN106471392B (zh) 2019-02-15
US10304218B2 (en) 2019-05-28
JPWO2016002084A1 (ja) 2017-04-27

Similar Documents

Publication Publication Date Title
CN105408940B (zh) 混合式(谱/非谱)成像探测器阵列和对应的处理电子设备
JP6256608B2 (ja) 画像再構成処理方法
EP3207406B1 (en) Pet detector scintillator arrangement with light sharing and depth of interaction estimation
US20090078876A1 (en) Method and system for using tissue-scattered coincidence photons for imaging
CN107684436B (zh) 在pet系统中确定入射光子的位置的方法
Mathews et al. Improving PET imaging for breast cancer using virtual pinhole PET half-ring insert
Pourmoghaddas et al. Scatter correction improves concordance in SPECT MPI with a dedicated cardiac SPECT solid-state camera
US8975587B2 (en) Positron CT apparatus and a reconstruction method
WO2014121072A2 (en) A spectral response effects (sre) compensation method for photon counting detectors (pcds)
WO2018163362A1 (ja) 散乱推定方法、散乱推定プログラム並びにそれを搭載したポジトロンct装置
JP2020060545A (ja) 陽電子放出撮像装置及び方法
JP7055606B2 (ja) 医用画像処理装置、医用画像診断装置及び医用画像処理方法
Efthimiou et al. Data-driven, energy-based method for estimation of scattered events in positron emission tomography
JP6761610B2 (ja) 吸収係数画像推定方法、吸収係数画像推定プログラム並びにそれを搭載したポジトロンct装置
Yang et al. Effect of CZT system characteristics on Compton scatter event recovery
US20230218243A1 (en) Medical image processing device, computer program, and nuclear medicine device
JP7302737B2 (ja) リストモード画像再構成方法および核医学診断装置
JP2015072125A (ja) 画像処理方法、画像再構成処理方法および断層画像撮影装置
Lee et al. A cubic gamma camera with an active collimator
US10354417B2 (en) Medical image processing apparatus and medical image diagnosis apparatus and medical image processing method
US20220319068A1 (en) Nuclear medicine diagnosis apparatus and nuclear medicine diagnosis method
JP2011153975A (ja) 核医学診断装置
JP2016018245A (ja) 画像再構成処理方法および画像処理方法
JP6206501B2 (ja) 断層画像処理方法およびそれにおける各工程を実行するための手段を備えた放射線放出型断層撮影装置
DiFilippo Instrumentation and principles of imaging: PET

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: 14896510

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2016530794

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 15320244

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14896510

Country of ref document: EP

Kind code of ref document: A1