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

Abstract

The invention discloses a dynamic PET image reconstruction method based on model-driven deep learning, which utilizes 3D space-time convolution to simultaneously extract the relevance of a time domain and a space domain of dynamic projection data, and front and back projection operators are integrated into a reconstruction network, so that the method has strong physical constraint and interpretability; the invention splits the dynamic PET image reconstruction problem into a plurality of cascaded reconstruction blocks, wherein each reconstruction block comprises a main network for updating a main image domain variable and a dual network for updating a dual measurement domain variable. The invention can obtain a high-quality dynamic PET tracer activity distribution image from the reconstruction of the ultra-low counting dynamic PET projection data, and solves the problems of poor interpretability and poor reconstruction effect of the current mainstream method.

Description

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

Technical Field

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

Background Art

Dynamic PET imaging can quantitatively characterize physiological parameters in biological tissues, and plays an indispensable role in tumor detection, cardiac disease characterization, and drug development. However, image reconstruction based on dynamic measurement data is very challenging due to the ill-posed nature of the PET reconstruction problem and the low count characteristics of single-frame images of dynamic PET data. Especially in the early time frame, the count of a single-frame image is very low, and the reconstructed image is severely affected by noise.

In addition, with the continuous improvement of imaging technology and detector level, ultrafast time frame imaging hardware technology can be realized, but the existing reconstruction algorithm cannot reconstruct high-quality ultra-low count projection data well. Traditional dynamic image reconstruction algorithms such as filtered back projection and maximum likelihood expectation maximum algorithm cannot model time domain signals, resulting in large differences in reconstructed images in different frames.

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 or 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 the reconstruction effect by introducing time priors on the basis of the kernel method. However, in practical applications, the performance is not very good for ultra-fast time frames and ultra-low count data. Part of the reason is that the prior information of this method is obtained from the data of a single reconstruction, and the advantages of data-driven methods are not combined. 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 a relatively good dynamic reconstruction effect by combining the traditional filtered back projection algorithm with a denoising neural network. However, this method has poor interpretability and generalization because it does not take into account the system matrix constraints related to the physical characteristics of the PET instrument.

It can be seen that the existing technologies either do not make good use of the correlation of dynamic projection data in the time domain, or lack physical constraints, resulting in unstable reconstruction results, and poor performance in the case of ultra-low counts, thus limiting the development of ultrafast time frame PET imaging technology. In order to obtain better reconstruction quality, the existing dynamic reconstruction methods often require projection data with a longer acquisition time. However, the long scanning time will inevitably introduce motion artifacts, further deteriorating the quality of the reconstructed image.

Summary of the invention

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, effectively utilize the spatiotemporal correlation of dynamic projection data, and combine the constraints of the physical projection matrix, and has strong interpretability.

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

(1) Using a detector to detect the biological tissue injected with radioactive drugs, and acquiring 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 comprising 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 the optimization problem into the corresponding saddle point problem;

(5) The main and dual networks are used to alternately update the main and dual variables to solve the above saddle point problem, 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 main 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.

Furthermore, 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.

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

R()表示正则项，I为探测器的总数量，T为时间帧的总数量，λ为正则项惩罚系数，Y_i,t表示动态正弦图投影数据Y中第i行第t列的元素值，

表示动态正弦图投影数据

中第i行第t列的元素值，

E()为期望函数。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 regularization term penalty coefficient, Yi _,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 the i-th row and t-th column,

E() is the expectation function.

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

Where: sup<> represents the supremum of a set, and L ^* (Y|X) represents the conjugate of L(Y|X).

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

表示第k次迭代的动态PET示踪剂活度分布图X ^k中第j行第t列的元素值，

分别表示第k-1次迭代的动态PET示踪剂活度分布图X ^k-1中第j行第t列的元素值，

表示第k次迭代的对偶变量h^k中第i行第t列的元素值，

表示第k-1次迭代的对偶变量h^k-1中第i行第t列的元素值，G^*表示系统响应矩阵G的共轭，k为大于0的自然数，j为自然数且1≤j≤N，N为PET示踪剂活度分布图的总像素点数量。Where: P() represents the main network, D() represents the dual network,

represents the value of the element in the jth row and tth column of the dynamic PET tracer activity distribution diagram ^Xk of the kth iteration,

They represent the element value of the jth row and tth column in the dynamic PET tracer activity distribution diagram Xk ^-1 of the k-1th iteration,

represents the value of the element in the i-th row and t-th column of the dual variable h ^k of the k-th iteration,

represents the element value of the i-th row and t-th column in the dual variable h ^k-1 of the k-1th iteration, G ^* represents the conjugate of the system response matrix G, k is a natural number greater than 0, j is a natural number and 1≤j≤N, and N is the total number of pixels in the PET tracer activity distribution map.

Furthermore, the input of the dual network includes the dual variable h ^k-1 , the dynamic PET tracer activity distribution map X ^{k -1} and the dynamic sinusoidal projection data Y, and the output is the iteratively updated dual variable h ^k , X ^k-1 is first subjected to a forward projection operation and then spliced with h ^k-1 and Y in the channel dimension, and then the spliced result is sequentially passed through four 3D spatiotemporal convolutional layers for spatiotemporal feature extraction, and finally the extracted features are spliced with h ^k-1 in the channel dimension and the output is h ^k .

Furthermore, the input of the main network includes the dual variable h ^k and the dynamic PET tracer activity distribution map X ^{k -1} , and the output is the iteratively updated dynamic PET tracer activity distribution map X ^k . h ^k is first back-projected and then concatenated with X ^k-1 in the channel dimension. The concatenated result is then sequentially passed through four 3D spatiotemporal convolutional layers for spatiotemporal feature extraction. Finally, the extracted features are concatenated with X k ^-1 in the channel dimension and the output is X ^k .

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

Furthermore, 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 in 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 the result and the label;

6.3 Based on the loss function, use the optimizer to iteratively update the model parameters through the gradient descent method until the loss function converges and the training is completed.

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.

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 dependency 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 similarity between PET images in different time frames, which is superior to the existing dynamic reconstruction methods based on deep learning.

The present invention proposes to use a principal-dual network to replace the proximal operator to update the principal variable and the dual variable, and expands 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 a strong learning representation ability. This is the first attempt at a model-based deep learning method in the field of dynamic PET image reconstruction.

The present invention performs very well on dynamic projection data with ultra-low counts. In practical applications, it can effectively solve the problem of long waiting time for patients during the dynamic acquisition process. In the experiment, we used both simulation data and scanning data of clinical mice to verify the method proposed by the present invention, and found that the present invention has a good reconstruction effect on dynamic data with ultra-low counts of only a few thousand in a single frame, indicating that the method proposed by the present invention will be particularly suitable for the application of human whole-body dynamic PET imaging and parametric PET imaging.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic flow chart of the steps of the dynamic PET image reconstruction method of the present invention.

FIG2 is a schematic diagram of the overall structure of the STPD-Net model of the present invention.

Figure 3 is a 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 image of the MLEM algorithm, the reconstruction image of the KEM-ST method, the reconstruction image of the LPD method, the reconstruction image of the FBPnet method, the reconstruction image of the present invention, and the true value label image, and from top to bottom are the third frame, the eighth frame, and the fifteenth frame.

DETAILED DESCRIPTION

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.

As shown in FIG1 , the dynamic PET image reconstruction method based on model-driven deep learning of the present invention comprises the following steps:

Training phase

(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.

(2) According to the principle of dynamic PET imaging, the measurement equation model is established:

Y＝G·X+R

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.

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

为泊松似然项，I表示总探测器的数量，T表示总时间帧的数量，

表示动态正弦图投影矩阵的均值，R()代表正则项，λ为正则项惩罚系数。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.

Using the properties of the dual problem, the problem is transformed into a saddle point problem:

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

Where: h represents the dual variable, k is the iteration index, G ^* represents the conjugate operator of the system matrix, that is, the back-projection operation, and the calculation is performed using the transpose of the system matrix.

The input of the dual network is the dual variable h ^k-1 , the dynamic PET image X ^k-1 and the projection data Y, and the output is the updated dual variable h ^k . The dynamic PET image X ^k-1 is first forward projected and then concatenated with the dual variable h ^k-1 and the projection data Y in the channel dimension. Then, it is subjected to four 3D spatiotemporal convolution layers for spatiotemporal feature extraction, and the 3D convolution kernel size is 3×3×3. Finally, it is concatenated with the dual variable h ^k-1 in the channel dimension and the updated dual variable h ^k is output.

The input of the main network is the dynamic PET image ^Xk-1 and the dual variable ^hk updated by the dual network. The dual variable ^hk updated by the dual network is first back-projected and then concatenated with the dynamic PET image ^Xk-1 in the channel dimension. Then it is extracted through four 3D spatiotemporal convolution layers and then concatenated with the dynamic PET image ^Xk-1 in the channel dimension to output the PET reconstructed image ^Xk .

(3) In the training phase, the dynamic PET tracer activity distribution map is used as a label and the dynamic sinusoidal data is used as an input to train the reconstruction model.

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 the initialization dual variables are both initialized to 0.

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.

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, and 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.

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.

Inference Phase

(1) Measure or simulate the dynamic sinusoidal projection data.

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

Below we conduct experiments based on simulated ultra-low count dynamic PET data to verify the effectiveness of this implementation; the data set contains 40 3D brain template data, each template is 128×128×40×18 in size, and a total of 18 frames are collected. The simulated tracer is ¹⁸ F-FDG, and the corresponding dynamic sinusoidal data are obtained 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.

STPD-Net was implemented using pytorch 1.7.0 and trained on an 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.

Figure 3 shows the reconstruction images of the present invention and other methods on different slices of low-count projection data, among which KEM-ST and FBPnet methods are the current mainstream methods. It can be seen from the figure that both the MLEM algorithm and the KEM-ST algorithm show very serious noise. The LPD method does not consider the association between data in different time frames, resulting in large differences in image structure between different frames; FBPnet has better control over noise, but the structural information cannot be well restored, and the accuracy needs to be improved. The method proposed by the present invention not only has the best recovery performance in structure, but also achieves the best results in noise suppression and accuracy.

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 familiar with 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. Improvements and modifications made by those skilled in the art to the present invention based on the disclosure of the present invention should be within the protection scope of the present invention.

Claims (10)

1. A dynamic PET image reconstruction method based on model-driven deep learning comprises the following steps:

(1) Detecting biological tissues injected with the radiopharmaceuticals by using a detector, and acquiring corresponding dynamic sinogram projection data Y;

(2) Reconstructing projection data Y of the dynamic sinogram to obtain a corresponding activity distribution map X of the dynamic PET tracer;

(3) Executing the steps (1) and (2) for multiple times to obtain a large number of samples, wherein each sample comprises dynamic sinogram projection data Y and a dynamic PET tracer activity distribution diagram X corresponding to the dynamic sinogram projection data Y, and further dividing all the samples into a training set, a verification set and a test set;

(4) Converting the dynamic reconstruction problem into a Poisson log-likelihood optimization problem with a regular term according to a dynamic PET measurement equation, and converting the optimization problem into a corresponding saddle point problem by using the property of a dual variable;

(5) Alternately updating a main variable and a dual variable by using a main dual network to solve the saddle point problem, thereby constructing an STPD-Net model for reconstructing the dynamic PET image, wherein the model is formed by cascading a plurality of reconstruction modules, and each reconstruction module is formed by connecting the main network and the dual network;

(6) Training the STPD-Net model by using Y in the training set sample as the input of the STPD-Net model and using X as a label, thereby obtaining a final dynamic PET image reconstruction model;

(7) And inputting Y in the test set sample into the dynamic PET image reconstruction model, and directly reconstructing and outputting a corresponding dynamic PET tracer activity distribution diagram.

2. The dynamic PET image reconstruction method according to claim 1, characterized in that: the expression of the dynamic PET measurement equation is as follows:

Y＝G·X+R

wherein: g is the system response matrix and R is the random and scattered noise terms.

3. The dynamic PET image reconstruction method according to claim 2, characterized in that: the poisson log-likelihood optimization problem expression with the regular term in the step (4) is as follows:

wherein: l (Y | X) is a Poisson likelihood term and

r () represents a regular term, I is the total number of detectors, T is the total number of time frames, λ is a regular term penalty factor, Y _i,t Represents the value of an element in the ith row and the tth column of the dynamic sinogram projection data Y @>

Representing dynamic sinogram projection data pick>

The value of the element in the ith row and the tth column,

e () is the desired function.

4. The dynamic PET image reconstruction method according to claim 3, characterized in that: the expression of the saddle point problem in the step (4) is as follows:

wherein: sup<>Denotes the set supremum, 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 the step (5), the expression for solving the saddle point problem by alternately updating the primary variable and the dual variable by using the primary-dual network is as follows:

wherein: p () denotes the master network, D () denotes the dual network,

dynamic PET tracer activity profile X representing the kth iteration ^k The value of an element in the jth row and the tth column is greater than>

Dynamic PET tracer activity profile X for the k-1 th iteration, respectively ^{k -1} The value of an element in the jth row and the tth column is greater than>

Dual variable h representing the kth iteration ^k The value of an element in the ith row and the tth column is greater or less>

Dual variable h representing the k-1 iteration ^k-1 In the ith row and the tth columnElement value of (1), G ^* And (3) representing the conjugation of a system response matrix G, wherein k is a natural number greater than 0, j is a natural number, j is greater than or equal to 1 and less than or equal to N, and N is the total pixel point number of the PET tracer activity distribution diagram.

6. The dynamic PET image reconstruction method according to claim 5, characterized in that: the input of the dual network comprises a dual variable h ^k-1 Dynamic PET tracer activity distribution diagram X ^k-1 And dynamic sinogram projection data Y output as dual variable h after iterative update ^k ，X ^k-1 After the forward projection operation, h ^k-1 Y splicing the channel dimensions, then sequentially performing space-time feature extraction on the spliced result through four 3D space-time convolution layers, and finally performing space-time feature extraction on the extracted features and h ^{k -1} Output is h after channel dimension splicing ^k 。

7. The dynamic PET image reconstruction method according to claim 5, characterized in that: the input of the main network comprises a dual variable h ^k And dynamic PET tracer activity profile X ^k-1 And the output is a dynamic PET tracer activity distribution diagram X after iterative update ^k ，h ^k After the back projection operation, the X-ray source is connected with the X-ray source ^k-1 Splicing channel dimensions, sequentially performing space-time feature extraction on the spliced result through four 3D space-time convolution layers, and finally performing space-time feature extraction on the extracted features and X ^k-1 The output is X after channel dimension splicing ^k 。

8. The dynamic PET image reconstruction method according to claim 6 or 7, characterized in that: the convolution kernel size adopted by the 3D space-time convolution layer is 3 multiplied by 3.

9. The dynamic PET image reconstruction method according to claim 1, characterized in that: the process of training the model in the step (6) is as follows:

6.1 initializing model parameters, including a bias vector and a weight matrix of each layer, a learning rate and an optimizer;

6.2, inputting the projection data Y of the dynamic sinogram in the training set sample into a model, outputting the forward propagation of the model to obtain a corresponding activity distribution map of the dynamic PET tracer, and calculating a loss function between the result and the label;

and 6.3, continuously and iteratively updating model parameters by using an optimizer according to the loss function through a gradient descent method until the loss function is converged, and finishing training.

10. The dynamic PET image reconstruction method according to claim 9, characterized in that: and after the training is finished, verifying the model by using a verification set sample, and taking the model which best expresses on the verification set as a final dynamic PET image reconstruction model.

Images