Scattering correction method and device, digital equipment and computer readable storage medium
By adopting the scattering correction method of the iterative algorithm of the maximum likelihood-desirable maximization iterative algorithm in PET imaging, the problems of low accuracy of scattering correction and complex calculation in the prior art are solved, efficient and accurate scattering correction are achieved, and the operation process is simplified.
Patent Information
- Application Number
- CN202311546008.3
- 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
The prior art has problems of low accuracy, complex calculation and low efficiency in scattering correction in PET imaging. In particular, the single scattering simulation method only considers single scattering, the double scattering simulation method and the Monte Carlo simulation method are complex in calculations, and it is difficult to obtain accurate activity images.
A scattering correction method based on the iterative algorithm based on the maximum likelihood-desirable maximization is adopted. By obtaining the correspondence function of the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, the number of unscattered single events corresponding to each response line is calculated, and the number of scattering conforming events is calculated based on this to achieve correction of the scattering event.
This method does not rely on activity images and attenuation images, simplifies the scattering correction process, improves the accuracy and efficiency of correction, and does not require adjustment of hyperparameter β, reducing operational complexity.
Smart Images

Figure CN120019791A_ABST
Abstract
Description
Technical Field
[0001] The present application relates to the field of data processing, and in particular, to a scatter correction method, apparatus, digital device, and computer-readable storage medium. Background Art
[0002] The working principle of Positron Emission Tomography (PET) 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 an 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 detector detects these two gamma photons at different times. If two scintillation crystals located on the Line of Response (LOR) in the detector detect two gamma photons 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 in which the time difference between two gamma photons generated by the same annihilation event reaching two scintillation crystals located on the response line 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 misidentified 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 localization 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 scatter only once in total) for each response line by selecting a scattering point in the input activity image and attenuation image, calculate the single-scattering events generated on all photon motion paths corresponding to this response line to obtain the single-scattering events on this response line. Finally, using only the tail data containing scattering events, the obtained single-scattering events are fitted by the Tail Fitting (TF) technique to obtain the stretching factor of the scattered 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, and 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 scattering 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, multi-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 multi-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 scattering events. However, for the simulation-based method to obtain the scattering 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, the prior art usually does not consider external radiation (i.e., the scattering 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] In view of this, there is an urgent need to provide a scatter correction method that is different from the simulation-based technique, does not depend on the activity image and attenuation image, and does not require an appropriate hyperparameter β.
[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] This application aims to provide a scatter correction method, device, digital device, and computer-readable storage medium to solve at least one problem existing in the prior art.
[0015] According to a first aspect of the present application, there is provided a scatter correction method, which includes: obtaining a correspondence function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; based on the correspondence function, using the maximum likelihood-expectation maximization iterative algorithm and the initial iteration value corresponding to the downsampled response line to obtain two one-dimensional gamma photon energy spectra corresponding to each downsampled response line; obtaining the number of unscattered single events corresponding to each downsampled response line based on the one-dimensional gamma photon energy spectrum; calculating the number of scattered coincidence events corresponding to each downsampled response line based on the total number of single events and the number of unscattered single events corresponding to each downsampled response line; performing upsampling processing 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, before obtaining the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, it further includes: obtaining the downsampled response line based on the detection data; and obtaining two one-dimensional detection energy spectra based on each downsampled response line.
[0017] In some embodiments, the downsampled response line is obtained based on the detection data, including: obtaining the coincidence response line based on the detection data; and performing downsampling processing on the coincidence response line to obtain the downsampled response line; wherein the detection data is obtained based on the target object.
[0018] In some embodiments, the one-dimensional detection energy spectrum is obtained based on each downsampled response line, including:
[0019] Splitting all the coincidence events corresponding to the downsampled response line into single events;
[0020] Dividing 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;
[0021] Based on the two parts of the divided single events, obtaining two one-dimensional detection energy spectra.
[0022] In some embodiments, the position information of the single event includes the position information of the downsampled detection module where the single event is detected.
[0023] In some embodiments, the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is:
[0024]
[0025] Among them, p represents the fuzzy response matrix, represents the expected value of the one-dimensional detection energy spectrum; X represents the one-dimensional gamma photon energy spectrum, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence event.
[0026] In some embodiments, the one-dimensional delayed energy spectrum corresponding to the random coincidence event is obtained by estimating the delayed coincidence event.
[0027] In some embodiments, obtaining the one-dimensional delayed energy spectrum includes: obtaining all the delayed coincidence events in the down-sampled response line based on the delayed coincidence window; splitting the delayed coincidence events into delayed single events;
[0028] Dividing 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;
[0029] Based on the two parts of the divided delayed single events, two one-dimensional delayed energy spectra are obtained;
[0030] Among them, the position information of the delayed single event includes the position information of the down-sampled detection module that detects the delayed single event.
[0031] In some embodiments, obtaining the fuzzy response matrix p includes: obtaining prior detection data based on the phantom; obtaining the prior down-sampled response line based on the prior detection data; dividing all the prior coincidence events corresponding to the prior down-sampled response line into prior prompt coincidence events and prior delayed coincidence events; obtaining the prior one-dimensional prompt energy spectrum based on the prior prompt coincidence events corresponding to the down-sampled response line; obtaining the prior one-dimensional delayed energy spectrum based on the prior delayed coincidence events corresponding to the down-sampled response line; obtaining the difference between the prior one-dimensional prompt energy spectrum and the prior one-dimensional delayed energy spectrum, normalizing the difference, and obtaining the one-dimensional intermediate energy spectrum; obtaining the fuzzy response matrix based on the one-dimensional intermediate energy spectrum.
[0032] In some embodiments, obtaining two one-dimensional gamma photon energy spectra corresponding to each down-sampled response line includes: obtaining the maximum likelihood-expectation maximization 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 based on the obtained initial iteration value through the maximum likelihood-expectation maximization iteration formula until the set number of iterations is reached, and obtaining the one-dimensional gamma photon energy spectrum corresponding to the one-dimensional detection energy spectrum.
[0033] In some embodiments, the maximum likelihood-expectation maximization iteration formula is:
[0034]
[0035] Among them, Yi In which, 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 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, p ij and p ib In, p is the blur response matrix corresponding to this downsampled response line, p ij is the value of the i-th row and j-th column of the blur response matrix corresponding to this downsampled response line, p ib is the value of the i-th row and b-th column of the blur 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.
[0036] In some embodiments, obtaining the initial iteration value includes: obtaining the global initial scattered photon energy spectrum; performing deblurring processing on the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value; performing stretching on the global intermediate initial iteration value to obtain the initial iteration values corresponding to each downsampled response line.
[0037] In some embodiments, obtaining the global initial scattered photon energy spectrum includes: obtaining the coincidence response lines based on the detection data; dividing all the coincidence events corresponding to all the coincidence response lines into prompt coincidence events and random coincidence events; obtaining the global one-dimensional prompt energy spectrum based on the prompt coincidence events; obtaining the 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 the response lines to obtain the global one-dimensional undelayed energy spectrum; obtaining the objective function of the global initial scattered photon energy spectrum based on the global one-dimensional undelayed energy spectrum.
[0038] In some embodiments, the objective function of the global initial scattered photon energy spectrum is:
[0039]
[0040] 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; 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.
[0041] In some embodiments, obtaining U i includes: obtaining prior detection data based on the prosthesis; 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 between the global prior prompt energy spectrum and the global prior delayed energy spectrum i .
[0042] In some embodiments, deblurring the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value, including: obtaining a deblurring 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 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; 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.
[0043] In some embodiments, the deblurring iteration formula is:
[0044]
[0045] where S i in S represents the global one-dimensional detection energy spectrum obtained by scanning the target object, i represents the i-th vertical bar region 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, 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 global blurred response matrix, P ib is the value of the i-th row and b-th column of the global blurred 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 obtained in the k-th iteration;
[0046]
[0047] where largeconstant represents a constant.
[0048] In some embodiments, a global fuzzy response matrix P is obtained. ij This includes: obtaining prior detection data based on a prosthesis; 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 a normalization process on this difference to obtain a global prior intermediate energy spectrum; and obtaining the global fuzzy response matrix based on the global prior intermediate energy spectrum.
[0049] In some embodiments, estimating the unscattered part of the global initial energy spectrum is performed through an estimation function, and the estimation function is:
[0050] X j init = δ × sum(PX init,sc ) × κ, j = 511 keV
[0051] where P is the global fuzzy response matrix; δ is a hyperparameter; κ is the ratio of unscattered photons to scattered photons;
[0052]
[0053] where sum() represents summation, and C and S respectively correspond to C i 、S i .
[0054] In some embodiments, the characterization function of the global intermediate initial iteration value is:
[0055]
[0056] where δ is a hyperparameter, P is the global fuzzy response matrix, and κ is the ratio of unscattered photons to scattered photons.
[0057] In some embodiments, stretching the intermediate initial iteration value to obtain the initial iteration value corresponding to each downsampled response line includes:
[0058] 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:
[0059]
[0060] where X 0 is the initial iteration value, sum() represents summation, P is the global fuzzy response matrix, X initis 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.
[0061] In some embodiments, obtaining the number of unscattered events corresponding to each downsampled response line based on the one-dimensional gamma photon energy spectrum includes: using the count value corresponding to the energy of the gamma photons that have not undergone scattering in the one-dimensional gamma photon energy spectrum as the number of unscattered single events; using the sum of the numbers of unscattered single events of the two one-dimensional gamma photon energy spectra corresponding to each downsampled response line as the number of unscattered events corresponding to the corresponding downsampled response line.
[0062] In some embodiments, calculating the number of scattered coincidence events corresponding to each downsampled response line is performed through the following formula:
[0063] SC = S tot - S unsc ,
[0064] or
[0065] or
[0066] where S tot is the total number of single events, and twice the total number of coincidence events corresponding to all response lines is the total number of single events corresponding to the downsampled response line; S unsc is the number of unscattered single events; SC is the number of scattered coincidence events.
[0067] In some embodiments, performing upsampling processing on the downsampled response line includes: 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.
[0068] 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:
[0069]
[0070] where 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; Indicates the weight of the downsampling crystal K for the upsampling crystal i; Indicates the weight of the downsampling crystal L for the upsampling crystal j; W tot Is a normalization factor, characterized as;
[0071]
[0072] Where N is the number of upsampling response lines included in one downsampling response line.
[0073] In some embodiments, calculating the number of scattered coincidence events corresponding to each upsampling response line includes: calculating the ratio of the total number of scattered coincidence events on each downsampling response line to the total number of coincidence events; calculating the number of scattered coincidence events corresponding to the corresponding upsampling response line according to the number of coincidence events on each upsampling response line included in each downsampling response line and the ratio.
[0074] According to the second aspect of the present application, a scatter correction method is provided. The scatter correction method includes: obtaining a correspondence function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; based on the correspondence function, using the maximum likelihood-expectation maximization iterative algorithm and the initial iterative value corresponding to the response line to obtain two one-dimensional gamma photon energy spectra corresponding to each response line; obtaining the number of unscattered single events corresponding to each response line based on the one-dimensional gamma photon energy spectrum; calculating the number of scattered coincidence events corresponding to each response line based on the total number of single events and the number of unscattered single events corresponding to each response line.
[0075] In some embodiments, the one-dimensional gamma photon energy spectrum is obtained based on each response line, including: obtaining the 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; setting the number of iterations, and performing iterations through the maximum likelihood-expectation maximization iterative algorithm iterative formula based on the obtained initial iterative value until the set number of iterations is reached, to obtain the one-dimensional gamma photon energy spectrum corresponding to the one-dimensional detection energy spectrum.
[0076] In some embodiments, obtaining the initial iterative 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 iterative value; stretching the global intermediate initial iterative value to obtain the initial iterative value corresponding to each response line.
[0077] In some embodiments, obtaining the 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 the global one-dimensional prompt energy spectrum based on the prompt coincidence events; obtaining the 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 the global one-dimensional undelayed energy spectrum; and obtaining the objective function of the global initial scattered photon energy spectrum based on the global one-dimensional undelayed energy spectrum.
[0078] In some embodiments, performing deblurring processing on the global initial scattered photon energy spectrum to obtain the global intermediate initial iteration value includes: obtaining the deblurring processing iteration formula based on the maximum likelihood-expectation maximization iteration algorithm; setting the 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.
[0079] In some embodiments, obtaining the number of unscattered events corresponding to each response line based on the one-dimensional gamma photon energy spectrum includes: using the count value corresponding to the energy of the gamma photon that has not undergone scattering in the one-dimensional gamma photon energy spectrum as the number of unscattered single events; and using the sum of the numbers of unscattered single events in the two one-dimensional gamma photon energy spectra corresponding to each response line as the number of unscattered events corresponding to the corresponding response line.
[0080] In some embodiments, calculating the number of scattered coincidence events corresponding to each response line is performed through the following formula:
[0081] SC = S tot - S unsc ,
[0082] or
[0083] or
[0084] where S tot is the total number of single events, and twice the total number of coincidence events corresponding to all response lines is the total number of single events corresponding to the downsampled response lines; S unsc is the number of unscattered single events; and SC is the number of scattered coincidence events.
[0085] In some embodiments, before obtaining the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, it further includes: obtaining response lines based on detection data; and obtaining two one-dimensional detection energy spectra based on each response line.
[0086] According to a third aspect of the present application, there is provided an image reconstruction method, which includes: 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; and performing image reconstruction by using an iterative image reconstruction algorithm based on the number of scattered coincidence events to obtain a reconstructed image.
[0087] According to a fourth aspect of the present application, there is provided a scattering correction device, which includes: a correspondence obtaining module configured to obtain a correspondence function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; a gamma photon energy spectrum obtaining module configured to obtain two one-dimensional gamma photon energy spectra corresponding to each downsampled response line based on the correspondence function by using a maximum likelihood-expectation maximization iterative algorithm and an initial iteration value corresponding to the downsampled response line; an unscattered single event number obtaining module configured to obtain the number of unscattered single events corresponding to each downsampled response line based on the one-dimensional gamma photon energy spectrum; a scattered coincidence event number obtaining module configured to calculate the number of scattered coincidence events corresponding to each downsampled response line based on the total number of single events and the number of unscattered single events corresponding to each downsampled response line; an upsampling module configured to perform upsampling processing on the downsampled response line to obtain an upsampled response line; and a scattered coincidence event number calculation module configured to calculate the number of scattered coincidence events corresponding to each upsampled response line.
[0088] According to a fifth aspect of the present application, there is provided a scattering correction device, which includes: a correspondence obtaining module configured to obtain a correspondence function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; a gamma photon energy spectrum obtaining module configured to obtain two one-dimensional gamma photon energy spectra corresponding to each response line based on the correspondence function by using a maximum likelihood-expectation maximization iterative algorithm and an initial iteration value corresponding to the response line; an unscattered single event number obtaining module configured to obtain the number of unscattered single events corresponding to each response line based on the one-dimensional gamma photon energy spectrum; and a scattered coincidence event number obtaining module configured to calculate the number of scattered coincidence events corresponding to each response line based on the total number of single events and the number of unscattered single events corresponding to each response line.
[0089] According to a sixth aspect of the present application, there is provided an image reconstruction device, which includes: a scattered coincidence event number obtaining module configured to obtain the number of scattered coincidence events corresponding to each response line based on the scattering correction device described in any one of the above embodiments; and a reconstruction module configured to perform image reconstruction by using an iterative image reconstruction algorithm based on the number of scattered coincidence events to obtain a reconstructed image.
[0090] According to a seventh aspect of the present application, there is provided a digital device, including: 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 implements the steps of the scatter correction method described in any one of the above embodiments.
[0091] According to an eighth aspect of the present application, there is provided a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the scatter correction method described in any one of the above embodiments.
[0092] Based on the above embodiments of the present application, the beneficial effects of the present application include: In the embodiments of the present application, scatter correction is performed based on the maximum likelihood-expectation maximization iteration algorithm (MLEM) algorithm, 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, simple to implement and highly accurate. BRIEF DESCRIPTION OF THE DRAWINGS
[0093] 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 to the present application. In the drawings:
[0094] Figure 1 Shows an exemplary flowchart of a scatter correction method according to an exemplary embodiment of the present application;
[0095] Figure 2 Shows a schematic diagram of crystal weight calculation in 4D linear interpolation according to an exemplary embodiment of the present application;
[0096] Figure 3 Shows an exemplary flowchart of a scatter correction method according to another exemplary embodiment of the present application;
[0097] Figure 4 Shows an exemplary flowchart of an image reconstruction method according to an exemplary embodiment of the present application;
[0098] Figure 5 Shows an exemplary module diagram of a scatter correction method according to an exemplary embodiment of the present application;
[0099] Figure 6 Shows an exemplary module diagram of a scatter correction method according to another exemplary embodiment of the present application;
[0100] Figure 7 Shows an exemplary module diagram of an image reconstruction device according to an exemplary embodiment of the present application. Detailed implementation manners
[0101] In the following, only some exemplary embodiments are briefly described. As can be recognized by those skilled in the art, 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 regarded as exemplary in nature rather than restrictive.
[0102] In the description of the present application, it should be understood that the orientation or positional relationship indicated by the terms "center", "longitudinal", "transverse", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", etc. is based on the orientation or positional relationship shown in the drawings, and is 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 therefore 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 indicating 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 "a plurality of" is two or more, unless otherwise specifically defined.
[0103] In the description of the present application, it should be noted that, unless otherwise clearly specified and defined, 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.
[0104] In the present application, unless otherwise clearly specified and defined, the fact that the first feature is "above" or "below" the second feature may include the direct contact between the first and second features, or may include the situation where the first and second features are not in direct contact but in contact through other features between them. Moreover, the fact that the first feature is "above", "over" and "on" the second feature includes that the first feature is directly above and obliquely above the second feature, or merely indicates that the relative height of the first feature in a certain dimension is higher than that of the second feature. The fact that the first feature is "below", "beneath" and "under" the second feature includes that the first feature is directly below and obliquely below the second feature, or merely indicates that the relative position of the first feature in a certain dimension is less than that of the second feature.
[0105] The following provides different embodiments or examples for implementing different structures of the present application. To simplify the present application, the components and settings of specific examples are described below. Of course, they are only examples and are not intended to limit the present application. The present 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, the present 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 the present application.
[0106] Some preferred embodiments of the present application will be 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 the present application.
[0107] Figure 1 is an exemplary flowchart of a scatter correction method shown in some embodiments of the present application. Refer to Figure 1 , the energy screening method 100 may include the following steps S110 - S160.
[0108] In the embodiments of the present application, before step S110, steps S0110 - S0120 are further included:
[0109] S0110, obtaining downsampled response lines based on the detection data.
[0110] In some embodiments of the present application, before step S0110, the following may further be included:
[0111] S001, obtaining detection data based on the target object.
[0112] In some embodiments, the target object may be a living object, including but not limited to humans, animals, etc.
[0113] In some embodiments, the detection data may be sampling data obtained based on the detection by a detector and then based on the multi - voltage threshold (MVT) method. For example, 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.
[0114] In some embodiments of the present application, step S0110 may further include the following steps:
[0115] S0111, obtaining coincidence response lines based on the detection data.
[0116] In some embodiments, the detection data is screened based on a time window and an energy window to obtain coincidence events, thereby obtaining coincidence response lines. Specifically, reference can be made to the prior art and will not be elaborated here.
[0117] S0112, downsample the coincidence response lines to obtain downsampled response lines.
[0118] 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 pair of a PET system. Since the number of coincidence events on a response line may be relatively small during normal scanning, it is difficult to obtain the corresponding energy spectrum based on the coincidence events corresponding to a single response line. Therefore, it is necessary to obtain downsampled response lines through downsampling processing. For example, originally, a detection module pair of a PET system included 3×3 = 9 coincidence response lines. After downsampling these 9 response lines, a downsampled response line merged from these 9 response lines was obtained. Specifically, the downsampling processing can be to merge the response lines included in a detection module pair, or to merge the response lines included in multiple adjacent (such as two axially, four axially and radially adjacent, etc.) detection module pairs.
[0119] In some embodiments, the downsampling processing includes merging multiple coincidence response lines that are parallel to each other or whose angles and spatial positions meet preset conditions into one downsampled response line, but is not limited thereto. The downsampling processing can specifically refer to the prior art and will not be elaborated here.
[0120] S0120, obtain two one-dimensional detection energy spectra based on each downsampled response line.
[0121] In some embodiments, the energy spectrum is represented as a one-dimensional energy histogram, where the abscissa represents energy and the ordinate represents gamma photon counts. 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 of these vertical bar regions is called a bin. Specifically, it can refer to the prior art and will not be elaborated here.
[0122] In some embodiments of the present application, step S0120 includes:
[0123] S0121, split all the coincidence events corresponding to the downsampled response line into single events.
[0124] It should be understood that a coincidence event contains two single events. Split each coincidence event corresponding to a downsampled response line into two corresponding single events. The single event contains position information, energy information, and time information.
[0125] S0122, 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.
[0126] In some embodiments, the position information of a single event includes the position information of the downsampling detection module that detects the single event. For example, one downsampling response line corresponds to two detection modules A and B. A part of the single event included in this downsampling response line corresponds to detection module A, and the other part corresponds to detection module B.
[0127] S0123, based on the two parts of the single event obtained by partitioning, obtain two one-dimensional detection energy spectra.
[0128] In some specific examples, according to the detection module information corresponding to the single event, two one-dimensional detection energy spectra corresponding to this downsampling response line can be obtained. Continuing with the above example, one one-dimensional detection energy spectrum can be obtained based on the single event corresponding to detection module A, and the other one-dimensional detection energy spectrum corresponding to this downsampling response line can be obtained based on the single event corresponding to detection module B.
[0129] The steps S110 - S160 of the embodiments of the present application are specifically as follows:
[0130] Step S110, obtain the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum.
[0131] In some embodiments, the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is:
[0132]
[0133] Wherein, represents the expected value of the one-dimensional detection energy spectrum obtained by scanning the target object corresponding to a single downsampling response line; X represents the one-dimensional gamma photon energy spectrum, p represents the fuzzy response matrix corresponding to a single downsampling response line, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence events included in a single response line.
[0134] In some embodiments, the one-dimensional delayed energy spectrum corresponding to the random coincidence events is obtained by estimating the delayed coincidence events.
[0135] In some specific embodiments, correspondingly to step S0120, obtaining the one-dimensional delayed energy spectrum includes:
[0136] S1111, based on the delayed coincidence window, obtain all the delayed coincidence events in the downsampling response line; split each delayed coincidence event into two corresponding delayed single events.
[0137] In some embodiments, the delay coincidence window includes a delay coincidence time window. Specifically, when performing delay processing on a signal, the obtained coincidence events are delay coincidence events, and the corresponding time window is called the delay coincidence time window. The delay coincidence time window is specifically obtained according to prior information. For example, it can be set to 2 ns - 10 ns. Specifically, how to obtain delay coincidence events can refer to the prior art and will not be elaborated here. It should be understood that a delay coincidence event includes two delayed single events. Each delay coincidence event included in a decimated response line is split into two corresponding delayed single events, and the delayed single events include position information, energy information, and time information.
[0138] S1112, 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.
[0139] In some specific embodiments, correspondingly to step S0122, the position information of the delayed single events includes the position information of the decimation detection module that detects the delayed single events.
[0140] S1113, based on the two divided parts of the delayed single events, obtain two one-dimensional delayed energy spectra.
[0141] In some specific examples, according to the detection module information corresponding to the delayed single events, two one-dimensional delayed energy spectra corresponding to the decimated response line can be obtained.
[0142] In some specific embodiments, obtaining the blurred response matrix p includes:
[0143] S1121, obtain prior detection data based on a phantom.
[0144] In some embodiments, the phantom can be a point source or a line source, but is not limited thereto.
[0145] S1122, obtain a prior decimated response line based on the prior detection data.
[0146] In some embodiments, correspondingly to step S0110, step S1122 includes:
[0147] S11221, obtain a prior coincidence response line based on the prior detection data;
[0148] S11222, perform decimation processing on the prior coincidence response line to obtain a prior decimated response line.
[0149] In some specific examples, how to obtain the prior coincidence response line and the prior decimated response line can refer to step S0110 and will not be elaborated here.
[0150] S1123. Divide all the prior coincidence events corresponding to the prior downsampled response lines into prior prompt coincidence events and prior delayed coincidence events.
[0151] In some embodiments, the prior prompt coincidence events are obtained based on a prompt coincidence time window.
[0152] In some specific examples, the prompt coincidence time window corresponding to the prior prompt coincidence events is set to 2 ns - 10 ns. When the prompt coincidence time window exceeds 10 ns, there may be a situation where 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 an 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 2 ns - 10 ns helps to eliminate the interference of the background signal in the subsequently obtained prior one-dimensional prompt coincidence energy spectrum.
[0153] In some embodiments, corresponding to step S1111, the delayed coincidence time window of the prior delayed coincidence events is set to 2 ns - 10 ns. Specifically, the range of the delayed coincidence time window and the prompt coincidence time window can be the same or different.
[0154] S1124. Obtain a prior one-dimensional prompt energy spectrum based on the prior prompt coincidence events corresponding to the downsampled response lines; obtain a prior one-dimensional delayed energy spectrum based on the prior delayed coincidence events corresponding to the downsampled response lines.
[0155] Regarding how to obtain the prior one-dimensional prompt energy spectrum and the prior one-dimensional delayed energy spectrum based on the downsampled response lines, reference can be specifically made to steps S0120 and S1111 - S1113, which will not be elaborated here.
[0156] S1125. Obtain the difference between the prior one-dimensional prompt energy spectrum and the prior one-dimensional delayed energy spectrum, perform normalization processing on this difference, and obtain a one-dimensional intermediate energy spectrum Y psf 。
[0157] In some embodiments, for the normalization processing of the difference between the prior one-dimensional prompt energy spectrum and the prior one-dimensional delayed energy spectrum, reference can be made to the prior art, which will not be elaborated here.
[0158] S1126. Based on the one-dimensional intermediate energy spectrum Y psf Obtain a blurring response matrix p.
[0159] In some embodiments, for obtaining the blurring 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.
[0160] Step S120: Based on the correspondence function, use the maximum likelihood - expectation maximization (MLEM) iterative algorithm and the initial iteration values corresponding to the downsampled response lines to obtain two one - dimensional gamma photon energy spectra corresponding to each downsampled response line.
[0161] In some embodiments, the one - dimensional gamma photon energy spectrum refers to the energy spectrum composed of the energy of gamma photons themselves.
[0162] In some embodiments, the maximum likelihood - expectation maximization (MLEM) iterative algorithm is obtained based on the standard maximum likelihood - expectation maximization algorithm in PET image reconstruction.
[0163] In some embodiments, step S120 includes:
[0164] S1210: Based on the correspondence 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.
[0165] In some embodiments, the iterative formula of the maximum likelihood - expectation maximization iterative algorithm is:
[0166]
[0167] where 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, b and j respectively represent the b - 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 blurring response matrix corresponding to this downsampled response line, p ij is the value of the i - th row and j - th column of the blurring response matrix corresponding to this downsampled response line, p ib is the value of the i - th row and b - th column of the blurring 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.
[0168] Since the iterative formula of the MLEM iterative algorithm has no constraint of a penalty term, the selection of the initial iteration value of Equation (2) is very important.
[0169] In some embodiments, obtaining the initial iteration value includes:
[0170] S1211. Obtain the global initial scattered photon energy spectrum.
[0171] 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.
[0172] In some embodiments, step S1211 includes:
[0173] S1211a. Obtain coincidence response lines based on the detection data.
[0174] In some specific embodiments, step S1211a can refer to step S0111, which will not be elaborated here.
[0175] S1211b. Divide all coincidence events corresponding to all coincidence response lines into prompt coincidence events and random coincidence events.
[0176] In some embodiments, how to specifically operate step S1211b can refer to the prior art, which will not be elaborated here.
[0177] S1211c. Obtain the global one-dimensional prompt energy spectrum based on the prompt coincidence events; obtain the global one-dimensional delayed energy spectrum based on the delayed coincidence events.
[0178] In some embodiments, the global one-dimensional prompt energy spectrum refers to the 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 the one-dimensional delayed energy spectrum obtained using all coincidence response lines corresponding to the detection data.
[0179] In some embodiments, the specific operation of S1211c can be analogously referred to step S1124, which will not be elaborated here.
[0180] S1211d. Subtract 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 the global one-dimensional non-delayed energy spectrum.
[0181] S1211e. Based on the global one-dimensional non-delayed energy spectrum, obtain the objective function of the global initial scattered photon energy spectrum.
[0182] In some embodiments, the objective function of the global initial scattered photon energy spectrum is:
[0183]
[0184] Among them, A i where A in A is the global initial scattered photon energy spectrum; i represents the i-th bar region in the energy spectrum; Ci where C in it is the global one - dimensional undelayed energy spectrum, obtained based on the target object, and U i where U in it is the global one - dimensional intermediate energy spectrum Y obtained by scanning 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 threshold and upper threshold of the set high - energy window (specifically set based on experience, usually 511 < hel < heu, and heu generally takes values such as 600, 650, or 700, etc.), and is the stretching coefficient.
[0185] In some embodiments, the acquisition of the global one - dimensional intermediate energy spectrum Y psf can refer to step S1125. Specifically, the one - dimensional intermediate energy spectrum obtained from all the response lines corresponding to the detection data acquired based on the target object is the global one - dimensional intermediate energy spectrum Y psf .
[0186] It should be understood that in the high - energy window of C i there are almost all unscattered photons, and the number of scattered photons is very small. By stretching U i with the total number of photons in the high - energy window, the unscattered photon energy spectrum in the complete energy window can be estimated. Then, subtracting the estimated unscattered photon energy spectrum (the part after the minus sign in Equation (3) above) from the global one - dimensional undelayed energy spectrum, the global initial scattered photon energy spectrum can be obtained.
[0187] 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 decimated response line is generally not enough. If only the number of photons on a single response line or a single decimated response line is used, S i will be affected by statistical noise. If directly obtaining the energy spectrum of a single response line or a single decimated response line 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 .
[0188] In some specific examples, as can be seen from the above, U i is the global one - dimensional intermediate energy spectrum Y obtained by scanning the entire PET system psf ; the acquisition of U i can refer to Y in step S1125 psfThe acquisition of which will not be elaborated here.
[0189] In some embodiments, the prompt energy spectrum obtained based on all response lines corresponding to prior detection data can be referred to as the global prior prompt energy spectrum, and the acquired delayed energy spectrum can be referred to as the global prior delayed energy spectrum.
[0190] S1212. Perform deblurring processing on the global initial scattered photon energy spectrum to obtain a global intermediate initial iteration value.
[0191] 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.
[0192] In some embodiments, step S1212 includes:
[0193] S1212a. Obtain a deblurring processing iteration formula based on the maximum likelihood - expectation maximization iteration algorithm.
[0194] In some specific embodiments, the deblurring processing iteration formula is:
[0195]
[0196] Wherein, S i in S 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, p ij and p ib in p is the blurred response matrix corresponding to this down - sampled response line, P ij is the value of the i - th row and j - th column of the global blurred response matrix, P ib is the value of the i - th row and b - th column of the global blurred 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;
[0197]
[0198] Wherein, largeconstant represents a very large constant, such as including but not limited to 1000000, etc.
[0199] In the embodiments of the present application, by introducing M j , the S in this formula can be made iOnly includes scattered photons, that is, after the blurring process, only the scattered part has a value, and the part corresponding to j < 511 keV is the scattered part.
[0200] In some embodiments, the global blurring response matrix refers to the blurring response matrix corresponding to all response lines based on the detection data.
[0201] In some specific examples, obtaining the global blurring response matrix P ij includes:
[0202] Obtaining prior detection data based on a phantom (such as a point source or a line source, etc.);
[0203] Obtaining all prior coincidence response lines based on the prior detection data;
[0204] Dividing all prior coincidence events corresponding to all prior coincidence response lines into global prior prompt coincidence events and global prior delayed coincidence events;
[0205] 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;
[0206] Obtaining the difference between the global prior prompt energy spectrum and the global prior delayed energy spectrum;
[0207] Normalizing the difference to obtain the global prior intermediate energy spectrum Y psf1 ;
[0208] Based on the global prior intermediate energy spectrum Y psf1 Obtaining the global blurring response matrix P ij 。
[0209] 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 blurring response matrix P psf1 based on the global prior intermediate energy spectrum Y ij , specifically, it can refer to the prior art and will not be elaborated here.
[0210] In some embodiments, the acquisition of the global prior intermediate energy spectrum Y psf1 can refer to step S1325. 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 。
[0211] S1212b, setting the prior initial iteration value and the number of iterations, and performing iterations through the de-blurring iteration formula until the number of iterations is reached to obtain the scattered part of the global initial energy spectrum.
[0212] In some specific examples, X can be set 0 to 1, and substituted into Equation (4) for iteration for several times, such as including but not limited to 50 - 500 times, to obtain the scattered part X init,sc of the global initial energy spectrum.
[0213] S1212c, estimate the unscattered part of the global initial energy spectrum based on the scattered part of the global initial energy spectrum.
[0214] In some embodiments, in the case of obtaining X init,sc , in order to estimate the unscattered part of the global initial energy spectrum, that is, 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 as follows:
[0215] X j init = δ × Sum(PX init,sc ) × κ, j = 511 keV (6)
[0216] where X j init 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 blurring response matrix, and P ij is the value of the i-th row and j-th column of the global blurring response matrix; δ is a hyperparameter.
[0217] In some specific examples, the ratio κ of unscattered photons to scattered photons can be specifically estimated by the following formula:
[0218]
[0219] where sum() represents summation, and C and S respectively correspond to C i and S i as described above.
[0220] S1212d, obtain 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.
[0221] In some specific embodiments, since gamma photons with an energy of 511 keV are mainly used in the PET system, 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 :
[0222]
[0223] Among them, δ is a hyperparameter, P is the global fuzzy response matrix, and κ is the ratio of unscattered photons to scattered photons.
[0224] In the embodiments of the present application, the S obtained by formula (3) i may have errors because formula (3) assumes that i there are no scattered photons in the high-energy window of C, but in fact i there may still be scattered photons in the high-energy window of C, resulting in i U being pulled up by the stretching coefficient and S i decreasing. Then, κ obtained based on formula (7) will be overestimated and can be balanced by the hyperparameter δ. Specifically, the influence of the overestimated κ is balanced by setting δ to a number less than 1. 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, i S is the global initial scattered photon energy spectrum, that is, the energy spectrum obtained based on all response lines detected for the target object. Response lines that do not pass through the target object are considered all scattered events. For response lines that pass through the target object, there are both scattered events and true events, that is, unscattered events. Therefore, the proportion of scattered events in response lines passing through the target object is lower than the proportion of scattered events in all response lines. A more accurate unscattered part can also be obtained by balancing with the hyperparameter δ.
[0225] S1213 stretches the global intermediate initial iteration value to obtain the initial iteration value corresponding to each downsampled response line.
[0226] In some embodiments, step S1213 includes:
[0227] Stretch the global intermediate initial iteration value based on the stretching function to obtain the initial iteration value corresponding to each downsampled response line.
[0228] In some specific examples, the stretching function is:
[0229]
[0230] Among them, 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 based on the target object scan 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.
[0231] In the embodiments of the present application, through stretching, the initial iteration value X 0 matching the number of single events included in the downsampled response line can be obtained.
[0232] S1220. Set the number of iterations. Iterate based on the obtained initial iteration value through the iteration formula of the maximum likelihood - expectation maximization iteration algorithm until the set number of iterations is reached, and obtain the one - dimensional gamma - photon energy spectrum corresponding to the one - dimensional detection energy spectrum.
[0233] In some specific embodiments, the number of iterations is specifically set as needed. For example, iterate 20 - 60 times. Usually, a better one - dimensional gamma - photon energy spectrum can be obtained within 60 iterations.
[0234] Step S130. Obtain the number of unscattered single events corresponding to each down - sampled response line based on the one - dimensional gamma - photon energy spectrum.
[0235] In some embodiments, use the count value corresponding to the energy of the gamma - photons that have not undergone scattering in the one - dimensional gamma - photon energy spectrum as the number of unscattered single events; use the sum of the number of unscattered single events of the two one - dimensional gamma - photon energy spectra corresponding to each down - sampled response line as the number of unscattered events corresponding to the corresponding down - sampled response line.
[0236] Step S140. Calculate the number of scattered coincidence events corresponding to each down - sampled response line based on the total number of single events and the number of unscattered single events corresponding to each down - sampled response line.
[0237] In some embodiments, step S160 is carried out through the following formula:
[0238] SC = S tot - S unsc
[0239] Or
[0240] Or
[0241] Where S tot is the total number of single events. Twice the total number of coincidence events corresponding to all response lines is the total number of single events corresponding to the down - sampled response line; S unsc is the number of unscattered single events; SC is the number of scattered coincidence events.
[0242] In some embodiments, step S140 includes:
[0243] Statistically count the total number of single events and the number of unscattered single events corresponding to each down - sampled response line.
[0244] In some specific examples, statistically counting the total number of single events corresponding to each down - sampled response line includes:
[0245] Statistically count all the response lines included in each down - sampled response line. Twice the total number of coincidence events corresponding to all response lines is the total number of single events corresponding to the down - sampled response line.
[0246] Step S150: Upsample the downsampled response lines to obtain upsampled response lines.
[0247] In some embodiments, upsampling the downsampled response lines includes:
[0248] Use bilinear interpolation method, 4D linear interpolation method or 5D linear interpolation method to process the downsampled response lines to obtain the upsampled response lines corresponding to each downsampled response line. The specific operation can refer to the prior art and will not be elaborated here.
[0249] Step S160: Calculate the number of scattered coincidence events corresponding to each upsampled response line.
[0250] In some embodiments, the scintillation crystals included in the detection module in the PET system are referred to as upsampled crystals, and the set of scintillation crystals corresponding to the downsampled response lines obtained by downsampling is referred to as downsampled crystals. For example, as Figure 2 shown, the downsampled crystal to which the upsampled crystal (the small squares within the dashed box) belongs is denoted as "downsampled central crystal c"; the downsampled crystal that is axially closest to the upsampled crystal is denoted as "downsampled axial adjacent crystal a"; the downsampled crystal that is radially closest to the upsampled crystal is denoted as "downsampled radial adjacent crystal t"; the downsampled crystal that is diagonally closest to the upsampled crystal is denoted as "downsampled opposite crystal ta". As Figure 2 shown, let the axial length of the downsampled crystal be D a , and the radial length be D t . An upsampled response line has two upsampled crystals, which are distinguished by superscripts "i" and "j" respectively. Taking the scintillation crystal with superscript "i" as an example, denote the radial and axial distances from the center of this upsampled crystal to the center of the "central crystal c" as Then the weights corresponding to the four downsampled crystals are as follows:
[0251] "Downsampled central crystal c":
[0252] "Downsampled axial adjacent crystal a":
[0253] "Downsampled radial adjacent crystal t":
[0254] "Downsampled opposite crystal ta": Similarly, for the four downsampled crystals corresponding to the upsampled crystal with superscript "j", there are also:
[0255]
[0256]
[0257]
[0258]
[0259] In some specific embodiments, when the 4D linear interpolation method is used to process the downsampled response lines, the following formula is used to calculate the number of scattered coincidence events corresponding to the upsampled response lines:
[0260]
[0261] 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 lines; 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 for the upsampled crystal i; represents the weight of the downsampled crystal L for the upsampled crystal j; W tot is a normalization factor, represented as;
[0262]
[0263] where N is the number of upsampled response lines included in a downsampled response line.
[0264] In some other embodiments, step S180 includes:
[0265] Calculating the ratio of the total number of the scattered coincidence events on each downsampled response line to the total number of coincidence events;
[0266] 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.
[0267] In the embodiments of the present application, the MLEM algorithm is used for scatter correction. This method does not require the hyperparameter δ, so there is no need to set the hyperparameter δ within a certain range. Compared with the existing one-step later iterative algorithm OSL-EDR (One Step Later eliminate detector response) and the optimized 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-sized phantoms require different δ. The embodiments of the present application use the MLEM algorithm for scatter correction, which is hardly affected by the hyperparameter δ and hardly needs to adjust the hyperparameter δ, making scatter correction more convenient.
[0268] 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, so that the accuracy of the scatter correction result is 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.
[0269] Figure 3 It is an exemplary flowchart of the scatter correction method 200 shown according to some embodiments of the present application.
[0270] Continue to refer to Figure 3 The scatter correction method 200 may include the following steps:
[0271] S210, obtain the correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum;
[0272] S220, based on the correspondence function, use the maximum likelihood-expectation maximization iterative algorithm and the initial iterative value corresponding to the response line to obtain two one-dimensional gamma photon energy spectra corresponding to each response line;
[0273] S230, obtain the number of unscattered single events corresponding to each response line based on the one-dimensional gamma photon energy spectrum;
[0274] S240, calculate the number of scattered coincidence events corresponding to each response line based on the total number of single events and the number of unscattered single events corresponding to each response line.
[0275] In some embodiments, before step S210, the method further includes the step of:
[0276] S0210, obtain the response line based on the detection data;
[0277] S0220, obtain two one-dimensional detection energy spectra based on each response line.
[0278] Different from the embodiment of the scatter correction method 100, after obtaining the response lines based on the detection data, the scatter correction method 100 performs downsampling processing on the response lines to obtain downsampled response lines, and then uses the maximum likelihood-expectation maximization iterative algorithm to process the downsampled response lines to obtain the number of scatter coincidence events corresponding to the downsampled response lines. After that, upsampling processing is performed on the downsampled response lines to obtain the number of scatter coincidence events of each response line corresponding to the downsampled response lines; in this embodiment, the maximum likelihood-expectation maximization iterative algorithm is directly used to process the response lines obtained based on the detection data to obtain the number of scatter coincidence events corresponding to each response line.
[0279] In the embodiments of the present application, the scatter correction method 200 can selectively combine the features of the scatter correction method 100 or other methods, and vice versa.
[0280] Figure 4 It is an exemplary flowchart of an image reconstruction method shown in some embodiments of the present application.
[0281] Continue to refer to Figure 4 , the image reconstruction method 300 may include the following steps:
[0282] S310, obtain the number of scatter coincidence events corresponding to each upsampled response line based on the scatter correction method described in the above embodiments;
[0283] S320, correct the scatter events based on the number of scatter coincidence events to obtain a reconstructed image.
[0284] In some specific embodiments, an iterative image reconstruction algorithm may be used for image reconstruction based on the number of scatter coincidence events to obtain a reconstructed image.
[0285] Step S320 may specifically refer to the prior art and will not be elaborated here.
[0286] Figure 5 It is an exemplary module diagram of a scatter correction device shown in some embodiments of the present application. As Figure 5 shown, the scatter correction device 400 may include:
[0287] A correspondence relationship acquisition module 410, configured to acquire a correspondence relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum;
[0288] The gamma photon energy spectrum acquisition module 420 is configured to obtain two one-dimensional gamma photon energy spectra corresponding to each downsampled response line based on a correspondence function, using the maximum likelihood-expectation maximization iterative algorithm and the initial iteration values corresponding to the downsampled response lines.
[0289] The unscattered single event number acquisition module 430 is configured to obtain the number of unscattered single events corresponding to each downsampled response line based on the one-dimensional gamma photon energy spectrum.
[0290] The scattered coincidence event number acquisition module 440 is configured to calculate the number of scattered coincidence events corresponding to each downsampled response line based on the total number of single events and the number of unscattered single events corresponding to each downsampled response line.
[0291] The upsampling module 450 is configured to perform upsampling processing on the downsampled response lines to obtain upsampled response lines.
[0292] The scattered coincidence event number calculation module 460 is configured to calculate the number of scattered coincidence events corresponding to each upsampled response line.
[0293] In some embodiments, the correspondence acquisition module 410 includes a delayed energy spectrum acquisition module configured to obtain a one-dimensional delayed energy spectrum; a blurred response matrix acquisition module configured to obtain prior detection data based on a phantom and obtain a global blurred response matrix and the blurred response matrices corresponding to each downsampled response line based on the prior detection data.
[0294] In some embodiments, the gamma photon energy spectrum acquisition module 420 includes a maximum likelihood-expectation maximization iterative algorithm iterative formula acquisition module configured to obtain the 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; an initial iteration value acquisition module configured to obtain the initial iteration values; and an iteration module configured to set the number of iterations and perform iterations based on the obtained initial iteration values through the maximum likelihood-expectation maximization iterative algorithm iterative formula until the set number of iterations is reached, and obtain the one-dimensional gamma photon energy spectrum corresponding to the one-dimensional detection energy spectrum.
[0295] In some specific embodiments, 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 deblurring processing 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 the initial iteration values corresponding to each downsampled response line.
[0296] In some embodiments, the scatter correction device 400 further includes a statistics module configured to count the total number of single events and the number of unscattered single events corresponding to each downsampled response line.
[0297] In an embodiment of the present application, the scatter correction device 400 can selectively incorporate the features of the scatter correction method 100 or 200 or other methods, and vice versa.
[0298] In an embodiment of the present application, the scatter correction device 400 can be used to implement the scatter correction method 100 or the methods described in other embodiments herein, and can selectively incorporate the features of the scatter correction method 200 or other methods, and vice versa.
[0299] Figure 6 is an exemplary block diagram of a scatter correction device according to some other embodiments of the present application. As Figure 6 shown, the scatter correction device 500 may include:
[0300] A correspondence acquisition module 510 configured to acquire a correspondence function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum;
[0301] A gamma photon energy spectrum acquisition module 520 configured to, based on the correspondence function, use the maximum likelihood-expectation maximization iterative algorithm and the initial iterative values corresponding to the lines of response to obtain two one-dimensional gamma photon energy spectra corresponding to each line of response;
[0302] An unscattered single event number acquisition module 530 configured to acquire the number of unscattered single events corresponding to each line of response based on the one-dimensional gamma photon energy spectrum;
[0303] A scattered coincidence event number acquisition module 540 configured to calculate the number of scattered coincidence events corresponding to each line of response based on the total number of single events and the number of unscattered single events corresponding to each line of response.
[0304] Different from the embodiment of the scatter correction device 400, after the scatter correction device 400 obtains the lines of response based on the detection data, it performs downsampling processing on the lines of response to obtain downsampled lines of response, and uses the maximum likelihood-expectation maximization iterative algorithm to process the downsampled lines of response to obtain the number of scattered coincidence events corresponding to the downsampled lines of response, and then performs upsampling processing on the downsampled lines of response to obtain the number of scattered coincidence events corresponding to each line of response; in this embodiment, the maximum likelihood-expectation maximization iterative algorithm is directly used to process the lines of response obtained based on the detection data to obtain the number of scattered coincidence events corresponding to each line of response.
[0305] In an embodiment 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 incorporate the features of the scatter correction device 400 or other devices, and can also selectively incorporate the features of the scatter correction method 100 or 200 or other methods, and vice versa.
[0306] Figure 7 is an exemplary block diagram of an image reconstruction apparatus according to some embodiments of the present application. As Figure 7 shown, the image reconstruction apparatus 600 may include:
[0307] a scattered coincidence event number acquisition module 610 configured to obtain the number of scattered coincidence events corresponding to each upsampled response line based on the scattered correction method described in any of the above embodiments;
[0308] a reconstruction module 620 configured to correct scattered events based on the number of scattered coincidence events and obtain a reconstructed image.
[0309] In the embodiments of the present application, the image reconstruction apparatus 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.
[0310] In some embodiments of the present application, the image reconstruction apparatus 600 may also include components or features of the scattered correction apparatuses 400 and 500 in a non - contradictory manner, and vice versa.
[0311] In some embodiments, the present application also provides a digital device, which includes: the apparatus described in any of the above embodiments.
[0312] 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 of the above embodiments.
[0313] 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 of the above embodiments. The computer program includes each program module / unit constituting the apparatus 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 embodiments. The computer program can also run on an electronic device as described in the embodiments of the present application.
[0314] 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 herein, 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.
[0315] Meanwhile, this application uses specific terms to describe the embodiments of this application. For example, "one embodiment", "an embodiment", and / or "some embodiments" mean a certain feature, structure, or characteristic related to at least one embodiment of this 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 this application does not necessarily refer to the same embodiment. In addition, certain features, structures, or characteristics in one or more embodiments of this application can be appropriately combined.
[0316] In addition, those skilled in the art can understand that various aspects of this application can be illustrated and described by several patentable types or situations, including any new and useful process, machine, product, or composition of matter, or any new and useful improvement thereof. Accordingly, various aspects of this 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 blocks", "modules", "engines", "units", "components", or "systems". In addition, various aspects of this application may be embodied as a computer product located in one or more computer-readable media, which includes computer-readable program code.
[0317] A computer storage medium may contain a propagated data signal containing computer program code, for example, 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.
[0318] The computer program code required for the operations of various parts of this 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. This 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).
[0319] In addition, unless clearly stated in the claims, the order of the processing elements and sequences, the use of numbers and letters, 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 considered useful embodiments of the invention are discussed through various examples in the above disclosure, it should be understood that such details only serve the purpose of illustration, and 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 an existing server or mobile device.
[0320] Similarly, it should be noted that, in order to simplify the expression of the disclosure of this application 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 object of this application are more than those mentioned in the claims. In fact, the features of the embodiments are less than all the features of the individual embodiments disclosed above.
[0321] In some embodiments, numbers are used to describe components and the quantity of attributes. It should be understood that such numbers used in the description of embodiments are, in some examples, modified by the modifiers "about", "approximately" or "substantially". Unless otherwise stated, "about", "approximately" or "substantially" indicate that the said numbers are allowed to vary by ±20%. Accordingly, in some embodiments, the numerical parameters used in the specification and claims are approximate values, which may vary according to the characteristics required by individual embodiments. In some embodiments, the numerical parameters should take into account the specified significant digits and adopt the method of retaining the general number of digits. Although the numerical ranges and parameters used in some embodiments of the present application to confirm the breadth of their scope are approximate values, in specific embodiments, such numerical settings are made as precise as possible within the feasible range.
[0322] For each patent, patent application, patent application publication, and other materials cited in the present application, such as articles, books, specifications, publications, documents, etc., their entire contents are hereby incorporated into the present application by reference. Except for the application history documents that are inconsistent with or conflict with the content of the present application, and also except for the documents that limit the broadest scope of the claims of the present application (currently or subsequently attached to the present application). It should be noted that if there are any inconsistencies or conflicts between the descriptions, definitions, and / or uses of terms in the attached materials of the present application and the content described in the present application, the descriptions, definitions, and / or uses of terms in the present application shall prevail.
[0323] Finally, it should be noted that the above are only example embodiments of the present application and are not intended 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 modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present application shall be included within the protection scope of the present application.
Claims
1. A scatter correction method, characterized in that: The scatter correction method comprises: Obtaining a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; Based on the corresponding relationship function, the maximum likelihood-expectation maximization iterative algorithm and the initial iteration value corresponding to the down-sampled response line are used to obtain two one-dimensional gamma photon energy spectra corresponding to each down-sampled response line; Based on the one-dimensional gamma photon energy spectrum, the number of unscattered single events corresponding to each down-sampled response line is obtained; Based on the total number of single events and the number of unscattered single events corresponding to each downsampled response line, the number of scattered coincident events corresponding to each downsampled response line is calculated; 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: Before obtaining the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, it also includes: Obtaining a downsampled response line based on the detection data; Two one-dimensional detection energy spectra are obtained based on each downsampled response line.
3. The scatter correction method according to claim 2, 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 the coincident response line to obtain a downsampled response line; Wherein, the detection data is acquired based on the target object.
4. The scatter correction method according to claim 2, characterized in that: The one-dimensional detection energy spectrum is obtained based on each down-sampling response line, including: Split all matching events corresponding to the downsampled response line into single events; The single event is divided into two parts according to the position information of the single event, 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.
5. The scatter correction method according to claim 4, characterized in that: The location information of the single event includes location information of the downsampling detection module that detects the single event.
6. The scatter correction method according to claim 1, characterized in that: The corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum is: Among them, p represents the fuzzy response matrix, represents the expected value of the one-dimensional detection energy spectrum; X represents the one-dimensional gamma photon energy spectrum, and r represents the one-dimensional delayed energy spectrum corresponding to the random coincidence event.
7. The scatter correction method according to claim 6, characterized in that: The one-dimensional delayed energy spectrum corresponding to the random coincidence event is obtained by delay coincidence event estimation.
8. The scatter correction method according to claim 7, 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 the delay coincidence events into 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.
9. The scatter correction method according to claim 6, 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.
10. The scatter correction method according to claim 1, characterized in that: Obtain two one-dimensional gamma photon energy spectra corresponding to each downsampled response line, including: 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 the maximum likelihood-expectation maximization iteration algorithm is iterated based on the obtained initial iteration value until the set number of iterations is reached to obtain the one-dimensional gamma photon energy spectrum corresponding to the one-dimensional detection energy spectrum.
11. The scatter correction method according to claim 10, 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 b-th and j-th vertical bar areas in the energy spectrum respectively, k is the number of iterations, p is the number of iterations, and 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.
12. The scatter correction method according to claim 10, 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.
13. The scatter correction method according to claim 12, 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.
14. The scatter correction method according to claim 13, characterized in that: The objective function of the global initial scattered photon energy spectrum is: Among them, S i S in the equation 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 scanning the entire PET system; 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.
15. The scatter correction method according to claim 14, 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 .
16. The scatter correction method according to claim 12, 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.
17. The scatter correction method according to claim 16, 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 b-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 and p ib The p in the equation is the fuzzy response matrix corresponding to the down-sampled response line, 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; The global one-dimensional gamma photon energy spectrum corresponding to the j-th vertical bar area 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 obtained at the k-th iteration; Among them, largeconstant represents the constant.
18. The scatter correction method according to claim 17, characterized in that: Get the global fuzzy response matrix P ij 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; 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.
19. The scatter correction method according to claim 16, characterized in that: The estimation of the unscattered part of the global initial energy spectrum is performed by an estimation function, which is: j = 511keV, 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 .
20. The scatter correction method according to claim 16, 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.
21. The scatter correction method according to claim 12, 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.
22. The scatter correction method according to claim 1, characterized in that: Based on the one-dimensional gamma photon energy spectrum, the number of unscattered events corresponding to each downsampled response line is obtained, including: The count value corresponding to the energy of unscattered gamma photons in the one-dimensional gamma photon energy spectrum is taken as the number of unscattered single events; the sum of the numbers of unscattered single events of the two one-dimensional gamma photon energy spectra corresponding to each down-sampling response line is taken as the number of unscattered events corresponding to the corresponding down-sampling response line.
23. The scatter correction method according to claim 1, characterized in that: Calculate the number of scattering coincidence events corresponding to each downsampled response line using the following formula: SC=S tot -S unsc , or or Among them, S tot is the total number of single events, and twice the total number of events corresponding to all response lines is the total number of single events corresponding to the downsampled response line; S unsc is the number of unscattered single events; SC is the number of scattered coincident events.
24. 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.
25. The scatter correction method according to claim 24, 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, characterized by; Where N is the number of up-sampling response lines contained in one down-sampling response line.
26. 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.
27. A scatter correction method, characterized in that: The scatter correction method comprises: Obtaining a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; Based on the corresponding relationship function, the maximum likelihood-expectation maximization iterative algorithm and the initial iteration value corresponding to the response line are used to obtain two one-dimensional gamma photon energy spectra corresponding to each response line; Based on the one-dimensional gamma photon energy spectrum, the number of unscattered single events corresponding to each response line is obtained; Based on the total number of single events and the number of unscattered single events corresponding to each response line, the number of scattered coincident events corresponding to each response line is calculated.
28. The scatter correction method according to claim 27, characterized in that: The one-dimensional gamma photon energy spectrum is obtained based on each response line, including: 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 the maximum likelihood-expectation maximization iteration algorithm is iterated based on the obtained initial iteration value until the set number of iterations is reached to obtain the one-dimensional gamma photon energy spectrum corresponding to the one-dimensional detection energy spectrum.
29. The scatter correction method according to claim 27, 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.
30. The scatter correction method according to claim 29, 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.
31. The scatter correction method according to claim 29, 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.
32. The scatter correction method according to claim 27, characterized in that: Based on the one-dimensional gamma photon energy spectrum, the number of unscattered events corresponding to each response line is obtained, including: The count value corresponding to the energy of unscattered gamma photons in the one-dimensional gamma photon energy spectrum is taken as the number of unscattered single events; the sum of the numbers of unscattered single events of the two one-dimensional gamma photon energy spectra corresponding to each response line is taken as the number of unscattered events corresponding to the corresponding response line.
33. The scatter correction method according to claim 27, characterized in that: Calculate the number of scattering coincidence events corresponding to each response line using the following formula: SC=S tot -S unsc , or or Among them, S tot is the total number of single events, and twice the total number of events corresponding to all response lines is the total number of single events corresponding to the downsampled response line; S unsc is the number of unscattered single events; SC is the number of scattered coincident events.
34. The scatter correction method according to claim 27, characterized in that: Before obtaining the corresponding relationship function between the one-dimensional detection energy spectrum and the one-dimensional gamma photon energy spectrum, it also includes: Acquiring a response line based on the detection data; Two one-dimensional detection spectra are obtained based on each response line.
35. 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 34; Based on the number of scattering coincidence events, an iterative image reconstruction algorithm is used to reconstruct the image and obtain a reconstructed image.
36. A scatter correction device, characterized in that: The scatter correction device comprises: A corresponding relationship acquisition module configured to acquire a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; A gamma photon energy spectrum acquisition module is configured to acquire two one-dimensional gamma photon energy spectra corresponding to each down-sampled response line based on a corresponding relationship function, using a maximum likelihood-expectation maximization iterative algorithm and an initial iteration value corresponding to the down-sampled response line; An unscattered single event number acquisition module, configured to acquire the number of unscattered single events corresponding to each down-sampled response line based on a one-dimensional gamma photon energy spectrum; A scattered coincidence event number acquisition module is configured to calculate the number of scattered coincidence events corresponding to each downsampled response line based on the total number of single events corresponding to each downsampled response line and the number of unscattered single events; 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.
37. A scatter correction device, characterized in that: The scatter correction device comprises: A corresponding relationship acquisition module configured to acquire a corresponding relationship function between a one-dimensional detection energy spectrum and a one-dimensional gamma photon energy spectrum; A gamma photon energy spectrum acquisition module is configured to acquire two one-dimensional gamma photon energy spectra corresponding to each response line based on a corresponding relationship function, using a maximum likelihood-expectation maximization iterative algorithm and an initial iteration value corresponding to the response line; An unscattered single event number acquisition module, configured to acquire the number of unscattered single events corresponding to each response line based on a one-dimensional gamma photon energy spectrum; The scattered coincidence event number acquisition module is configured to calculate the scattered coincidence event number corresponding to each response line based on the total number of single events corresponding to each response line and the number of unscattered single events.
38. 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 each response line based on the scattering correction device according to claim 36 or 37; 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.
39. 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 when the computer program is executed by the processor, the steps of the scatter correction method according to any one of claims 1 to 34 are implemented.
40. 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 scatter correction method according to any one of claims 1 to 34 are implemented.
Citation Information
Patent Citations
PET random coincidence correction
CN106233336A
Scattering correction method and device, imaging system and computer readable storage medium
CN113506355A
Methods and apparatuses for reconstructing incident energy spectrum for a detector
US20180038970A1
Energy-Based Scatter Correction for PET Sinograms
US20210059629A1
System and method for medical imaging
US20230360794A1