A method for attenuation correction of PET system based on flow model

Abstract

The present invention discloses a PET system attenuation correction method based on a flow model. The flow model adopted by the method is a completely reversible model. The forward and reverse mappings share the same parameters. The structure itself is a consistency constraint. The model uses the spatial correlation of adjacent slices and adopts a multi-slice input and single-slice output structure. The model is composed of multiple reversible blocks, each of which is composed of a boosted affine coupling layer and a reversible 1×1 convolutional layer, and uses several small u‑nets to learn the transformation parameters of the boosted affine coupling layer. The present invention avoids additional CT or MR scanning, saves scanning costs for patients, and reduces the damage of CT radiation to patients. Compared with similar methods, a higher quality attenuation-corrected PET image can be obtained.

Description

Technical Field

[0001] The invention belongs to the technical field of biomedical image analysis, and in particular relates to a PET system attenuation correction method based on a flow model. Background Art

[0002] Attenuation correction is an important step in quantitative positron emission tomography (PET) analysis, which can improve the contrast of images and provide valuable information for doctors' diagnosis. In clinical practice, additional computed tomography (CT) or magnetic resonance (MR) scans are required to perform attenuation correction on PET images, which increases the cost and time of patient examinations. In addition, CT imaging is radioactive and harmful to the human body, and is more harmful to pediatric patients. Although MR imaging is not radioactive, it does not contain attenuation information; therefore, the anatomical structure of MR images is mainly used, but this often leads to inaccurate attenuation correction results, especially in the lungs.

[0003] Recently, some deep learning-based methods have been used to solve the attenuation correction problem without CT. Some researchers first obtain synthetic CT images (sCT) from MR images or non-attenuation corrected (NAC) PET images, and then use sCT for attenuation correction. However, there are some differences between the sCT images obtained by this method and the real CT images. When the erroneous sCT images are used for attenuation correction, the errors will be further transferred to the attenuation-corrected PET images.

[0004] Other researchers have obtained attenuation correction images directly from NAC images and achieved good results. For example, Yang et al. [Sinogram-based attenuation correction in PET / CT] and Shiri et al. [Direct attenuation correction of brain PET images using only emission data via a deep convolutional encoder-decoder (Deep-DAC)] only let the neural network learn the conversion from NAC images to AC images from the data, without considering the consistency constraint, that is, to obtain NAC images from AC images in reverse. Adding consistency constraints can make the neural network converge to a more reliable solution and reduce the common overfitting problem of deep learning; Dong et al. added consistency constraints using CycleGAN in [Deep learning-based attenuation correction in the absence of structural information for whole-body positron emission tomography imaging], but the forward and reverse generators have independent parameters and are not exactly inverse functions of each other; in addition, it uses block operations, which slows down the operation speed and will appear checkerboard when splicing image blocks.

[0005] The flow model is a fully reversible model, its forward and reverse mapping share the same parameters, and the structure itself is a consistency constraint. Once the forward mapping is trained, the reverse mapping is determined accordingly, which has the advantage of solving the PET attenuation correction problem. Summary of the invention

[0006] In view of the above, the present invention provides a PET system attenuation correction method based on a flow model, which can avoid additional CT or MR scans, save scanning costs for patients, reduce the damage of CT radiation to patients, and obtain higher quality attenuation-corrected PET images.

[0007] A PET system attenuation correction method based on a flow model comprises the following steps:

[0008] (1) Collect CT data of the scanned object, then inject the tracer into the scanned object, and use the PET system to scan and collect PET sinogram data;

[0009] (2) Convert the CT image into an attenuation image at 511 KeV energy;

[0010] (3) reconstructing an unattenuation-corrected PET image based on the sinogram data, and then using the PET image to calculate an unattenuation-corrected SUV image;

[0011] (4) reconstructing an attenuation-corrected PET image based on the sinogram data and the CT attenuation image, and using the PET image to calculate an attenuation-corrected SUV image;

[0012] (5) Establishing a flow model for attenuation correction, taking the unattenuated SUV image as input and the attenuated SUV image as the true value label, and training the flow model;

[0013] (6) The trained flow model can be used to reconstruct the sinogram data to obtain an unattenuation-corrected PET image and perform attenuation correction.

[0014] Furthermore, the scanning mode of the PET system in step (1) can be either static scanning or dynamic scanning.

[0015] Furthermore, in step (2), the CT image may be converted into an attenuation image by using a stretching-based method, a segmentation-based method, or a stretching-segmentation hybrid method.

[0016] Furthermore, in step (3), the un-attenuated PET image reconstructed using the sinogram data may be reconstructed using an analytical reconstruction method similar to filtered back projection, or an iterative reconstruction method.

[0017] Furthermore, in step (3), the unattenuated corrected SUV image is calculated by the following formula:

[0018]

[0019] Among them: SUV i is the SUV value of the i-th voxel in the unattenuated SUV image, C PET,i is the radioactivity concentration value of the i-th voxel in the unattenuated PET image, Weight is the weight of the scanned object, and Dose is the injected dose.

[0020] Furthermore, in the process of reconstructing the attenuation-corrected PET image based on the sinogram data and the CT attenuation image in step (4), the CT attenuation image needs to be registered with the PET image.

[0021] Furthermore, the flow model is composed of a plurality of reversible modules cascaded together, and the reversible module is composed of a reversible 1×1 convolutional layer connected with a lifted affine coupling layer.

[0022] Furthermore, the 1×1 convolutional layer is used to implement a permutation operation, that is, to disrupt the channel order of the input image. The matrix W is multiplied by the input image and the output is provided to the affine coupling layer. The expression of the matrix W is as follows:

[0023] W=PL(U+diag(v))

[0024] Where: P is a permutation matrix, L is a lower triangular matrix whose diagonal elements are all 1, U is an upper triangular matrix whose diagonal elements are all 0, and v is a set of numerically learnable vectors.

[0025] Furthermore, the affine coupling layer splits the input image into two parts M1 and M2 by channel, and then calculates N1 and N2 by the following formula, and finally splices N1 and N2 by channel and outputs them;

[0026] N1=M1+r(M2)

[0027] N2=M2⊙exp(s(M1))+t(M1)

[0028] Where: ⊙ represents the Hadamard product, r(), s() and t() are functions with independent parameters, which can be implemented using a neural network.

[0029] Furthermore, the specific process of training the flow model in step (5) is as follows:

[0030] 5.1 Initialize the model parameters, including randomly initializing the matrix W in the 1×1 convolution layer and initializing the parameters in the affine coupling layer using the Xavier normal distribution;

[0031] 5.2 Input the unattenuated corrected SUV images of three adjacent slices into the model, forward propagate to obtain three frames of attenuation corrected SUV images, and select the middle frame, i.e. the attenuation corrected SUV image corresponding to the middle slice, as the final output;

[0032] 5.3 Calculate the loss function between the SUV image output by the model and the SUV image after attenuation correction;

[0033] 5.4 Based on the loss function, the model parameters are continuously iterated and updated using the gradient descent method until the loss function converges and the training is completed.

[0034] Furthermore, the loss function is expressed as follows:

[0035]

[0036] Among them: loss is the loss function, The model outputs the attenuation-corrected SUV image corresponding to the middle slice, y1 is the true value label corresponding to the middle slice, y0 and y2 are the true value labels corresponding to the previous slice and the next slice, respectively, λ is the regularization coefficient, ‖‖1 and ‖‖2 represent the 1-norm and 2-norm, respectively.

[0037] Furthermore, the step 5.2 adopts a training mode of multiple frame input and single frame output.

[0038] The present invention proposes a reversible flow model to directly obtain AC PET SUV images from NAC PET SUV images. The model utilizes the spatial correlation of adjacent slices and adopts a multi-slice input form; the model consists of 8 reversible modules, each of which consists of an enhanced affine coupling layer and a reversible 1×1 convolutional layer. The present invention utilizes the structural reversibility of the flow model and adds consistency constraints to the proposed attenuation correction framework, which can alleviate the overfitting problem of deep learning to a certain extent and improve the accuracy of the results. In addition, the present invention does not need to scan additional CT or MR images for attenuation correction of PET images, which can not only reduce the scanning cost and time of patients, reduce the radiation damage of CT to patients, but also increase the patient throughput of the hospital and improve the diagnostic efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a schematic diagram of the implementation process of the method of the present invention.

[0040] Figure 2 It is a schematic diagram of the overall structure of the flow model of the present invention.

[0041] Figure 3 It is a schematic structural diagram of the affine coupling layer of the present invention.

[0042] Figure 4 This is a detailed structural diagram of small u-net.

[0043] Figure 5 For 13 The NAC PET SUV image of N-Ammonia obtains AC PETSUV images of 6 time frames under different methods; the 1st to 5th rows are NAC PET SUV images, labels of AC PET SUV images, AC PET SUV images obtained by AC-flow, AC PET SUV images obtained by Deep-DAE, and AC PET SUV images obtained by DCNN; the 1st to 6th columns are the 1st frame, the 4th frame, the 7th frame, the 10th frame, the 13th frame, and the 16th frame. DETAILED DESCRIPTION

[0044] In order to describe the present invention more specifically, the technical solution of the present invention is described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0045] like Figure 1 As shown, the attenuation correction method of the PET system based on the flow model of the present invention is generally divided into two stages:

[0046] Training phase

[0047] (1) Collect CT data of the scanned object.

[0048] (2) Inject the tracer into the scanned object and collect the raw data of PET, i.e., the PET sinogram. The scanning mode of the PET system can be either static scanning or dynamic scanning. According to the Poisson distribution of photons, establish an imaging model:

[0049]

[0050] in: represents the expectation of the sinogram of the attenuation factor, and is the estimated attenuation coefficient, i.e., the attenuation image; η ij is the probability that a pair of photons emitted by the jth voxel is detected by the i-th bin of the sinogram; l ij represents the intersection length of the projection line i and the jth voxel; represents the estimated radioactivity concentration distribution; s represents the detection efficiency; σ represents the scattered sinogram; ρ represents the smoothed random sinogram; and J represents the total number of voxel points.

[0051] (3) Converting the CT image into an attenuation image at 511 KeV energy; there are many methods, including but not limited to stretching-based methods, segmentation-based methods, and stretching-segmentation hybrid methods.

[0052] (4) Without using the attenuation image, directly reconstruct the non-attenuation corrected PET image (NAC PET); the reconstruction method can be either an analytical reconstruction method similar to filtered back projection or an iterative reconstruction method. According to the weight of the scanned object and the injection dose, the standard uptake value (SUV) is calculated to obtain the NAC PET SUV image. The calculation formula of SUV is as follows:

[0053]

[0054] Where: C PET It represents the radioactivity concentration of PET image, Weight represents the patient's weight, and Dose represents the injected dose.

[0055] (5) Register the attenuation map and the PET image, and use the CT-based attenuation image to reconstruct the PET image after attenuation correction (AC PET); the reconstruction method can be either an analytical reconstruction method similar to filtered backprojection or an iterative reconstruction method. Calculate the SUV value according to the body weight and injection dose of the scanned object to obtain the AC PET SUV image.

[0056] (6) Determine the detailed structure of the flow model for attenuation correction, adopting a multi-frame input and single-frame output mode, as Figure 2 shown. The input of the model is the NAC PET SUV images of three adjacent slices, and the output is the AC PET SUV image of the middle frame. This model is called AC-flow and consists of several reversible modules, and each reversible module consists of an enhanced affine coupling layer and a reversible 1×1 convolutional layer.

[0057] The structure of the affine coupling layer is as Figure 3 shown. For the affine coupling layer, given the input image where C, H, and W represent the total number of channels, the height of the image, and the width of the image respectively, then the output image can be obtained by the following two formulas:

[0058] n 1:c = m 1:c + r(m c+1:C )

[0059] n c+1:C = m c+1:C ⊙ exp(s(m 1:c )) + t(m 1:c )

[0060] where: c < C, ⊙ represents the Hadamard product, r(·), s(·), t(·) represent arbitrary functions with independent parameters, and they can be represented using neural networks. In this embodiment, 3 small u-nets are used to independently learn the functions r(·), s(·), t(·), Figure 4 shows the detailed structure of the small u-net. The inverse transformation of the affine coupling layer can be expressed by the following formula:

[0061] m c+1:C = (n c+1:C - t(n 1:c )) ⊙ exp(-s(n 1:c ))

[0062] m 1:c = n 1:c - r(m c+1:C )

[0063] From this formula, the input can be losslessly inferred from the output of the affine coupling layer.

[0064] The reversible 1×1 convolutional layer is used to implement the permutation operation, that is, to disrupt the channel order of the input image; given the input and the weight matrix Then output It can be expressed as:

[0065] n i,j =Wm i,j

[0066] The corresponding inverse transform is:

[0067] m i,j =W -1 n i,j

[0068] In order to reduce the amount of calculation, LU matrix decomposition is used to parameterize W:

[0069] W=PL(U+diag(v))

[0070] Where: P is a permutation matrix, L is a lower triangular matrix with all 1s on the diagonal, U is an upper triangular matrix with all 0s on the diagonal, and v is a vector.

[0071] (7) Using NAC PET SUV images as input and CT-based attenuation-corrected AC PET SUV images as labels, a flow model for attenuation correction is trained.

[0072] Initialize the flow model. For the reversible 1×1 convolutional layer, first sample a random rotation matrix W to initialize the parameters, then calculate the corresponding value of P (remains unchanged during training) and the corresponding initial values ​​of L, U, and v, and use the Xavier normal distribution to initialize other modules. During training, L is always constrained to be a lower triangular matrix with all diagonals equal to 1, and the lower triangular elements other than the diagonal are learnable; U is always constrained to be an upper triangular matrix with all diagonals equal to 0, and the upper triangular elements other than the diagonal are learnable; the sign of the elements of the v vector remains unchanged (consistent with the sign at initialization), and the absolute value of the elements of the v vector is learnable.

[0073] The NAC PET SUV images of three adjacent slices are input into the flow model, forward propagated, and three frames of images are obtained. The middle frame is selected as the final output, that is, the AC PET SUV image of the middle slice. The loss function between the output AC PET SUV image of the middle slice and the CT attenuation corrected AC PET SUV image is calculated:

[0074]

[0075] in: represents the AC PET SUV image of the middle slice of the output, y1 represents the label corresponding to the middle slice, is the core loss function, which is used to constrain the output AC PET SUV image and label of the middle slice to be as close as possible; y0 and y2 represent the labels of the previous slice and the next slice, respectively. and is a consistency loss function, which is used to constrain the intermediate slice image and the adjacent images to be as close as possible; λ is a regularization parameter, which is used to adjust the ratio between the core loss function and the consistency loss function.

[0076] Calculate the gradient of the loss function with respect to each variable, and use the Adam optimizer to update the learnable parameters in the flow model until the value of the loss function is small enough or the number of epochs reaches the upper limit.

[0077] Inference Phase

[0078] (1) Inject a tracer into the scanned object and collect the raw PET data.

[0079] (2) The NAC PET image is directly reconstructed without using the attenuation image.

[0080] (3) Calculate the SUV value based on the scanned subject’s weight and injection dose to obtain the NAC PET SUV image.

[0081] (4) Taking the NAC PET SUV image as input, directly output the AC PET SUV image.

[0082] In the following, we conduct experiments based on the cardiac gating data of patients to verify the effectiveness of this implementation.

[0083] The patient dataset used contains 36 patients, and the tracer injected is 13 N-Ammonia, the acquisition method is cardiac gated imaging, the scanning device is SIEMENS Biography 16PET / CT, and the scanning time is 20 minutes; the data of the first 3 minutes are discarded, and the remaining 17 minutes are constructed into NAC PET images and CT-AC PET images. The reconstructed image contains 16 gated time frames, each frame consists of 65 slices, and the size is 168×168. The standardized uptake value (SUV) image of each patient is calculated, and the NACPET SUV image and CT-AC PET SUV image are directly used as input and label without normalization operation.

[0084] The dataset was randomly divided into a training set (25 patients), a validation set (6 patients), and a test set (6 patients). The AC-flow proposed in the present invention was implemented using Pytorch 1.2.0 and trained on an Ubuntu 18.04 LTS server with TITAN RTX. Adam was used as the optimizer with an initial learning rate of 0.0005, an exponential decay factor of 0.9, and a batch size of 16. The training set was augmented online by rotation, scaling, and flipping. A total of 50 epochs were trained, and the test set was evaluated using the epoch with the best result on the validation set.

[0085] We evaluate the experimental results from two aspects: the whole image and the heart region. For the whole image, we use the SSIM and PSNR indicators; for the heart region, Carimas 2.10 is used to segment the 3D left ventricular myocardium and generate 3D masks for the patients in the test set to evaluate the performance of the proposed method; for the heart region, the evaluation indicators used are relative error (RE) and absolute error (AE).

[0086] The sizes of c1, c2, and c3 in the AC-flow of the present invention are 32, 64, and 96 respectively, and the size of the regularization parameter λ in the loss is 0.1. The results of AC-flow are compared with the most advanced Deep-DAC and DCNN. Both methods are encoding-decoding convolutional neural networks to achieve the conversion from NAC PET images to AC PET images; DCNN uses multi-slice input, while Deep-DAC uses single-slice input. For fair comparison, online data enhancement is performed on the two methods, and the batch size is set to 16. The model with the best performance on the validation set is selected to evaluate the performance on the test set.

[0087] Table 1 shows the evaluation indicators of different methods. The SSIM and PSNR of the images obtained by the AC-flow proposed in the present invention are the highest among the three methods, and the RE and AE in the heart area are the lowest. From the indicators, DCNN has an advantage over Deep-DAC, which shows that multi-slice input is better than single-slice input; although the same multi-slice input is used, AC-flow is still better than DCNN, which may benefit from the reversible structure of AC-flow; the reversibility of AC-flow constitutes its own constraint, which enables the model to converge to a more reliable solution and alleviate the problem of overfitting. Figure 5 The ACPET SUV images of the left ventricular region of a certain patient in different frames are shown. For the left ventricular myocardium, the AC images of DCNN and AC-flow are closer to the label than Deep-DAC, which is consistent with the indicators in Table 1. The experimental results show that the method of the present invention has advantages in solving PET attenuation correction images.

[0088] Table 1

[0089]

[0090] The above description of the embodiments is to facilitate the understanding and application of the present invention by those skilled in the art. It is obvious that those skilled in the art can easily make various modifications to the above embodiments and apply the general principles described herein to other embodiments without creative work. Therefore, the present invention is not limited to the above embodiments, and improvements and modifications made to the present invention by those skilled in the art based on the disclosure of the present invention should be within the scope of protection of the present invention.

Claims

1. A PET system attenuation correction method based on a flow model, comprising the following steps: (1) Collect CT data of the scanned object, then inject the tracer into the scanned object, and use the PET system to scan and collect PET sinogram data; (2) Convert the CT image into an attenuation image at 511 KeV energy; (3) reconstructing an unattenuation-corrected PET image based on the sinogram data, and then using the PET image to calculate an unattenuation-corrected SUV image; (4) reconstructing an attenuation-corrected PET image based on the sinogram data and the CT attenuation image, and using the PET image to calculate an attenuation-corrected SUV image; (5) Establishing a flow model for attenuation correction, taking the unattenuated SUV image as input and the attenuated SUV image as the true value label, and training the flow model; (6) The trained flow model can be used to reconstruct the sinogram data to obtain an unattenuation-corrected PET image and perform attenuation correction.

2. The PET system attenuation correction method according to claim 1, characterized in that: The scanning mode of the PET system in step (1) can be either static scanning or dynamic scanning.

3. The PET system attenuation correction method according to claim 1, characterized in that: In the step (2), the CT image may be converted into an attenuation image by using a stretching-based method, a segmentation-based method or a stretching-segmentation hybrid method.

4. The PET system attenuation correction method according to claim 1, characterized in that: In the step (3), the unattenuated corrected SUV image is calculated by the following formula; Among them: SUV i is the SUV value of the i-th voxel in the unattenuated SUV image, C PET,i is the radioactivity concentration value of the i-th voxel in the unattenuated PET image, Weight is the weight of the scanned object, and Dose is the injected dose.

5. The PET system attenuation correction method according to claim 1, characterized in that: In the step (3), the un-attenuation-corrected PET image is reconstructed using the sinogram data by using an analytical reconstruction method similar to filtered back projection, or an iterative reconstruction method. In the step (4), the attenuation-corrected PET image is reconstructed based on the sinogram data and the CT attenuation image, and the CT attenuation image needs to be registered with the PET image.

6. The PET system attenuation correction method according to claim 1, characterized in that: The flow model is composed of a plurality of reversible modules cascaded together, and the reversible module is composed of a reversible 1×1 convolutional layer connected with a lifted affine coupling layer.

7. The PET system attenuation correction method according to claim 6, characterized in that: The 1×1 convolutional layer is used to implement the permutation operation, that is, to disrupt the channel order of the input image. The matrix W is multiplied by the input image and the output is provided to the affine coupling layer. The expression of the matrix W is as follows: W=PL(U+diag(v)) Where: P is a permutation matrix, L is a lower triangular matrix whose diagonal elements are all 1, U is an upper triangular matrix whose diagonal elements are all 0, and v is a set of numerically learnable vectors.

8. The PET system attenuation correction method according to claim 6, characterized in that: The affine coupling layer splits the input image into two parts M1 and M2 by channel, and then calculates N1 and N2 by the following formula, and finally splices N1 and N2 by channel and outputs them; N1=M1+r(M2) N2=M2⊙exp(s(M1))+t(M1) Where: ⊙ represents the Hadamard product, r(), s() and t() are functions with independent parameters, which can be implemented using a neural network.

9. The PET system attenuation correction method according to claim 1, characterized in that: The specific process of training the flow model in step (5) is as follows: 5.1 Initialize the model parameters, including randomly initializing the matrix W in the 1×1 convolution layer and initializing the parameters in the affine coupling layer using the Xavier normal distribution; 5.2 Input the unattenuated corrected SUV images of three adjacent slices into the model, forward propagate to obtain three frames of attenuation corrected SUV images, and select the middle frame, i.e. the attenuation corrected SUV image corresponding to the middle slice, as the final output; 5.3 Calculate the loss function between the SUV image output by the model and the SUV image after attenuation correction; 5.4 Based on the loss function, the model parameters are continuously updated iteratively using the gradient descent method until the loss function converges and the training is completed.

10. The PET system attenuation correction method according to claim 9, characterized in that: The expression of the loss function is as follows: Among them: loss is the loss function, The model outputs the SUV image after attenuation correction corresponding to the middle slice, y1 is the true value label corresponding to the middle slice, y0 and y2 are the true value labels corresponding to the previous slice and the next slice, respectively, λ is the regularization coefficient, ‖ ‖1 and ‖ ‖2 represent the 1-norm and 2-norm, respectively.

Images