Scattering correction method and device, digital equipment and computer readable storage medium
Through the maximum likelihood-desired maximization iterative algorithm combined with the moment estimation method, the problems of low accuracy of scattering correction and complex calculation in the prior art are solved, and efficient and accurate scattering correction effect is achieved.
Patent Information
- Application Number
- CN202311546053.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2023-11-20
- Publication Date
- 2025-05-20
- Estimated Expiration
- 2043-11-20
AI Technical Summary
In the existing signal sampling technology, simulation-based scattering correction methods have problems such as low accuracy, complex calculation and low efficiency, especially in the case of multi-scattering.
The maximum likelihood-expected maximization iterative algorithm is used to combine the moment estimation method to obtain the downsampled response line based on the detection data, and the objective function is estimated to obtain the number of scattered events through two-dimensional instant matching energy histogram, delay matching energy histogram, probability density function of scattered photons, and probability density function of unscattered photons.
Scatter correction that does not depend on activity images and attenuation images is realized, the calculation process is simplified, the accuracy and efficiency of correction results are improved, and the dependence on hyperparameter β is reduced.
Smart Images

Figure CN120019792A_ABST
Abstract
Description
Technical Field
[0001] The present application relates to the technical field of signal sampling, and particularly relates to a scatter correction method, device, digital device, and computer-readable storage medium. Background Art
[0002] The working principle of positron emission tomography (PET for short) is to label a radioactive nuclide onto a compound that can participate in the blood flow or metabolic process of a living tissue, and inject the compound into a living body. The positrons emitted by the radioactive nuclide in the living body then combine with the negative electrons in the living body, resulting in the annihilation event of an electron pair and generating two gamma photons with equal energy and opposite directions. Since the flight directions of these two gamma photons are different, the time when the detector detects these two gamma photons is also different. If two scintillation crystals located on the line of response (LOR for short) in the detector detect two gamma photons respectively within a specified coincidence time window (for example, 0 to 15 nanoseconds), the event of detecting these two gamma photons can be called a coincidence event.
[0003] Coincidence events generally include true coincidence events, scatter coincidence events, and random coincidence events. Among them, a true coincidence event refers to an event where the time difference between two gamma photons generated by the same annihilation event reaching two scintillation crystals located on the line of response is within the coincidence time window. A random coincidence event is a false coincidence event. In a random coincidence event, the two detected gamma photons come from different annihilation events but are misinterpreted as two gamma photons that occur "simultaneously" within the coincidence time window. A scatter coincidence event refers to the following event: for two gamma photons generated by the same annihilation event detected, one of the gamma photons changes its flight direction due to physical effects such as Compton scattering and / or Rayleigh scattering during flight.
[0004] Among these three types of coincidence events, since the data collected for random coincidence events and scatter coincidence events may be incorrect, this may affect the resolution, contrast, and positioning accuracy of PET imaging. Therefore, it is crucial to correct the collected coincidence events.
[0005] Currently, the more commonly used scatter correction methods mainly include multi-energy window technology, convolution / deconvolution technology, and simulation-based technology, etc. Among these technologies, the most accurate and widely used is the simulation-based technology.
[0006] Simulation-based techniques generally include single-scattering simulation methods, double-scattering simulation methods, and Monte Carlo simulation methods. Among them, the single-scattering simulation method is to obtain the photon motion path of single scattering (i.e., a pair of gamma photons undergoes a total of one scattering) by selecting a scattering point in the input activity image and attenuation image for each response line. By calculating the single-scattering events generated on all photon motion paths corresponding to this response line, the single-scattering events on this response line are obtained. Finally, using only the tail data containing scattering events, the single-scattering events obtained are fitted by the Tail Fitting (TF) technique to obtain the stretching factor of the scatter coincidence events, and this stretching factor is applied to all data to achieve scatter correction. The double-scattering simulation method is usually used in combination with the single-scattering simulation method. The main difference between it and the single-scattering simulation method is that the double-scattering simulation method determines the photon motion path through two scattering points. The Monte Carlo simulation method mainly determines the total scatter distribution by simulating the motion of each pair of gamma photons determined based on the input activity image and attenuation image.
[0007] In the process of implementing this application, the inventors found that the simulation-based techniques have at least the following problems:
[0008] (1) The single-scattering simulation method only considers the case of single scattering. However, in practice, multiple scattering may occur, which may lead to low accuracy of the scatter correction result, thereby resulting in low quality of the reconstructed image. Although the double-scattering simulation method and the Monte Carlo simulation method consider the case of multiple scattering, the calculation process is complex and the amount of calculation is very large, which will slow down the image reconstruction process and the efficiency is not high.
[0009] (2) Before reconstructing the activity image, it is necessary to first obtain the scatter events. However, for the simulation-based method to obtain scatter events, it is necessary to know the activity image first, and it is difficult to obtain an accurate activity image in the prior art.
[0010] (3) Since it is difficult to obtain the activity image outside the field of view, therefore, the prior art usually does not consider external radiation (i.e., the scatter events generated by gamma photon pairs located outside the field of view), which may lead to low accuracy of the correction result, thereby affecting the quality of the reconstructed image.
[0011] (4) The currently adopted OSL-EDR and OT-EDR algorithms need to select appropriate hyperparameters β to obtain better results. The optimal β requires some pre-experiments to adjust, and different sizes of phantoms require different β, which brings many inconveniences to the use of the OSL-EDR and OT-EDR algorithms.
[0012] (5) Currently, the moment estimation method is used to obtain the estimated objective function corresponding to each downsampled response line. In the estimated objective function, the probability density function of unscattered photons is defaulted to be the same, and the probability density function of scattered photons is also defaulted to be the same. As a result, the number of scattered coincidence events corresponding to the obtained downsampled response line may not be accurate. In view of this, there is an urgent need to provide a scattering correction method that is different from the simulation-based technology and does not depend on the activity image and the attenuation image.
[0013] The content in the background art section is only the technology known to the inventor and is not regarded as representing the prior art in this field. Summary of the Invention
[0014] The present application aims to provide a scattering correction method, device, digital device, and computer-readable storage medium to solve at least one problem existing in the prior art.
[0015] According to the first aspect of the present application, a scattering correction method is provided. The scattering correction method includes: obtaining a downsampled response line based on detection data; obtaining an estimated objective function corresponding to each downsampled response line by using the moment estimation method based on the maximum likelihood-expectation maximization iterative algorithm, where the estimated objective function is obtained based on the two-dimensional prompt coincidence energy histogram, two-dimensional delayed coincidence energy histogram, probability density function of scattered photons, and probability density function of unscattered photons corresponding to each downsampled response line; solving the estimated objective function to obtain the number of scattered coincidence events corresponding to each downsampled response line; performing an upsampling process on the downsampled response line to obtain an upsampled response line; and calculating the number of scattered coincidence events corresponding to each upsampled response line.
[0016] In some embodiments, obtaining the downsampled response line based on detection data includes: obtaining a coincidence response line based on detection data; and performing a downsampling process on the coincidence response line to obtain a downsampled response line.
[0017] In some embodiments, the detection data is obtained based on a target object.
[0018] In some embodiments, obtaining an estimated objective function corresponding to each downsampled response line by using the moment estimation method based on the maximum likelihood-expectation maximization iterative algorithm includes: obtaining a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, and two one-dimensional detection energy spectra corresponding to each downsampled response line based on each coincidence event corresponding to each downsampled response line, and obtaining the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding downsampled response line by using the maximum likelihood-expectation maximization iterative algorithm based on the two one-dimensional detection energy spectra.
[0019] In some embodiments, based on the coincidence events corresponding to each downsampled response line, two one-dimensional detection energy spectra are obtained, including: splitting all the coincidence events corresponding to the downsampled response line into single events; dividing the single events into two parts according to the position information of the single events, with each part corresponding to a one-dimensional detection energy spectrum; and obtaining two one-dimensional detection energy spectra based on the two divided parts of single events.
[0020] In some embodiments, based on the two one-dimensional detection energy spectra, the probability density functions of scattered photons and unscattered photons corresponding to the corresponding downsampled response line are obtained by using the maximum likelihood-expectation maximization iterative algorithm, including: obtaining the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum; based on the two obtained one-dimensional detection energy spectra and the correspondence function, using the maximum likelihood-expectation maximization iterative algorithm and the obtained initial iteration value corresponding to the corresponding downsampled response line to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra; and estimating the probability density functions of two scattered photons and two unscattered photons corresponding to the corresponding downsampled response line based on the two one-dimensional gamma photon energy spectra.
[0021] In some embodiments, the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is:
[0022]
[0023] where, represents the expected value of the one-dimensional detection energy spectrum obtained based on the scanning of the target object corresponding to a single downsampled response line; X represents the one-dimensional gamma photon energy spectrum, p represents the blurring response matrix corresponding to a single downsampled response line, and r represents the one-dimensional delay energy spectrum corresponding to the random coincidence events included in a single response line.
[0024] In some embodiments, the one-dimensional delay energy spectrum corresponding to the random coincidence events is estimated by delayed coincidence events.
[0025] In some embodiments, obtaining the one-dimensional delay energy spectrum includes: obtaining all the delayed coincidence events in the downsampled response line based on the delayed coincidence window; splitting each delayed coincidence event into two corresponding delayed single events; dividing the delayed single events into two parts according to the position information of the delayed single events, with each part corresponding to a one-dimensional delay energy spectrum; and obtaining two one-dimensional delay energy spectra based on the two divided parts of delayed single events; where the position information of the delayed single events includes the position information of the downsampled detection module where the delayed single event is detected.
[0026] In some embodiments, obtaining the blurred response matrix p includes: obtaining prior detection data based on a prosthesis; obtaining prior downsampled response lines based on the prior detection data; dividing all prior coincidence events corresponding to the prior downsampled response lines into prior prompt coincidence events and prior delayed coincidence events; obtaining a prior one-dimensional prompt energy spectrum based on the prior prompt coincidence events corresponding to the downsampled response lines; obtaining a prior one-dimensional delayed energy spectrum based on the prior delayed coincidence events corresponding to the downsampled response lines; obtaining the difference between the prior one-dimensional prompt energy spectrum and the prior one-dimensional delayed energy spectrum, normalizing the difference to obtain a one-dimensional intermediate energy spectrum; and obtaining the blurred response matrix based on the one-dimensional intermediate energy spectrum.
[0027] In some embodiments, obtaining the two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra includes: obtaining a maximum likelihood-expectation maximization (MLE-EM) iteration formula based on the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum; setting the number of iterations, and iterating through the MLE-EM iteration formula based on the obtained initial iteration value until the set number of iterations is reached to obtain the two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
[0028] In some embodiments, the MLE-EM iteration formula is:
[0029]
[0030] where Y i in Y represents the one-dimensional detection energy spectrum obtained based on scanning a target object corresponding to a single downsampled response line, i represents the i-th vertical bar region in the energy spectrum, b and j respectively represent the i-th and j-th vertical bar regions in the energy spectrum, k is the number of iterations, p ij and p ib in p is the blurred response matrix corresponding to this downsampled response line, p ij is the value of the i-th row and j-th column of the blurred response matrix corresponding to this downsampled response line, p ib is the value of the i-th row and b-th column of the blurred response matrix corresponding to this downsampled response line; r i represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in this downsampled response line; is the one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region of the energy spectrum obtained in the (k + 1)-th iteration; respectively represent the one-dimensional gamma photon energy spectra corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration.
[0031] In some embodiments, obtaining an initial iteration value includes: obtaining a global initial scattered photon energy spectrum; performing deblurring processing on the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value; stretching the global intermediate initial iteration value to obtain the initial iteration value corresponding to each downsampled response line.
[0032] In some embodiments, obtaining a global initial scattered photon energy spectrum includes: obtaining coincidence response lines based on detection data; dividing all coincidence events corresponding to all coincidence response lines into prompt coincidence events and random coincidence events; obtaining a global one-dimensional prompt energy spectrum based on the prompt coincidence events; obtaining a global one-dimensional delayed energy spectrum based on the delayed coincidence events; subtracting the global one-dimensional delayed energy spectrum corresponding to the delayed coincidence events from the global one-dimensional prompt energy spectrum corresponding to the prompt coincidence events on all response lines to obtain a global one-dimensional undelayed energy spectrum; obtaining an objective function of the global initial scattered photon energy spectrum based on the global one-dimensional undelayed energy spectrum.
[0033] In some embodiments, the objective function of the global initial scattered photon energy spectrum is:
[0034]
[0035] wherein, S i in S is the global initial scattered photon energy spectrum; i represents the i-th vertical bar region in the energy spectrum, C i in C is the global one-dimensional undelayed energy spectrum, obtained based on the target object, U i in U is the global one-dimensional intermediate energy spectrum obtained by scanning the entire PET system, and the global one-dimensional intermediate energy spectrum is obtained based on the phantom; hel and heu are the lower threshold and upper threshold of the set high-energy window, is the stretching coefficient.
[0036] In some embodiments, obtaining U i includes: obtaining prior detection data based on the phantom; obtaining all prior coincidence response lines based on the prior detection data; dividing all prior coincidence events corresponding to all prior coincidence response lines into global prior prompt coincidence events and global prior delayed coincidence events; obtaining a global prior one-dimensional prompt energy spectrum based on the global prior prompt coincidence events; obtaining a global prior one-dimensional delayed energy spectrum based on the global prior delayed coincidence events; obtaining the difference U i .
[0037] In some embodiments, deblurring the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value includes: obtaining a deblurring iteration formula based on the maximum likelihood-expectation maximization iteration algorithm; setting a prior initial iteration value and iteration conditions, and performing iteration through the deblurring iteration formula until the iteration conditions are met to obtain the scattered part of the global initial energy spectrum; estimating the unscattered part of the global initial energy spectrum based on the scattered part of the global initial energy spectrum; and obtaining a global intermediate initial iteration value based on the scattered part of the global initial energy spectrum and the unscattered part of the global initial energy spectrum.
[0038] In some embodiments, the deblurring iteration formula is:
[0039]
[0040] where S i in S represents the global one-dimensional detection energy spectrum obtained based on scanning the target object, i represents the i-th vertical bar region in the energy spectrum, b and j respectively represent the i-th and j-th vertical bar regions in the energy spectrum, k is the number of iterations, P ij , P ib in P is the global blurring response matrix, P ij is the value of the i-th row and j-th column of the global blurring response matrix, P ib is the value of the i-th row and b-th column of the global blurring response matrix; r i represents the global one-dimensional delayed energy spectrum corresponding to all random coincidence events; is the global one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region of the energy spectrum obtained in the (k + 1)-th iteration; respectively represent the global one-dimensional gamma photon energy spectra corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration;
[0041]
[0042] where largeconstant represents a constant.
[0043] In some embodiments, obtaining the global blurring response matrix P includes: obtaining prior detection data based on a phantom; obtaining all prior coincidence response lines based on the prior detection data; dividing all prior coincidence events corresponding to all prior coincidence response lines into global prior prompt coincidence events and global prior delayed coincidence events; obtaining a global prior one-dimensional prompt energy spectrum based on the global prior prompt coincidence events; obtaining a global prior one-dimensional delayed energy spectrum based on the global prior delayed coincidence events; obtaining the difference between the global prior prompt energy spectrum and the global prior delayed energy spectrum; performing normalization processing on the difference to obtain a global prior intermediate energy spectrum; and obtaining the global blurring response matrix based on the global prior intermediate energy spectrum.
[0044] In some embodiments, estimating the unscattered part of the global initial energy spectrum is performed through an estimation function, and the estimation function is:
[0045]
[0046] where P is the global blur response matrix; δ is a hyperparameter; κ is the ratio of unscattered photons to scattered photons;
[0047]
[0048] where sum() represents summation, and C and S respectively correspond to C i and S i .
[0049] In some embodiments, the characterization function of the global intermediate initial iteration value is:
[0050]
[0051] where δ is a hyperparameter, P is the global blur response matrix, and κ is the ratio of unscattered photons to scattered photons.
[0052] In some embodiments, stretching the intermediate initial iteration value to obtain the initial iteration value corresponding to each downsampled response line includes:
[0053] Stretching the global intermediate initial iteration value based on a stretching function to obtain the initial iteration value corresponding to each downsampled response line, and the stretching function is:
[0054]
[0055] where X 0 is the initial iteration value, sum() represents summation, P is the global blur response matrix, and X init is the global intermediate initial iteration value; Y is the one-dimensional detection energy spectrum obtained based on the scanning of the target object corresponding to a single downsampled response line, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in a single response line.
[0056] In some embodiments, estimating the probability density functions of two scattered photons and the probability density functions of two unscattered photons corresponding to the corresponding downsampled response lines based on the two one-dimensional gamma photon energy spectra includes: dividing the one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part; performing restoration and deblurring processing on the scattered part and the unscattered part of the two one-dimensional gamma photon energy spectra to obtain the two scattered parts and the two unscattered parts corresponding to the corresponding downsampled response lines of the two one-dimensional detection energy spectra; performing normalization processing on the two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra, and estimating the probability density functions of the two scattered photons and the probability density functions of the two unscattered photons corresponding to the corresponding downsampled response lines based on the obtained normalization values.
[0057] In some embodiments, the scattered part is characterized as:
[0058]
[0059] The unscattered part is characterized as:
[0060]
[0061] Wherein, is the scattered part, is the unscattered part.
[0062] In some embodiments, dividing the one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part includes: dividing the one-dimensional gamma photon energy spectrum based on the energy value when the gamma photon does not scatter, the part less than the energy value is the scattered part, and the part equal to the energy value is the unscattered part.
[0063] In some embodiments, the estimation objective function is:
[0064]
[0065] Wherein, σ 0 respectively represent the number of true coincidence events and do not require subsequent processing; σ 1 +σ 2 +σ 3 represents the number of scattered coincidence events AC, which is an unknown quantity to be solved; the calculation methods of the coefficients α 1,m , α 2,n , β 1,m , β 2,n are:
[0066]
[0067]
[0068] Wherein, and is the probability density function of unscattered photons corresponding to a detection module in the downsampled response line, and is the probability density function of scattered photons corresponding to another detection module in the downsampled response line;
[0069] wherein, the calculated σ 1 +σ 2 +σ 3 value is the number of scattered coincidence events corresponding to the corresponding downsampled response line.
[0070] In some embodiments, upsampling processing is performed on the downsampled response line, including: using a bilinear interpolation method, a 4D linear interpolation method, or a 5D linear interpolation method to process the downsampled response line to obtain the upsampled response line corresponding to each downsampled response line.
[0071] In some embodiments, when using the 4D linear interpolation method to process the downsampled response line, the following formula is used to calculate the number of scattered coincidence events corresponding to the upsampled response line:
[0072]
[0073] wherein, represents the number of scattered coincidence events of the upsampled response line formed by the upsampled crystals i and j after 4D linear interpolation; S i represents the set of the downsampled central crystal, the downsampled axially adjacent crystals, the downsampled radially adjacent crystals, and the downsampled relative crystals corresponding to the upsampled crystal i; K and L represent the numbers of the crystals corresponding to the downsampled response line; represents the number of scattered coincidence events on the downsampled response line formed by the downsampled crystals K and L; represents the weight of the downsampled crystal K with respect to the upsampled crystal i; represents the weight of the downsampled crystal L with respect to the upsampled crystal j; W tot is a normalization factor, represented as;
[0074]
[0075] where N is the number of upsampled response lines included in a downsampled response line.
[0076] In some embodiments, calculating the number of scattered coincidence events corresponding to each upsampled response line includes: calculating the ratio of the total number of the scattered coincidence events on each downsampled response line to the total number of coincidence events; calculating the number of scattered coincidence events corresponding to the corresponding upsampled response line according to the number of coincidence events on each upsampled response line included in each downsampled response line and the ratio.
[0077] According to a second aspect of the present application, a scatter correction method is provided. The scatter correction method includes: obtaining response lines based on detection data; obtaining an estimated objective function corresponding to each response line by using the moment estimation method based on the maximum likelihood-expectation maximization (MLE-EM) iterative algorithm, where the estimated objective function is obtained based on a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each response line; solving the estimated objective function to obtain the number of scattered coincidence events corresponding to each response line; and calculating the number of scattered coincidence events corresponding to each response line.
[0078] In some embodiments, obtaining an estimated objective function corresponding to each response line by using the moment estimation method based on the MLE-EM iterative algorithm includes: obtaining a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, and two one-dimensional detection energy spectra corresponding to each response line based on each coincidence event corresponding to each response line, and obtaining a probability density function of scattered photons and a probability density function of unscattered photons corresponding to the corresponding response line by using the MLE-EM iterative algorithm based on the two one-dimensional detection energy spectra.
[0079] In some embodiments, obtaining two one-dimensional detection energy spectra based on each coincidence event corresponding to each response line includes: splitting all coincidence events corresponding to the response line into single events; dividing the single events into two parts according to the position information of the single events, with each part corresponding to a one-dimensional detection energy spectrum; and obtaining two one-dimensional detection energy spectra based on the two divided parts of single events.
[0080] In some embodiments, obtaining a probability density function of scattered photons and a probability density function of unscattered photons corresponding to the corresponding response line by using the MLE-EM iterative algorithm based on the two one-dimensional detection energy spectra includes: obtaining a corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum; obtaining two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra by using the MLE-EM iterative algorithm and the obtained initial iteration value corresponding to the corresponding response line based on the two obtained one-dimensional detection energy spectra and the corresponding relationship function; and estimating a probability density function of two scattered photons and a probability density function of two unscattered photons corresponding to the corresponding response line based on the two one-dimensional gamma photon energy spectra.
[0081] In some embodiments, obtaining two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra includes: obtaining an MLE-EM iterative formula based on the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum; setting the number of iterations, and performing iterations to the set number of iterations by using the MLE-EM iterative formula based on the obtained initial iteration value to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
[0082] In some embodiments, obtaining an initial iteration value includes: obtaining a global initial scattered photon energy spectrum; performing deblurring processing on the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value; and stretching the global intermediate initial iteration value to obtain the initial iteration value corresponding to each response line.
[0083] In some embodiments, obtaining a global initial scattered photon energy spectrum includes: obtaining coincidence response lines based on detection data; dividing all coincidence events corresponding to all coincidence response lines into prompt coincidence events and random coincidence events; obtaining a global one-dimensional prompt energy spectrum based on the prompt coincidence events; obtaining a global one-dimensional delayed energy spectrum based on the delayed coincidence events; subtracting the global one-dimensional delayed energy spectrum corresponding to the delayed coincidence events from the global one-dimensional prompt energy spectrum corresponding to the prompt coincidence events on all response lines to obtain a global one-dimensional non-delayed energy spectrum; and obtaining an objective function of the initial scattered photon energy spectrum based on the global one-dimensional non-delayed energy spectrum.
[0084] In some embodiments, performing deblurring processing on the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value includes: obtaining a deblurring processing iteration formula based on the maximum likelihood-expectation maximization iteration algorithm; setting a prior initial iteration value and the number of iterations, and performing iterations through the deblurring processing iteration formula until the number of iterations is reached to obtain the scattered part of the global initial energy spectrum; estimating the unscattered part of the global initial energy spectrum based on the scattered part of the global initial energy spectrum; and obtaining the global intermediate initial iteration value based on the scattered part of the global initial energy spectrum and the unscattered part of the global initial energy spectrum.
[0085] In some embodiments, estimating the probability density functions of two scattered photons and the probability density functions of two unscattered photons corresponding to a corresponding response line based on the two one-dimensional gamma photon energy spectra includes: dividing the one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part; performing restoration and blurring response processing on the scattered part and the unscattered part of the two one-dimensional gamma photon energy spectra to obtain the two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra corresponding to the corresponding response line; performing normalization processing on the two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra, and estimating the probability density functions of two scattered photons and the probability density functions of two unscattered photons corresponding to the corresponding response line based on the obtained normalization values.
[0086] In some embodiments, dividing the one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part includes: dividing the one-dimensional gamma photon energy spectrum based on the energy value when the gamma photon does not undergo scattering, where the part less than the energy value is the scattered part, and the part equal to the energy value is the unscattered part.
[0087] According to a third aspect of the present application, there is provided an image reconstruction method, the image reconstruction method comprising: obtaining the number of scattered coincidence events corresponding to each response line based on the scattering correction method described in any one of the above embodiments; correcting the scattering events based on the number of scattered coincidence events to obtain a reconstructed image.
[0088] According to a fourth aspect of the present application, there is provided a scattering correction device, the scattering correction device comprising: a downsampled response line obtaining module configured to obtain downsampled response lines based on detection data; an estimation module configured to obtain an estimated objective function corresponding to each downsampled response line based on the maximum likelihood-expectation maximization iteration algorithm and using the moment estimation method, the estimated objective function being obtained based on a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each downsampled response line; a scattered coincidence event number obtaining module configured to solve the estimated objective function to obtain the number of scattered coincidence events corresponding to each downsampled response line; an upsampling module configured to perform upsampling processing on the downsampled response lines to obtain upsampled response lines; and a scattered coincidence event number calculation module configured to calculate the number of scattered coincidence events corresponding to each upsampled response line.
[0089] According to a fifth aspect of the present application, there is provided a scattering correction device, the scattering correction device comprising: a response line obtaining module configured to obtain response lines based on detection data; an estimation module configured to obtain an estimated objective function corresponding to each response line based on the maximum likelihood-expectation maximization iteration algorithm and using the moment estimation method, the estimated objective function being obtained based on a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each response line; and a scattered coincidence event number obtaining module configured to solve the estimated objective function to obtain the number of scattered coincidence events corresponding to each response line.
[0090] According to a sixth aspect of the present application, there is provided an image reconstruction device, the image reconstruction device comprising: a scattered coincidence event number obtaining module configured to obtain the number of scattered coincidence events corresponding to response lines based on the scattering correction device described in any one of the above embodiments; and a reconstruction module configured to perform image reconstruction based on the number of scattered coincidence events using an iterative image reconstruction algorithm to obtain a reconstructed image.
[0091] According to a seventh aspect of the present application, there is provided a digital device, comprising: a memory, a processor, and a computer program stored on the memory and executable on the processor, the computer program, when executed by the processor, implementing the steps of the method described in any one of the above embodiments.
[0092] According to an eighth aspect of the present application, there is provided a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the method described in any one of the above embodiments are implemented.
[0093] Based on the above embodiments of the present application, the beneficial effects of the present application include one or a combination of more than one of the following effects:
[0094] In some embodiments, a solution that combines the maximum likelihood-expectation maximization iterative algorithm with the moment estimation method is adopted. Based on the maximum likelihood-expectation maximization iterative algorithm, the moment estimation method is used to obtain the estimated objective function corresponding to each downsampled response line. On the one hand, the maximum likelihood-expectation maximization iterative algorithm (MLEM) algorithm is used for scatter correction, which does not depend on the activity image and the attenuation image, making the scatter correction easy to implement. And this method does not require the hyperparameter β, so there is no need to set the hyperparameter β within a certain range, that is, it is hardly affected by the hyperparameter β and almost no adjustment of the hyperparameter β is required, making the scatter correction more convenient. On the other hand, the unscattered photon probability density function and in the estimated objective function are not necessarily the same, and the scattered photon probability density function and are not necessarily the same either. Instead, it is specifically determined according to the two detection modules corresponding to the downsampled response line, and the number of scattered coincidence events corresponding to the downsampled response line obtained is more accurate. BRIEF DESCRIPTION OF THE DRAWINGS
[0095] The embodiments of the present application will be described in detail below with reference to the drawings. Here, the drawings forming a part of the present application are used to provide a further understanding of the present application. The exemplary embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation of the present application. In the drawings:
[0096] Figure 1 An exemplary flowchart showing a scatter correction method according to an exemplary embodiment of the present application;
[0097] Figure 2 A schematic diagram showing crystal weight calculation in 4D linear interpolation according to an exemplary embodiment of the present application;
[0098] Figure 3 An exemplary flowchart showing a scatter correction method according to another exemplary embodiment of the present application;
[0099] Figure 4 An exemplary flowchart showing an image reconstruction method according to an exemplary embodiment of the present application;
[0100] Figure 5 An exemplary module diagram showing a scatter correction method according to an exemplary embodiment of the present application;
[0101] Figure 6Exemplary block diagram showing a scatter correction method according to another exemplary embodiment of the present application;
[0102] Figure 7 Exemplary block diagram showing an image reconstruction apparatus according to an exemplary embodiment of the present application. Detailed Description
[0103] In the following, only some exemplary embodiments are briefly described. As those skilled in the art will recognize, the described embodiments can be modified in various different ways without departing from the spirit or scope of the present application. Therefore, the drawings and the description are to be regarded as illustrative in nature and not restrictive.
[0104] In the description of the present application, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", etc. indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present application and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and thus should not be construed as a limitation of the present application. In addition, the terms "first" and "second" are only used for description and cannot be understood as indicating or implying relative importance or implicitly specifying the number of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more similar features. In the description of the present application, the meaning of "plurality" is two or more unless otherwise specifically defined.
[0105] In the description of the present application, it should be noted that, unless otherwise clearly specified and limited, the terms "mounted", "connected", and "coupled" should be understood in a broad sense. For example, it may be a fixed connection, a detachable connection, or an integral connection; it may be a mechanical connection, an electrical connection, or a communication connection with each other; it may be directly connected, or indirectly connected through an intermediate medium, and it may be the internal communication of two elements or the interaction relationship between two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present application can be understood according to specific circumstances.
[0106] In this application, unless otherwise clearly specified and defined, the first feature being "on" or "under" the second feature may include direct contact between the first and second features, or may include the first and second features not being in direct contact but being in contact through additional features therebetween. Moreover, the first feature being "above", "over" and "on top of" the second feature includes the first feature being directly above and obliquely above the second feature, or merely indicating that the relative height of the first feature in a certain dimension is higher than that of the second feature. The first feature being "under", "beneath" and "underneath" the second feature includes the first feature being directly below and obliquely below the second feature, or merely indicating that the relative position of the first feature in a certain dimension is lower than that of the second feature.
[0107] The following provides different embodiments or examples for implementing different structures of this application. To simplify this application, the components and settings of specific examples are described below. Of course, they are only examples and are not intended to limit this application. This application may repeat the reference numerals in different examples, and this repetition is for the purpose of simplification and clarity, and does not itself indicate the relationship between the various embodiments and / or settings described. In addition, this application provides examples of various specific processes and materials, but those skilled in the art can apply other processes and / or use alternative materials according to the inspiration of this application.
[0108] Some preferred embodiments of this application are described below with reference to the accompanying drawings. It should be noted that the following description is for illustrative purposes and is not intended to limit the protection scope of this application.
[0109] Figure 1 is an exemplary flowchart of a scatter correction method shown in some embodiments of this application.
[0110] Continue to refer to Figure 1 , the scatter correction method 100 may include the following steps.
[0111] S110, obtaining downsampled response lines based on the detection data.
[0112] In some embodiments of this application, before step S110, it may further include:
[0113] Step S0111, obtaining detection data based on the target object.
[0114] In some embodiments, the target object may be a living object, including but not limited to humans, animals, etc.
[0115] In some embodiments, the detection data may be sampling data obtained based on the detection of a detector and then based on multi-voltage threshold (MVT) sampling. Specifically, the detection data includes voltage threshold time pairs. Specifically, the specific design of the detector can refer to the prior art and will not be elaborated here.
[0116] In some embodiments of the present application, step S110 may further include the following steps:
[0117] S111. Obtain a coincidence response line based on the detection data.
[0118] In some embodiments, the detection data is screened based on a time window and an energy window to obtain coincidence events, thereby obtaining a coincidence response line. Specifically, reference may be made to the prior art and details are not described herein.
[0119] S112. Perform downsampling processing on the coincidence response line to obtain a downsampled response line.
[0120] Those skilled in the art can understand that a response line refers to the connection line of a pair of scintillation crystals in a detection module of a PET system. Since the number of coincidence events on a response line during normal scanning may be relatively small, it is difficult to obtain the corresponding energy spectrum based on the coincidence events corresponding to one response line. It is necessary to perform downsampling processing to obtain a downsampled response line. For example, originally, a detection module pair of a PET system included 3×3 = 9 coincidence response lines. After downsampling the 9 response lines, a downsampled response line formed by merging the 9 response lines was obtained. Specifically, the downsampling processing may be to merge the response lines included in a detection module pair, or may also merge the response lines included in multiple adjacent (such as two axially, 4 axially and radially adjacent, etc.) detection module pairs.
[0121] In some embodiments, the downsampling processing includes merging multiple coincidence response lines that are parallel to each other or whose included angle and spatial position meet preset conditions into one downsampled response line, but is not limited thereto. The downsampling processing can specifically refer to the prior art and details are not described herein.
[0122] S120. Based on the maximum likelihood-expectation maximization iterative algorithm, use the moment estimation method to obtain the estimated objective function corresponding to each downsampled response line, and the estimated objective function is obtained based on the two-dimensional prompt coincidence energy histogram, two-dimensional delayed coincidence energy histogram, probability density function of scattered photons, and probability density function of unscattered photons corresponding to each downsampled response line.
[0123] In some embodiments, the maximum likelihood-expectation maximization iterative algorithm MLEM (maximum likelihood-expectation maximization, abbreviated as MLEM) is obtained based on the standard maximum likelihood-expectation maximization algorithm in PET image reconstruction. Specifically, reference may be made to the prior art and details are not described herein.
[0124] In some embodiments, step S120 includes:
[0125] Based on each coincidence event corresponding to each downsampled response line, obtain a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, and two one-dimensional detection energy spectra corresponding to each downsampled response line. Based on the two one-dimensional detection energy spectra, use the maximum likelihood-expectation maximization iterative algorithm to obtain the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding downsampled response line.
[0126] In some embodiments, the energy spectrum is characterized as a one-dimensional energy histogram, where the abscissa represents energy and the ordinate represents the gamma photon count. The one-dimensional detection energy spectrum refers to the energy spectrum formed by the energy obtained when gamma photons are detected by the detector. It should be understood that the energy spectrum histogram includes multiple vertical bar regions, and each vertical bar region is called a bin. For details, reference can be made to the prior art and will not be elaborated here.
[0127] In some embodiments, based on each coincidence event corresponding to each downsampled response line, obtaining two one-dimensional detection energy spectra includes:
[0128] S1211, split all the coincidence events corresponding to the downsampled response line into single events.
[0129] It should be understood that one coincidence event contains two single events. Split each coincidence event included in a downsampled response line into two corresponding single events. The single event contains position information, energy information, and time information.
[0130] S1212, divide the single events into two parts according to the position information of the single events, and each part corresponds to a one-dimensional detection energy spectrum.
[0131] In some embodiments, the position information of the single event includes the position information of the downsampled detection module that detects the single event. For example, a downsampled response line corresponds to two detection modules A and B. One part of the single events included in this downsampled response line corresponds to detection module A, and the other part corresponds to detection module B.
[0132] S1213, based on the two parts of single events divided, obtain two one-dimensional detection energy spectra.
[0133] In some specific examples, according to the detection module information corresponding to the single event, two one-dimensional detection energy spectra corresponding to the downsampled response line can be obtained. Continuing with the above example, one one-dimensional detection energy spectrum can be obtained based on the single events corresponding to detection module A, and the other one-dimensional detection energy spectrum corresponding to the corresponding downsampled response line can be obtained based on the single events corresponding to detection module B.
[0134] In some specific embodiments, based on the two one-dimensional detection energy spectra, the maximum likelihood-expectation maximization iterative algorithm is used to obtain the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding downsampled response lines, including:
[0135] S1221. Obtain the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum.
[0136] In some embodiments, the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is:
[0137]
[0138] Wherein, represents the expected value of the one-dimensional detection energy spectrum obtained based on the scanning of the target object corresponding to a single downsampled response line; X represents the one-dimensional gamma photon energy spectrum, p represents the blurring response matrix corresponding to a single downsampled response line, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in a single response line.
[0139] In some embodiments, the one-dimensional delayed energy spectrum corresponding to the random coincidence event is obtained by estimating the delayed coincidence event.
[0140] In some specific embodiments, obtaining the one-dimensional delayed energy spectrum includes:
[0141] S1221a. Based on the delayed coincidence window, obtain all the delayed coincidence events in the downsampled response line; split each delayed coincidence event into two corresponding delayed single events.
[0142] In some embodiments, the delayed coincidence window includes a delayed coincidence time window. Specifically, when the signal is delayed, the corresponding coincidence event obtained is the delayed coincidence event, and the corresponding time window is called the delayed coincidence time window. The delayed coincidence time window is specifically obtained according to prior information. For example, it can be set to 2 ns - 10 ns. Specifically, regarding how to obtain the delayed coincidence event specifically, reference can be made to the prior art and will not be elaborated here.
[0143] It should be understood that a delayed coincidence event includes two delayed single events. Split each delayed coincidence event included in a downsampled response line into two corresponding delayed single events. The delayed single event includes position information, energy information, and time information.
[0144] S1221b. Divide the delayed single events into two parts according to the position information of the delayed single events, and each part corresponds to a one-dimensional delayed energy spectrum.
[0145] In some specific embodiments, correspondingly to step S132, the position information of the delayed single event includes the position information of the downsampled detection module that detects the delayed single event.
[0146] S1221c. Based on the two - part delayed single event obtained by partitioning, two one - dimensional delayed energy spectra are acquired; wherein, the position information of the delayed single event includes the position information of the down - sampled detection module that detects the delayed single event.
[0147] In some specific examples, according to the information of the down - sampled detection module corresponding to the delayed single event, two one - dimensional delayed energy spectra corresponding to the down - sampled response line can be obtained.
[0148] In some specific embodiments, acquiring the fuzzy response matrix p includes:
[0149] S1221a'. Based on the phantom, prior detection data is acquired.
[0150] In some embodiments, the phantom can be a point source or a line source, but is not limited thereto.
[0151] S1221b'. Based on the prior detection data, prior down - sampled response lines are acquired.
[0152] In some embodiments, corresponding to step S110, step S1221b' includes:
[0153] Based on the prior detection data, prior coincidence response lines are acquired;
[0154] The prior coincidence response lines are down - sampled to obtain prior down - sampled response lines.
[0155] In some specific examples, how to acquire the prior coincidence response lines and the prior down - sampled response lines can refer to step S110, which will not be elaborated here.
[0156] S1221c'. All prior coincidence events corresponding to the prior down - sampled response lines are divided into prior prompt coincidence events and prior delayed coincidence events.
[0157] In some embodiments, the prior prompt coincidence events are acquired based on the prompt coincidence time window.
[0158] In some specific examples, the prompt coincidence time window corresponding to the prior prompt coincidence events is set to 2ns - 10ns. When the prompt coincidence time window exceeds 10ns, it may occur that the background signal in the subsequently obtained prior one - dimensional prompt coincidence energy spectrum is more than the background signal in the delayed coincidence energy spectrum, resulting in the inability to cancel out, generating corresponding error interference and affecting the accuracy of the scatter correction result. Setting the prompt coincidence time window within the range of 2ns - 10ns helps to eliminate the interference of the background signal in the subsequently obtained prior one - dimensional prompt coincidence energy spectrum.
[0159] In some embodiments, corresponding to step S1221a, the delay coincidence time window of the prior delay coincidence event is set to 2 ns - 10 ns. Specifically, the range of the delay coincidence time window and the immediate coincidence time window may be the same or different.
[0160] S1221d’, obtain the prior one-dimensional immediate energy spectrum based on the prior immediate coincidence events corresponding to the downsampled response lines; obtain the prior one-dimensional delay energy spectrum based on the prior delay coincidence events corresponding to the downsampled response lines.
[0161] Regarding how to obtain the prior one-dimensional immediate energy spectrum and the prior one-dimensional delay energy spectrum based on the downsampled response lines, reference can be specifically made to steps S1211 - S1213 and S1221a - S1221c, which will not be elaborated here.
[0162] S1221e’, obtain the difference between the prior one-dimensional immediate energy spectrum and the prior one-dimensional delay energy spectrum, perform normalization processing on this difference, and obtain the one-dimensional intermediate energy spectrum Y psf 。
[0163] In some embodiments, for the normalization processing of the difference between the prior one-dimensional immediate energy spectrum and the prior one-dimensional delay energy spectrum, reference can be made to the prior art, which will not be elaborated here.
[0164] S1221f’, based on the one-dimensional intermediate energy spectrum Y psf Obtain the blurred response matrix moment p.
[0165] In some embodiments, for obtaining the blurred response matrix p based on the one-dimensional intermediate energy spectrum Y psf Reference can be made to the prior art, which will not be elaborated here.
[0166] S1222, based on the two obtained one-dimensional detection energy spectra and the corresponding relationship function, use the maximum likelihood - expectation maximization iterative algorithm and the initial iteration values corresponding to the obtained downsampled response lines to obtain the two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
[0167] In some embodiments, the one-dimensional gamma photon energy spectrum refers to the energy spectrum composed of the energy of the gamma photons themselves.
[0168] In some embodiments, step S1222 includes:
[0169] S122210, based on the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, obtain the iterative formula of the maximum likelihood - expectation maximization iterative algorithm.
[0170] In some embodiments, the iterative formula of the maximum likelihood - expectation maximization iterative algorithm is:
[0171]
[0172] Among them, Y i in Y represents the one-dimensional detection energy spectrum obtained by scanning the target object corresponding to a single downsampled response line. i represents the i-th vertical bar region (i.e., "bin") in the energy spectrum, and b and j respectively represent the b-th and j-th vertical bar regions in the energy spectrum. k is the number of iterations, and p ij and p ib in p is the blurred response matrix corresponding to this downsampled response line. p ij is the value of the i-th row and j-th column of the blurred response matrix corresponding to this downsampled response line. p ib is the value of the i-th row and b-th column of the blurred response matrix corresponding to this downsampled response line, which is the blurred response matrix corresponding to this downsampled response line; r i represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in this downsampled response line; is the one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region obtained in the (k + 1)-th iteration; respectively represent the one-dimensional gamma photon energy spectra corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration.
[0173] Since the iteration formula of the MLEM iteration algorithm has no constraint of a penalty term, the selection of the initial iteration value of Equation (2) is very important.
[0174] In some embodiments, obtaining the initial iteration value includes:
[0175] S122211, obtaining the global initial scattered photon energy spectrum.
[0176] In some embodiments, "global" refers to an object obtained using all coincidence response lines corresponding to the detection data. For example, the global initial scattered photon energy spectrum refers to the initial scattered photon energy spectrum obtained using all coincidence response lines corresponding to the detection data.
[0177] In some embodiments, step S132211 includes:
[0178] S122211a, obtaining the coincidence response lines based on the detection data.
[0179] In some specific embodiments, step S122211a may refer to step S112, which will not be elaborated here.
[0180] S122211b, dividing all coincidence events corresponding to all coincidence response lines into prompt coincidence events and random coincidence events.
[0181] In some embodiments, how to specifically operate step S122211b can refer to the prior art, which will not be elaborated here.
[0182] S122211c obtains a global one-dimensional prompt energy spectrum based on prompt coincidence events and a global one-dimensional delayed energy spectrum based on delayed coincidence events.
[0183] In some embodiments, the specific operation of S122211c can be analogously referred to step 1221d’, which will not be elaborated here.
[0184] S122211d subtracts the global one-dimensional delayed energy spectrum corresponding to the delayed coincidence events from the global one-dimensional prompt energy spectrum corresponding to the prompt coincidence events on all response lines to obtain a global one-dimensional undelayed energy spectrum.
[0185] In some embodiments, the global one-dimensional prompt energy spectrum refers to a one-dimensional prompt energy spectrum obtained using all coincidence response lines corresponding to the detection data, and the global one-dimensional delayed energy spectrum refers to a one-dimensional delayed energy spectrum obtained using all coincidence response lines corresponding to the detection data.
[0186] S122211e obtains an objective function for the global initial scattered photon energy spectrum based on the global one-dimensional undelayed energy spectrum.
[0187] In some specific embodiments, the objective function for the global initial scattered photon energy spectrum is:
[0188]
[0189] wherein, S i in S is the global initial scattered photon energy spectrum; i represents the i-th vertical bar region in the energy spectrum; C i in C is the global one-dimensional undelayed energy spectrum, obtained based on the target object, U i in U is the global one-dimensional intermediate energy spectrum Y obtained from the scan of the entire PET system psf ; the global one-dimensional intermediate energy spectrum Y psf is obtained based on the phantom; hel and heu are the lower and upper thresholds of a set high-energy window (specifically set based on experience, usually 511 < hel < heu, and heu is generally taken as 600, 650 or 700, etc.), is the stretching coefficient.
[0190] In some embodiments, the acquisition of the global one-dimensional intermediate energy spectrum Y psf can refer to step S1221e’. Specifically, the one-dimensional intermediate energy spectrum obtained using all response lines corresponding to the detection data obtained based on the target object is the global one-dimensional intermediate energy spectrum Y psf .
[0191] It should be understood that almost all photons in the high-energy window of C i are unscattered photons, and the number of scattered photons is very small. By the total number of photons in the high-energy window, U iBy stretching, the energy spectrum of unscattered photons in the complete energy window can be estimated. Then, subtracting the estimated energy spectrum of unscattered photons from the global one-dimensional undelayed energy spectrum can obtain the global initial scattered photon energy spectrum.
[0192] It should be understood that the prerequisite for the success of Equation (3) is that C i contains a sufficient number of photons. According to prior information, the number of photons on all coincidence response lines is sufficient, while the number of photons on a single response line or a single downsampled response line is generally insufficient. If only the number of photons on a single response line or a single downsampled response line is used, S i will be affected by statistical noise. If the energy spectrum of a single response line or a single downsampled response line is directly obtained based on Equation (3), due to being severely affected by statistical noise, the obtained S i has very poor accuracy. In the embodiments of the present application, since the global one-dimensional undelayed energy spectrum C i is obtained based on all coincidence response lines corresponding to the detection data, the number of photons in C i is usually sufficient, which can improve the accuracy of the obtained S i .
[0193] In some specific examples, as can be seen from the above, U i is the global one-dimensional intermediate energy spectrum Y psf obtained by scanning the entire PET system. To obtain U i , reference can be made to the acquisition of the energy spectrum Y in the steps psf .
[0194] In some embodiments, the prompt energy spectrum obtained based on all response lines corresponding to prior detection data can be called the global prior prompt energy spectrum, and the obtained delayed energy spectrum can be called the global prior delayed energy spectrum.
[0195] S122212. Perform deblurring processing on the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value.
[0196] In some embodiments, the global intermediate initial iteration value refers to the intermediate initial iteration value obtained based on all response lines corresponding to the detection data.
[0197] In some embodiments, step S122212 includes:
[0198] S122212a. Obtain the deblurring processing iteration formula based on the maximum likelihood-expectation maximization iteration algorithm.
[0199] In some specific embodiments, the deblurring processing iteration formula is:
[0200]
[0201] Among them, Si where S in it represents the global one-dimensional detection energy spectrum obtained by scanning the target object, i represents the i-th vertical bar region (i.e., "bin") in the energy spectrum, b and j respectively represent the b-th and j-th vertical bar regions in the energy spectrum, k is the number of iterations, and P ij , P ib in it, P is the global fuzzy response matrix, and P ij is the value of the i-th row and j-th column of the global fuzzy response matrix, and P ib is the value of the i-th row and b-th column of the global fuzzy response matrix; r i represents the global one-dimensional delayed energy spectrum corresponding to all random coincidence events; is the global one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region obtained in the (k + 1)-th iteration; respectively represent the global one-dimensional gamma photon energy spectra corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration;
[0202]
[0203] where largeconstant represents a very large constant, such as including but not limited to 1000000, etc.
[0204] In the embodiments of the present application, by introducing M into formula (4) j , the S in this formula i can only include scattered photons, that is, after the blurring process of this formula, only the scattered part has a value, and the part corresponding to j < 511 keV is the scattered part.
[0205] In some embodiments, the global fuzzy response matrix refers to the fuzzy response matrix based on all response lines corresponding to the detection data.
[0206] In some specific examples, obtaining the global fuzzy response matrix P ij includes:
[0207] obtaining prior detection data based on a phantom (such as a point source or a line source, etc.);
[0208] obtaining all prior coincidence response lines based on the prior detection data;
[0209] dividing all prior coincidence events corresponding to all prior coincidence response lines into global prior prompt coincidence events and global prior delayed coincidence events;
[0210] obtaining a global prior one-dimensional prompt energy spectrum based on the global prior prompt coincidence events; obtaining a global prior one-dimensional delayed energy spectrum based on the global prior delayed coincidence events;
[0211] obtaining the difference between the global prior prompt energy spectrum and the global prior delayed energy spectrum;
[0212] Normalize the difference to obtain the global prior intermediate energy spectrum Y psf1 ;
[0213] Based on the global prior intermediate energy spectrum Y psf1 Obtain the global fuzzy response matrix P ij .
[0214] In some embodiments, the obtaining of the global prior intermediate energy spectrum Y psf1 can refer to S1221e'. Specifically, the one-dimensional intermediate energy spectrum obtained from all response lines corresponding to the prior detection data obtained based on the phantom is the global prior intermediate energy spectrum Y psf1 .
[0215] In some specific embodiments, how to normalize the difference to obtain the global prior intermediate energy spectrum Y psf1 ; and how to obtain the global fuzzy response matrix P psf1 based on the global prior intermediate energy spectrum Y ij can specifically refer to the prior art and will not be elaborated here.
[0216] S122212b, set the prior initial iteration value and the number of iterations, and perform iterations through the deblurring iteration formula until the number of iterations is reached to obtain the scattered part of the global initial energy spectrum.
[0217] In some specific examples, X 0 can be set to 1 and substituted into Equation (4) for several iterations, such as including but not limited to 50 - 500 times, to obtain the scattered part X init,sc of the global initial energy spectrum.
[0218] S122212c, estimate the unscattered part of the global initial energy spectrum based on the scattered part of the global initial energy spectrum.
[0219] In some embodiments, when X init,sc is obtained, in order to estimate the unscattered part of the global initial energy spectrum, that is, in order to obtain the value of the global initial energy spectrum at j = 511 keV, the unscattered part of the global initial energy spectrum can be estimated based on the estimation function, and the estimation function is:
[0220]
[0221] where represents the estimation function segment of the global intermediate initial iteration value function corresponding to j = 511 keV, κ is the ratio of unscattered photons to scattered photons; here P is the global fuzzy response matrix, and P ij is the value of the i-th row and j-th column of the global fuzzy response matrix; δ is a hyperparameter.
[0222] In some specific examples, the ratio κ of unscattered photons to scattered photons can be specifically estimated by the following formula:
[0223]
[0224] where sum() represents summation, and C and S are the C i and S i .
[0225] S122212d, based on the scattered part of the global initial energy spectrum and the unscattered part of the global initial energy spectrum, obtain the global intermediate initial iteration value.
[0226] In some specific embodiments, since gamma photons are mainly used in the PET system and the energy of gamma photons is 511 keV, therefore, when the energy is greater than 511 keV, the value of the corresponding global initial energy spectrum is 0. Based on this, the global intermediate initial iteration value function is X init :
[0227]
[0228] where δ is a hyperparameter, P has the same meaning as P ij above, both are the global blur response matrix, κ is the ratio of unscattered photons to scattered photons; δ is a hyperparameter.
[0229] In the embodiments of the present application, the S obtained by formula (3) i may have errors. The reason is that formula (3) assumes that there are no scattered photons in the high-energy window of C i , but in fact, there may still be scattered photons in the high-energy window of C i , resulting in U i being pulled up by the stretching coefficient, S i decreasing, then κ obtained based on formula (7) will be overestimated, and it can be balanced by the hyperparameter δ. Specifically, by setting δ to a number less than 1 to balance the overestimation effect of κ. It is found through experiments that good results can be obtained when δ takes values between 0.3 and 0.5, and basically no adjustment is required. On the other hand, S i is the global initial scattered photon energy spectrum, that is, the energy spectrum obtained based on all response lines detected by the target object. Response lines that do not pass through the target object are considered all scattering events. Response lines that pass through the target object have both scattering events and true events, that is, unscattered events. Therefore, the proportion of scattering events in the response lines passing through the target object is lower than the proportion of scattering events in all response lines. A more accurate unscattered part can also be obtained by balancing with the hyperparameter δ.
[0230] S122213, stretch the global intermediate initial iteration value to obtain the initial iteration value corresponding to each downsampled response line.
[0231] In some embodiments, step S122212 includes:
[0232] Stretching the global intermediate initial iteration value based on a stretching function to obtain the initial iteration value corresponding to each downsampled response line.
[0233] In some specific examples, the stretching function is:
[0234]
[0235] where X 0 is the initial iteration value, sum() represents summation, P is the global fuzzy response matrix, and X init is the global intermediate initial iteration value; Y is the one-dimensional detection energy spectrum obtained by scanning the target object corresponding to a single downsampled response line, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in a single response line.
[0236] In the embodiments of the present application, through stretching, the initial iteration value X matching the number of single events included in the downsampled response line can be obtained 0 .
[0237] S122220. Set the number of iterations, and based on the obtained initial iteration value, perform iterations through the iteration formula of the maximum likelihood-expectation maximization iteration algorithm until the set number of iterations is reached, to obtain the two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
[0238] In some specific embodiments, the number of iterations can be set according to specific needs, for example, iterating 20 - 60 times. Generally, a better one-dimensional gamma photon energy spectrum can be obtained within 60 iterations. If the number of iterations is too large, the effect will deteriorate.
[0239] S1223. Estimate the probability density functions of the two scattered photons and the probability density functions of the two unscattered photons corresponding to the corresponding downsampled response lines based on the two one-dimensional gamma photon energy spectra.
[0240] In some embodiments, step S1223 includes:
[0241] Step S12231. Divide the one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part.
[0242] In some embodiments, the obtained one-dimensional gamma photon energy spectrum can be segmented based on the energy value when the gamma photon does not scatter. The part less than this energy value is the scattered part, and the part equal to this energy value is the unscattered part.
[0243] In some specific examples, based on the energy of 511 keV of the gamma photon that does not scatter, the obtained one-dimensional gamma photon energy spectrum, i.e., X, can be divided into a scattered part X sc and an unscattered part Xus , as follows:
[0244]
[0245] Step S12232, perform restoration and deblurring processing on the scattered part X and the unscattered part X of the two one-dimensional gamma photon energy spectra to obtain the two scattered parts and the two unscattered parts of the corresponding two one-dimensional detection energy spectra of the downsampled response lines. sc and the unscattered part X us to obtain the two scattered parts and the two unscattered parts of the corresponding two one-dimensional detection energy spectra of the downsampled response lines.
[0246] In some embodiments, the scattered part and the unscattered part of the one-dimensional detection energy spectrum can be calculated based on X and X obtained from Equation (7): sc and X us The scattered part and the unscattered part of the one-dimensional detection energy spectrum are calculated as follows:
[0247] S = PX sc ;
[0248] U = PX us ;
[0249] where S represents the scattered part of the one-dimensional detection energy spectrum, U represents the unscattered part of the one-dimensional detection energy spectrum, and P is the blurring response matrix corresponding to the corresponding downsampled response line.
[0250] Step S12233, perform normalization processing on the two scattered parts and the two unscattered parts of the two one-dimensional detection energy spectra, and the obtained normalized values are used to estimate the probability density functions of the two scattered photons and the two unscattered photons corresponding to the corresponding downsampled response lines.
[0251] In some embodiments, the two-dimensional prompt coincidence energy histogram and the two-dimensional delayed coincidence energy histogram corresponding to each downsampled response line can be obtained with reference to the prior art and will not be elaborated here.
[0252] In some embodiments, after obtaining the two-dimensional prompt coincidence energy histogram, the two-dimensional delayed coincidence energy histogram, the probability density function of the scattered photons, and the probability density function of the unscattered photons corresponding to each downsampled response line, the estimation objective function can be obtained.
[0253] In some specific embodiments, the estimation objective function is:
[0254]
[0255] where σ 0 respectively represent the number of true coincidence events and do not require subsequent processing; σ 1 + σ 2 + σ 3 represents the number of scattered coincidence events SC, which is an unknown quantity to be solved; the coefficient α 1,m , α2,n , β 1,m , β 2,n The calculation method of is as follows:
[0256]
[0257]
[0258] Wherein, and are the probability density functions of unscattered photons corresponding to a detection module in the downsampled response line, and are the probability density functions of scattered photons corresponding to another detection module in the downsampled response line.
[0259] S130, solve the estimated objective function to obtain the number of scattered coincidence events corresponding to each downsampled response line.
[0260] In some embodiments, the value of σ obtained by calculating formula (11) 1 +σ 2 +σ 3 is the number of scattered coincidence events corresponding to the corresponding downsampled response line.
[0261] In the embodiments of the present application, the probability density functions of unscattered photons and in the estimated objective function are not necessarily the same, and the probability density functions of scattered photons and are not necessarily the same either, but are specifically determined according to the two detection modules corresponding to the downsampled response line. Based on this, the number of scattered coincidence events corresponding to the downsampled response line obtained is more accurate.
[0262] S140, perform upsampling processing on the downsampled response line to obtain an upsampled response line.
[0263] In some embodiments, performing upsampling processing on the downsampled response line includes:
[0264] Use bilinear interpolation method, 4D linear interpolation method or 5D linear interpolation method to process the downsampled response line to obtain the upsampled response line corresponding to each downsampled response line. The specific operation can refer to the prior art and will not be elaborated here.
[0265] S150, calculate the number of scattered coincidence events corresponding to each upsampled response line.
[0266] In some embodiments, the scintillation crystals included in the detection module in the PET system are called upsampled crystals, and the set of scintillation crystals corresponding to the downsampled response line obtained by downsampling processing is called downsampled crystals. For example, as Figure 2As shown, the downsampling crystal to which the upsampling crystal belongs is denoted as "downsampling central crystal c"; the downsampling crystal that is axially closest to the upsampling crystal is denoted as "downsampling axially adjacent crystal a"; the downsampling crystal that is radially closest to the upsampling crystal is denoted as "downsampling radially adjacent crystal t"; and the downsampling crystal that is diagonally closest to the upsampling crystal is denoted as "downsampling opposite crystal ta".
[0267] As Figure 2 shown, let the axial length of the downsampling crystal be D a , and the radial length be D t . An upsampling response line has two upsampling crystals, which are distinguished by superscripts "i" and "j" respectively. Taking the scintillation crystal with superscript "i" as an example, the radial and axial distances from the center of this upsampling crystal to the center of "central crystal c" are denoted as Then the weights corresponding to the four downsampling crystals are as follows:
[0268] "Downsampling central crystal c":
[0269] "Downsampling axially adjacent crystal a":
[0270] "Downsampling radially adjacent crystal t":
[0271] "Downsampling opposite crystal ta":
[0272] Similarly, for the four downsampling crystals corresponding to the upsampling crystal with superscript "j", there are also:
[0273]
[0274]
[0275]
[0276]
[0277] In some specific embodiments, when using the 4D linear interpolation method to process the downsampling response line, the following formula is used to calculate the number of scattered coincidence events corresponding to the upsampling response line:
[0278]
[0279] where represents the number of scattered coincidence events of the upsampling response line formed by upsampling crystals i and j after 4D linear interpolation; S iDenote the set of downsampled central crystals, axially adjacent downsampled crystals, radially adjacent downsampled crystals, and relatively downsampled crystals corresponding to the upsampled crystal i; K and L denote the numbers of the crystals corresponding to the downsampled response lines; Characterize the number of scatter coincidence events on the downsampled response line formed by the downsampled crystals K and L; Denote the weight of the downsampled crystal K for the upsampled crystal i; Denote the weight of the downsampled crystal L for the upsampled crystal j; W tot Is a normalization factor, characterized as;
[0280]
[0281] Where N is the number of upsampled response lines included in one downsampled response line.
[0282] In some other embodiments, step S150 includes:
[0283] Calculate the ratio of the total number of the scatter coincidence events on each downsampled response line to the total number of coincidence events;
[0284] Calculate the number of scatter coincidence events corresponding to the corresponding upsampled response line according to the number of coincidence events on each upsampled response line included in each downsampled response line and the ratio.
[0285] In some embodiments, a solution that combines the maximum likelihood-expectation maximization iterative algorithm with the moment estimation method is adopted. Based on the maximum likelihood-expectation maximization iterative algorithm, the moment estimation method is used to obtain the estimated objective function corresponding to each downsampled response line. On the one hand, the maximum likelihood-expectation maximization iterative algorithm (MLEM) is used for scatter correction. It does not depend on the activity image and the attenuation image, making the scatter correction easy to implement. And this method does not require the hyperparameter β, so there is no need to set the hyperparameter β within a certain range, that is, it is hardly affected by the hyperparameter β. Compared with the existing one-step later iterative algorithm OSL-EDR (OneStep Later eliminate detector response) and the optimization transfer iterative algorithm (Optimization transfer eliminate detector response, abbreviated as OT-EDR), both OSL-EDR and OT-EDR need to select an appropriate hyperparameter β to obtain better results, and the optimal β needs some pre-experiments to adjust, and different sizes of phantoms require different β. The embodiment of the present application uses the MLEM algorithm for scatter correction and hardly needs to adjust the hyperparameter β, making the scatter correction more convenient. On the other hand, the unscattered photon probability density function and in the estimated objective function are not necessarily the same, and the scattered photon probability density function and are not necessarily the same, but are specifically determined according to the two detection modules corresponding to the downsampled response line, and the number of scattered coincidence events corresponding to the downsampled response line obtained is more accurate.
[0286] As can be seen from the above, according to some embodiments of the present application, various scattering situations are considered, such as external scattering such as random scattering and scattering events that do not pass through the target object, making the accuracy of the scatter correction result relatively high, thereby improving the quality of the subsequent reconstructed image. In addition, compared with the double-scattering simulation method and the Monte Carlo simulation method in the prior art, the scatter correction method according to the present application has the advantages of simple calculation process and small calculation amount.
[0287] Figure 3 It is an exemplary flowchart of the scatter correction method 200 shown in some embodiments of the present application.
[0288] Continue to refer to Figure 3 The scatter correction method 200 may include the following steps:
[0289] S210, obtaining response lines based on detection data;
[0290] S220. Based on the maximum likelihood - expectation maximization iterative algorithm, use the moment estimation method to obtain the estimated objective function corresponding to each response line, where the estimated objective function is obtained based on the two - dimensional prompt coincidence energy histogram, two - dimensional delayed coincidence energy histogram, probability density function of scattered photons, and probability density function of unscattered photons corresponding to each response line;
[0291] S230. Solve the estimated objective function to obtain the number of scattered coincidence events corresponding to each response line.
[0292] Different from the embodiment of the scattering correction method 100, after obtaining the response lines based on the detection data, the scattering correction method 100 performs down - sampling processing on the response lines to obtain down - sampled response lines, and then based on the maximum likelihood - expectation maximization iterative algorithm and combined with the matrix estimation method, after obtaining the number of scattered coincidence events corresponding to the down - sampled response lines, it performs up - sampling processing on the down - sampled response lines to obtain the number of scattered coincidence events corresponding to each response line of the down - sampled response lines; in this embodiment, the maximum likelihood - expectation maximization iterative algorithm and the matrix estimation method are directly used to process the response lines obtained based on the detection data to obtain the number of scattered coincidence events corresponding to each response line.
[0293] In the embodiments of the present application, the scattering correction method 200 can selectively combine the features of the scattering correction method 100 or other methods, and vice versa.
[0294] Figure 4 It is an exemplary flowchart of an image reconstruction method shown in some embodiments of the present application.
[0295] Continue to refer to Figure 4 , the image reconstruction method 300 may include the following steps:
[0296] S310. Obtain the number of scattered coincidence events corresponding to each response line based on the scattering correction method described in the above embodiments;
[0297] S320. Correct the scattering events based on the number of scattered coincidence events to obtain a reconstructed image.
[0298] In some specific embodiments, based on the number of scattered coincidence events, an iterative image reconstruction algorithm can be used for image reconstruction to obtain a reconstructed image.
[0299] Step S320 can specifically refer to the prior art and will not be elaborated here.
[0300] Figure 5 It is an exemplary module diagram of a scattering correction device shown in some embodiments of the present application. As Figure 5 shown, the scattering correction device 400 may include:
[0301] A down-sampling response line acquisition module 410, configured to acquire down-sampling response lines based on detection data;
[0302] An estimation module 420, configured to obtain an estimated objective function corresponding to each down-sampling response line based on the maximum likelihood-expectation maximization iterative algorithm and using the moment estimation method, where the estimated objective function is obtained based on the two-dimensional prompt coincidence energy histogram, two-dimensional delayed coincidence energy histogram, probability density function of scattered photons, and probability density function of unscattered photons corresponding to each down-sampling response line;
[0303] A scattered coincidence event number acquisition module 430, configured to solve the estimated objective function to obtain the number of scattered coincidence events corresponding to each down-sampling response line;
[0304] An up-sampling module 440, configured to perform up-sampling processing on the down-sampling response lines to obtain up-sampling response lines;
[0305] A scattered coincidence event number calculation module 450, configured to calculate the number of scattered coincidence events corresponding to each up-sampling response line.
[0306] In some embodiments, the estimation module 420 includes:
[0307] An energy spectrum acquisition module, configured to obtain a two-dimensional prompt coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, and two one-dimensional detection energy spectra corresponding to each down-sampling response line based on each coincidence event corresponding to each down-sampling response line;
[0308] A probability density acquisition module, configured to obtain the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding down-sampling response line based on the two one-dimensional detection energy spectra using the maximum likelihood-expectation maximization iterative algorithm.
[0309] In some embodiments, the estimation module 420 further includes:
[0310] A correspondence acquisition module, configured to obtain the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum;
[0311] A gamma photon energy spectrum acquisition module, configured to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra based on the two obtained one-dimensional detection energy spectra, the correspondence function, using the maximum likelihood-expectation maximization iterative algorithm, and the initial iteration value corresponding to the corresponding down-sampling response line;
[0312] A probability density estimation module, configured to estimate the probability density function of two scattered photons and the probability density function of two unscattered photons corresponding to the corresponding down-sampling response line based on the two one-dimensional gamma photon energy spectra.
[0313] In some specific embodiments, the probability density estimation module includes:
[0314] A partitioning module configured to partition a one-dimensional gamma photon energy spectrum into a scattered part and an unscattered part;
[0315] A restoring blurred response module configured to perform a restoring blurred response process on the scattered part and the unscattered part of two one-dimensional gamma photon energy spectra to obtain two scattered parts and two unscattered parts corresponding to two one-dimensional detection energy spectra corresponding to the corresponding downsampled response lines;
[0316] A density estimation module configured to perform a normalization process on the two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra, and estimate the probability density functions of two scattered photons and the probability density functions of two unscattered photons corresponding to the corresponding downsampled response lines based on the obtained normalized values.
[0317] In some specific embodiments, the gamma photon energy spectrum acquisition module includes a maximum likelihood-expectation maximization iterative algorithm iterative formula acquisition module configured to obtain a maximum likelihood-expectation maximization iterative algorithm iterative formula based on the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum;
[0318] An initial iteration value acquisition module configured to obtain an initial iteration value;
[0319] An iteration module configured to set iteration conditions and perform iterations based on the obtained initial iteration value through the maximum likelihood-expectation maximization iterative algorithm iterative formula until the set iteration conditions are reached, to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
[0320] In some specific examples, the initial iteration value acquisition module includes a global initial scattered photon energy spectrum configured to obtain a global initial scattered photon energy spectrum; an intermediate initial iteration value acquisition module configured to perform a deblurring process on the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value; and a stretching module configured to stretch the global intermediate initial iteration value to obtain initial iteration values corresponding to each downsampled response line.
[0321] In the embodiments of the present application, the scatter correction device 400 can selectively combine the features of the scatter correction method 100 or 200 or other methods, and vice versa.
[0322] Figure 6 is an exemplary module diagram of a scatter correction device according to some embodiments of the present application. As Figure 6 shown, the scatter correction device 500 may include:
[0323] A response line acquisition module 510 configured to obtain response lines based on detection data;
[0324] An estimation module 520 configured to obtain an estimated objective function corresponding to each response line by using a moment estimation method based on a maximum likelihood - expectation maximization iterative algorithm, where the estimated objective function is obtained based on a two - dimensional prompt - coincidence energy histogram, a two - dimensional delayed - coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each response line;
[0325] A scattered - coincidence event number acquisition module 530 configured to solve the estimated objective function to obtain the number of scattered - coincidence events corresponding to each response line.
[0326] Different from the embodiment of the scatter correction device 400, after the scatter correction device 400 obtains response lines based on detection data, it performs downsampling processing on the response lines to obtain downsampled response lines, and then based on the maximum likelihood - expectation maximization iterative algorithm and combined with a matrix estimation method, after obtaining the number of scattered - coincidence events corresponding to the downsampled response lines, it performs upsampling processing on the downsampled response lines to obtain the number of scattered - coincidence events corresponding to each response line of the downsampled response lines; in this embodiment, the maximum likelihood - expectation maximization iterative algorithm and combined with a matrix estimation method are used to directly process the response lines obtained based on detection data to obtain the number of scattered - coincidence events corresponding to each response line.
[0327] In the embodiments of the present application, the scatter correction device 500 can be used to implement the scatter correction method 200 or the methods described in other embodiments herein, and can selectively combine the features of the scatter correction device 400 or other devices, and can also selectively combine the features of the scatter correction method 100 or 200 or other methods, and vice versa.
[0328] Figure 7 It is an exemplary module diagram of an image reconstruction device shown in some embodiments of the present application. As Figure 7 shown, the image reconstruction device 600 may include:
[0329] A scattered - coincidence event number acquisition module 610 configured to obtain the number of scattered - coincidence events corresponding to each response line based on the scatter correction method described in any of the above embodiments;
[0330] A reconstruction module 620 configured to correct scatter events based on the number of scattered - coincidence events to obtain a reconstructed image.
[0331] In the embodiments of the present application, the image reconstruction device 600 can be used to implement the image reconstruction method 300 or the methods described in other embodiments herein, and can selectively combine the features of the image reconstruction method 300 or other methods, and vice versa.
[0332] In some embodiments of the present application, the image reconstruction device 600 may also include components or features of the scatter correction devices 400 and 500 in a non - contradictory manner, and vice versa.
[0333] In some embodiments, the present application also provides a digital device, which includes: the device described in any one of the above - mentioned embodiments.
[0334] In some embodiments, the present application also provides a digital device, which may include a memory, a processor, and a computer program stored on the memory and executable on the processor. When the computer program is executed by the processor, it can implement the steps of the method described in any one of the above - mentioned embodiments.
[0335] Although not shown, in some embodiments, a computer - readable storage medium is also provided, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the method described in any one of the above - mentioned embodiments. The computer program includes each program module / unit that constitutes the device according to the embodiments of the present application. When the computer program constituted by each program module / unit is executed, it can implement the functions corresponding to each step in the method described in the above - mentioned embodiments. The computer program can also run on an electronic device as described in the embodiments of the present application.
[0336] The basic concepts have been described herein. Obviously, for those skilled in the art, the above - detailed disclosure is only an example and does not constitute a limitation to the present application. Although not explicitly stated here, those skilled in the art may make various modifications, improvements, and corrections to the present application. Such modifications, improvements, and corrections are proposed in the present application, so such modifications, improvements, and corrections still belong to the spirit and scope of the exemplary embodiments of the present application.
[0337] Meanwhile, the present application uses specific terms to describe the embodiments of the present application. Such as "one embodiment", "an embodiment", and / or "some embodiments" mean a certain feature, structure, or characteristic related to at least one embodiment of the present application. Therefore, it should be emphasized and noted that the "one embodiment" or "an embodiment" or "an alternative embodiment" mentioned twice or more at different positions in the present application does not necessarily refer to the same embodiment. In addition, certain features, structures, or characteristics in one or more embodiments of the present application can be appropriately combined.
[0338] In addition, those skilled in the art can understand that various aspects of the present application can be illustrated and described by several patentable types or situations, including any new and useful process, machine, product, or combination of substances, or any new and useful improvement thereof. Accordingly, various aspects of the present application can be executed entirely by hardware, entirely by software (including firmware, resident software, microcode, etc.), or by a combination of hardware and software. The above hardware or software can all be referred to as "data block", "module", "engine", "unit", "component", or "system". In addition, various aspects of the present application may be embodied as a computer product located in one or more computer-readable media, which includes computer-readable program code.
[0339] A computer storage medium may contain a propagated data signal containing computer program code, such as on a baseband or as part of a carrier wave. This propagated signal may have various forms of manifestation, including electromagnetic form, optical form, etc., or a suitable combination thereof. A computer storage medium can be any computer-readable medium other than a computer-readable storage medium, which can be connected to an instruction execution system, apparatus, or device to implement communication, propagation, or transmission for use of a program. The program code located on the computer storage medium can be propagated through any suitable medium, including radio, cable, fiber optic cable, RF, or similar media, or any combination of the above media.
[0340] The computer program code required for the operation of each part of the present application can be written in any one or more programming languages, including object-oriented programming languages such as Java, Scala, Smalltalk, Eiffel, JADE, Emerald, C++, C#, VB.NET, Python, etc., conventional procedural programming languages such as C language, Visual Basic, Fortran 2003, Perl, COBOL 2002, PHP, ABAP, dynamic programming languages such as Python, Ruby, and Groovy, or other programming languages, etc. The program code can run entirely on the user's computer, or run as an independent software package on the user's computer, or run partially on the user's computer and partially on a remote computer, or run entirely on a remote computer or server. In the latter case, the remote computer can be connected to the user's computer through any network form, such as a local area network (LAN) or a wide area network (WAN), or connected to an external computer (for example, through the Internet), or in a cloud computing environment, or used as a service such as software as a service (SaaS).
[0341] In addition, unless explicitly stated in the claims, the order of the processing elements and sequences, the use of numerical and alphabetical characters, or the use of other names in this application are not used to limit the order of the processes and methods of this application. Although some currently useful embodiments of the invention are discussed through various examples in the above disclosure, it should be understood that such details are for illustrative purposes only. The appended claims are not limited to the disclosed embodiments. On the contrary, the claims are intended to cover all modifications and equivalent combinations that conform to the essence and scope of the embodiments of this application. For example, although the system components described above can be implemented by hardware devices, they can also be implemented only through software solutions, such as installing the described system on existing servers or mobile devices.
[0342] Similarly, it should be noted that, in order to simplify the description of this application disclosure and thus help the understanding of one or more embodiments of the invention, in the foregoing description of the embodiments of this application, sometimes multiple features are merged into one embodiment, drawing, or description thereof. However, this disclosure method does not mean that the features required by the subject matter of this application are more than those mentioned in the claims. In fact, the features of the embodiments are fewer than all the features of the single embodiments disclosed above.
[0343] In some embodiments, numbers describing the components and attribute quantities are used. It should be understood that such numbers used in the description of the embodiments are modified by the modifiers "about", "approximately", or "substantially" in some examples. Unless otherwise stated, "about", "approximately", or "substantially" indicate that the stated numbers allow a variation of ±20%. Accordingly, in some embodiments, the numerical parameters used in the specification and claims are approximate values, and these approximate values may change according to the characteristics required by individual embodiments. In some embodiments, the numerical parameters should consider the specified significant digits and adopt the method of retaining general digits. Although the numerical ranges and parameters used in some embodiments of this application to confirm the breadth of their scope are approximate values, in specific embodiments, such numerical settings are as precise as possible within the feasible range.
[0344] For each patent, patent application, patent application publication, and other materials cited in this application, such as articles, books, specifications, publications, documents, etc., their entire contents are hereby incorporated into this application by reference. This excludes the application history files that are inconsistent with or conflict with the content of this application, and also excludes the files that limit the broadest scope of the claims of this application (currently or subsequently appended to this application). It should be noted that if there are inconsistencies or conflicts between the descriptions, definitions, and / or uses of terms in the attached materials of this application and the content described in this application, the descriptions, definitions, and / or uses of terms in this application shall prevail.
[0345] Finally, it should be noted that the above are only exemplary embodiments of the present application and are not used to limit the present application. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or perform equivalent replacements for some of the technical features. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.
Claims
1. A scatter correction method, characterized in that: The scatter correction method comprises: Obtaining a downsampled response line based on the detection data; Based on the maximum likelihood-expectation maximization iterative algorithm, the estimated objective function corresponding to each downsampled response line is obtained using the moment estimation method, wherein the estimated objective function is obtained based on the two-dimensional instantaneous coincidence energy histogram, the two-dimensional delayed coincidence energy histogram, the probability density function of scattered photons, and the probability density function of unscattered photons corresponding to each downsampled response line; Solving the estimation objective function to obtain the number of scattering coincidence events corresponding to each downsampled response line; Performing upsampling processing on the downsampled response line to obtain an upsampled response line; Count the number of scatter coincidence events corresponding to each upsampled response line.
2. The scatter correction method according to claim 1, characterized in that: The down-sampling response line is obtained based on the detection data, including: Acquiring a coincidence response line based on the detection data; Downsampling is performed on the coincident response line to obtain a downsampled response line.
3. The scatter correction method according to claim 2, characterized in that: The detection data is acquired based on the target object.
4. The scatter correction method according to claim 1, characterized in that: Based on the maximum likelihood-expectation maximization iterative algorithm, the moment estimation method is used to obtain the estimated objective function corresponding to each downsampled response line, including: Based on each coincidence event corresponding to each down-sampling response line, a two-dimensional immediate coincidence energy histogram, a two-dimensional delayed coincidence energy histogram and two one-dimensional detection energy spectra corresponding to each down-sampling response line are obtained. Based on the two one-dimensional detection energy spectra, a maximum likelihood-expectation maximization iterative algorithm is used to obtain the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding down-sampling response line.
5. The scatter correction method according to claim 4, characterized in that: Based on each coincidence event corresponding to each down-sampled response line, two one-dimensional detection energy spectra are obtained, including: Split all matching events corresponding to the downsampled response line into single events; The single event is divided into two parts according to its position information, and each part corresponds to a one-dimensional detection energy spectrum; Based on the divided two-part single event, two one-dimensional detection energy spectra are obtained.
6. The scatter correction method according to claim 4, characterized in that: Based on the two one-dimensional detection energy spectra, a maximum likelihood-expectation maximization iterative algorithm is used to obtain a probability density function of scattered photons and a probability density function of unscattered photons corresponding to the corresponding downsampled response lines, including: Obtaining a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; Based on the two acquired one-dimensional detection energy spectra and the corresponding relationship function, a maximum likelihood-expectation maximization iterative algorithm and the acquired initial iteration values corresponding to the corresponding down-sampling response lines are used to acquire two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra; The probability density functions of two scattered photons and two unscattered photons corresponding to the corresponding down-sampled response lines are estimated based on the two one-dimensional gamma photon energy spectra.
7. The scatter correction method according to claim 6, characterized in that: The corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is: in, represents the expected value of the one-dimensional detection energy spectrum corresponding to a single down-sampled response line based on the scanning of the target object; X represents the one-dimensional gamma photon energy spectrum, p represents the fuzzy response matrix corresponding to the single down-sampled response line, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in the single response line.
8. The scatter correction method according to claim 7, characterized in that: The one-dimensional delayed energy spectrum corresponding to the random coincidence event is obtained by delay coincidence event estimation.
9. The scatter correction method according to claim 8, characterized in that: Obtaining a one-dimensional delay spectrum includes: Acquire all delay coincidence events in the downsampled response line based on the delay coincidence window; split each delay coincidence event into two corresponding delay single events; The delayed single event is divided into two parts according to the position information of the delayed single event, and each part corresponds to a one-dimensional delayed energy spectrum; Based on the divided two-part delayed single event, two one-dimensional delayed energy spectra are obtained; The location information of the delayed single event includes location information of a downsampling detection module that detects the delayed single event.
10. The scatter correction method according to claim 7, characterized in that: Obtaining the fuzzy response matrix p includes: Acquiring prior detection data based on the prosthesis; Obtaining a priori down-sampling response lines based on priori detection data; Divide all a priori coincidence events corresponding to the a priori downsampling response line into a priori immediate coincidence events and a priori delayed coincidence events; Based on the prior instant coincidence event corresponding to the down-sampled response line, a priori one-dimensional instant energy spectrum is obtained; based on the priori delayed coincidence event corresponding to the down-sampled response line, a priori one-dimensional delayed energy spectrum is obtained; Obtaining a difference between a priori one-dimensional instantaneous energy spectrum and a priori one-dimensional delayed energy spectrum, performing normalization processing on the difference, and obtaining a one-dimensional intermediate energy spectrum; Obtain the fuzzy response matrix based on the one-dimensional intermediate energy spectrum.
11. The scatter correction method according to claim 6, characterized in that: Obtaining two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra, comprising: Based on the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, the iterative formula of the maximum likelihood-expectation maximization iterative algorithm is obtained; The number of iterations is set, and based on the obtained initial iteration value, it is iterated to the set number of iterations through the maximum likelihood-expectation maximization iterative algorithm iterative formula to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
12. The scatter correction method according to claim 11, characterized in that: The maximum likelihood-expectation maximization iterative algorithm iteration formula is: Among them, Y i Where Y represents the one-dimensional detection energy spectrum corresponding to a single down-sampled response line based on scanning the target object, i represents the i-th vertical bar area in the energy spectrum, b and j represent the i-th and j-th vertical bar areas in the energy spectrum respectively, k is the number of iterations, p ij and p ib The p in the equation is the fuzzy response matrix corresponding to the down-sampled response line. ij is the value of the i-th row and j-th column of the fuzzy response matrix corresponding to the downsampled response line, p ib is the value of the i-th row and b-th column of the fuzzy response matrix corresponding to the downsampled response line; r i Indicates the one-dimensional delayed energy spectrum corresponding to the random coincidence event included in the down-sampled response line; The one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region of the energy spectrum obtained in the k+1-th iteration; Respectively represent the one-dimensional gamma photon energy spectra corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration.
13. The scatter correction method according to claim 11, characterized in that: Get the initial iteration value, including: Obtain the global initial scattered photon energy spectrum; Defuzzifying the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value; The global intermediate initial iteration value is stretched to obtain the initial iteration value corresponding to each down-sampling response line.
14. The scatter correction method according to claim 13, characterized in that: Get the global initial scattered photon energy spectrum, including: Acquiring a coincidence response line based on the detection data; Divide all the coincidence events corresponding to all the coincidence response lines into immediate coincidence events and random coincidence events; Based on the instant coincidence event, a global one-dimensional instant energy spectrum is obtained; based on the delayed coincidence event, a global one-dimensional delayed energy spectrum is obtained; The global one-dimensional instantaneous energy spectrum corresponding to the instantaneous coincidence event on all response lines is subtracted from the global one-dimensional delayed energy spectrum corresponding to the delayed coincidence event to obtain the global one-dimensional undelayed energy spectrum; Based on the global one-dimensional undelayed energy spectrum, an objective function is used to obtain the global initial scattered photon energy spectrum.
15. The scatter correction method according to claim 14, characterized in that: The objective function of the global initial scattered photon energy spectrum is: Among them, S i The S in is the global initial scattered photon energy spectrum; i represents the i-th vertical bar area in the energy spectrum, C i C is the global one-dimensional undelayed energy spectrum, obtained based on the target object, and U i U in the figure is the global one-dimensional intermediate energy spectrum obtained by the entire PET system scan, and the global one-dimensional intermediate energy spectrum is obtained based on the phantom; hel and heu are the lower and upper thresholds of the set high-energy window, is the stretch coefficient.
16. The scatter correction method according to claim 15, characterized in that: Get U i include: Acquiring prior detection data based on the prosthesis; Acquire all a priori coincident response lines based on a priori detection data; Divide all prior coincidence events corresponding to all prior coincidence response lines into global prior immediate coincidence events and global prior delayed coincidence events; Based on the global prior instant coincidence event, a global prior one-dimensional instant energy spectrum is obtained; based on the global prior delayed coincidence event, a global prior one-dimensional delayed energy spectrum is obtained; Get the difference U between the global prior instantaneous energy spectrum and the global prior delayed energy spectrum i .
17. The scatter correction method according to claim 13, characterized in that: Defuzzify the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value, including: Obtaining the defuzzification iterative formula based on the maximum likelihood-expectation maximization iterative algorithm; The a priori initial iteration value and iteration condition are set, and the defuzzification processing iteration formula is used for iteration until the iteration condition is met to obtain the scattered part of the global initial energy spectrum; estimating an unscattered portion of the global initial energy spectrum based on the scattered portion of the global initial energy spectrum; A global intermediate initial iteration value is obtained based on a scattered part of the global initial energy spectrum and an unscattered part of the global initial energy spectrum.
18. The scatter correction method according to claim 17, characterized in that: The defuzzification iterative formula is: Among them, S i The S in the above expression represents the global one-dimensional detection energy spectrum obtained by scanning the target object, i represents the i-th vertical bar area in the energy spectrum, b and j represent the i-th and j-th vertical bar areas in the energy spectrum respectively, k is the number of iterations, and P is the number of iterations. ij , P ib P in it is the global fuzzy response matrix, P ij is the value of the i-th row and j-th column of the global fuzzy response matrix, P ib is the value of the i-th row and b-th column of the global fuzzy response matrix; r i Represents the global one-dimensional delayed energy spectrum corresponding to all random coincidence events; is the global one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar region of the energy spectrum obtained in the k+1-th iteration; Respectively represent the global one-dimensional gamma photon energy spectrum corresponding to the b-th and j-th vertical bar regions of the energy spectrum obtained in the k-th iteration; Among them, largeconstant represents the constant.
19. The scatter correction method according to claim 18, characterized in that: Obtaining the global fuzzy response matrix P includes: Acquiring prior detection data based on the prosthesis; Acquire all a priori coincident response lines based on a priori detection data; Divide all prior coincidence events corresponding to all prior coincidence response lines into global prior immediate coincidence events and global prior delayed coincidence events; Based on the global prior instant coincidence event, a global prior one-dimensional instant energy spectrum is obtained; based on the global prior delayed coincidence event, a global prior one-dimensional delayed energy spectrum is obtained; Obtain the difference between the global a priori instantaneous energy spectrum and the global a priori delayed energy spectrum; The difference is normalized to obtain a global a priori intermediate energy spectrum; The global fuzzy response matrix is obtained based on the global prior intermediate energy spectrum.
20. The scatter correction method according to claim 17, characterized in that: The estimation of the unscattered part of the global initial energy spectrum is performed by an estimation function, which is: Where P is the global fuzzy response matrix; δ is a hyperparameter; κ is the ratio of unscattered photons to scattered photons; Among them, sum() means sum, C and S correspond to C i , S i .
21. The scatter correction method according to claim 17, characterized in that: The characterization function of the global intermediate initial iteration value is: Among them, δ is a hyperparameter, P is the global blur response matrix, and κ is the ratio of unscattered photons to scattered photons.
22. The energy screening method according to claim 13, characterized in that: The intermediate initial iteration values are stretched to obtain the initial iteration values corresponding to each downsampling response line, including: The global intermediate initial iteration value is stretched based on the stretching function to obtain the initial iteration value corresponding to each down-sampled response line, and the stretching function is: Among them, X 0 is the initial iteration value, sum() represents summation, P is the global fuzzy response matrix, X init is the global intermediate initial iteration value; Y is the one-dimensional detection energy spectrum corresponding to a single down-sampling response line obtained based on the target object scanning, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in a single response line.
23. The scatter correction method according to claim 6, characterized in that: Estimating the probability density functions of two scattered photons and two unscattered photons corresponding to the corresponding downsampled response lines based on the two one-dimensional gamma photon energy spectra includes: The one-dimensional gamma photon energy spectrum is divided into a scattered part and an unscattered part; Performing fuzzy response recovery processing on the scattered parts and unscattered parts of the two one-dimensional gamma photon energy spectra to obtain two scattered parts and two unscattered parts corresponding to the two one-dimensional detection energy spectra corresponding to the corresponding down-sampled response lines; The two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra are normalized, and the probability density functions of the two scattered photons and the two unscattered photons corresponding to the corresponding downsampled response lines are estimated based on the obtained normalized values.
24. The scatter correction method according to claim 23, characterized in that: The scattering part is characterized by: The unscattered portion is characterized by: in, is the scattering part, is the unscattered part.
25. The scatter correction method according to claim 23, characterized in that: The one-dimensional gamma photon energy spectrum is divided into scattered and unscattered parts, including: The one-dimensional gamma photon energy spectrum is divided based on the energy value of the gamma photon when it is not scattered. The part less than the energy value is the scattered part, and the part equal to the energy value is the unscattered part.
26. The scatter correction method according to claim 1, characterized in that: The estimation objective function is: Among them, σ0 represents the number of true coincidence events, which does not require subsequent processing; σ1+σ2+σ3 represents the number of scattered coincidence events SC, which is the unknown quantity to be calculated; the coefficient α 1,m ,α 2,n ,β 1,m ,β 2,n The calculation method is: in, and is the probability density function of unscattered photons corresponding to a detection module in the downsampled response line, and is the probability density function of scattered photons corresponding to another detection module in the down-sampling response line; The calculated σ1+σ2+σ3 value is the number of scattering coincidence events corresponding to the corresponding downsampled response line.
27. The scatter correction method according to claim 1, characterized in that: Perform upsampling on the downsampled response line, including: The down-sampled response lines are processed by a bilinear interpolation method, a 4D linear interpolation method or a 5D linear interpolation method to obtain up-sampled response lines corresponding to the down-sampled response lines.
28. The scatter correction method according to claim 27, characterized in that: When the down-sampled response line is processed by the 4D linear interpolation method, the following formula is used to calculate the number of scattering coincidence events corresponding to the up-sampled response line: in, Characterizes the number of scattering coincidence events of the upsampled response line composed of upsampled crystals i and j after 4D linear interpolation; S i represents the set of downsampling central crystal, downsampling axially adjacent crystals, downsampling radially adjacent crystals and downsampling relative crystals corresponding to the upsampling crystal i; K and L represent the numbers of the crystals corresponding to the downsampling response lines; Characterize the number of scattering coincidence events on the down-sampling response line formed by the down-sampling crystals K and L; represents the weight of the downsampling crystal K for the upsampling crystal i; W represents the weight of the downsampling crystal L for the upsampling crystal j; tot is the normalization factor, which is represented by: Where N is the number of up-sampling response lines contained in one down-sampling response line.
29. The scatter correction method according to claim 1, characterized in that: Calculate the number of scattering coincidence events corresponding to each upsampled response line, including: Calculating the ratio of the total number of the scattering coincident events to the total number of coincident events on each of the down-sampled response lines; The number of scattering coincident events corresponding to the corresponding up-sampled response line is calculated according to the number of coincident events on each up-sampled response line included in each down-sampled response line and the ratio.
30. A scatter correction method, characterized in that: The scatter correction method comprises: Acquiring a response line based on the detection data; Based on the maximum likelihood-expectation maximization iterative algorithm, the estimated objective function corresponding to each response line is obtained by using the moment estimation method, and the estimated objective function is obtained based on the two-dimensional immediate coincidence energy histogram, the two-dimensional delayed coincidence energy histogram, the probability density function of scattered photons and the probability density function of unscattered photons corresponding to each response line; Solving the estimation objective function to obtain the number of scattering coincidence events corresponding to each response line; Calculate the number of scattering coincidence events corresponding to each response line.
31. The scatter correction method according to claim 30, characterized in that: Based on the maximum likelihood-expectation maximization iterative algorithm, the moment estimation method is used to obtain the estimated objective function corresponding to each response line, including: Based on each coincidence event corresponding to each response line, a two-dimensional immediate coincidence energy histogram, a two-dimensional delayed coincidence energy histogram and two one-dimensional detection energy spectra corresponding to each response line are obtained. Based on the two one-dimensional detection energy spectra, a maximum likelihood-expectation maximization iterative algorithm is used to obtain the probability density function of scattered photons and the probability density function of unscattered photons corresponding to the corresponding response line.
32. The scatter correction method according to claim 31, characterized in that: Based on each coincidence event corresponding to each response line, two one-dimensional detection energy spectra are obtained, including: Split all matching events corresponding to the response line into single events; The single event is divided into two parts according to its position information, and each part corresponds to a one-dimensional detection energy spectrum; Based on the divided two-part single event, two one-dimensional detection energy spectra are obtained.
33. The scatter correction method according to claim 31, characterized in that: Based on the two one-dimensional detection energy spectra, a maximum likelihood-expectation maximization iterative algorithm is used to obtain a probability density function of scattered photons and a probability density function of unscattered photons corresponding to the corresponding response lines, including: Obtaining a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; Based on the two acquired one-dimensional detection energy spectra and the corresponding relationship function, two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra are acquired using a maximum likelihood-expectation maximization iterative algorithm and the initial iteration values corresponding to the acquired corresponding response lines; The probability density functions of two scattered photons and two unscattered photons corresponding to the corresponding response lines are estimated based on the two one-dimensional gamma photon energy spectra.
34. The scatter correction method according to claim 33, characterized in that: Obtaining two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra, comprising: Based on the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, the iterative formula of the maximum likelihood-expectation maximization iterative algorithm is obtained; The number of iterations is set, and based on the obtained initial iteration value, it is iterated to the set number of iterations through the maximum likelihood-expectation maximization iterative algorithm iterative formula to obtain two one-dimensional gamma photon energy spectra corresponding to the two one-dimensional detection energy spectra.
35. The scatter correction method according to claim 33 or 34, characterized in that: Get the initial iteration value, including: Obtain the global initial scattered photon energy spectrum; Defuzzifying the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value; The global intermediate initial iteration value is stretched to obtain the initial iteration value corresponding to each response line.
36. The scatter correction method according to claim 35, characterized in that: Get the global initial scattered photon energy spectrum, including: Acquiring a coincidence response line based on the detection data; Divide all the coincidence events corresponding to all the coincidence response lines into immediate coincidence events and random coincidence events; Based on the instant coincidence event, a global one-dimensional instant energy spectrum is obtained; based on the delayed coincidence event, a global one-dimensional delayed energy spectrum is obtained; The global one-dimensional instantaneous energy spectrum corresponding to the instantaneous coincidence event on all response lines is subtracted from the global one-dimensional delayed energy spectrum corresponding to the delayed coincidence event to obtain the global one-dimensional undelayed energy spectrum; The objective function of obtaining the energy spectrum of the initial scattered photons is based on the global one-dimensional undelayed energy spectrum.
37. The scatter correction method according to claim 35, characterized in that: Defuzzify the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value, including: Obtaining the defuzzification iterative formula based on the maximum likelihood-expectation maximization iterative algorithm; The a priori initial iteration value and the number of iterations are set, and the defuzzification processing iteration formula is iterated until the number of iterations is reached to obtain the scattered part of the global initial energy spectrum; estimating an unscattered portion of the global initial energy spectrum based on the scattered portion of the global initial energy spectrum; A global intermediate initial iteration value is obtained based on a scattered part of the global initial energy spectrum and an unscattered part of the global initial energy spectrum.
38. The scatter correction method according to claim 33, characterized in that: Estimating the probability density functions of two scattered photons and two unscattered photons corresponding to the corresponding response lines based on the two one-dimensional gamma photon energy spectra includes: The one-dimensional gamma photon energy spectrum is divided into a scattered part and an unscattered part; Performing fuzzy response recovery processing on the scattered parts and unscattered parts of the two one-dimensional gamma photon energy spectra to obtain two scattered parts and two unscattered parts corresponding to the two one-dimensional detection energy spectra corresponding to the corresponding response lines; The two scattered parts and the two unscattered parts corresponding to the two one-dimensional detection energy spectra are normalized, and the probability density functions of the two scattered photons and the two unscattered photons corresponding to the corresponding response lines are estimated based on the obtained normalized values.
39. The scatter correction method according to claim 38, characterized in that: The one-dimensional gamma photon energy spectrum is divided into scattered and unscattered parts, including: The one-dimensional gamma photon energy spectrum is divided based on the energy value of the gamma photon when it is not scattered. The part less than the energy value is the scattered part, and the part equal to the energy value is the unscattered part.
40. An image reconstruction method, characterized in that: The image reconstruction method comprises: Obtaining the number of scattering coincidence events corresponding to each response line based on the scattering correction method described in any one of claims 1 to 39; The scattering events are corrected based on the number of scattering coincidence events to obtain a reconstructed image.
41. A scatter correction device, characterized in that: The scatter correction device comprises: A down-sampling response line acquisition module configured to acquire a down-sampling response line based on the detection data; An estimation module is configured to obtain an estimation target function corresponding to each downsampled response line using a moment estimation method based on a maximum likelihood-expectation maximization iterative algorithm, wherein the estimation target function is obtained based on a two-dimensional immediate coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each downsampled response line; A scattering coincidence event number acquisition module, configured to solve the estimation objective function and acquire the scattering coincidence event number corresponding to each down-sampled response line; An up-sampling module configured to perform up-sampling processing on the down-sampled response line to obtain an up-sampled response line; The scattered coincidence event number calculation module is configured to calculate the scattered coincidence event number corresponding to each up-sampled response line.
42. A scatter correction device, characterized in that: The scatter correction device comprises: a response line acquisition module configured to acquire a response line based on the detection data; An estimation module is configured to obtain an estimation target function corresponding to each response line using a moment estimation method based on a maximum likelihood-expectation maximization iterative algorithm, wherein the estimation target function is obtained based on a two-dimensional immediate coincidence energy histogram, a two-dimensional delayed coincidence energy histogram, a probability density function of scattered photons, and a probability density function of unscattered photons corresponding to each response line; The scattering coincidence event number acquisition module is configured to solve the estimation objective function and acquire the scattering coincidence event number corresponding to each response line.
43. An image reconstruction device, characterized in that: The image reconstruction device comprises: A scattering coincidence event number acquisition module, configured to acquire the number of scattering coincidence events corresponding to the response line based on the scattering correction device according to claim 41 or 42; The reconstruction module is configured to perform image reconstruction based on the number of scattering coincidence events using an iterative image reconstruction algorithm to obtain a reconstructed image.
44. A digitizing device, characterized in that: include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the computer program implements the steps of the method according to any one of claims 1 to 40 when executed by the processor.
45. A computer-readable storage medium, characterized in that The storage medium stores a computer program, and when the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 40 are implemented.
Citation Information
Patent Citations
Monte Carlo scatter removal correction method
CN112949156A
Scattering correction method and device, imaging system and computer readable storage medium
CN113506355A
Image reconstruction method, apparatus, system, and computer-readable storage medium
WO2023035361A1