A method for scatter correction of PET images

Abstract

The present invention relates to a method for scatter correction of PET images, comprising: for the detection data of the scatter region in the PET detection data, obtaining a third scatter distribution of the detection data; implementing fitting and recombination processing based on the Mask mask matrix, the second scatter distribution and the third scatter distribution of the known scatter region to obtain scaling factors corresponding to three direction variables of angle, radial direction and axial direction; constructing a cost function of the scaling factors and iteratively obtaining the scaling factors as the optimal solution; multiplying the scaling factors of the optimal solution by the second scatter distribution to obtain a fourth scatter distribution as the scatter distribution after correction of the scatter region; performing spline interpolation on the fourth scatter distribution to obtain a surface of the scatter distribution within the detection region, and downsampling the surface to obtain a continuous first scatter distribution required for PET reconstruction. The above method realizes obtaining a more accurate scatter distribution without other modality images, and has higher correction accuracy, which helps to improve the image quality.

Description

Technical Field

[0001] The present invention relates to the field of medical imaging, and particularly to an image scatter correction method in a positron emission computed tomography system. Background Art

[0002] Positron emission tomography (PET) is a high-end nuclear medicine imaging diagnostic device. In actual operation, radioactive nuclides (such as 18 F, 11 C, etc.) are used to label metabolic substances and injected into the human body. Then, the PET system is used to perform functional metabolic imaging on the patient to reflect the situation of life metabolic activities, so as to achieve the purpose of diagnosis. Currently, commercially available positron emission tomography (PET) is usually integrated with other modality imaging systems, such as computed tomography (CT) or magnetic resonance imaging (MRI), to achieve the purpose of simultaneously imaging the anatomical structure of the patient. This can accurately locate the distribution imaging of PET nuclides and improve the accuracy of lesion localization. Finally, functional imaging and anatomical imaging are fused in the same machine, combining the advantages of dual-modal imaging, enabling a clear understanding of the overall condition of the whole body at a glance, achieving the purpose of early detection of lesions and diagnosis of diseases, and having more advantages in guiding the diagnosis and treatment of tumors, heart, and brain diseases.

[0003] During the acquisition process of a positron emission tomography (PET) system, photons may undergo Compton scattering with human tissues before reaching the detector and change their flight directions. Due to the limited energy resolution of the detector, these scattered events are wrongly recorded as true coincidence events, confusing the position information of the nuclides, and thus generating scatter artifacts in the image, seriously affecting the image quality. Especially during three-dimensional data acquisition, the number of scattered coincidences may reach 30%-60% of the total count, which makes scatter correction one of the key links in PET reconstruction.

[0004] However, existing scatter correction largely depends on CT or MR images. If there are obvious movements of the patient between CT or MR scans and PET scans, the patient has a large body weight, or there is an obvious truncation, it will lead to obvious errors in the scatter distribution and significant artifacts in the image. Summary of the Invention

[0005] (I) Technical Problems to be Solved

[0006] In view of the above-mentioned disadvantages and deficiencies of the prior art, the present invention provides a method for scatter correction of PET images, which realizes obtaining a more accurate scatter distribution without other modality images, and at the same time realizes retaining the accuracy of scatter distribution correction as much as possible on the premise of reducing noise.

[0007] (II) Technical solution

[0008] In order to achieve the above object, the main technical solutions adopted by the present invention include:

[0009] In a first aspect, an embodiment of the present invention provides a method for scatter correction of PET images, which includes:

[0010] 101. For the detection data in the scatter region of the collected PET detection data, obtain a third scatter distribution S3 after preprocessing the detection data;

[0011] 102. Based on a pre-established Mask mask matrix of the scatter region, a second scatter distribution S2, and the third scatter distribution S3, perform fitting and recombination processing to obtain a scaling factor SF corresponding to three direction variables of angle, radial, and axial;

[0012] The second scatter distribution S2 is a one-dimensional scatter distribution that is pre-estimated and only includes single-scattering events;

[0013] 103. Construct a cost function of the scaling factor SF, and obtain the scaling factor as the optimal solution through an iterative method;

[0014] 104. Multiply the scaling factor of the optimal solution by the second scatter distribution S2 to obtain a fourth scatter distribution S4 as the scatter distribution after correction of the scatter region;

[0015] 105. Perform spline interpolation on the fourth scatter distribution S4 to obtain a surface of the scatter distribution within the detection region, and perform downsampling on the surface to obtain a continuous first scatter distribution S1 required for PET reconstruction.

[0016] Optionally, the 101 includes:

[0017] The PET detection data includes: true counts, random counts, and scatter counts;

[0018] The scatter region is the region obtained by removing the region where the scanned object is located from the region to which the PET detection data belongs;

[0019] Perform random correction on the detection data of the scatter region to obtain scatter counts

[0020] Formula 1:

[0021] Perform Gaussian filtering on the said scattering count to obtain a third scattering distribution S3;

[0022] Formula 2: S3 = G * y R

[0023] In Formula 2, * represents convolution; G represents the Gaussian kernel function;

[0024] R = [r1, r2, … r i …, r N T represents the average value of random noise, i represents the response line, y = [y1, y2, … y i …, y N T represents the detected data, and N represents the size of the sinogram in sinogram-based reconstruction.

[0025] Optionally, the 102 includes:

[0026] 102-1. Establish a mask matrix Mask for the scattering region using the threshold method / RANSAC algorithm / Canny operator / Roberts operator;

[0027] Formula 3:

[0028] where, △m is a preset threshold;

[0029] 102-2. Based on the mask matrix Mask, estimate the second scattering distribution S2 using the SSS method / Monte Carlo simulation method, S2 = [s21, s22 …, s2 N T ; and

[0030] 102-3. Recombine the mask matrix Mask and the second scattering distribution S2 in three directions: angular, radial, and axial to obtain Mask' and S2';

[0031]

[0032]

[0033] 102-4. Fit the second scattering distribution S2' and the third scattering distribution S3 in the angular, radial, and axial directions to obtain a fitted scaling factor SF;

[0034]

[0035] where, TN, RN, and PN represent the number of angular samples, radial samples, and axial samples in the recombined data.

[0036] ​​​Optionally, the 103 includes:

[0037] Based on the recombined Mask', S2' and S3, construct the cost function formula four of the scaling factor, and obtain the fitting factor f in the scaling factor when the L2-norm in the cost function is minimized through an iterative solution method (tn,rn,pn) ;

[0038] Formula four:

[0039] And Δ|f (tn,rn,pn) | ≤ Δd;

[0040] Where f (tn,rn,pn) is the value of the fitting factor in the scaling factor SF at the position (tn, rn, pn), mask (tn,rn,pn) is the value of the scattering region mask matrix at the position (tn, rn, pn), s2 (tn,rn,pn) is the value of the second scattering distribution S2 at the position (tn, rn, pn), s3 (tn,rn,pn) is the value of the third scattering distribution S3 at the position (tn, rn, pn), and Δd = 5.

[0041] Optionally, the 104 includes:

[0042]

[0043] Formula five: s4 (tn,rn,pn) = s2 (tn,rn,pn) × f (tn,rn,pn) , tn = 1…TN; rn = 1…RN; pn = 1…PN.

[0044] Optionally, the 105 includes:

[0045] Perform B-spline surface interpolation on the fourth scattering distribution in the angular, radial, and axial directions;

[0046] Specifically, upsample the variables in the three directions of angle, radial, and axial. Then, the sampling sequence of the fourth scattering distribution in the angular direction is defined as [1,…,tn,…,TN], and the corresponding sampling points are [u1,…,u tn ,…,u TN ; the radial sampling sequence is defined as [1,…,rn,…,RN], and the corresponding sampling points are [v1,…,v rn ,…,v RN ; the axial sampling sequence is defined as [1,…,pn,…,PN], and the corresponding sampling points are [w1,…,w pn ,…,w PN ;

[0047] Fit each sampling curve to obtain the fitting function corresponding to each sampling curve;

[0048] For the fitting surface function of the fourth scattering distribution within the given sampling range of the spline, it is Equation Six:

[0049] Equation Six:

[0050] u tn ≤u≤u tn+1 ,v rn ≤v≤v rn+1 ,w pn ≤w≤w pn+1 ;

[0051] where, B tn,i (u) represents the i-th order B-spline basis function in the angular direction within the range [u tn ,u tn+1 , B rn,,j (v) represents the j-th order B-spline basis function in the radial direction within the range [v rn ,v rn+1 , B pn,k (w) represents the k-th order B-spline basis function in the axial direction within the range [w pn ,w pn+1 ;

[0052] P i,j,k is the spline control point of the basis function within the ranges [u tn ,u tn+1 ,[v rn ,v rn+1 ,[w pn ,w pn+1 ;

[0053] I, J, K are the maximum orders of the pre-determined spline basis functions; and, the cardinality B tn,i( u) of the n-th i-th order B-spline is expressed using the De Boor-Cox recurrence formula as Equation Seven:

[0054] Equation Seven:

[0055] The 0-th order basis function is expressed as Equation Eight:

[0056] Equation Eight:

[0057] Optionally, the 105 includes:

[0058] The steps for solving the control point P of the surface function equation C(u, v, w) include:

[0059] First, fix u tn, for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the angular direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the TN parametric curves respectively;

[0060] Second, fix v rn , for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the radial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the RN parametric curves respectively;

[0061] Third, fix w pn , for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the axial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the PN parametric curves respectively;

[0062] Through the first to the third, obtain the basis B tn,i (u), B rn,j (u), B pn,k (u), and the control points P. Finally, downsample the surface C(u, v, w) to obtain the first scattering distribution S1, the angular direction sampling number TN, the radial sampling number RN, and the axial sampling number PN.

[0063] In the second aspect, an embodiment of the present invention provides a PET image reconstruction method, including:

[0064] Based on the PET detection data collected by the PET system, use the method described in any of the above first aspects to obtain the first scattering distribution,

[0065] Input the detection data and the first scattering distribution into a pre-established reconstruction function to obtain the reconstructed PET image.

[0066] In the third aspect, a PET system includes: a memory and a processor; computer program instructions are stored in the memory, and the processor executes the computer program instructions stored in the memory to specifically execute the PET image reconstruction method described in the above second aspect.

[0067] (III) Beneficial effects

[0068] In the present invention, the existing method is used to estimate the second scattering distribution, and combined with the measured data after random correction, three-dimensional surface fitting is performed in the scattering region, so that the scattering correction applied to the reconstructed image is more accurate and effectively improves the image quality.

[0069] In specific processing, three-dimensional surface fitting and interpolation are performed on the scattering distribution estimation obtained by the single-photon scattering simulation algorithm, which more accurately considers the scattering distribution differences of the object at different angles, different radial positions, and different axial positions. Compared with the traditional scattering correction algorithm, it has higher correction accuracy and helps to improve the image quality.

[0070] In addition, in order to avoid the influence of the second-modal image on the scattering distribution, when the second-modal image does not match the PET image, the PET acquisition data is directly used to estimate the scattering spatial distribution, which improves the robustness of the system.

[0071] The method of the present invention does not rely on other-modal images to obtain the position distribution of the detected object in the projection space like the traditional method, so it has stronger applicability.

[0072] In the present invention, the object distribution information is obtained by applying the PET detection data itself, and a more accurate scattering distribution can be obtained without other-modal images, obtaining higher-quality images while significantly reducing the radiation dose to the patient. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a schematic flow chart of the scattering correction method for PET images provided by an embodiment of the present invention;

[0074] Figure 2 It is a schematic diagram of the PET reconstruction image obtained by calculating the scattering distribution using the traditional algorithm;

[0075] Figure 3 is based on Figure 1 The method shown is a schematic diagram of the PET reconstruction image obtained by calculating the scattering distribution. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0076] In order to better explain the present invention for easy understanding, the present invention will be described in detail below with reference to the drawings through specific embodiments.

[0077] There is currently provided a single-scatter simulation correction (SSS) method, which is widely used for scatter correction in PET reconstruction. The SSS method simulates the scatter distribution by calculating the probability that a coincidence gamma photon experiences a single scattering event before being detected. Since there are multiple scattering components in the true detection data, a single SSS cannot accurately determine the proportional relationship between the scattering component and the total coincidence detected. The relative amount of the scattering contribution is usually determined by a sinogram radial tail fitting method or a method based on Monte Carlo simulation. When using the tail-fitting-based method for scatter estimation, in order to reduce noise, for each axial position, the sinogram of the scatter distribution simulated by SSS and the sinogram of the detection data corrected by random are usually summed along the angular and radial directions in the scatter region, and then linear correction is performed to ensure that the scatter component values of the two are equal. The purpose of this method is to sacrifice the scatter distribution accuracy to reduce the data noise in the calibration process, losing the information of the angular and radial direction distributions during the summation process, and unable to perform accurate correction when there are large differences in the angular and radial direction distributions of the simulated scatter. In addition, the SSS scatter correction depends to a large extent on CT or MR images, and obvious patient movement, large patient body weight, or obvious truncation between CT or MR scans and PET scans will lead to obvious errors in the scatter distribution, resulting in significant artifacts in the image and affecting the doctor's diagnosis.

[0078] That is to say, during the PET image reconstruction process, the SSS method is usually used to estimate the scatter distribution. However, this distribution cannot quantitatively reflect the actual scatter count. Therefore, the tail-fitting method is usually required to estimate the scatter correction scaling factor to ensure that the estimated scatter distribution after correction conforms to the actual scatter. The calculation of this scaling factor is usually estimated using the total count in the scatter region of each axial layer. Only one factor can be obtained for each layer, ignoring the differences in different angular and radial positions, resulting in inaccurate estimation of the scatter distribution and affecting the image quality.

[0079] In order to solve the problem of inaccurate scatter estimation, considering the differences in the actual scatter distribution, this application performs a three-dimensional surface fitting on the scatter distribution estimate obtained by SSS and the measurement result after random correction to achieve the purpose of noise reduction and accurate fitting.

[0080] Furthermore, in order to minimize the influence of CT or MR, when there is a significant mismatch between the distributions of CT or MR and PET, the PET data itself is used to estimate the scatter distribution, avoiding the scatter estimation error caused by CT or MR, and thus obtaining a stable high-quality image.

[0081] Embodiment 1

[0082] As Figure 1As shown in the figure, this embodiment provides a method for scatter correction of PET images. The method of this embodiment can be implemented on any electronic device, preferably in a computing device associated with a PET detector. The method of this embodiment may include the following steps:

[0083] 101. For the detection data in the scatter region of the acquired PET detection data, obtain the third scatter distribution S3 after preprocessing the detection data.

[0084] In this embodiment, the PET detection data includes: true counts, random counts, and scatter counts; the detection data in the scatter region includes: random counts and scatter counts.

[0085] The scatter region is the region obtained by removing the region where the scanned object is located from the region to which the PET detection data belongs. The scatter region in the embodiments of the present invention is the region of interest for scatter. In some of the following examples, the scatter region is used, and in some examples, the region of interest for scatter is used, both of which represent the same region without the scanned object.

[0086] The preprocessing of the detection data may include:

[0087] Perform random correction on the detection data in the scatter region to obtain scatter counts y R ;

[0088]

[0089] Perform Gaussian filtering on the scatter counts y R to obtain the third scatter distribution S3;

[0090] S3 = G * y R where * represents convolution, G represents the Gaussian kernel function, R = [r1, r2,..., r N T represents the average value of random noise, i represents the response line, y = [y1, y2,..., y N T represents the detected data, and N represents the size of the sinogram in the sinogram-based reconstruction.

[0091] 102. Based on the pre-established Mask mask matrix of the scatter region, the second scatter distribution S2, and the third scatter distribution S3, perform fitting and recombination processing to obtain the scaling factor SF corresponding to the three direction variables of angle, radial, and axial.

[0092] The second scatter distribution S2 in this embodiment is a one-dimensional scatter distribution including only single-scatter events estimated by the SSS estimation method / Monte Carlo simulation algorithm. The SSS estimation method / Monte Carlo simulation algorithm is an existing method, which is not limited in this embodiment and will not be elaborated further.​​

[0093] In this embodiment, the second scattering distribution S2 is not limited to the above two ways of acquisition. Any method capable of acquiring a one-dimensional scattering distribution can be used. This step is only for subsequent processing using the existing second scattering distribution.

[0094] It should be noted that the scattering region is the edge information of the scanned object. In this embodiment, the image edge extraction technology in the field of machine vision is applied to obtain the edge information of the scanned object in the PET detection space. Usually, there is a relatively obvious boundary between the scanned object and the region outside the scanned object in the projection space to which the PET detection data belongs. In this embodiment, the boundary in the projection space to which the PET detection data belongs can be directly used to determine the mask matrix Mask of the scattering region of interest.

[0095] In actual scatter correction, a certain redundancy is required for delimiting the scattering region to improve the robustness of the system and avoid unnecessary true coincidences within the selected region causing scatter correction errors. Therefore, the selection of the scattering region can tolerate the low resolution and high noise of PET data relative to CT or MR, can meet the requirements of scatter correction, and can avoid the errors caused by the above-mentioned mismatch with the second modality imaging position.

[0096] 103. Construct a cost function for the scaling factor SF and obtain the scaling factor as the optimal solution through an iterative method.

[0097] 104. Multiply the scaling factor of the optimal solution by the second scattering distribution S2 to obtain the fourth scattering distribution S4 as the scatter distribution after scatter region correction.

[0098] 105. Perform spline interpolation on the fourth scattering distribution S4 to obtain the surface of the scatter distribution within the detection region, and downsample the surface to obtain the continuous first scattering distribution S1 required for PET reconstruction.

[0099] In the above method, using the PET detection data itself to obtain the object distribution information can obtain a more accurate scatter distribution without other modality images, and can obtain a higher-quality image while significantly reducing the radiation dose of the patient.

[0100] In the above embodiment, a three-dimensional surface fitting can be performed on the scatter distribution estimated by the existing SSS method to make the scatter correction more accurate and significantly improve the image quality. During the processing, the influence of the second modality image on the scatter distribution is effectively avoided, that is, without considering the problem of whether the second modality image matches the PET image, and directly using the PET detection data to estimate the scatter spatial distribution to obtain the first scatter distribution.

[0101] Embodiment 2

[0102] To better illustrate the method of the embodiments of the present invention, the method of the embodiments of the present invention will be described in detail below.

[0103] Before describing the method of the embodiments of the present invention, the PET reconstruction process will be described first, because reconstruction requires scatter correction.

[0104] Currently, the PET acquisition process can be modeled by the following formula (1):

[0105]

[0106] In formula (1), y = [y1, y2,..., y N T represents the detected data, N represents the dimension of the detected data. For listmode reconstruction, N is the number of detected events; for sinogram-based reconstruction, N is the size of the sinogram; if the acquisition events contain time-of-flight (TOF) information, N should also include the dimension of the time-of-flight discretization. After obtaining the first scatter distribution S1 described below, subsequent listmode reconstruction or sinogram reconstruction can be performed.

[0107] x = [x1, x2,... x j ,..., x M T represents the unknown PET radioactivity concentration distribution image, and M represents the size of the PET image discrete space. A = [A ij is the system matrix, which mathematically expresses the probability that the point source j at a spatial position in the PET system is detected by the line of response (LOR) i, and reflects the physical characteristics of the PET detection system.

[0108] R = [r1, r2,..., r i ,..., r N T represents the average value of random noise, and S1 = [s11, s12,..., s1 i ,..., s1 N T represents the average value of scatter noise, which is defined as the first scatter distribution.

[0109] In traditional PET image reconstruction, the SSS method is used to obtain the first scatter distribution S1. Since the existing first scatter distribution S1 depends on CT / MR images (i.e., the second modality images), and only the differences in one-dimensional directions can be considered, ignoring the differences in different angular and radial directions, the estimation of the scatter distribution is inaccurate, affecting the image quality.

[0110] ​​​​For example, the Mask mask matrix (i.e., the scatter mask vector) defined artificially in the traditional method is: Mask = [mask1, mask2, …, mask N T ;

[0111]

[0112] The traditional definition method of the scatter mask is to project the image (CT or MR) of the second modality, and then perform threshold segmentation or edge extraction on the projection distribution of the second modality image to obtain the contour distribution of the scanned object in the detector space. The area outside the contour is defined as the scatter area / scatter region of interest, and is calibrated with the Mask mask matrix.

[0113] Therefore, if there is an obvious mismatch between the image of the second modality and PET during the PET system scanning process, it is inaccurate to obtain the scatter region of interest Mask using the second modality image at this time.

[0114] Therefore, the scatter correction method in this embodiment is the process of obtaining an accurate and continuous first scatter distribution S1. When obtaining the scatter mask, the second modality image is not required, and the differences in three-dimensional directions are considered during the obtaining process. The specific steps are as follows:

[0115] The first step: The detection data of the scatter region of interest includes: scatter counts and random counts. Therefore, random correction is performed on the detection data in the scatter region of interest, and the remaining counts after random correction are the scatter counts, expressed as

[0116]

[0117] Since PET detection data usually has relatively large noise, and the ideal scatter distribution belongs to a low-frequency signal, therefore, the detection data of the scatter region of interest also needs to be filtered for high-frequency signals, such as using Gaussian filtering to remove high-frequency signals.

[0118] That is, Gaussian filtering is performed on after the random correction formula (3) to obtain the third scatter distribution S3:

[0119] S3 = G * y R (4)

[0120] Among them, * in formula (4) represents convolution, and G represents the Gaussian kernel function.

[0121] It can be understood that the random correction and Gaussian filtering in this embodiment are both existing methods, and this embodiment does not limit them.

[0122] ​Step 2: In this embodiment, the mask matrix Mask of the scatter region of interest is obtained by using the threshold method:

[0123]

[0124] where △m is the threshold, which is usually selected according to experience or can be calculated by using an adaptive algorithm.

[0125] In this embodiment, the method for obtaining the mask matrix Mask includes, but is not limited to, the threshold method. It can also be obtained by using various methods such as the RANSAC algorithm, Canny operator, Roberts operator, and machine learning to obtain edge information. The scatter region of interest is the edge information of the scanned object. That is, in this embodiment, the image edge extraction technology in the field of machine vision is applied to obtain the edge position information of the object in the PET detection space.

[0126] In the traditional method, the above formula (5) can also be used as the corresponding mask matrix Mask.

[0127] Generally, there is an obvious boundary between the scanned object and the area outside the scanned object in the projection space to which the PET detection data belongs. Therefore, in this embodiment, the boundary in the projection space to which the PET detection data belongs can be directly used to determine the mask matrix Mask of the scatter region of interest.

[0128] Step 3: Since the second scatter distribution S2 estimated by the existing SSS estimation / Monte Carlo simulation algorithm only considers single-scatter events, the second scatter distribution S2 obtained by the SSS method and only containing single-scatter events is S2 = [s21, s22…, s2 N T .

[0129] Therefore, the second scatter distribution S2 estimated by the existing SSS / Monte Carlo simulation algorithm is inaccurate, and in this embodiment, further optimization processing can be performed based on the second scatter distribution S2.

[0130] Specifically, the second scatter distribution S2 is fitted with the third scatter distribution S3 representing the scatter counts in the detection data obtained in the first step to obtain the fitting scaling factor SF.

[0131] SF: The fitting scaling factor SF = [f1, f2,…, f N T .

[0132] For example, a single-scatter event can be an event that has undergone one reflection. There are single-scatter events and multi-scatter events in the scatter analysis region of interest, and all possible situations are considered in the calculation in this embodiment.

[0133] ​​Step 4: Recombine the above one-dimensional sine diagram and split it into three direction variables: angular sampling, radial sampling, and axial sampling.

[0134] It can be understood that the sine diagram is a data arrangement method for collecting data and is commonly used in the reconstruction of PET. The reconstruction formula in the aforementioned formula (1) is in one-dimensional form. Therefore, in the following, a three-dimensional fitting of the scattering distribution is required. Thus, the one-dimensional data, i.e., the one-dimensional sine diagram, needs to be recombined into a three-dimensional form. Here, it should be noted that because there is a recombination process, superscripts are added to the representation symbols of each number, but the arrays / elements inside remain unchanged.

[0135] Specifically, the second scattering distribution S2 can be recombined into S2':

[0136]

[0137] The scaling factor SF' can be obtained by fitting the recombined second scattering distribution S2' with the third scattering distribution S3 representing the scattered counts in the detection data obtained in the first step;

[0138] The recombined scaling factor SF' can be:

[0139]

[0140] At the same time, the masked matrix Mask' after recombination can be:

[0141]

[0142] In the above formula, TN represents the number of angular samplings in the sine diagram / recombined data, RN represents the number of radial samplings in the sine diagram / recombined data, and PN represents the number of axial samplings in the sine diagram / recombined data. The sine diagram here arranges the collected listmode data according to the angle, radial, and axial directions, and the rearranged and recombined data is the sine diagram.

[0143] N = TN·RN·PN;

[0144] That is, the number of angular samplings, the number of radial samplings, and the number of axial samplings are multiplied to obtain the size N of the sine diagram.

[0145] In the traditional SSS processing, the fitting is to sum the second scattering distribution S2 along the three directions of angle, radial, and axial. It loses the information of the scattering distribution differences in the three dimensions of angle, radial, and axial, resulting in inaccurate estimation of the scattering distribution. Therefore, in this embodiment, the differences in each dimension are taken into account, and recombination and scattering distribution fitting are performed in three directions to obtain a more accurate third scattering distribution.

[0146] Step 5: Construct the cost function as in formula (6). The cost function satisfies the minimum of L2-norm and the SQP algorithm is used to solve it to obtain the scaling factor of the optimal solution.

[0147] That is to say, the scattering distribution within the entire scattering region of interest is fitted along three directions to obtain an accurate scaling factor SF'. The optimal solution in the scaling factor SF' is solved. L2-norm is an existing method for constructing a cost function and is used to calculate the minimum of the sum of squares.

[0148]

[0149] where f (tn,rn,pn) takes the value of the element (fitting factor) in the scaling factor at the position (tn, rn, pn), mask (tn,rn,pn) takes the value of the mask matrix at the position (tn, rn, pn), s2 (tn,rn,pn) is the value of the second scattering distribution S2 at the position (tn, rn, pn), s3 (tn,rn,pn) takes the value of the third scattering distribution S3 at the position (tn, rn, pn).

[0150] In this embodiment, in order to ensure the robustness of the fitting process and avoid overfitting, the variable range of the scaling factor is restricted. Based on the fact that there is not a particularly large deviation between the third scattering distribution S3 containing multiple scatterings and the first scattering distribution S1 of single scattering. The variation range is based on experience. Without loss of generality, it can be selected to limit Δd = 5.

[0151] The scaling factor f of the scattering region of interest is obtained through formula (6) (tn,rn,pn) is a non-linear constrained optimization problem without an explicit analytical solution. Existing algorithms are used for iterative solution. For example, but not limited to, the SQP (Sequential Quadratic Programming) algorithm can be used.

[0152] Step 6: Multiply the scaling factor of the optimal solution by the second scattering distribution S2 to obtain the fourth scattering distribution S4 within the scattering region of interest.

[0153]

[0154] s4 (tn,rn,pn) = s2 (tn,rn,pn) × f (tn,rn,pn) , tn = 1…TN; rn = 1…RN; pn = 1…PN (7)

[0155] The fourth scattering distribution S4 obtained in this step belongs to the calibrated scattering distribution of the scattering region of interest. At this time, the scattering distribution of the scanned object region has not been obtained yet.

[0156] That is to say, the fourth scattering distribution S4 only corrects the scattering distribution in the region of interest of scattering. At this time, the scattering distribution inside the detected object cannot be accurately described. To ensure that the scattering distribution in the entire detection space is smooth and continuous, spline interpolation needs to be performed on the fourth scattering distribution S4 to obtain the accurate and continuous first scattering distribution S1.

[0157] In addition, due to the large noise in PET detection data, if the fourth scattering distribution S4 is directly used for fitting to obtain the first scattering distribution S1, there will be large oscillations in the fitting results. Therefore, it is necessary to use the interpolation method to establish a smooth and continuous surface, and then sample to obtain the first scattering distribution S1.

[0158] The spline interpolation in this embodiment can adopt the method of piecewise polynomials, which solves the continuity between the curved surfaces and makes its shape controllable. The B-spline uses the Cox recurrence formula for stable, convenient and reliable calculation. Therefore, B-spline surface interpolation is performed on the fourth scattering distribution S4 in three directions.

[0159] It should be noted that B-spline fitting is a common existing spline fitting method. The B-spline (B-spline) is a special representation form of spline curves in numerical analysis and is the abbreviation of the basis spline. Because the B-spline is stable and reliable and ensures the continuity of spline interpolation, meeting the condition that the scattering distribution has continuity, the B-spline interpolation method is used in this embodiment for fitting the scattering in the region of interest of scattering and the region of the scanned object.

[0160] The B-spline surface is formed by constructing multiple spline curves in different directions multiple times. To improve the accuracy of spline interpolation, the variables in the three directions of angle, radial and axial are first upsampled, which are respectively represented as U = [u1, u2, …, u TN , V = [v1, v2, …, v RN , W = [w1, w2, …, w PN . Then the sampling sequence in the angle direction of S4 is defined as [1, …, tn, …, TN], and the corresponding sampling points are [u1, …, u tn , …, u TN ; the radial sampling sequence is defined as [1, …, rn, …, RN], and the corresponding sampling points are [v1, …, v rn , …, v RN ; the axial sampling sequence is defined as [1, …, pn, …, PN], and the corresponding sampling points are [w1, …, w pn , …, w PN . Then, the fitting function corresponding to each segment of spline curve is obtained by fitting respectively.

[0161] It can be understood that B-spline interpolation is a piecewise interpolation algorithm, sampling multiple points between two known points, and these points are represented by u, v, w. utn ≤ u ≤ u tn+1 where u tn , u tn+1 respectively represent the two endpoints of the data segment that needs to be interpolated, and the surface function of this segment is expressed as C tn,rn,pn (u, v, w). Without loss of generality, for the spline sampling range u tn ≤ u ≤ u tn+1 , v rn ≤ v ≤ v rn+1 , w pn ≤ w ≤ w pn+1 , the fitting surface function of the fourth scattering distribution S4 is expressed as C tn,rn,pn (u, v, w):

[0162]

[0163] where B tn,i (u) represents the i-th order B-spline basis function in the angular direction within the range [u tn , u tn+1 , B rn,,j (v) represents the j-th order B-spline basis function in the radial direction within the range [v rn , v rn+1 , and B pn,k (w) represents the k-th order B-spline basis function in the axial direction within the range [w pn , w pn+1 .

[0164] At this time, the cardinality B tn,i (u) of the n-th i-th order B-spline can be expressed by the De Boor-Cox recurrence formula as follows:

[0165]

[0166] The 0-th order basis function is expressed as:

[0167]

[0168] I, J, K are the maximum orders of the spline basis functions and can be freely selected according to actual needs.

[0169] P i,j,k is the spline control point of the basis function within the ranges [u tn , u tn+1 , [v rn , v rn+1 , [w pn , w pn+1 . It can be understood that the interpolation coefficient points calculated in B-spline interpolation are called spline control points, which are used to control the shape of the spline.

[0170] In the entire variable space, the three-dimensional spline control point P is expressed as:

[0171]

[0172] The spline fitting function is expressed as:

[0173]

[0174] u tn ≤u≤u tn+1 ,v rn ≤v≤v rn+1 ,w pn ≤w≤w

[0175] Steps to find the control points P of the surface equation C(u, v, w) after fitting:

[0176] 1) Fix u tn , and for S4 (tn = 1...TN; rn = 1...RN; pn = 1...PN), calculate the partial derivative vectors along the angular direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the TN parametric curves respectively.

[0177] 2) Fix v rn , and for S4 (tn = 1...TN; rn = 1...RN; pn = 1...PN), calculate the partial derivative vectors along the radial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the RN parametric curves respectively.

[0178] 4) Fix w pn , and for S4 (tn = 1...TN; rn = 1...RN; pn = 1...PN), calculate the partial derivative vectors along the axial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the PN parametric curves respectively.

[0179] Through the above steps, the basis functions B tn,i (u), B rn,j (u), B pn,k (u), and the control points P can be obtained. Finally, downsample the surface C(u, v, w) to obtain the calculated smooth and accurate first scattering distribution S1;

[0180] The number of samples in the angular direction TN, the number of samples in the radial direction RN, and the number of samples in the axial direction PN in the first scattering distribution S1.

[0181] The interpolation method in this step uses the B-spline uniform parameter method. Other B-spline algorithms, such as the centripetal parameter method and the accumulated chord length method, are also within the scope of protection. At the same time, interpolation can also use other mature algorithms such as natural splines, Bezier splines, smooth uniform splines, or non-uniform splines, etc., which are also within the scope of protection.

[0182] In this embodiment, the fitting factor f in the scaling factor of the region of interest for scattering is solved by a three-dimensional fitting method (tn,rn,pn) , accurately evaluating the distribution differences of scattering at different positions in the angular, radial, and axial directions, while reducing the influence of noise and having a wider range of practical applications.

[0183] The calculation process of the above method has no restrictions on the definitions of fitting and interpolation algorithms, and is applicable to any fitting and interpolation methods, such as polynomial fitting, least squares fitting, spline interpolation, Newton interpolation, etc.

[0184] With the improvement of PET detection technology, especially the improvement of time-of-flight resolution, the signal-to-noise ratio of the acquired data will be greatly improved. Therefore, the scattering correction method in this embodiment can retain the accuracy of scattering distribution correction as much as possible on the premise of reducing noise.

[0185] Embodiment III

[0186] The embodiment of the present invention provides a PET image reconstruction method, which can be implemented on any electronic device. The method includes:

[0187] Based on the PET detection data collected by the PET system, the first scattering distribution is obtained by using the method described in any of the above embodiments.

[0188] The detection data and the first scattering distribution are input into a pre-established reconstruction function (such as the above formula (1)) to obtain the reconstructed PET image.

[0189] In this embodiment, the SSS method is used to estimate the scattering distribution, and three-dimensional surface fitting is performed on the measured data after random correction in the scattering region, making the scattering correction more accurate and significantly improving the quality of the reconstructed image. In the above reconstruction process, dividing the mask matrix does not require the use of a second modality image, effectively removing the influence of the second modality image on the scattering distribution.

[0190] Figure 2 For the PET reconstructed image obtained by using the traditional algorithm to calculate the scattering distribution, obvious artifacts appear due to the inability to accurately calculate the differences in the angular direction, radial direction, and axial direction. Figure 3 For the PET image obtained by using the method of Embodiment 1 or Embodiment 2 to calculate the scattering distribution, under the same reconstruction method and parameters, the scattering distribution is more accurate, there are no obvious artifacts, and the image quality is better.

[0191] In addition, the embodiment of the present invention also provides a PET system, which is characterized by including a memory and a processor; computer program instructions are stored in the memory, and the processor executes the computer program instructions stored in the memory to specifically execute the PET image reconstruction method described in any of the above embodiments.

[0192] It should be noted that in the claims, any reference signs placed between parentheses shall not be construed as limiting the claim. The word "comprising" does not exclude the presence of elements or steps not listed in a claim. The word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements. The present invention can be implemented by means of hardware including several different elements and by means of a suitably programmed computer. In a claim listing several means, several of these means can be embodied by the same hardware item. The use of the terms first, second, third, etc. is for convenience only and does not denote any order. These terms can be construed as part of the name of the element.

[0193] In addition, it should be noted that in the description of this specification, the descriptions of terms such as "one embodiment", "some embodiments", "embodiment", "example", "specific example" or "some examples", etc. mean that the specific features, structures, materials or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic descriptions of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner. In addition, without conflict, those skilled in the art can combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0194] Although the preferred embodiments of the present invention have been described, those skilled in the art can make additional changes and modifications after learning the basic creative concept. Therefore, the claims should be construed to include the preferred embodiments as well as all changes and modifications falling within the scope of the present invention.

[0195] Obviously, those skilled in the art can make various modifications and variations to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, the present invention should also include these modifications and variations.

Claims

1. A method for scatter correction of PET images, characterized in that, Including:

101. For the detection data of the scatter region in the collected PET detection data, obtain the third scatter distribution S3 after preprocessing the detection data; 102. Based on the pre-established Mask mask matrix of the scatter region, the second scatter distribution S2, and the third scatter distribution S3, perform fitting and recombination processing to obtain the scaling factor SF corresponding to the three direction variables of angle, radial, and axial; specifically, perform recombination on the mask matrix Mask and the second scatter distribution S2 in the three directions of angle, radial, and axial to obtain Mask' and S2'; perform fitting processing on the second scatter distribution S2' and the third scatter distribution S3 in the angle, radial, and axial directions to obtain the fitted scaling factor SF; The second scatter distribution S2 is a one-dimensional scatter distribution that is pre-estimated and only includes single scatter events; 103. Construct the cost function of the scaling factor SF and obtain the scaling factor as the optimal solution through an iterative method; Specifically, based on the recombined Mask', S2' and S3, the cost function formula four of the scaling factor is constructed, and the fitting factor f in the scaling factor when the L2-norm in the cost function is minimized is obtained through an iterative solution method (tn,rn,pn) ; Formula Four: and Δ|f (tn,rn,pn) | ≤ Δd; where f (tn,rn,pn) is the value of the fitting factor in the scaling factor SF at the position (tn, rn, pn), mask (tn,rn,pn) is the value of the scattering region mask matrix at the position (tn, rn, pn), s2 (tn,rn,pn) is the value of the second scattering distribution S2 at the position (tn, rn, pn), s3 (tn,rn,pn) is the value of the third scattering distribution S3 at the position (tn, rn, pn), △d = 5; 104. Multiply the scaling factor of the optimal solution by the second scatter distribution S2 to obtain the fourth scatter distribution S4 as the scatter distribution after correction of the scatter region; 105. Perform spline interpolation on the fourth scatter distribution S4 to obtain the surface of the scatter distribution within the detection region, and perform downsampling on the surface to obtain the continuous first scatter distribution S1 required for PET reconstruction.

2. The scatter correction method according to claim 1, wherein The 101 includes: The PET detection data includes: true counts, random counts, and scatter counts; The scatter region is the region obtained by removing the region where the scanned object is located from the region to which the PET detection data belongs; Randomly correct the detection data of the scattering region to obtain the scattering count ; Formula 1: i = 1,...N Perform Gaussian filtering on the said scattering count to obtain a third scattering distribution S3; Formula 2: S3 = G * y R In formula two, * represents convolution; G represents the Gaussian kernel function; R = [r1, r2, … r i …, r N T represents the average value of random noise, i represents the response line, y = [y1, y2, … y i …, y N T represents the detected data, and N represents the size of the sinogram in the sinogram-based reconstruction.​​ 3. The scatter correction method according to claim 1, wherein The 102 includes: 102-1. Use the threshold method / RANSAC algorithm / Canny operator / Roberts operator to establish the mask matrix Mask of the scatter region; Formula III: Where, △m is a preset threshold; 102-2. Based on the mask matrix Mask, estimate the second scattering distribution S2 through the SSS method / Monte Carlo simulation method, where S2 = [s21, s22…, s2 N T ; and​ 102-3. Perform recombination on the mask matrix Mask and the second scatter distribution S2 in the three directions of angle, radial, and axial to obtain Mask' and S2'; 102-4. Perform fitting processing on the second scatter distribution S2' and the third scatter distribution S3 in the angle, radial, and axial directions to obtain the fitted scaling factor SF; Where, TN, RN, and PN represent the number of angle sampling, the number of radial sampling, and the number of axial sampling in the recombination data.

4. The scatter correction method according to claim 3, wherein The 104 includes: Formula Five: s4 (tn,rn,pn) = s2 (tn,rn,pn) × f (tn,rn,pn) , where tn = 1…TN; rn = 1…RN; pn = 1…PN.

5. The scatter correction method according to claim 4, wherein The 105 includes: Perform B-spline surface interpolation on the fourth scatter distribution in the angle, radial, and axial directions; Specifically, upsampling is performed on the variables in the three directions of angle, radial direction, and axial direction. The sampling sequence of the fourth scattering distribution in the angle direction is defined as [1, …, tn, …, TN], and the corresponding sampling points are [u1, …, u tn , …, u TN ; the radial sampling sequence is defined as [1, …, rn, …, RN], and the corresponding sampling points are [v1, …, v rn , …, v RN ; the axial sampling sequence is defined as [1, …, pn, …, PN], and the corresponding sampling points are [w1, …, w pn , …, w PN ; Fit each sampled curve to obtain the fitting function corresponding to each sampled curve; The fitting surface function of the fourth scatter distribution for the given sampling range of the spline is formula six: u tn u ≤ u ≤ tn+1 , v rn v ≤ v ≤ rn+1 , w pn w ≤ w ≤ pn+1 ; Among them, B tn,i (u) represents the i-th order B-spline basis function within the range of [u tn , u tn+1 , and B rn,j (v) represents the j-th order B-spline basis function within the range of [v rn , v rn+1 , and B pn,k (w) represents the k-th order B-spline basis function within the range of [w pn , w pn+1 . P i,j,k is the spline control point of the basis function within the range [[u tn , u tn+1 , [v rn , v rn+1 , [w pn , w pn+1 ; I, J, and K are the maximum orders of the pre-determined spline basis functions; and, the basis B of the i-th B-spline at the n-th knot t tn,i( u) is expressed using the De Boor-Cox recurrence formula as Formula Seven: Formula VII: The 0th order basis function is expressed as formula eight: Formula VIII:

6. The scatter correction method according to claim 5, wherein The 105 includes: The steps of solving the control point P of the surface function equation C(u, v, w) include: First, fix u tn , for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the angular direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of TN parametric curves respectively; Second, fix v rn , for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the radial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of RN parameter curves respectively; Third, fix w pn , for S4 (tn = 1... TN; rn = 1... RN; pn = 1... PN), calculate the partial derivative vector along the axial direction to determine the boundary constraint conditions, and use the chasing method to find the control vertices of the PN parametric curves respectively; By the first to the third, the cardinal B of the B-spline is obtained tn,i (u), B rn,j (u), B pn,k (u), and the control points P. Finally, the surface C(u, v, w) is downsampled to obtain the first scattering distribution S1, the angular direction sampling number TN, the radial sampling number RN, and the axial sampling number PN.

7. A PET image reconstruction method, characterized in that, Including: Based on the PET detection data collected by the PET system, use the method described in any one of claims 1 to 6 above to obtain the first scatter distribution, Input the detection data and the first scatter distribution into a pre-established reconstruction function to obtain the reconstructed PET image.

8. A PET system, characterized in that, Including: A memory and a processor; computer program instructions are stored in the memory, and the processor executes the computer program instructions stored in the memory to specifically execute the PET image reconstruction method described in claim 7 above.

Images