A dynamic PET image reconstruction method based on model-driven deep learning

Abstract

This paper discloses a dynamic PET image reconstruction method based on model-driven deep learning. This method utilizes 3D spatiotemporal convolution to simultaneously extract correlations in both the temporal and spatial domains of dynamic projection data. Forward and backward projection operators are integrated into the reconstruction network, resulting in strong physical constraints and interpretability. The method splits the dynamic PET image reconstruction problem into several cascaded reconstruction blocks, each of which contains a primary network for updating primary image domain variables and a dual network for updating dual measurement domain variables. This method can reconstruct high-quality dynamic PET tracer activity distribution images from ultra-low count dynamic PET projection data, addressing the issues of current mainstream methods with limited interpretability and poor reconstruction quality.

Description

Technical Field

[0001] The present invention belongs to the field of PET imaging technology, and in particular relates to a dynamic PET image reconstruction method based on model-driven deep learning. Background Art

[0002] Dynamic PET imaging can quantify physiological parameters within biological tissues and plays an indispensable role in tumor detection, cardiac disease characterization, and drug development. However, image reconstruction from dynamic measurement data is extremely challenging due to the ill-posed nature of the PET reconstruction problem and the low counts of single frames of dynamic PET data. Especially in early time frames, the counts of single frames are very low, and the reconstructed images are severely affected by noise.

[0003] Furthermore, with the continuous advancement of imaging technology and detector capabilities, ultrafast time-frame imaging is now feasible in hardware. However, existing reconstruction algorithms are not well suited to high-quality reconstruction of ultra-low-count projection data. Traditional dynamic image reconstruction algorithms, such as filtered back projection and maximum likelihood expectation maximization, cannot model time-domain signals, resulting in significant differences in reconstructed images across different frames.

[0004] With the development of deep learning technology, a number of new deep learning-based solutions have emerged in the field of PET image reconstruction, such as direct learning and model-based learning. Reference [G. Wang and J. Qi, “PET image reconstruction using kernel method,” IEEE transactions on medical imaging, vol. 34, no. 1, pp. 61–71, 2014] improves reconstruction performance by introducing a temporal prior based on the kernel method. However, in practical applications, the performance is not very good for ultra-fast time frames and ultra-low count data. This is partly because the prior information of this method is obtained from single-shot reconstruction data and does not incorporate the advantages of data-driven methods. Reference [B. Wang and H. Liu, “FBP-Net for direct reconstruction of dynamic PET images,” Physics in Medicine & Biology, vol. 65, no. 23, p. 235008, 2020] achieves relatively good dynamic reconstruction results by combining a traditional filtered back projection algorithm with a denoising neural network. However, this method does not consider the system matrix constraints related to the physical characteristics of the PET instrument, resulting in poor interpretability and generalization.

[0005] This suggests that existing technologies either fail to fully exploit the temporal correlation of dynamic projection data or lack physical constraints, resulting in unstable reconstruction results and poor performance at ultra-low counts, thus limiting the development of ultrafast time-frame PET imaging. To achieve superior reconstruction quality, existing dynamic reconstruction methods often require projection data with longer acquisition times. However, these long scan times inevitably introduce motion artifacts, further degrading the quality of the reconstructed image. Summary of the Invention

[0006] In view of the above, the present invention provides a dynamic PET image reconstruction method based on model-driven deep learning, which can reconstruct high-quality dynamic PET tracer activity distribution maps from ultra-low count dynamic projection data. It can effectively utilize the spatiotemporal correlation of dynamic projection data, and combines the constraints of the physical projection matrix, and has strong interpretability.

[0007] A dynamic PET image reconstruction method based on model-driven deep learning includes the following steps:

[0008] (1) Using a detector to detect the biological tissue injected with radioactive drugs, and collecting the corresponding dynamic sinusoidal projection data Y;

[0009] (2) reconstructing the dynamic sinusoidal projection data Y to obtain the corresponding dynamic PET tracer activity distribution map X;

[0010] (3) performing steps (1) and (2) multiple times to obtain a large number of samples, each sample including dynamic sinusoidal projection data Y and its corresponding dynamic PET tracer activity distribution map X, and then dividing all samples into a training set, a validation set, and a test set;

[0011] (4) According to the dynamic PET measurement equation, the dynamic reconstruction problem is transformed into a Poisson log-likelihood optimization problem with a regularization term, and the properties of the dual variables are used to transform this optimization problem into the corresponding saddle point problem;

[0012] (5) The saddle point problem is solved by alternately updating the primary and dual variables using the primary-dual network, thereby constructing the STPD-Net model for dynamic PET image reconstruction. The model is composed of a cascade of several reconstruction modules, each of which is composed of a primary network and a dual network.

[0013] (6) Using Y in the training set samples as the input of the STPD-Net model and X as the label, the STPD-Net model is trained to obtain the final dynamic PET image reconstruction model;

[0014] (7) Input Y in the test set sample into the dynamic PET image reconstruction model, and the corresponding dynamic PET tracer activity distribution map can be directly reconstructed and output.

[0015] Furthermore, the expression of the dynamic PET measurement equation is as follows:

[0016] Y=G·X+R

[0017] Where: G is the system response matrix, and R is the random and scattered noise term.

[0018] Furthermore, the Poisson log-likelihood optimization problem with regularization term in step (4) is expressed as follows:

[0019]

[0020] Where: L(Y|X) is the Poisson likelihood term and R() represents the regularization term, I is the total number of detectors, T is the total number of time frames, λ is the penalty coefficient of the regularization term, and Y i,t Represents the element value of the i-th row and t-th column in the dynamic sinusoidal projection data Y, Represents dynamic sinogram projection data The value of the element in row i and column t, E() is the expectation function.

[0021] Furthermore, the expression of the saddle point problem in step (4) is as follows:

[0022]

[0023] Among them: sup<> represents the supremum of the set, L * (Y|X) represents the conjugate of L(Y|X).

[0024] Furthermore, in step (5), the expression for solving the saddle point problem by alternately updating the primary and dual variables using the primary-dual network is as follows:

[0025]

[0026]

[0027] Where: P() represents the main network, D() represents the dual network, represents the dynamic PET tracer activity distribution map X at the kth iteration k The value of the element in row j and column t, They represent the dynamic PET tracer activity distribution diagram X of the k-1th iteration respectively. k-1 The value of the element in row j and column t, represents the dual variable h of the kth iterationk The value of the element in row i and column t, represents the dual variable h of the k-1th iteration k-1 The element value of the i-th row and t-th column in G * represents the conjugate of the system response matrix G, k is a natural number greater than 0, j is a natural number with 1≤j≤N, and N is the total number of pixels in the PET tracer activity distribution map.

[0028] Furthermore, the input of the dual network includes the dual variable h k-1 , Dynamic PET tracer activity distribution diagram X k -1 And the dynamic sinusoidal projection data Y, the output is the iteratively updated dual variable h k , X k-1 After the forward projection operation, h k-1 And Y are spliced ​​in the channel dimension, and then the spliced ​​results are sequentially passed through four 3D spatiotemporal convolution layers to extract spatiotemporal features, and finally the extracted features are combined with h k-1 After channel dimension splicing, the output is h k .

[0029] Furthermore, the input of the main network includes the dual variable h k and dynamic PET tracer activity distribution map X k -1 The output is the iteratively updated dynamic PET tracer activity distribution map X k , h k After the reverse projection operation, k-1 The channel dimension is spliced, and the spliced ​​result is sequentially passed through four 3D spatiotemporal convolution layers to extract spatiotemporal features. Finally, the extracted features are combined with X k-1 After channel dimension splicing, the output is X k .

[0030] Furthermore, the convolution kernel size used in the 3D spatiotemporal convolution layer is 3×3×3.

[0031] Furthermore, the process of training the model in step (6) is as follows:

[0032] 6.1 Initialize model parameters, including the bias vector and weight matrix of each layer, learning rate, and optimizer;

[0033] 6.2 Input the dynamic sinusoidal projection data Y from the training set samples into the model, forward propagate the model output to obtain the corresponding dynamic PET tracer activity distribution map, and calculate the loss function between this result and the label;

[0034] 6.3 Based on the loss function, the optimizer is used to iteratively update the model parameters using the gradient descent method until the loss function converges and the training is completed.

[0035] Furthermore, after the training is completed, the model is verified using the validation set samples, and the model with the best performance on the validation set is used as the final dynamic PET image reconstruction model.

[0036] The present invention proposes to use 3D spatiotemporal convolution to simultaneously extract the spatiotemporal correlation of dynamic PET projection data. Compared with other reconstruction methods based on 2D convolutional neural networks, the present invention can well model the dependencies between different time frames of dynamic projection data, has a good reconstruction and restoration effect on single-frame images, and also has a good grasp of the structural similarities between PET images in different time frames, which is superior to existing dynamic reconstruction methods based on deep learning.

[0037] The present invention proposes to use a principal-dual network to replace the proximal operator to update the principal and dual variables, and expand the principal-dual hybrid gradient descent algorithm into a model-based deep neural network. It ensures a certain degree of interpretability in mathematical derivation while having strong learning representation capabilities. This is the first attempt at a model-based deep learning method in the field of dynamic PET image reconstruction.

[0038] The present invention performs well on ultra-low count dynamic projection data. In practical applications, it can effectively solve the problem of long patient waiting times during the dynamic acquisition process. In experiments, we used both simulated data and clinical mouse scanning data to verify the method proposed by the present invention. We found that the present invention has a good reconstruction effect on single-frame dynamic data with only an ultra-low count of several thousand, indicating that the method proposed by the present invention is particularly suitable for applications in human whole-body dynamic PET imaging and parametric PET imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 Schematic diagram of the steps of the dynamic PET image reconstruction method of the present invention.

[0040] Figure 2 Schematic diagram of the overall structure of the STPD-Net model of the present invention.

[0041] Figure 3 Schematic diagram comparing the reconstruction results of different methods on different slices of low-count dynamic PET projection data; from left to right are the reconstruction images of the MLEM algorithm, the reconstruction images of the KEM-ST method, the reconstruction images of the LPD method, the reconstruction images of the FBPnet method, the reconstruction images of the present invention, and the true value label images; from top to bottom are the third frame, the eighth frame, and the fifteenth frame. DETAILED DESCRIPTION

[0042] In order to describe the present invention more specifically, the technical solution of the present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0043] like Figure 1 As shown, the dynamic PET image reconstruction method based on model-driven deep learning of the present invention includes the following steps:

[0044] Training phase

[0045] (1) Radioactive drugs are injected into the target tissue and detected, dynamic sinusoidal projection data Y is collected, and the corresponding PET tracer activity distribution map X is reconstructed.

[0046] (2) Based on the principle of dynamic PET imaging, a measurement equation model is established:

[0047] Y=G·X+R

[0048] Where: G is the system response matrix, R is the random and scattered noise term, and the system response matrix is ​​calculated using the ray simulation method.

[0049] The above inverse imaging problem is solved using the Poisson log-likelihood optimization method with a regularization term:

[0050]

[0051] in: is the Poisson likelihood term, I represents the total number of detectors, T represents the total number of time frames, represents the mean of the dynamic sinusoidal projection matrix, R() represents the regularization term, and λ is the penalty coefficient of the regularization term.

[0052] Using the properties of the dual problem, the problem is transformed into the form of a saddle point problem:

[0053]

[0054] The saddle point problem is solved by using the STPD-Net method of alternating updating of the primary and dual variables, as shown in Figure 2 As shown in Figure 1, STPD-Net is composed of several cascaded reconstruction modules. Each reconstruction module contains a main network to update the main variable X and a dual network D to update the introduced dual variable:

[0055]

[0056]

[0057] Where: h represents the dual variable, k is the iteration index, G * The conjugate operator of the system matrix, i.e., the back-projection operation, is calculated using the transpose of the system matrix.

[0058] The input of the dual network is the dual variable h k-1 , dynamic PET image X k-1 and projection data Y, the output is the updated dual variable h k , dynamic PET image X k-1 After the forward projection operation and the dual variable h k-1 , the projected data Y is spliced ​​in the channel dimension, and then the spatiotemporal features are extracted through four 3D spatiotemporal convolution layers, with a 3D convolution kernel size of 3×3×3; finally, the dual variable h k-1 After channel dimension splicing, the output is the updated dual variable h k .

[0059] The main network input is a dynamic PET image X k-1 and the dual variable h after the dual network is updated k , the dual variable h after the dual network is updated k After back-projection operation and dynamic PET image X k-1 The channel dimension is spliced, and then the spatiotemporal convolution is performed through four 3D spatiotemporal convolution layers, and then combined with the dynamic PET image X k-1 After channel dimension splicing, the output is the PET reconstructed image X k .

[0060] (3) In the training phase, the reconstruction model is trained using the dynamic PET tracer activity distribution map as the label and the dynamic sinusoidal data as the input.

[0061] First, the model parameters are initialized, including network parameter initialization and reconstructed image initialization. The model parameters are initialized using the Kaiming method, and the initialization image and initialization dual variables are both initialized to 0.

[0062] Then the dynamic sinusoidal projection data Y in the training sample is input into the STPD-Net model, and the output result of the last iterative module is obtained by forward propagation as the result of model reconstruction.

[0063] Then, the MSE loss function between the model output and the dynamic PET tracer activity distribution map and the gradient of the loss function with respect to each variable are calculated. The Adam optimizer is used to update all learnable parameters in the model until the value of the loss function remains basically unchanged, at which point the training is completed.

[0064] Finally, the validation set samples were used for verification, and the model with the best performance in the validation set was used as the final dynamic PET image reconstruction model.

[0065] Inference stage

[0066] (1) Measure or simulate dynamic sinusoidal projection data.

[0067] (2) The dynamic sinusoidal data and the initialization image are used as input, and the trained reconstruction model directly outputs the dynamic PET tracer activity distribution map.

[0068] Below we conduct experiments based on simulated ultra-low count dynamic PET data to verify the effectiveness of this implementation. The dataset contains 40 3D brain template data, each template has a size of 128×128×40×18, and a total of 18 frames were collected. The simulated tracer is 18 F-FDG, obtains the corresponding dynamic sinusoidal data by simulating projection and adding noise. The count of a single-frame sinusoidal graph is about 1e4, of which 33 samples are used as training data, 2 samples are used as validation set data, and the remaining 5 samples are used for testing.

[0069] STPD-Net was implemented using PyTorch 1.7.0 and trained on a Ubuntu server host with a TITAN-X. The optimizer was Adam, the initial learning rate was 0.0001, the batch size was 4, and a total of 200 epochs were trained. The best performing model on the validation set was used for testing.

[0070] Figure 3 The images show reconstructions of different slices of low-count projection data using the present invention and other methods. KEM-ST and FBPnet are currently the mainstream methods. The images show that both the MLEM and KEM-ST algorithms exhibit significant noise. The LPD method, because it fails to consider the correlation between data in different time frames, results in significant differences in image structure between frames. FBPnet has good noise control, but poor structural recovery, and its accuracy needs improvement. The method proposed in this paper not only achieves the best structural recovery, but also achieves the best results in noise suppression and accuracy.

[0071] The above description of the embodiments is intended to facilitate understanding and application of the present invention by those skilled in the art. It is apparent that those skilled in the art can readily make various modifications to the above embodiments and apply the general principles described herein to other embodiments without requiring creative effort. Therefore, the present invention is not limited to the above embodiments. Any improvements or modifications made by those skilled in the art based on the disclosure of the present invention should fall within the scope of protection of the present invention.

Claims

1. A dynamic PET image reconstruction method based on model-driven deep learning, comprising the following steps: (1) Using a detector to detect the biological tissue injected with radioactive drugs, and collecting the corresponding dynamic sinusoidal projection data Y; (2) reconstructing the dynamic sinusoidal projection data Y to obtain the corresponding dynamic PET tracer activity distribution map X; (3) performing steps (1) and (2) multiple times to obtain a large number of samples, each sample including dynamic sinusoidal projection data Y and its corresponding dynamic PET tracer activity distribution map X, and then dividing all samples into a training set, a validation set, and a test set; (4) According to the dynamic PET measurement equation, the dynamic reconstruction problem is transformed into a Poisson log-likelihood optimization problem with a regularization term, and the properties of the dual variables are used to transform this optimization problem into the corresponding saddle point problem; (5) The saddle point problem is solved by alternately updating the primary and dual variables using the primary-dual network, thereby constructing the STPD-Net model for dynamic PET image reconstruction. The model is composed of a cascade of several reconstruction modules, each of which is composed of a primary network and a dual network. (6) Using Y in the training set samples as the input of the STPD-Net model and X as the label, the STPD-Net model is trained to obtain the final dynamic PET image reconstruction model; (7) Input Y in the test set sample into the dynamic PET image reconstruction model, and the corresponding dynamic PET tracer activity distribution map can be directly reconstructed and output.

2. The dynamic PET image reconstruction method according to claim 1, wherein: The expression of the dynamic PET measurement equation is as follows: Y=G·X+R Where: G is the system response matrix, and R is the random and scattered noise term.

3. The dynamic PET image reconstruction method according to claim 2, wherein: The Poisson log-likelihood optimization problem with regularization in step (4) is expressed as follows: Where: L(Y|X) is the Poisson likelihood term and R() represents the regularization term, I is the total number of detectors, T is the total number of time frames, λ is the penalty coefficient of the regularization term, and Y i,t Represents the element value of the i-th row and t-th column in the dynamic sinusoidal projection data Y, Represents dynamic sinogram projection data The value of the element in row i and column t, E() is the expectation function.

4. The dynamic PET image reconstruction method according to claim 3, wherein: The expression of the saddle point problem in step (4) is as follows: Among them: sup<> represents the supremum of the set, L * (Y|X) represents the conjugate of L(Y|X).

5. The dynamic PET image reconstruction method according to claim 4, characterized in that: In step (5), the expression for solving the saddle point problem by alternately updating the primary and dual variables using the primary-dual network is as follows: Where: P() represents the main network, D() represents the dual network, represents the dynamic PET tracer activity distribution map X at the kth iteration k The value of the element in row j and column t, They represent the dynamic PET tracer activity distribution diagram X of the k-1th iteration respectively. k -1 The value of the element in row j and column t, represents the dual variable h of the kth iteration k The value of the element in row i and column t, represents the dual variable h of the k-1th iteration k-1 The element value of the i-th row and t-th column in G * represents the conjugate of the system response matrix G, k is a natural number greater than 0, j is a natural number with 1≤j≤N, and N is the total number of pixels in the PET tracer activity distribution map.

6. The dynamic PET image reconstruction method according to claim 5, characterized in that: The input of the dual network includes the dual variable h k-1 , Dynamic PET tracer activity distribution diagram X k-1 And the dynamic sinusoidal projection data Y, the output is the iteratively updated dual variable h k , X k-1 After the forward projection operation, h k-1 And Y are spliced ​​in the channel dimension, and then the spliced ​​results are sequentially passed through four 3D spatiotemporal convolution layers to extract spatiotemporal features, and finally the extracted features are combined with h k -1 After channel dimension splicing, the output is h k .

7. The dynamic PET image reconstruction method according to claim 5, characterized in that: The input of the main network includes the dual variable h k and dynamic PET tracer activity distribution map X k-1 The output is the iteratively updated dynamic PET tracer activity distribution map X k , h k After the reverse projection operation, k-1 The channel dimension is spliced, and the spliced ​​result is sequentially passed through four 3D spatiotemporal convolution layers to extract spatiotemporal features. Finally, the extracted features are combined with X k-1 After channel dimension splicing, the output is X k .

8. The dynamic PET image reconstruction method according to claim 6 or 7, characterized in that: The convolution kernel size used in the 3D spatiotemporal convolution layer is 3×3×3.

9. The dynamic PET image reconstruction method according to claim 1, wherein: The process of training the model in step (6) is as follows: 6.1 Initialize model parameters, including the bias vector and weight matrix of each layer, learning rate, and optimizer; 6.2 Input the dynamic sinusoidal projection data Y from the training set samples into the model, forward propagate the model output to obtain the corresponding dynamic PET tracer activity distribution map, and calculate the loss function between this result and the label; 6.3 Based on the loss function, the optimizer is used to iteratively update the model parameters using the gradient descent method until the loss function converges and the training is completed.

10. The dynamic PET image reconstruction method according to claim 9, characterized in that: After training is completed, the model is verified using the validation set samples, and the model with the best performance on the validation set is used as the final dynamic PET image reconstruction model.

Images