, with the ability to observe molecular-level activities inside the tissue through the injection of specific radioactive tracers. Though PET has high sensitivity compared with other imaging modalities, its image resolution and signal to noise ratio (SNR) are still low due to various physical degradation factors and low coincident-photon counts detected. Improving PET image quality is essential, especially in applications like small lesion detection, brain imaging and longitudinal studies. Over the past decades, multiple advances have been made in PET system instrumentation, such as exploiting time of flight (TOF) information, enabling depth of interaction capability  and extending the solid angle coverage [6, 7].
With the wide adoption of iterative reconstruction in clinical scanners, more accurate point spread function (PSF) modeling can be used to take various degradation factors into consideration . In addition, various post processing approaches and iterative reconstruction methods have been developed by making use of local patch statistics, prior anatomical or temporal information. Denoising methods, such as the HYPR processing , non-local mean denoising [10, 11] and guided image filtering 
have been developed and show better performance in bias-variance tradeoff or partial volume correction than the conventional Gaussian filtering. In regularized image reconstruction, entropy or mutual information based methods[13, 14, 15], segmentation based methods[16, 17], and gradient based methods [18, 19] have been developed by penalizing the difference between the reconstructed image and the prior information in specific domains. The Bowsher’s method  adjusts the weight of the penalty based on similarity metrics calculated from prior images. Methods based on sparse representations [21, 22, 23, 24, 25, 26], have also shown better image qualities in both static and dynamic reconstructions. Most of the aforementioned methods require prior information from the same patient which is not always available due to instrumentation limitation or long scanning time, which may hamper the practical application of these methods. Recently a new method is developed to use information in longitudinal scans , but can only be applied to specific studies.
In this paper, we explore the potential of using existing inter-patient information via deep neural network to improve PET image reconstruction. Over the past several years, deep neural networks have been widely and successfully applied to computer vision tasks, such as image segmentation , object detection 
and image super resolution, due to the availability of large data sets, advances in optimization algorithms and emerging of effective network structures. Recently, it has been applied to medical imaging, such as image denoising and artifact reduction, using convolutional neural network (CNN) [31, 32, 33, 34] or generative adversarial network (GAN) . It showed comparable or superior results to the iterative reconstruction but at a faster speed. In this paper, we propose a new framework to integrate deep CNN in PET image reconstruction. The network structure is a combination of U-net structure  and the residual network . Different from existing CNN based image denoising methods, we use a CNN trained with iterative reconstructions of low-counts data as the input and high-counts reconstructions as the label to represent the unknown PET image to be reconstructed. Rather than feeding a noisy image into the CNN, we use the CNN to define the feasible set of valid PET images. To our knowledge, this is the first of its kind in the applications of neural network in medical imaging. The solution is formulated as the solution of a constrained optimization problem and sought by using the alternating direction method of multipliers (ADMM) algorithm . The proposed method is validated using both simulation and hybrid real data.
The main contributions of this paper include (1) using dynamic data of prior patients to train a network for PET denoising and (2) proposing to incorporate the neural network into the iterative reconstruction framework and demonstrating better performance than the denoising approach. This paper is organized as follows. Section 2 introduces the theory and optimization algorithm. Section 3 describes the Monte Carlo simulations and real data used in the evaluation. Experimental results are shown in Section 4, followed by discussions in Section 5. Finally conclusions are drawn in Section 6.
Ii-a PET data model
In PET image reconstruction, the measured data
can be modeled as a collection of independent Poisson random variables and its meanis related to the unknown image through an affine transform
is the detection probability matrix, withdenoting the probability of photons originating from voxel being detected by detector . denotes the expectation of scattered events, and denotes the expectation of random coincidences. is the number of lines of response (LOR) and is the number of pixels in image space. The log-likelihood function can be written as
The maximum likelihood estimate of the unknown imagecan be found by
Ii-B Representing PET images using neural network
Previously, the kernel method  used a kernel representation to represent the image , through which the prior temporal or anatomical information can be embedded into the kernel matrix . Inspired by this idea, here we represent the unknown image as
where represents the neural network and denotes the input to the neural network. Through this representation, inter-patient information and intra-patient information can be included into the iterative reconstruction framework through pre-training the neural network using existing data.
Our network implemented in this work is based on the U-net structure 
and also includes the batch normalization layer. The overall network architecture is summarized in Fig.2. It consists of repetitive applications of 1) x
convolutional layer, 2) batch normalization layer, 3) ReLU layer, 4) convolutional layer with stride 2 for down-sampling, 5) transposed convolutional layer with stride 2 for up-sampling, and 6) identity mapping layer that adds the left-side feature layer to the right-side. In our implementation, there are three major modifications compared to the original U-net:
using convolutional layer with stride 2 to down-sample the image instead of using max pooling layer, to construct a fully convolutional network;
directly adding the left side feature to the right side instead of concatenating, to reduce the number of training parameters;
connecting the input to the output, to construct a residual network .
The left-hand side of the architecture aims to compress the input path layer by layer, an “encoder” part, while the right-hand side is to expand the path, a “decode” part. This neural network has 19 convolutional layers in total and the largest feature size is 512. To reduce computational cost, the network denoises PET image one slice at a time. However, the input layer has five channels to provide information from 4 neighboring axial slices for effective noise removal. We have found that if the input did not contain the neighboring axial information, there will be artifacts in the axial direction. The network is trained with reconstructed images from low counts data as the input and the images reconstructed from high counts data as the label.
The maximum likelihood estimate of the unknown image can be calculated as
The objective function in (7) is difficult to solve due to the nonlinearity of the neural network representation. Here we transfer it to the constrained format as below
We use the Augmented Lagrangian format for the constrained optimization problem in (8) as
which can be solved by the ADMM algorithm iteratively in three steps
where and is calculated by
It can be verified that the constructed surrogate function fulfills the following two conditions:
The final update equation for pixel after maximizing (17) is
Subproblem (11) is a non-linear least square problem. In order to solve it, we need to compute the gradient of the objective function with respect to the input . As it is difficult to calculate the Jacobian matrix or Hessian matrix of the objective function with respect to the input in current network platform, we use a first-order method as follows
where is the step size. In our implementation, was chosen so that the objective function in subproblem (11) can be monotonic decreasing. For our neural network, the input have five channels, which include four neighboring axial slices. Therefore, equation (19) should be modified to include the first-order gradients from the other four neighboring slices. The final update equation is changed to
where is the number of channels, , and
is the spatial input size. In order to accelerate the convergence speed, Nesterov momentum method was used in subproblem (11) . In our implementation, we run one iteration for subproblem (10) and then run five iterations for subproblem (11). As subproblem (11) is a non-linear problem, it is very easy to be trapped into a local minimum and it is thus essential to assign a good initial for . In our implementation, we first ran MLEM for 30 iterations and used its CNN output as the initial for . The overall algorithm flowchart is presented in Algorithm 1.
Ii-D Implementation details and reference methods
The neural network was implemented using TensorFlow 1.0 on a NVIDIA GTX 1080Ti. The network input size isand the output size is . During training, Adam algorithm  was used as the optimizer and the cost function was calculated as the L2 norm between the network outputs and the label images. The first-order gradient used in subproblem (11) was implemented using the tf.gradient function in TensorFlow.
We compared the proposed methods with the post-reconstruction Gaussian filtering and a penalized reconstruction. For the penalized reconstruction, the fair penalty was used with the form
The fair penalty approaches the L-1 penalty when and is similar to the quadratic penalty when . In our implementation, was set to be of the mean image intensity in order to have the edge preserving capability. MAP EM algorithm was used in the penalized reconstruction . In order to accelerate the convergence, 10 iterations of MLEM algorithm was used for “warming up” before running the MAP EM algorithm.
Iii Experimental setup
Iii-a Simulation study
The computer simulation modeled the geometry of a GE 690 scanner . The scanner consists of LYSO crystals forming a ring of diameter of 81 cm with an axial field of view (FOV) of 157 mm. The crystal size for this scanner is . Nineteen XCAT phantoms with different organ sizes and genders were employed in the simulation . Apart from the major organs, thirty hot spheres of diameters ranging from 12.8 mm to 22.4 mm were inserted into eighteen phantoms as lung lesions for the training images. For the test image, five lesions with diameter 12.8 mm were inserted. The time activity curves (TAC) of different tissues were generated mimicking an FDG scan using a two-tissue-compartment model with an analytic blood input function . In order to simulate the population difference, each kinetic parameter was modelled as a Gaussian variable with coefficient of variation equal to 0.1. The mean values of the kinetic parameters are presented in Table I [45, 46]. The TACs using the mean kinetic parameters are shown in Fig. 4. The system matrix used in the data generation and image reconstruction was computed by using the multi-ray tracing method , which modeled the inter-crystal photon penetration. The image matrix size is 128 128 49 and the voxel size is 3.27 3.27 3.27 . Noise-free sinogram data were generated by forward-projecting the ground-truth images using the system matrix and the attenuation map. Poisson noise was then introduced to the noise-free data by setting the total count level to be equivalent to an 1-hour FDG scan with 5 mCi injection. Uniform random and scatter events were simulated and accounted for 60% of the noise free prompt data in all time frames to match those observed in real data-sets. During image reconstruction, all the correction factors were assumed to be known exactly.
To generate the training data, forty-minutes data from 20 min to 60 min post injection were combined into a high count sinogram and reconstructed as the label image for training. The high count data was down-sampled to 1/10th of the counts and reconstructed as the input image. In order to account for different noise levels, images reconstructed at iteration 20, 40, 60 using ML EM algorithm were all used in the training phase. In total ( of slices per phantom) ( of phantoms) (
of different iterations) training data pairs were generated. Different rotations and translations were applied to each training pair to enable larger data capacity for the training. The training data set was separated randomly into 45 batches for every epoch. In total 1000 epochs were run. Three training pair examples are shown in Fig.6.
During the evaluation, 20 low-counts realizations of the testing phantom, generated by pooling the last 40 min data together and resampling with a 1/10 ratio, were reconstructed using different methods. For quantitative comparison, contrast recovery (CR) vs. the standard deviation (STD) curves were plotted. The CR was computed from the lung lesion regions as
where is the number of realizations, is the average uptake of all the lung lesions in the test phantom. The background STD was computed as
where is the average of the background ROI means over realizations, and is the total number of background ROIs chosen.
Iii-B Hybrid real data
Six patient data sets of one hour FDG dynamic scan acquired on a GE 690 scanner with 5 mCi dose injection were employed in this study. Training and validation data were generated in the same way as that in the simulation. The system matrix used in the reconstruction is the same as the one used in the simulation. Normalization, attenuation correction, randoms and scatters were generated using the manufacturer software and included in image reconstruction. Five patient data sets were used in the training and the last one was left for validation. As no ground truth exist in the real data-sets, 27 lesions were inserted in the training data and 5 in the testing data to generate the hybrid real data-sets for quantitative analysis. The diameters of the lesions inserted into the training data sets range from 12.8 mm to 22.4 mm and the diameter for the lesions inserted in the testing data is 12.8 mm. The intensity of all the lesions were simulated as a Gaussian random variable with coefficient of variation equal to 0.2 to simulate the population difference. In order to increase the training samples, for each patient data set we have generated five low-dose realizations from the high-counts data. Training pairs of iteration 20, 40, 60 were also included to account for different noise levels. In total ( of slices per data set) ( of patients) ( of different iterations) ( of realizations) training data sets were generated with different rotations and translations. Three pairs of the training examples are shown in Fig. 8.
Twenty realizations of the low dose data sets were resampled from the testing data and reconstructed to evaluate the noise performance. Forty-seven background ROIs were chosen in the liver region to calculate the STD as presented in (23). For lesion quantification, images with and without the inserted lesion were reconstructed and the difference was taken to obtain the lesion only image and compared with the ground truth. The lesion contrast recovery was calculated as in (22).
Iv-a Simulation results
Fig. 12 shows three orthogonal slices of the reconstructed images using different methods. From the image appearance, we can see that the proposed iterative CNN method can generate images with a higher lung lesion uptake and reveal more vessel details in the lung region as compared with the CNN denoising method. This is beneficial as the CNN method is criticized for over-smoothing and losing small structures due to the L2 norm used as the cost function. Both CNN approaches are better than the traditional Gaussian post filtering method as the images have less noise but also keep all the detailed features, such as the thin myocardium regions. The penalized reconstruction result has a high lesion uptake, but also has some noise spots showing up in different regions. These observations are consistent with the quantitative results shown in Fig. 14. In terms of the bias-variance trade-off, the proposed iterative CNN method has the best performance among all methods.
Iv-B Real data results
Fig. 18 shows three orthogonal slices of the reconstructed images using the real data-set by different methods. We can see that the uptake of the inserted lesion in the iterative CNN method is higher than the CNN denoising method, same conclusion as in the simulation study. In addition, the iterative CNN method produced the clearest image details in the spinal regions compared with all other methods. The penalized reconstruction can preserve lesion uptake and reduce image noise, but can also present cartoon-like patterns, especially in the high uptake regions. The results using CNN methods are more natural with no obvious artifacts. The quantitative results are presented in Fig. 20. From the figure, we can see that about two-fold STD reduction can be achieved by the CNN methods, compared with the Gaussian filtering method.
Many prior works have used CNN in CT or MRI denoising. Here we use CNN as the image representation and embedded it into PET iterative reconstruction framework, where no prior arts exist. Compared with the CNN denoising approach, the proposed iterative CNN method has a constraint from the measured data, which can help recover some small features that are removed or annihilated by the image denoising methods. Higher contrast recovery of the lesions shown in both simulation and real data sets demonstrate this benefit.
Previously the kernel method has been applied successfully in both static and dynamic PET image reconstructions. When using the kernel method, we need to explicitly specify the basis function when constructing the kernel matrix. This is not needed for CNN and the whole network representation is more data-driven. The biggest advantage of the proposed method is that more generalized prior information, such as the inter-patient scanning information, can be included in the image representation. In addition, when the prior information is from multiple resources, such as both the temporal and anatomical information, it is hard to specify how to combine those information in the kernel method. For neural network, we can use multiple input channels to aggregate the information and let the network decide the optimum combination in the training phase.
As for the optimization process of the proposed iterative CNN method, the most challenging part is Subproblem (11) as it is a non-linear problem. As the computation of the Jacobian matrix is difficult due to the platform limitation, currently we choose a first-order method with Nesterov momentum to solve it. However, it is easy to get trapped in local minimums. In our experiment, we found that if the initial value of is a uniform image, the result is very poor. In our proposed solution, we used the EM results after 30 iterations as the input, which can make the results more stable. Better optimization methods and more effective initial choosing strategies need further investigations.
The network structure used in this study is the modified U-net structure, which is a fully convolutional network. One drawback of CNN is that it will remove some of the small structures in the final output. Though our proposed iterative framework can overcome this issue, better network structures, which can preserve more features, can make our proposed iterative framework work better. For example, our proposed approach can be also fit for GAN. After the generator network is trained through GAN, it can be included into the iterative framework based on the proposed method. Besides, though the data model used here is PET, it can also be used in CT or MRI reconstruction framework.
In this work, we proposed an iterative PET image reconstruction framework by using convolutional neural network representation. Both simulated XCAT data and real data sets were used in the evaluation. Quantitative results show that the proposed iterative CNN method performs better than the CNN denoising method as well as the Gaussian filter and penalized reconstruction methods regarding contrast recovery vs. noise trade-offs. Future work will focus on more clinical data sets evaluation.
-  T. Beyer, D. W. Townsend, T. Brun et al., “A combined PET/CT scanner for clinical oncology,” Journal of nuclear medicine, vol. 41, no. 8, pp. 1369–1379, 2000.
-  R. N. Gunn, M. Slifstein, G. E. Searle et al., “Quantitative imaging of protein targets in the human brain with PET,” Physics in medicine and biology, vol. 60, no. 22, p. R363, 2015.
-  J. Machac, “Cardiac positron emission tomography imaging,” Seminars in nuclear medicine, vol. 35, no. 1, pp. 17–36, 2005.
-  J. S. Karp, S. Surti, M. E. Daube-Witherspoon et al., “Benefit of time-of-flight in PET: experimental and clinical results,” Journal of Nuclear Medicine, vol. 49, no. 3, pp. 462–470, 2008.
-  Y. Yang, J. Qi, Y. Wu et al., “Depth of interaction calibration for PET detectors with dual-ended readout by PSAPDs,” Physics in medicine and biology, vol. 54, no. 2, p. 433, 2009.
-  J. K. Poon, M. L. Dahlbom, W. W. Moses et al., “Optimal whole-body PET scanner configurations for different volumes of LSO scintillator: a simulation study,” Physics in medicine and biology, vol. 57, no. 13, p. 4077, 2012.
-  K. Gong, S. Majewski, P. E. Kinahan et al., “Designing a compact high performance brain pet scanner—simulation study,” Physics in medicine and biology, vol. 61, no. 10, p. 3681, 2016.
-  K. Gong, J. Zhou, M. Tohme et al., “Sinogram blurring matrix estimation from point sources measurements with rank-one approximation for fully 3D PET,” IEEE Transactions on Medical Imaging, 2017.
-  B. T. Christian, N. T. Vandehey, J. M. Floberg et al., “Dynamic pet denoising with hypr processing,” Journal of Nuclear Medicine, vol. 51, no. 7, pp. 1147–1154, 2010.
-  J. Dutta, R. M. Leahy, and Q. Li, “Non-local means denoising of dynamic pet images,” PloS one, vol. 8, no. 12, p. e81390, 2013.
-  C. Chan, R. Fulton, R. Barnett et al., “Postreconstruction nonlocal means filtering of whole-body pet with an anatomical prior,” IEEE Transactions on medical imaging, vol. 33, no. 3, pp. 636–650, 2014.
-  J. Yan, J. C.-S. Lim, and D. W. Townsend, “Mri-guided brain pet image filtering and partial volume correction,” Physics in medicine and biology, vol. 60, no. 3, p. 961, 2015.
-  J. Nuyts, “The use of mutual information and joint entropy for anatomical priors in emission tomography,” in 2007 IEEE Nuclear Science Symposium Conference Record, vol. 6. IEEE, 2007, pp. 4149–4154.
-  J. Tang and A. Rahmim, “Bayesian PET image reconstruction incorporating anato-functional joint entropy,” Physics in Medicine and Biology, vol. 54, no. 23, p. 7063, 2009.
-  S. Somayajula, C. Panagiotou, A. Rangarajan et al., “PET image reconstruction using information theoretic anatomical priors,” IEEE Transactions on Medical Imaging, vol. 30, no. 3, pp. 537–549, 2011.
-  C. Comtat, P. E. Kinahan, J. A. Fessler et al., “Clinically feasible reconstruction of 3D whole-body PET/CT data using blurred anatomical labels,” Physics in Medicine and Biology, vol. 47, no. 1, p. 1, 2001.
-  K. Baete, J. Nuyts, W. Van Paesschen et al., “Anatomical-based FDG-PET reconstruction for the detection of hypo-metabolic regions in epilepsy,” IEEE Transactions on Medical Imaging, vol. 23, no. 4, pp. 510–519, 2004.
-  M. J. Ehrhardt, P. Markiewicz, M. Liljeroth et al., “PET reconstruction with an anatomical MRI prior using parallel level sets,” IEEE Transactions on Medical Imaging, vol. 35, no. 9, pp. 2189–2199, 2016.
-  F. Knoll, M. Holler, T. Koesters et al., “Joint MR-PET reconstruction using a multi-channel image regularizer,” IEEE Transactions on Medical Imaging, 2016.
-  J. E. Bowsher, H. Yuan, L. W. Hedlund et al., “Utilizing MRI information to estimate F18-FDG distributions in rat flank tumors,” in Nuclear Science Symposium Conference Record, 2004 IEEE, vol. 4. IEEE, 2004, pp. 2488–2492.
-  S. Chen, H. Liu, P. Shi et al., “Sparse representation and dictionary learning penalized image reconstruction for positron emission tomography,” Physics in medicine and biology, vol. 60, no. 2, p. 807, 2015.
-  M. S. Tahaei and A. J. Reader, “Patch-based image reconstruction for PET using prior-image derived dictionaries,” Physics in Medicine and Biology, vol. 61, no. 18, p. 6833, 2016.
-  J. Tang, B. Yang, Y. Wang et al., “Sparsity-constrained PET image reconstruction with learned dictionaries,” Physics in Medicine and Biology, vol. 61, no. 17, p. 6347, 2016.
-  G. Wang and J. Qi, “PET image reconstruction using kernel method,” IEEE Transactions on Medical Imaging, vol. 34, no. 1, pp. 61–71, 2015.
-  W. Hutchcroft, G. Wang, K. T. Chen et al., “Anatomically-aided PET reconstruction using the kernel method,” Physics in Medicine and Biology, vol. 61, no. 18, p. 6668, 2016.
-  P. Novosad and A. J. Reader, “MR-guided dynamic PET reconstruction with the kernel method and spectral temporal basis functions,” Physics in Medicine and Biology, vol. 61, no. 12, p. 4624, 2016.
-  S. Ellis and A. J. Reader, “Simultaneous maximum a posteriori longitudinal pet image reconstruction,” Physics in medicine and biology, vol. 62, no. 17, pp. 6963–6979, 2017.
-  O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2015, pp. 234–241.
-  S. Ren, K. He, R. Girshick et al., “Faster r-cnn: Towards real-time object detection with region proposal networks,” in Advances in neural information processing systems, 2015, pp. 91–99.
-  C. Dong, C. C. Loy, K. He et al., “Image super-resolution using deep convolutional networks,” IEEE transactions on pattern analysis and machine intelligence, vol. 38, no. 2, pp. 295–307, 2016.
S. Wang, Z. Su, L. Ying et al.
, “Accelerating magnetic resonance imaging via deep learning,” inBiomedical Imaging (ISBI), 2016 IEEE 13th International Symposium on. IEEE, 2016, pp. 514–517.
-  E. Kang, J. Min, and J. C. Ye, “A deep convolutional neural network using directional wavelets for low-dose x-ray ct reconstruction,” arXiv preprint arXiv:1610.09736, 2016.
-  H. Chen, Y. Zhang, W. Zhang et al., “Low-dose ct via convolutional neural network,” Biomedical optics express, vol. 8, no. 2, pp. 679–694, 2017.
-  D. Wu, K. Kim, G. E. Fakhri et al., “A cascaded convolutional nerual network for x-ray low-dose ct image denoising,” arXiv preprint arXiv:1705.04267, 2017.
-  J. M. Wolterink, T. Leiner, M. A. Viergever et al., “Generative adversarial networks for noise reduction in low-dose ct,” IEEE Transactions on Medical Imaging, 2017.
K. He, X. Zhang, S. Ren et al., “Deep residual learning for image
Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
S. Boyd, N. Parikh, E. Chu et al., “Distributed optimization and
statistical learning via the alternating direction method of multipliers,”
Foundations and Trends® in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
-  J. Qi, R. M. Leahy, S. R. Cherry et al., “High-resolution 3D Bayesian image reconstruction using the microPET small-animal scanner,” Physics in medicine and biology, vol. 43, no. 4, p. 1001, 1998.
-  G. Wang and J. Qi, “Penalized likelihood pet image reconstruction using patch-based edge-preserving regularization,” IEEE transactions on medical imaging, vol. 31, no. 12, pp. 2194–2204, 2012.
-  Y. Nesterov, “A method of solving a convex programming problem with convergence rate o (1/k2),” in Soviet Mathematics Doklady, vol. 27, no. 2, 1983, pp. 372–376.
-  D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
-  V. Bettinardi, L. Presotto, E. Rapisarda et al., “Physical performance of the new hybrid pet/ct discovery-690,” Medical physics, vol. 38, no. 10, pp. 5394–5411, 2011.
-  W. Segars, G. Sturgeon, S. Mendonca et al., “4d xcat phantom for multimodality imaging research,” Medical physics, vol. 37, no. 9, pp. 4902–4915, 2010.
-  D. Feng, K.-P. Wong, C.-M. Wu et al., “A technique for extracting physiological parameters and the required input function simultaneously from PET image measurements: Theory and simulation study,” IEEE Transactions on Information Technology in Biomedicine, vol. 1, no. 4, pp. 243–254, 1997.
-  H. Qiao and J. Bai, “Dynamic simulation of fdg-pet image based on vhp datasets,” in Complex Medical Engineering (CME), 2011 IEEE/ICME International Conference on. IEEE, 2011, pp. 154–158.
-  N. A. Karakatsanis, M. A. Lodge, A. K. Tahari et al., “Dynamic whole-body pet parametric imaging: I. concept, acquisition protocol optimization and clinical application,” Physics in medicine and biology, vol. 58, no. 20, p. 7391, 2013.
-  J. Zhou and J. Qi, “Fast and efficient fully 3D PET image reconstruction using sparse system matrix factorization with GPU acceleration,” Physics in medicine and biology, vol. 56, no. 20, p. 6739, 2011.