CN110428384B - Method for acquiring correction information for attenuation correction of PET images of respiration or heart - Google Patents
Method for acquiring correction information for attenuation correction of PET images of respiration or heart Download PDFInfo
- Publication number
- CN110428384B CN110428384B CN201910731364.XA CN201910731364A CN110428384B CN 110428384 B CN110428384 B CN 110428384B CN 201910731364 A CN201910731364 A CN 201910731364A CN 110428384 B CN110428384 B CN 110428384B
- Authority
- CN
- China
- Prior art keywords
- pet
- phase
- image
- attenuation coefficient
- distribution
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 89
- 238000012937 correction Methods 0.000 title claims abstract description 55
- 230000029058 respiratory gaseous exchange Effects 0.000 title claims description 18
- 238000009826 distribution Methods 0.000 claims abstract description 167
- 239000011159 matrix material Substances 0.000 claims abstract description 85
- 238000001514 detection method Methods 0.000 claims abstract description 60
- 230000009466 transformation Effects 0.000 claims abstract description 46
- 230000033001 locomotion Effects 0.000 claims abstract description 39
- 230000000737 periodic effect Effects 0.000 claims abstract description 10
- 238000005457 optimization Methods 0.000 claims abstract description 9
- 238000013501 data transformation Methods 0.000 claims abstract description 4
- 230000000694 effects Effects 0.000 claims description 55
- 230000008569 process Effects 0.000 claims description 22
- 238000012545 processing Methods 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 8
- 238000012804 iterative process Methods 0.000 claims description 7
- 238000005259 measurement Methods 0.000 claims description 6
- 230000004044 response Effects 0.000 claims description 6
- 238000013507 mapping Methods 0.000 claims description 4
- 239000000523 sample Substances 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 2
- 238000002600 positron emission tomography Methods 0.000 abstract description 200
- 238000003384 imaging method Methods 0.000 description 14
- 210000001519 tissue Anatomy 0.000 description 11
- 238000004364 calculation method Methods 0.000 description 10
- 230000000747 cardiac effect Effects 0.000 description 7
- 230000000241 respiratory effect Effects 0.000 description 7
- 210000000056 organ Anatomy 0.000 description 6
- 238000003745 diagnosis Methods 0.000 description 5
- 230000003902 lesion Effects 0.000 description 5
- 238000012636 positron electron tomography Methods 0.000 description 5
- 210000004072 lung Anatomy 0.000 description 4
- 239000002184 metal Substances 0.000 description 4
- 210000004872 soft tissue Anatomy 0.000 description 4
- 238000013473 artificial intelligence Methods 0.000 description 3
- 230000002503 metabolic effect Effects 0.000 description 3
- 238000004445 quantitative analysis Methods 0.000 description 3
- 206010028980 Neoplasm Diseases 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 210000003484 anatomy Anatomy 0.000 description 2
- 238000002591 computed tomography Methods 0.000 description 2
- 238000013170 computed tomography imaging Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 230000004927 fusion Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 208000014644 Brain disease Diseases 0.000 description 1
- 238000012879 PET imaging Methods 0.000 description 1
- 210000001015 abdomen Anatomy 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000010247 heart contraction Effects 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 230000010349 pulsation Effects 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 229940121896 radiopharmaceutical Drugs 0.000 description 1
- 239000012217 radiopharmaceutical Substances 0.000 description 1
- 230000002799 radiopharmaceutical effect Effects 0.000 description 1
- 238000001959 radiotherapy Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 210000003781 tooth socket Anatomy 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/80—Geometric correction
-
- 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/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- 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
-
- 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]
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30048—Heart; Cardiac
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Nuclear Medicine (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention discloses a correction information acquisition method for carrying out attenuation correction on a PET (positron emission tomography) image of a breath or a heart, which comprises the following steps: acquiring PET detection data with flight time information and a gating signal reflecting periodic motion of a breath or a heart in the data; acquiring combined data of each of B phases in a period; obtaining combined data for each phase
And
obtaining a transformation matrix of combined data transformation of any two phases in B phases, and establishing a transformation matrix sum
And a transformation matrix sum
A second relation of (1); establishing an objective function of a PET radioactivity distribution estimated value and a linear attenuation coefficient image which comprise a basic phase; performing iterative optimization on the objective function to obtain x meeting the requirement of maximizing the objective function0And mu0An estimate of (d). All PET data are used in the PET image reconstruction, and the signal to noise ratio of the reconstructed image is improved due to the fact that data of a certain phase position are not used.
Description
Technical Field
The invention relates to the field of medical imaging, in particular to a method and a system for acquiring correction information for attenuation correction of PET (positron emission tomography) images of a breath or a heart in a Positron Emission Tomography (PET) system.
Background
Positron Emission tomography (pet) (positron Emission tomography) is a high-end nuclear medicine image diagnostic device. In practice using radionuclides (e.g. of the type18F、11C, etc.) mark the metabolic substance and inject the nuclide into the human body, and then the PET system is used for carrying out functional metabolic imaging on the patient to reflect the condition of the life metabolic activity. PET scanning is usually matched with other modality scanning (such as CT, MRI, etc.) to obtain anatomical structure imaging of a patient, so that on one hand, nuclide distribution can be accurately positioned, and the accuracy of lesion positioning is improved; on the other hand, the tissue density distribution of the patient can be correspondingly obtained, the attenuation property (attenuation coefficient) of the tissue to the ray is calculated, then the attenuation property is applied to image reconstruction, attenuation correction is carried out on PET data, finally, an image of the actual tissue radioactivity distribution is obtained, and the distribution of the radiopharmaceutical in the patient body is accurately and quantitatively analyzed. The two images are finally fused in the same machine, and the advantages of functional imaging and anatomical imaging are compatible, so that the overall situation of the whole body is understood clearly, the aims of early focus finding and disease diagnosis are fulfilled, and the method has more advantages in guiding diagnosis and treatment of tumors, heart diseases and brain diseases.
In practical applications, the complete modality scan can be completed in a short time. Taking a CT scan as an example, the obtained CT image is a snapshot at almost a certain moment. However, PET scanning is slow and typically takes several minutes per position, making it impossible to complete the data acquisition while the patient is holding his breath. Cardiac activity and respiratory motion can negatively impact PET imaging: on one hand, the PET acquisition can superpose data acquired when the moving focus is at different positions, and the corresponding imaging reflects an average effect of the focus position, which inevitably brings about the reduction of resolution, so that the accuracy of image blurring and quantitative analysis of the SUV value (standardized uptake value) is reduced; on the other hand, under the influence of respiration and cardiac motion, the PET image and other modality images have a certain degree of mismatching on the position and phase of the same lesion, which causes the registration and fusion of the PET image and other modality images to generate deviation (particularly near the diaphragm with the maximum respiratory motion amplitude), and the instant other modality images are used for attenuation correction of the PET image with the average phase, so that errors inevitably occur, local artifacts are generated on the PET image, and the quantitative analysis is influenced.
In order to remove the influence caused by the heart beating or the correction of the respiratory motion, the conventional method uses a gating technique, that is, a motion signal of the heart or the respiration is extracted through an external gating acquisition device or through PET data, the acquired PET data is classified and analyzed according to the period phase by taking the motion signal as a reference, and the PET data in the same phase is reconstructed to obtain a reconstructed image corresponding to each phase. Because the focus is approximately kept still in the same phase interval in the motion period, the gating technology reduces motion artifacts, improves image quality, improves the accuracy of focus positioning and the accuracy of SUV value measurement, can obtain the motion information of the focus three-dimensional space image changing along with time, and has important significance for clinical diagnosis, staging and delineation of a focus radiotherapy target area.
Although the PET reconstructed image obtained by the gating method can effectively remove the influence of respiratory motion and heart pulsation, the method has the following disadvantages in practical application:
1. each phase is reconstructed separately after PET data acquisition, so much less data is available during reconstruction relative to a static scan, which results in a low signal-to-noise ratio of the final gated image. Thus, even if the blurring effect of the object motion on the lesion is removed, the added noise will greatly reduce the detectability of small lesions.
2. Because other modality images correspond to a certain phase position, at most only one phase in the PET multi-phase images obtained by the gating method can be matched with the anatomical structure images, which still causes the deviation (particularly near the diaphragm with the largest motion amplitude) generated by registration and fusion of the PET images and other modality images, and the error necessarily occurs when the average PET image is subjected to attenuation correction by using the instant multi-modality images, so that local artifacts can be generated on the PET images, the quantitative analysis is influenced, and the accurate diagnosis of the tumors in the chest and abdomen and the formulation of a treatment plan can be influenced.
3. The scan range of PET will typically be larger than the scan range of other modalities, such as CT or MRI. Other modality imaging is likely to fail to provide a sufficiently large scan range when scanning a relatively heavy patient, which can result in attenuation images being truncated. The application of such incomplete attenuation information in PET reconstruction also produces attenuation artifacts.
4. When PET is imaged in conjunction with other phantoms, satisfactory attenuation-corrected images, such as PET/MR imaging, are sometimes not obtained. Compared to CT imaging, MR mainly utilizes magnetic spin imaging, not tissue density distribution imaging, and therefore does not directly provide accurate information about tissue attenuation properties. At present, the algorithm for attenuation correction by means of MR imaging is complex in application and low in precision, and attenuation artifacts are easily generated. In addition, the bed and the MR coils cannot be displayed on the MR image, which also has an effect on the subsequent attenuation correction.
5. The gated acquisition method ultimately provides a series of reconstructed images corresponding to different phases within the motion cycle, which may reflect changes in the location of the lesion within the motion cycle. The physician observes the gated images in an animated form, which requires more software functions to be developed than in a static scan, such as: adjusting the animation frame rate, optimizing settings improves animation display quality, outputting animations, etc., which adds additional software development costs. In practical application, doctors usually browse images of all phases in an animation mode, and select a proper frame to be put into a report according to display effect. On one hand, the workload of a doctor is increased, on the other hand, only one frame of image is finally shown in a diagnosis report, other reconstruction results are not used, a large amount of information redundancy exists, and data and calculation time are wasted.
Disclosure of Invention
To solve at least one of the above problems, the present invention provides a correction information acquisition method for attenuation correction of a PET image of a respiration or a heart, a PET activity distribution image reconstruction method, and a PET system.
In order to achieve the purpose, the invention adopts the main technical scheme that:
in a first aspect, the present invention provides a correction information acquisition method for attenuation correction of a PET image of a respiration or a heart, including:
s1, acquiring PET detection data with flight time information during PET system scanning, and acquiring a gating signal which is associated with the PET detection data and reflects periodic motion of a breath or a heart;
s2, acquiring combined data of each phase in B phases in one period based on the gating signals and the PET detection data; the combined data is a data combination of the PET detection data with the same phase in each period;
s3 obtaining PET radioactivity distribution estimation value of combined data of each phase
And linear attenuation coefficient image estimation
S4, obtaining a transformation matrix of combined data transformation of any two phases in the B phases, and establishing the transformation matrix and the transformation matrix
And a conversion matrix and the first relation of
A second relation of (1);
s5, acquiring a PET radioactivity distribution estimated value including a basic phase and an objective function of a linear attenuation coefficient image based on the first relational expression and the second relational expression;
s6, performing iterative optimization on the objective function by adopting an iterative mode to obtain x meeting the requirement of maximizing the objective function0And mu0An estimate of (d).
Optionally, the method further comprises:
s7, according to the x0And mu0By means of a phase transformation matrix, PET radioactivity distribution estimation values and linear attenuation coefficient images of other phases in a period are obtained.
Optionally, the acquiring a gating signal reflecting respiratory or cardiac periodic motion associated with the PET detection data in step S1 includes:
s11, acquiring the gating signals by adopting an external gating device during the acquisition of the PET detection data;
or,
s11a, extracting the gating signal based on the PET detection data.
Optionally, the step S2 includes:
s21, dividing the PET detection data by adopting the gating signals to obtain the PET detection data in each period;
s22, segmenting the PET detection data in each period to obtain B phases of PET detection data in each period, and combining the PET detection data in the same phase of all periods to obtain combined data of each phase;
b is a positive integer greater than 1.
Optionally, the step S3 includes:
performing joint distribution estimation of the PET activity distribution image and the linear attenuation coefficient image on the combined data of each phase to obtain a PET radioactivity distribution estimation value corresponding to each phase
And linear attenuation coefficient image estimation
Optionally, the source image is a base phase 0 phase, and the first relation is:
each pixel in the phase image after non-rigid change is taken as a mapping sum of all pixels in the source image, and the formula (2) is expressed as a formula (3) in a matrix form:
wherein，
Is the jth pixel value, M, of the radioactivity image of the b-th phasexAs a transformation matrix corresponding to the PET radioactivity x,
is a transformation matrix of phase 0 to phase b,
is an identity matrix, J denotes the size of the discrete space of the PET image, b denotes the sequential number of phases,
and
PET radioactivity distribution for phase 0 and phase b; ξ represents another variable, independent of variable J, with a range of values 1-J,
the ξ -th pixel value representing the PET radioactivity image of phase 0;
the second relation is:
equation (4) is expressed in matrix form:
Mμfor the transformation matrix corresponding to the linear attenuation coefficient mu,
is a transformation matrix of phase 0 to phase b,
is an identity matrix; k is the size of the discrete space of the linear attenuation coefficient, eta represents another variable independent of the variable K and has a value range of 1-K,
and an η -th pixel value of the linear attenuation coefficient image representing the 0-th phase.
Optionally, the step S5 includes:
s51, linear attenuation coefficient image estimation based on the first relational expression, conversion matrix and corresponding to each phase
The method is applied to a modeling relational expression in the process of acquiring PET detection data in advance;
i＝1,...,N；t＝1,...,T；b＝1,…,B
s52, obtaining a log-likelihood function of the PET detection data based on the modeling relational expression of the PET detection data obtained based on the PET detection data obeying Poisson distribution;
wherein equation (7) is optimized as:
s53, carrying out iterative solution on the log-likelihood function to obtain PET radioactivity distribution corresponding to a basic phase; in the optimization process, the linear attenuation coefficient distribution mu is firstly maintained0As a constant, activity distribution x for PET0To maximize the objective function; reselecting to maintain PET activity distribution x0Is constant, for an unknown linear attenuation coefficient distribution mu0To maximize the objective function, and alternately perform the maximization operation to obtain x satisfying the requirement of the maximization objective function0And mu0An estimated value of (d);
wherein n is the number of iterations; the detection data is divided into three dimensions of a gating phase b, a flight time t and a sinogram coordinate i; j and xi represent PET radioactivity distribution image coordinates, independent of each other; k represents linear attenuation coefficient image coordinates;
y＝[y1t,y2t,…,yNT]Trepresents the detected data, i.e. photon pairs, N represents the size of the sinogram of the detected data, and T represents the dimension of TOF;
x＝[x1,x2,…,xJ]Trepresenting a PET activity distribution image; j represents the size of discrete space of the PET image;
μ＝[μ1,μ2,…,μK]Texpressing unknown linear attenuation coefficient distribution, and expressing K as the size of a discrete space of a linear attenuation coefficient image; a ═ Aijt]The probability that a spatial position point source j is detected by a response line LOR i and the flight time TOF is t in the PET system is described in a mathematical form as a system matrix, and the physical characteristics of the system are reflected;
l＝[lik]the linear attenuation coefficient matrix represents the track crossing length when the LOR i passes through the space position point source k; r ═ r1t,r2t,…,rNT]TMeans representing random noise and scattering noise; r isbA step of representing the b-th phaseMean value of machine and scattering noise, ybProbe data representing the b-th phase;
Mxas a transformation matrix corresponding to the PET radioactivity x,
is a transformation matrix of phase 0 to phase b,
is an identity matrix, J denotes the size of the discrete space of the PET image, b denotes the sequential number of phases,
and
PET radioactivity distribution for phase 0 and phase b; ξ represents another variable, independent of variable J, with a range of values 1-J,
the ξ -th pixel value representing the PET radioactivity image of phase 0;
Mμfor the transformation matrix corresponding to the linear attenuation coefficient mu,
is a transformation matrix of phase 0 to phase b,
is an identity matrix; k is the size of the linear attenuation coefficient discrete space,
and
is the linear attenuation coefficient distribution of phase 0 and phase b; eta and beta represent another variable independent of the variable K, over a range of 1-K
And an η -th pixel value of the linear attenuation coefficient image representing the 0-th phase.
Optionally, in the iterative process, dividing the PET detection data into a plurality of subsets according to the angle direction of the detection response line, and performing iterative processing on a log-likelihood function for each subset; within the iteration cycle, there is an inner loop for each subset, the inner loop of each subset using only the measurement data contained in the subset;
alternatively, the method further comprises:
s7, according to the x0And mu0By means of a phase transformation matrix, PET radioactivity distribution estimation values and linear attenuation coefficient images of other phases in a period are obtained.
In a second aspect, the present invention also provides a method for reconstructing a PET activity distribution image, including:
acquiring output values of PET radioactivity distribution and linear attenuation coefficient distribution images corresponding to the basic phase by adopting any one of the methods;
applying the output values of the PET activity distribution images and the linear attenuation coefficient distribution images corresponding to the basic phases to PET activity distribution image reconstruction scanned by a PET system to obtain reconstructed images corresponding to the basic phases;
or,
acquiring output values of the PET radioactivity distribution x and the linear attenuation coefficient distribution image u corresponding to any phase by adopting any method;
the method is applied to the reconstruction of the PET activity distribution image scanned by the PET system according to the output value of the PET activity distribution x corresponding to any phase, and a reconstructed image corresponding to a certain phase is obtained, or the obtained reconstructed images of all phases are combined to obtain reconstructed images of all phases.
In a third aspect, the present invention further provides a PET system, which includes an image acquisition and processing device;
the image acquisition and processing device carries out image reconstruction by using the method.
The invention has the beneficial effects that:
the invention provides a reconstruction method for respiratory and cardiac motion correction. First, the acquired data is sorted by phase according to the respiratory or cardiac motion waveform. Secondly, a linear attenuation coefficient image matched with the self phase is iteratively extracted from the PET data with Time Of Flight (TOF) (time Of flight) information at each phase. And finally, using the PET data of all phases and the corresponding attenuation information in the reconstruction process to iteratively obtain a PET reconstruction image. Therefore, on one hand, all PET data are used in the reconstruction process instead of data of a certain phase, the signal to noise ratio of the image is greatly improved compared with the original gating data reconstruction method, and on the other hand, the linear attenuation coefficient image is directly derived from the PET data with flight time information, so that the registration of attenuation correction and the PET image is very ideal, and the attenuation artifact problem of the PET image is effectively solved.
Compared with the traditional gated image, the motion correction method provided by the invention has the following advantages: the signal to noise ratio of the image is greatly improved, the diagnosis of doctors is facilitated, and small focuses are prevented from being annihilated in noise;
because the attenuation correction information in the reconstruction process is from the PET data, when the images of other modalities of PET are not matched due to respiration or heartbeat and patient movement, the attenuation correction can still be correctly carried out on the images, the image quality is improved, and more accurate images are provided for analysis and application of doctors; attenuation correction is directly carried out by using PET acquisition information, and the method is not limited by other imaging modes (such as PET/MRI) which are not easy to extract attenuation information, so that the attenuation correction can be conveniently carried out; the algorithm is applied without the problem of attenuation image truncation, and a doctor can scan a large and heavy patient conveniently.
Drawings
FIG. 1 is a schematic flow diagram of the method of the present invention;
fig. 2a to 2c are schematic diagrams illustrating the comparison between the method of the present invention and the reconstructed image according to the prior art.
Detailed Description
For the purpose of better explaining the present invention and to facilitate understanding, the present invention will be described in detail by way of specific embodiments with reference to the accompanying drawings.
In order to better understand the scheme of the embodiment of the invention, the following outlines the scheme of the embodiment of the invention.
During PET system acquisition, Time Of Flight (TOF) information Of a photon pair is usually acquired, that is, the Time difference between two photons in the photon pair reaching a detector ring is measured, and the approximate position Of an annihilation event on a coincidence attenuation curve is estimated according to the speed Of light.
The flight time information is applied to the reconstruction process of the PET image (also called as PET radioactivity distribution image/PET activity distribution image), so that the signal to noise ratio of the image can be obviously improved, and higher image quality can be obtained. In addition, if the PET image is not matched with the attenuation image, the flight time information is introduced in the reconstruction process, and the attenuation artifact in the PET image can be effectively reduced.
It can be seen that the PET data with time-of-flight information, itself, contains attenuation information. The invention can effectively extract a Linear attenuation coefficient distribution image (Linear attenuation coefficient image) from TOF information of PET data in an iterative manner, generates an attenuation correction parameter (attenuation correction factor) by projecting the Linear attenuation coefficient distribution image to a detection data space, and applies the attenuation correction parameter to PET image reconstruction in real time.
Further, since the linear attenuation coefficient image is directly derived from the PET data with time-of-flight information, the attenuation correction can be accurately matched with the PET image. In the algorithm application process, in order to accurately evaluate the missing part of the truncated attenuation image or eliminate the influence of obvious artifacts in the attenuation image and enable the PET result to be more accurate, the invention provides a reconstruction method for correcting respiration and heart motion. All acquired data is used in the PET reconstruction process and matched attenuation information is extracted from the data with time-of-flight information, thereby improving the signal-to-noise ratio of the image and eliminating attenuation artifacts. .
It should be noted that, in the present application, the PET image may be a partial image including respiration, or the PET image may be a partial image including a heart, or the PET image may include a part or all of respiration or a heart.
In the acquisition process of current PET systems, the PET acquisition process can be modeled as the following equation (1):
in formula (1), y ═ y1t,y2t,…,yNT]TRepresenting detected data, i.e., photon pairs, N represents the size of the sinogram of the detected data (the sinogram is used to characterize the detected data space), and T represents the dimension of TOF.
x＝[x1,x2,…,xJ]TRepresenting an unknown PET image, i.e. a PET activity distribution image; j is expressed as the size of the discrete space of the PET image.
μ＝[μ1,μ2,…,μK]TAnd the unknown linear attenuation coefficient distribution is represented, the dimension of the attenuation coefficient is independent of the flight time, and K is represented as the size of a discrete space of a linear attenuation coefficient image.
A＝[Aijt]For the system matrix, the system matrix can be understood as a mathematical expression of the probability that a spatial position point source j is detected by a line of response (LOR) i in the PET system and TOF is t, which reflects the physical characteristics of the PET system.
l＝[lik]And the linear attenuation coefficient matrix represents the track crossing length of the LOR i when the LOR i passes through the space position point source k.
r＝[r1t,r2t,…,rNT]TMean values of random noise and scattering noise are indicated.
Therefore, in the embodiment of the present invention, the subsequent calculation processing is performed based on the above formula (1), specifically referring to the first embodiment.
Example one
The present application proposes a correction information acquisition method for attenuation correction of PET images of the respiration or heart. Referring to fig. 1, the method comprises the following specific steps:
s1, acquiring PET detection data with flight time information during the scanning of the PET system, and acquiring a gating signal which reflects the periodic movement of the breath or the heart and is related to the PET detection data.
For example, the gating signal may be extracted by an external gating device to reflect the periodic motion of the breath or heart in this embodiment.
In another embodiment, the PET data itself may be utilized to extract the gating signal to reflect the periodic motion of the breath or heart.
S2, acquiring combined data of each phase in B phases in one period based on the gating signals and the PET detection data; the combined data is a data combination of the PET detection data with the same phase in each period;
s3 obtaining PET radioactivity distribution estimation value of combined data of each phase
And linear attenuation coefficient image estimation
b＝1,…B。
In the specific implementation process, the acquired PET data are segmented and recombined into B phases according to cycles, and joint distribution estimation of a PET activity distribution image and a linear attenuation coefficient image is performed on the data of each phase by using the following embodiment to obtain a corresponding PET activity distribution image estimation of each phase
And linear attenuation coefficient image estimation
Since the attenuation information is directly extracted from the PET data, the gated image of each phase has correct attenuation correction, but only approximate 1/B data is used for reconstructing the gated image of each phase, and the gated PET activity distribution image and the linear attenuation coefficient image are very noisy.
S4, obtaining a transformation matrix of combined data transformation of any two phases in the B phases, and establishing the transformation matrix and the transformation matrix
And a conversion matrix and the first relation of
The second relation of (1).
For better understanding, in the present application, the PET activity distribution images corresponding to the gated phases have a correlation with each other, reflecting the periodic variation of radioactivity distribution in the human body affected by respiration or cardiac motion, the gated images of different phases satisfy a non-rigid transformation with each other without loss of generality, and the source image is the base phase 0 phase, then:
each pixel of the non-rigidly transformed phase image can be regarded as the mapping sum of all pixels in the source image. Equation (2) can be written in matrix form:
Mxfor transforming matrices, for calibratingNon-rigidly varying weights of different phase PET images relative to each other.
Is a transformation matrix of phase 0 to phase b, wherein
Is an identity matrix.
Further, the distribution images of the linear attenuation coefficients corresponding to the gating phases have correlation with each other, which reflects the periodic variation of tissue density distribution in the human body influenced by respiration or heart motion, the gating images of different phases satisfy non-rigid transformation with each other, without loss of generality, and the source image is the base phase 0 phase, then:
each pixel of the non-rigidly transformed phase image can be regarded as the mapping sum of all pixels in the source image. Equation (2) can be written in matrix form:
Mμthe transformation matrix is used to calibrate the non-rigid variation weights of the different phase linear attenuation coefficient images with respect to each other.
Is a transformation matrix of phase 0 to phase b, wherein
Is an identity matrix.
It is emphasized that the PET activity distribution image transformation matrix MxAnd linear attenuation coefficient distribution image conversion matrix MμThe base phases of (a) are not required to be the same, but for convenience the same base phases are typically used.
In this embodiment, the transformation matrix is obtained by a registration algorithm, such as demons method, based on images of different phases. The images of different phases can be obtained by using a gating technique in step S2, or can be obtained by a hardware method, such as 4D-CT, so that high-precision images of different phases can be obtained. The selection of the transformation matrix parameters should follow the following principle: namely, each phase image obtained after conversion is closest to the actually acquired image. The transformation matrix solution can therefore be transformed into a feasible extremum solving problem. The calculation method of the transformation matrix is the prior art, and is not described in detail in the application.
S5, acquiring a PET radioactivity distribution estimated value including a basic phase and an objective function of a linear attenuation coefficient image based on the first relational expression and the second relational expression;
s6, performing iterative optimization on the objective function by adopting an iterative mode to obtain x meeting the requirement of maximizing the objective function0And mu0An estimate of (d).
For the above step S5 and step S6, the following (the following substeps, not shown in the figure) are explained:
s51, applying the transformation matrix to equation (1), for phase b:
i＝1,...,N；t＝1,...,T；b＝1,…,B
s52, the PET detection data obey Poisson distribution, unknowns are a PET activity distribution image x and a linear attenuation coefficient distribution image mu, and a log-likelihood function of the detection data is expressed as follows:
the detection data is divided into three dimensions, gated phase b, time of flight t and sinogram coordinates i.
S53, substituting equation (6) into equation (7), ignoring terms that are not related to unknowns, the log-likelihood function can be written as:
as can be seen from the formula (8), the unknown number in the log-likelihood function is the PET activity distribution image of the fundamental phase
And linear attenuation coefficient distribution image
Since the motion is modeled in the function using a transformation matrix, the estimation is performed by optimizing the log-likelihood function
And
the device is not influenced by respiration or heart movement;
and
are estimated from the PET detection data and are perfectly matched to each other, thus eliminating motion artifacts; the whole data is used in the estimation process, so the signal-to-noise ratio of the image quality is greatly improved compared with the traditional gating image.
S54, obtaining unknown number x0And mu0The estimate of (2) needs to be optimized for its jointly distributed objective function (8). Due to the log-likelihood function described above for the unknown x0And mu0The method is a very complex function, and the formula (8) is difficult to directly obtain an analytic solution, so that an iterative algorithm is required to gradually approximate an optimal solution so as to maximize a log-likelihood function. For the feasibility of algorithm implementation, the optimization process needs to be simplified: firstly, useMaintaining a linear attenuation coefficient distribution mu0As a constant, activity distribution x for PET0To maximize the objective function; reselecting to maintain PET activity distribution x0Is constant, for an unknown linear attenuation coefficient distribution mu0The maximum objective function is realized, so that the maximum operation is alternately performed, attenuation correction is continuously corrected to approximate to the real attenuation condition, and finally, x meeting the requirement of the maximum objective function is obtained0And mu0An estimate of (d).
S55, keeping linear attenuation coefficient distribution mu0As a constant, activity distribution x for PET0To maximize the objective function, the iterative calculation formula is:
and n is the iteration number.
S56, maintaining PET activity distribution x0Is constant, for a linear attenuation coefficient distribution mu0To maximize the objective function, the iterative calculation formula is:
and n is the iteration number.
S57, alternately iterating the step S55 and the step S56 to ensure the PET activity distribution image x0And linear attenuation coefficient image mu0The estimated values of (c) converge.
In particular, in the iterative process, the measurement data may be divided into several subsets according to the angular direction of the probe response Line (LOR), and the subset division follows the principle that the angular distribution is uniformly symmetrical. When the iterative algorithm is implemented, an inner loop aiming at the subset is added within the iteration loop, and only the measurement data contained in the corresponding subset is used in each calculation. The entire measurement data is used in its entirety at the end of the subset cycle. The method for dividing according to the subsets can greatly improve the calculation efficiency of the algorithm on the premise of not increasing the calculation amount of the algorithm.
In practical application, the original log-likelihood function (8) can be adjusted by carrying out priori knowledge constraint on the linear attenuation coefficient to generate a new objective function with a priori knowledge item. However, the optimization idea is not greatly different from the actual algorithm, and the redundancy is not repeated for the moment, so that the accuracy of the final quantification is ensured.
Due to the influence of breathing or heart movement, the non-rigid movement of organs in the human body occurs near the organs, so that the conversion matrix for describing the non-rigid movement between PET phase images and between linear attenuation coefficient phase images is very sparse, and the calculation time required for converting the phase images is negligible.
In the implementation process of the invention, the transformation matrix can be continuously updated to be closer to the true non-rigid deformation by the iterative updating mode of the step S5, the noise influence is reduced, and the PET activity distribution reconstruction image and the linear attenuation coefficient reconstruction image corresponding to the basic phase are estimated.
For better illustration, the method of the present application, the method further includes the following step S7:
s7: after the basic phase PET activity distribution reconstructed image is obtained, the PET activity distribution reconstructed images of other phases can be obtained through the phase transformation matrix by using the formula (2).
Example two
The embodiment of the invention provides a correction information acquisition method for performing attenuation correction on a PET activity distribution image, which comprises the following steps:
q0, acquiring PET detection data with time-of-flight information and other modality images while the PET system is scanning.
For example, other modality images may include: CT images or MR images.
Q1, modeling the PET detection data based on the PET detection data obeying Poisson distribution to obtain a log-likelihood function L (x, mu, y) of a formula (p 1);
formula (II)
Wherein y ═ y1t,y2t,…,yNT]TRepresenting the detection data, N representing the size of the sinogram of the detection data, and T representing the dimension of the time of flight TOF; x ═ x1,x2,…,xJ]TRepresenting the distribution of unknown PET radioactivity, J is expressed as the size of discrete space of the PET image; mu ═ mu1,μ2,.....,μK]TExpressing unknown linear attenuation coefficient distribution, and expressing K as the size of a discrete space of a linear attenuation coefficient image; a ═ Aijt]Is a system matrix; l ═ lik]Is a linear attenuation coefficient matrix, r ═ r1t,r2t,…,rNT]TMeans representing random noise and scattering noise;
is an attenuation correction parameter;
q2, obtaining the linear attenuation coefficient distribution image mu according to the other mode images0。
For example: when other mode images are CT images, the data of the CT images are converted into photon linear attenuation coefficient distribution images under 511KeV energy by a bilinear method to obtain linear attenuation coefficient distribution images mu0；
When the other modality image is an MR image, the linear attenuation coefficient distribution image mu0Is a theoretical linear attenuation coefficient value directly assigned according to prior knowledge.
Taking the PET/MR imaging system as an example, the MR image is segmented for different regions (such as soft tissue, fat, lung, air, etc.), and then corresponding theoretical linear attenuation coefficient values are directly given (such as selecting soft tissue region to assign 0.0975cm value)-1The fat region was assigned a value of 0.0864cm-1Assigned lung region of 0.0224cm-1The air zone is assigned a value of 0).
Q3, if artifact or truncation exists in the linear attenuation coefficient distribution image or the linear attenuation coefficient distribution image is not matched with the PET image, acquiring a complete linear attenuation coefficient distribution R (mu);
wherein the PET image is generated directly from the PET detection data.
It can be understood that if the attenuation coefficient distribution image has no obvious artifact and is matched with the PET image, the linear attenuation coefficient distribution can be directly applied for attenuation correction, and if the linear attenuation coefficient distribution image has obvious artifact, truncation or mismatch, the complete linear attenuation coefficient distribution R (μ) needs to be obtained.
In particular, the region of the linear attenuation coefficient distribution image in the CT image or the MR image, which is free from the artifact region and matches the PET image, is automatically determined using a threshold value, or is determined using an artificial intelligence recognition technique.
In this embodiment, R (μ) ═ gW + μ (E-S) ═ μ0SW + μ (E-S); formula (p6)
Wherein E represents a K-order unit matrix, and W is a weight value matrix of prior attenuation coefficient distribution and is used for adjusting the weight of the prior attenuation coefficient distribution in the iterative process; g is a defined region of linear attenuation coefficient distribution without artifacts, g ═ mu0S；
And S is a mask matrix of a region which is matched with the PET image and has no artifact in the linear attenuation coefficient distribution image determined according to the prior condition.
It will be appreciated that,
the K-order diagonal matrix W is:
the weights of the required prior attenuation coefficient distributions are different for different tissue organs or regions wkThe assignments are different. W is the same weight of the required a priori attenuation coefficient distribution for different tissue organs or regionskFor equal weight values or using scalar quantitiesThe numerical values indicated.
Q4 image mu based on linear attenuation coefficient distribution0Optimizing a log-likelihood function L (x, mu, y) by adopting an iterative algorithm, updating a mu value of each iteration by adopting R (mu) in the optimization process, and acquiring an estimated value of x and mu serving as correction information when an iteration termination condition is met;
the R (mu) is a mask matrix and mu which are determined according to prior conditions and are free of artifact areas in the linear attenuation coefficient distribution image and matched with the PET image0μ for each iteration is determined for updating μ for each iteration.
In the present embodiment, in order to better understand the step Q4, the step Q4 is explained as follows.
Q41, holding initial value μ0Fixing, and optimizing a log-likelihood function L (x, mu, y) by using an MLEM iterative reconstruction algorithm, namely obtaining a first estimated value of an unknown number x by the following formula (p 4);
q42, keeping the first estimation value of x fixed, optimizing a log-likelihood function L (x, mu, y) by using an MLTR algorithm, namely obtaining the first estimation value of an unknown number mu through the following formula (p 5);
q43, obtaining R (mu) by adopting a formula (p6) according to the first estimation value of the unknown mu, and updating the first estimation value of mu by adopting R (mu) to obtain an updated mu value;
q44, based on the updated μ value, repeating the above-mentioned processes from step Q41 to step Q43, and when the iteration termination condition is satisfied, taking the final unknown x and the estimated value of μ as the final output values.
In practical application, in the implementation process of the invention, the linear attenuation coefficient distribution μ is kept as a constant, the PET activity distribution x is obtained by adopting the MLEM iterative reconstruction algorithm to maximize the objective function (such as the formula p4), then the PET activity distribution x is selected and kept as a constant, the objective function is maximized for the unknown linear attenuation coefficient distribution μ (such as the formula p5), and the complete artifact-free linear attenuation coefficient distribution R (μ) is calculated according to the obtained μ. And the operation is performed alternately, attenuation correction is continuously corrected to approximate to the real attenuation condition, and finally the estimated values of x and mu meeting the requirement of the maximum objective function are obtained.
The method of the embodiment extracts the characteristic tissue of the object in the iterative process, and introduces the prior knowledge to adjust the iterative process, so that the iterative result approaches to an ideal value, and the completeness and accuracy of the final attenuation image are ensured.
For a better understanding of the above Q2, the process of Q2 is described in detail below.
Since PET scanning is usually used in matching with other modality imaging, the linear attenuation coefficient distribution calculated based on other modality images is defined as μ0(μ0As well as the initial linear attenuation coefficient distribution of the MLEM algorithm).
Taking a PET/CT imaging system as an example, a high signal-to-noise ratio image obtained by the CT system can be utilized to convert a CT value into a photon linear attenuation coefficient distribution image under 511KeV energy through a bilinear method, and mu is obtained at the moment0；
Taking the PET/MR imaging system as an example, the MR image is segmented for different regions (such as soft tissue, fat, lung, air, etc.), and then corresponding theoretical linear attenuation coefficient values are directly given (such as selecting soft tissue region to assign 0.0975cm value)-1The fat region was assigned a value of 0.0864cm-1Assigned lung region of 0.0224cm-1The air region is assigned a value of 0), in which case mu is obtained0。
If the attenuation coefficient distribution image has no obvious artifact and is matched with the PET image, the linear attenuation coefficient distribution can be directly applied to attenuation correction, and when the linear attenuation coefficient distribution image has obvious artifact, truncation or mismatching, an artifact-free region matched with the PET image is selected for different artifacts and the region mask matrix is generated
The definition is as follows:
for the selection of the region with obvious artifacts and matching with the PET, the region can be directly and manually drawn, or different selection drawing methods such as automatic threshold drawing, artificial intelligence identification and the like can be used (for example, a metal artifact region is drawn by depending on a threshold in a CT image, a tissue organ with artifacts is extracted by depending on an image segmentation technology in an MR image, and the like).
For the above R (μ), R (μ) ═ gW + μ (E-S) ═ μ0SW+μ(E-S) (p6)
Determining a region of the linear attenuation coefficient distribution without artifacts as g, where g ═ mu0S, the distribution of this region (region on other modality image) has no obvious artifact and matches with the PET image, but the integrity of the linear attenuation coefficient distribution cannot be guaranteed, so that the linear attenuation coefficient distribution g cannot be directly utilized for attenuation correction. A complete linear attenuation coefficient distribution R (μ) needs to be obtained by a formula (p6), that is, the attenuation coefficient distribution μ in the iterative computation process and the prior attenuation coefficient distribution g are weighted, so that the missing part of the incomplete attenuation coefficient distribution g is obtained by μ weighting.
E represents an identity matrix of order K,
w is a weight value matrix of prior attenuation coefficient distribution and is used for adjusting the weight of the prior attenuation coefficient distribution in the iterative process. The K diagonal matrix W may be defined as:
the weights of the required prior attenuation coefficient distributions are different for different tissue organs or regions wkCan be assigned according to different conditions(ii) a For the same integral weight of the prior attenuation coefficient distribution, w can be setkW may be replaced with a scalar quantity, provided that the weight values are the same.
In this embodiment, the complete linear attenuation coefficient distribution is calculated by not limited to the weighting method of the above formula (p6), and the complete linear attenuation coefficient distribution may be calculated by using different detection fusion methods such as an SIFT/SURF automatic detection and splicing algorithm, an artificial intelligence detection and identification algorithm, and the like.
EXAMPLE III
The invention also provides a PET image reconstruction method, which comprises the following steps:
m01, obtaining output values of the PET radioactivity distribution x and the linear attenuation coefficient distribution image μ corresponding to the fundamental phase by using any one of the correction information obtaining methods described in the above embodiments;
m02, applying the output values of the PET activity distribution x and the linear attenuation coefficient distribution image mu corresponding to the basic phase in the PET activity distribution image reconstruction scanned by the PET system to obtain the reconstructed PET image corresponding to the basic phase.
In another embodiment, the output values of the PET activity distribution x and the linear attenuation coefficient distribution image u corresponding to any phase in the period can be obtained by using any one of the methods described in the above embodiments, and then the method is applied to PET activity distribution image reconstruction scanned by a PET system to obtain a PET image of each phase.
In addition, the PET images of each phase may be combined to obtain a reconstructed PET image of one cycle.
It should be noted that the output values of x and μ obtained in the first embodiment are represented by an array, which is used to represent the values of the pixels in the PET radioactivity distribution and the linear attenuation coefficient distribution.
In practical application, the PET activity distribution image reconstruction method carries out image reconstruction on a single bed, and then PET activity distribution images of the whole scanning space are spliced;
or splicing the PET detection data to be reconstructed of all the beds, and performing image reconstruction on the spliced detection data by adopting a PET activity distribution image reconstruction method to obtain a PET activity distribution image of the whole scanning space.
That is to say, when multiple beds are collected, each bed can be selected to perform attenuation correction calculation and simultaneously reconstruct to obtain a PET activity distribution image of each bed, and then the PET activity distribution images are spliced together; alternatively, the whole detection data can be spliced together, attenuation correction is performed in the whole scanning space at one time, and simultaneously, a PET activity distribution image of the whole scanning space is reconstructed.
The method of the embodiment can be used for the heavy patients or the patients with certain parts (such as arms, hands and the like) of the body exceeding the scanning visual field of other modes in certain special cases, the attenuation images are truncated, the attenuation correction can still be carried out on the PET images, the complete PET multi-mode images are provided, the image quality is improved, and more accurate images are provided for the analysis and application of doctors.
In addition, for the patient with the attenuation image artifact, such as a PET/CT scanning patient with a cardiac pacemaker or a metal tooth socket in the body, the CT image has an obvious metal artifact, so that accurate attenuation correction can be performed, and the influence of the metal artifact is eliminated; when the PET multi-modality images are not matched due to respiration or heartbeat and patient movement, accurate attenuation correction can be performed on the PET images.
For better comparison, fig. 2a is a PET image with motion artifact correction using the algorithm of the present invention. By way of comparison, fig. 2b is a conventional PET image with mismatched CT images to produce blurred motion artifacts, and fig. 2c is a conventional PET image with gating techniques to remove blurred motion artifacts but with a low signal-to-noise ratio.
Furthermore, the invention also provides a PET system, which comprises an image acquisition and processing device;
the image acquisition and processing device performs image reconstruction by using the method described in the third embodiment.
It should be understood that the above description of specific embodiments of the present invention is only for the purpose of illustrating the technical lines and features of the present invention, and is intended to enable those skilled in the art to understand the contents of the present invention and to implement the present invention, but the present invention is not limited to the above specific embodiments. It is intended that all such changes and modifications as fall within the scope of the appended claims be embraced therein.
Claims (9)
1. A correction information acquisition method for attenuation correction of a PET image of a respiration or a heart, characterized by comprising:
s1, acquiring PET detection data with flight time information during PET system scanning, and acquiring a gating signal which is associated with the PET detection data and reflects periodic motion of a breath or a heart;
s2, acquiring combined data of each phase in B phases in one period based on the gating signals and the PET detection data; the combined data is a data combination of the PET detection data with the same phase in each period;
s3 obtaining PET radioactivity distribution estimation value of combined data of each phase
And linear attenuation coefficient image estimation
S4, obtaining a transformation matrix of combined data transformation of any two phases in the B phases, and establishing the transformation matrix and the transformation matrix
And a conversion matrix and the first relation of
A second relation of (1);
s5, acquiring a PET radioactivity distribution estimated value including a basic phase and an objective function of a linear attenuation coefficient image based on the first relational expression and the second relational expression;
s6, adopting an iteration mode to aim at the targetPerforming iterative optimization on the scalar function to obtain x meeting the requirement of the maximized objective function0And mu0An estimated value of (d);
the step S5 includes:
s51, linear attenuation coefficient image estimation based on the first relational expression, conversion matrix and corresponding to each phase
The method is applied to a modeling relational expression in the process of acquiring PET detection data in advance;
s52, obtaining a log-likelihood function of the PET detection data based on the PET detection data obeying Poisson distribution and the modeling relational expression of the PET detection data;
wherein equation (7) is optimized as:
s53, carrying out iterative solution on the log-likelihood function to obtain PET radioactivity distribution corresponding to a basic phase; in the optimization process, the linear attenuation coefficient distribution mu is firstly maintained0As a constant, activity distribution x for PET0To maximize the objective function; reselecting to maintain PET activity distribution x0Is constant, for an unknown linear attenuation coefficient distribution mu0To maximize the objective function, and alternately perform the maximization operation to obtain x satisfying the requirement of the maximization objective function0And mu0An estimated value of (d);
wherein n is the number of iterations; the detection data is divided into three dimensions of a gating phase b, a flight time t and a sinogram coordinate i; j and xi represent PET radioactivity distribution image coordinates, independent of each other; k represents linear attenuation coefficient image coordinates;
y＝[y1t,y2t,…,yNT]Trepresents the detected data, i.e. photon pairs, N represents the size of the sinogram of the detected data, and T represents the dimension of TOF;
x＝[x1,x2,…,xJ]Trepresenting a PET activity distribution image; j represents the size of discrete space of the PET image;
μ＝[μ1,μ2,…,μK]Texpressing unknown linear attenuation coefficient distribution, and expressing K as the size of a discrete space of a linear attenuation coefficient image; a ═ Aijt]The probability that a spatial position point source j is detected by a response line LORi and the time of flight TOF is t in the PET system is described as a system matrix, and the physical characteristics of the system are reflected;
l＝[lik]the linear attenuation coefficient matrix represents the track crossing length of the LORi when the LORi passes through a space position point source k; r ═ r1t,r2t,…,rNT]TMeans representing random noise and scattering noise; r isbMean value of random noise and scattered noise, y, representing the b-th phasebProbe data representing the b-th phase;
Mxas a transformation matrix corresponding to the PET radioactivity x,
is a transformation matrix of phase 0 to phase b,
is an identity matrix, J denotes the size of the discrete space of the PET image, b denotes the sequential number of phases,
and
PET radioactivity distribution for phase 0 and phase b; ξ represents another variable, independent of variable J, with a range of values 1-J,
the ξ -th pixel value representing the PET radioactivity image of phase 0;
Mμfor the transformation matrix corresponding to the linear attenuation coefficient mu,
is a transformation matrix of phase 0 to phase b,
is an identity matrix; k is the size of the linear attenuation coefficient discrete space,
and
is the linear attenuation coefficient distribution of phase 0 and phase b; η and β represent another variable, independent of the variable K, with a value in the range 1-K,
and an η -th pixel value of the linear attenuation coefficient image representing the 0-th phase.
2. The method of claim 1, further comprising:
s7, according to the x0And mu0By means of a phase transformation matrix, PET radioactivity distribution estimation values and linear attenuation coefficient images of other phases in a period are obtained.
3. The method according to claim 1, wherein the acquiring of the gating signal reflecting the periodic motion of the respiration or heart associated with the PET detection data in step S1 comprises:
s11, acquiring the gating signals by adopting an external gating device during the acquisition of the PET detection data;
or,
s11a, extracting the gating signal based on the PET detection data.
4. The method according to claim 1, wherein the step S2 includes:
s21, dividing the PET detection data by adopting the gating signals to obtain the PET detection data in each period;
s22, segmenting the PET detection data in each period to obtain B phases of PET detection data in each period, and combining the PET detection data in the same phase of all periods to obtain combined data of each phase;
b is a positive integer greater than 1.
5. The method according to claim 1, wherein the step S3 includes:
performing joint distribution estimation of the PET activity distribution image and the linear attenuation coefficient image on the combined data of each phase to obtain a PET radioactivity distribution estimation value corresponding to each phase
And linear attenuation coefficient image estimation
6. The method of claim 1,
the source image is a base phase 0 phase, and the first relation is:
each pixel in the phase image after non-rigid change is taken as a mapping sum of all pixels in the source image, and the formula (2) is expressed as a formula (3) in a matrix form:
wherein,
is the jth pixel value, M, of the radioactivity image of the b-th phasexAs a transformation matrix corresponding to the PET radioactivity x,
is a transformation matrix of phase 0 to phase b,
is an identity matrix, J denotes the size of the discrete space of the PET image, b denotes the sequential number of phases,
and
PET radioactivity distribution for phase 0 and phase b; ξ represents another variable, independent of variable J, with a range of values 1-J,
the ξ -th pixel value representing the PET radioactivity image of phase 0;
the second relation is:
equation (4) is expressed in matrix form:
Mμfor the transformation matrix corresponding to the linear attenuation coefficient mu,
is a transformation matrix of phase 0 to phase b,
is an identity matrix; k is the size of the discrete space of the linear attenuation coefficient, eta represents another variable independent of the variable K and has a value range of 1-K,
and an η -th pixel value of the linear attenuation coefficient image representing the 0-th phase.
7. The method according to claim 1, characterized in that in the iterative process, the PET detection data is divided into a plurality of subsets according to the angle direction of the detection response line, and the iterative processing of the log-likelihood function is performed for each subset; within the iteration cycle, there is an inner loop for each subset, the inner loop of each subset being for the measurement data contained by the subset;
alternatively, the method further comprises:
s7, according to the x0And mu0By means of a phase transformation matrix, PET radioactivity distribution estimation values and linear attenuation coefficient images of other phases in a period are obtained.
8. A method of PET activity distribution image reconstruction, comprising:
acquiring output values of PET radioactivity distribution and linear attenuation coefficient distribution images corresponding to a basic phase by adopting the method of any one of the claims 1 and 3 to 7;
applying the output values of the PET activity distribution images and the linear attenuation coefficient distribution images corresponding to the basic phases to PET activity distribution image reconstruction scanned by a PET system to obtain reconstructed images corresponding to the basic phases;
or,
acquiring output values of the PET radioactivity distribution x and the linear attenuation coefficient distribution image u corresponding to any phase by using the method of any one of the claims 2 to 7;
the method is applied to the reconstruction of the PET activity distribution image scanned by the PET system according to the output value of the PET activity distribution x corresponding to any phase, and a reconstructed image corresponding to a certain phase is obtained, or the obtained reconstructed images of all phases are combined to obtain reconstructed images of all phases.
9. A PET system is characterized by comprising an image acquisition and processing device;
the image acquisition processing device performs image reconstruction using the method of claim 8.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910731364.XA CN110428384B (en) | 2019-08-08 | 2019-08-08 | Method for acquiring correction information for attenuation correction of PET images of respiration or heart |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910731364.XA CN110428384B (en) | 2019-08-08 | 2019-08-08 | Method for acquiring correction information for attenuation correction of PET images of respiration or heart |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110428384A CN110428384A (en) | 2019-11-08 |
| CN110428384B true CN110428384B (en) | 2021-11-16 |
Family
ID=68413401
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910731364.XA Active CN110428384B (en) | 2019-08-08 | 2019-08-08 | Method for acquiring correction information for attenuation correction of PET images of respiration or heart |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110428384B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112842271B (en) * | 2021-01-11 | 2021-12-17 | 武汉理工大学 | Physiological signal separation and extraction system and method based on optical fiber sensing |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103282941A (en) * | 2011-01-05 | 2013-09-04 | 皇家飞利浦电子股份有限公司 | Method and apparatus to detect and correct motion in list-ode pet data with a gated signal |
| CN105705097A (en) * | 2013-11-08 | 2016-06-22 | 皇家飞利浦有限公司 | Empirical beam hardening correction for differential phase contrast CT |
| CN106548473A (en) * | 2016-11-07 | 2017-03-29 | 赛诺联合医疗科技（北京）有限公司 | A kind of method and device for building phase image |
| CN106618628A (en) * | 2017-01-24 | 2017-05-10 | 昆明理工大学 | Breathing movement gating correction and attenuation correction method based on PET/CT imaging |
| CN107610198A (en) * | 2017-09-20 | 2018-01-19 | 赛诺联合医疗科技（北京）有限公司 | PET image attenuation correction method and device |
| CN107638188A (en) * | 2017-09-28 | 2018-01-30 | 江苏赛诺格兰医疗科技有限公司 | Image attenuation bearing calibration and device |
| CN109961419A (en) * | 2019-03-26 | 2019-07-02 | 江苏赛诺格兰医疗科技有限公司 | The correction information acquiring method of correction for attenuation is carried out to PET activity distributed image |
| CN109978966A (en) * | 2019-03-21 | 2019-07-05 | 江苏赛诺格兰医疗科技有限公司 | The correction information acquiring method of correction for attenuation is carried out to PET activity distributed image |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP6170549B2 (en) * | 2012-05-21 | 2017-07-26 | コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． | Rapid scatter estimation in PET reconstruction |
-
2019
- 2019-08-08 CN CN201910731364.XA patent/CN110428384B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103282941A (en) * | 2011-01-05 | 2013-09-04 | 皇家飞利浦电子股份有限公司 | Method and apparatus to detect and correct motion in list-ode pet data with a gated signal |
| CN105705097A (en) * | 2013-11-08 | 2016-06-22 | 皇家飞利浦有限公司 | Empirical beam hardening correction for differential phase contrast CT |
| CN106548473A (en) * | 2016-11-07 | 2017-03-29 | 赛诺联合医疗科技（北京）有限公司 | A kind of method and device for building phase image |
| CN106618628A (en) * | 2017-01-24 | 2017-05-10 | 昆明理工大学 | Breathing movement gating correction and attenuation correction method based on PET/CT imaging |
| CN107610198A (en) * | 2017-09-20 | 2018-01-19 | 赛诺联合医疗科技（北京）有限公司 | PET image attenuation correction method and device |
| CN107638188A (en) * | 2017-09-28 | 2018-01-30 | 江苏赛诺格兰医疗科技有限公司 | Image attenuation bearing calibration and device |
| CN109978966A (en) * | 2019-03-21 | 2019-07-05 | 江苏赛诺格兰医疗科技有限公司 | The correction information acquiring method of correction for attenuation is carried out to PET activity distributed image |
| CN109961419A (en) * | 2019-03-26 | 2019-07-02 | 江苏赛诺格兰医疗科技有限公司 | The correction information acquiring method of correction for attenuation is carried out to PET activity distributed image |
Non-Patent Citations (2)
| Title |
|---|
| PET/CT呼吸门控显像对肺部结节的应用研究;武志芳;《中国博士学位论文全文数据库 医药卫生科技辑》;20110815(第8期);E063-6 * |
| Respiratory Motion Compensation for PET/CT with Motion Information Derived from Matched Attenuation Corrected Gated PET Data;Yihuan Lu 等;《Journal of Nuclear Medicine》;20180209;1-26 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110428384A (en) | 2019-11-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11020077B2 (en) | Simultaneous CT-MRI image reconstruction | |
| CN109009200B (en) | System and method for positron emission tomography image reconstruction | |
| US8774481B2 (en) | Atlas-assisted synthetic computed tomography using deformable image registration | |
| Li et al. | Motion correction for improved target localization with on-board cone-beam computed tomography | |
| US8411915B2 (en) | Motion compensation in functional imaging | |
| CN106491151B (en) | PET image acquisition method and system | |
| JP5734664B2 (en) | Image Restoration Method Using Dilution Constraint Correction | |
| US7117026B2 (en) | Physiological model based non-rigid image registration | |
| CN109961419B (en) | Correction information acquisition method for attenuation correction of PET activity distribution image | |
| CN102067176B (en) | Radiological imaging incorporating local motion monitoring, correction, and assessment | |
| Kyme et al. | Motion estimation and correction in SPECT, PET and CT | |
| CN109978966B (en) | Correction information acquisition method for attenuation correction of PET activity distribution image | |
| CN104700438B (en) | Image rebuilding method and device | |
| Chan et al. | Non-rigid event-by-event continuous respiratory motion compensated list-mode reconstruction for PET | |
| CN110458779B (en) | Method for acquiring correction information for attenuation correction of PET images of respiration or heart | |
| US20110082368A1 (en) | Reconstruction of dynamical cardiac spect for measuring tracer uptake and redistribution | |
| CN114387364A (en) | Linear attenuation coefficient acquisition method and reconstruction method for PET image reconstruction | |
| Li et al. | Deep learning based joint PET image reconstruction and motion estimation | |
| CN110428384B (en) | Method for acquiring correction information for attenuation correction of PET images of respiration or heart | |
| CN115439572A (en) | Attenuation correction coefficient image acquisition method and PET image reconstruction method | |
| US10417793B2 (en) | System and method for data-consistency preparation and image reconstruction | |
| CN106548473B (en) | A kind of method and device constructing phase image | |
| Ozsahin et al. | Monte Carlo simulation of PET/MR scanner and assessment of motion correction strategies | |
| Li et al. | 3D and 4D medical image registration combined with image segmentation and visualization | |
| Wang et al. | Motion correction strategies for enhancing whole-body PET imaging |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |