Low Dose PET Image Reconstruction with Total Variation Using Alternating Direction Method

In this paper, a total variation (TV) minimization strategy is proposed to overcome the problem of sparse spatial resolution and large amounts of noise in low dose positron emission tomography (PET) imaging reconstruction. Two types of objective function were established based on two statistical models of measured PET data, least-square (LS) TV for the Gaussian distribution and Poisson-TV for the Poisson distribution. To efficiently obtain high quality reconstructed images, the alternating direction method (ADM) is used to solve these objective functions. As compared with the iterative shrinkage/thresholding (IST) based algorithms, the proposed ADM can make full use of the TV constraint and its convergence rate is faster. The performance of the proposed approach is validated through comparisons with the expectation-maximization (EM) method using synthetic and experimental biological data. In the comparisons, the results of both LS-TV and Poisson-TV are taken into consideration to find which models are more suitable for PET imaging, in particular low-dose PET. To evaluate the results quantitatively, we computed bias, variance, and the contrast recovery coefficient (CRC) and drew profiles of the reconstructed images produced by the different methods. The results show that both Poisson-TV and LS-TV can provide a high visual quality at a low dose level. The bias and variance of the proposed LS-TV and Poisson-TV methods are 20% to 74% less at all counting levels than those of the EM method. Poisson-TV gives the best performance in terms of high-accuracy reconstruction with the lowest bias and variance as compared to the ground truth (14.3% less bias and 21.9% less variance). In contrast, LS-TV gives the best performance in terms of the high contrast of the reconstruction with the highest CRC.


Introduction
Positron emission tomography (PET) is a nuclear image modality that can produce 3D functional images of biological processes inside the human body [1]. PET has now become an indispensable tool in cardiac/brain research and cancer diagnosis/treatment. However, the reconstruction of low-dose PET images has remained a challenge because of the large amount PLOS ONE | DOI: 10 of noise and the sparse spatial resolution [2]. To meet the challenge, researchers have proposed a number of different iterative statistical methods to reconstruct the PET images based on different statistical assumptions of PET measurements, such as maximum likelihood (ML), expectation maximization (EM) [3][4][5], maximum a posteriori (MAP) [6][7][8], and penalized weighted least-squares (PWLS) [9,10]. Nevertheless, these iterative statistical methods still have several drawbacks. For example they cannot easily handle low signal-to-noise ratio (SNR) data. Furthermore, after the subtraction of random events, the corrections for scanner sensitivity and dead time, and the corrections for attenuation and scatter, the actual statistical property of measured data is quite complex and it does not exactly follow the Poisson distribution [11,12]. The nature of PET images is an additional factor to consider in PET image reconstruction [13]. One important PET image feature is the edge. The edge information usually represents the sharp variation in images, such as the object boundaries, which is very important and useful for clinical diagnosis, e.g., of tumors [14]. Given the nature of PET images, researchers have proposed total variation (TV) based methods for both the image space and the projection space in PET image reconstruction [15]. TV was incorporated to provide edge-preserving guidance for the reconstruction [15], and it is well known that it suppresses noise effectively while preserving sharp edges [16][17][18]. Based on the complexity of the statistical property of PET data, there should be several types of TV-constraint PET imaging models for different statistical assumptions. However, there have been no discussions or research studies on this point in the PET imaging community. In contrast, two types of TV incorporated models have been developed in the image processing community: Gaussian-TV [19] and Poisson-TV [20]. The Gaussian-TV model includes an energy functional to control the type of smoothness of solutions. The Poisson-TV model includes a Poisson data fidelity term to meet the statistical assumption. Both can be used for PET image reconstruction theoretically. However, it remains unknown which one is the better choice for various PET application areas. For example, the low bias and variance of the reconstructed image is more important for kinetic estimation, whereas high contrast is required for tumor diagnosis. In addition, a complication of the TV model is the strong nonlinearity in the data fidelity, and therefore, problems or issues arise in the computation of minimizers [20].
There are three classes of existing solvers for the TV problem. Algorithms in the first class are based on smoothing the TV term, since TV is non-smooth, which causes the main difficulty in solving the TV problem. A number of methods based on smoothing have been proposed, one of which is the second-order cone program (SOCP) [21,22]. The SOCP solver reformulates TV minimization as a second order cone programming problem [22], which can be solved by interior-point algorithms. The SOCP solver can easily be adapted to various convex TV models with distinct terms and constraints with high accuracy. However, the speed of SOCP is very slow, because it embeds the interior-point algorithm and directly solves a linear system at each iteration. The second class of algorithms for TV problems comprises those based on the iterative shrinkage/thresholding (IST) algorithms, which have been proposed by several researchers in the last few years [23][24][25][26]. IST is able to minimize the TV constraint term with some non-quadratic and non-smooth regularization terms. The convergence rate of IST algorithms heavily relies on the linear observation operator, and the convergence rate of algorithms in this class is not sufficiently fast. Furthermore, IST-based algorithms for TV deconvolution require that a de-noising subproblem be solved at each iteration and they cannot take advantage of problem structures. The third class of algorithms for TV problems comprises those based on seeking the minimizer or maximizer of the original constrained problem by a sequence of unconstrained subproblems. Methods belonging to this class add a quadratic penalty term instead of the normal constraint term in the objective function. The penalty term is the square of the constraint violation with the multiplier. Because of its simplicity and intuitive appeal, this approach is widely used. The well-known augmented Lagrangian method (ALM) [27,28] belongs to this class. In the augmented Lagrangian method, the Lagrangian multiplier is introduced and is estimated at each iteration in the objective function. However, it requires that multipliers go to infinity to guarantee the convergence, which may cause the ill-condition problem numerically.
To overcome the issue that the solvers mentioned above are weak in terms of efficiency and robustness, in this study an alternating direction method (ADM) was applied [29,30] to solve the TV problem. The proposed ADM constitutes implementable variants of the classical ALM for optimization problems with separable structures and linear constraints. In this method, the TV regularization term is split into two terms with the aid of a new slack variable so that an alternating minimization scheme can be coupled to minimize the approximate objective function. This split and the use of the alternating minimization scheme not only accelerate the convergence rate of the solution, but also result in improved accuracy, as well as in robustness of the reconstructed results to the noise in the data set.
We used Monte Carlo simulated data, phantom data, and real patient data to validate the performance of the proposed algorithm as perceived both quantitatively and visually. Experimental results show that the proposed algorithm is highly effective in preserving sharp image edges and more details while eliminating staircase artifacts. In addition, the performances of the Poisson-TV and LS-TV methods on PET data at different counting levels are also evaluated in this work.
The rest of this paper is organized as follows. In Section 2, first the PET imaging model is reviewed, and then, two objective functions and the corresponding method to solve them are suggested. In Section 3, the experimental setup, results, and comparisons of the existing methods and the proposed methods are presented. Section 4 constitutes the conclusion.

PET imaging model
PET acquired data are organized in a series of parallel slices that can be reconstructed independently. Every slice of raw data collected by a PET scanner constitutes a list of coincidence events representing near-simultaneous detection of annihilation photons by a pair of detectors. Each coincidence event represents a line in space connecting the two detectors along which the positron emission occurred (the line of response (LOR)). The raw data from PET are organized in a sinogram.
Therefore, PET image reconstruction problems are specific cases of the following general inverse problem: find an estimate of radioactive activity map u from a measurement b by In the process of PET imaging reconstruction, u is the reconstruction image and A is the system matrix that describes the tomographic geometry and the physical factors.

LS-TV
Problem formulation. Assuming piecewise constant behavior of PET images, we introduce total variation (TV) regularization into PET reconstruction.
The problem is formulated as min u TVðuÞ; s:t: where TV(u) is the PET reconstruction with the TV defined regularization; it is defined as ∑ i kD i uk, the sum of the discrete gradient of activity map u of every pixel i. Solution. It is difficult to directly obtain the solution of Eq (2). Therefore, we introduce a new auxiliary variable, w. At each pixel, an auxiliary variable w i is introduced into the term k.k. The purpose of this process is to transfer D i u out of k.k. The optimization problem of Eq (2) is clearly equivalent to To deal with the constraints, we transform the constrained Problem (3) to an equivalent unconstrained problem using an augmented Lagrangian function [31]. The corresponding augmented Lagrangian function of Problem (3) is where υ i , β i , λ, and μ are the multipliers of the four penalty terms. The first term of Eq (4) is the regularization term, and the remaining terms are penalty terms. The second and fourth terms are linear parts, whereas the third term and fifth terms are quadratic parts. These parts ensure the accuracy and robustness of the TV constraint. In order to make the result of every term a number rather than a matrix, transposition of υ and λ is utilized in Eq (4). To solve Eq (4) efficiently, the ADM, which was originally proposed to handle parabolic and elliptic differential equations, was embedded here. This algorithm is a variant of the classical ALM. When the classical ALM approaches the solutions of the original Problem (4), the nice separable structure emerging from Eq (4) in both the objective function and the constraint is weaker. This drawback, however, can be completely overcome by the ADM. In the ADM, the solution of Eq (4) is transformed to solve three subproblems at each iteration for variables u k and w i,k and parameters λ and υ.
In ADM, let u k and w i,k represent the true minimizers of Eq (4) at the kth iteration. w i,k+1 can be attained by To solve this minimization problem, Eq (5) can be separated as two subproblems for variables u and w based on ADM. We first fix the value of u and calculate the solution of w. Therefore, only the terms containing w are useful. The corresponding subproblem can be expressed as the following problem: For given β > 0, the minimizer of Eq (6) is given by the 2D shrinkage-like formula [32]. w i,k+1 can be calculated by With w i,k+1 , we can achieve u k+1 by The constant terms do not influence the minimum, and thus, this subproblem is equivalent to the problem Its gradient is By enforcing d k (u) = 0, we can obtain the exact minimizer of Eq (9) directly. However, this calculation is too costly to implement numerically, in particular when the matrix is large. To obtain u k+1 more efficiently, the steepest descent method with an appropriate step length is used iteratively by applying the recurrence formulation: where α k is the step length. Each iteration of the steepest descent method demands that the gradient be updated to update the estimation value of u k+1 . Therefore, the step length should be chosen carefully to obtain an accurate solution. It remains to choose α. It is suggested that the aggressive manner [33][34][35] be used to choose the step length for the steepest descent method, which is called the BB (Barzilai and Borwein) step or BB method. This method is applied to choose α: Parameters. There are several parameters in the algorithm. Among these, β, υ i , and λ are initialized as 1, 1, 1, based on [36]. μ is the most important parameter, since it determines the weight of the data fitting term. Therefore, to achieve the best performance, the value of μ should be set according to the noise level in observation b. For example, the higher the noise level is, the smaller μ should be. μ was manually tuned in this study. υ i and λ should be updated provided that Eq (4) is minimized at each iteration. According to the formula proposed by Hestenes and Powell [28,[37][38][39][40], the update formulas of multipliers follow The program terminates after a certain number of iterations (300 in this study), or when the relative stopping criterion (based on empirical estimations) is reached: Poisson-TV Problem formulation. Because of the issue of the low SNR and the Poisson distribution of PET measurements, we build the objective function as min u TVðuÞ þ m Solution. To solve the objective function in Eq (16), the ADM is also used to solve some subproblems at each iteration to approach a solution of Eq (16). We add two auxiliary variables, Z S and S, into the function, as in Eq (3). Then, we construct the augmented Lagrangian objective function for Eq (16): In a more compact form, we have We can obtain u by where the u sub-problem is solved by an EM-like method. We rewrite the objective function with items containing only u and take the expectation step with respect to the unobservable variables w i j and calculate the surrogate function.
We minimize this surrogate function F by zeroing its derivative with respect to S jm : Thus, x kþ1 jm is the root of Finally, the updating rule for x kþ1 jm is Moreover, Z S in Eq (19) is updated as Convergence analysis The analysis of the convergence properties of the proposed TV regulation process is presented in this section. β and β S strongly affect convergence in both the Poisson-TV and LS-TV methods. The selection of these parameters also affects the convergence rate of the proposed method. The convergence of the sequence (w and u) in the Poisson-TV and (u and S) in the LS-TV was proven in the convergence analysis section in [32,41]. It has been proven that the proposed method can converge to the final solutions from any initial points. To demonstrate the effectiveness of the ADM, we also provide a comparison of the convergence speed of the proposed method and IST method (a detailed description of this is algorithm is given in [23]) in the results section.

Ethics Statement
Before the Results section, we state that no human and animal research was involved in this work.

Results
In this section, to validate and evaluate the proposed method, the results of simulation and clinical experiments are presented. In the Monte Carlo simulation experiments, a Hoffman brain phantom and Zubal phantom were used and we simulated projections using GATE [42]. All the tests were performed using MATLAB on a PC with an Intel i7-3770 3.40 GHz and 8 GB RAM. In the clinical experiment, we applied our method to a typical PET scanner. All the scans were executed using a Hamamatsu PET scanner. To obtain a good visual effect of our method, we compared the EM [43], Poisson-TV, and LS-TV methods.

Monte Carlo simulations
The scanner simulated in GATE was the Hamamatsu SHR-22000 scanner, which consists of 32 rings with 24576 (3072/ring) BGO crystals of 2.8mm Ã 6.95mm Ã 30mm, 768 PMTs, and an 838 mm detector ring diameter. The activity maps of the Hoffman brain phantom (Fig 1(a)) and Zubal phantom (Fig 1(b)) were used as the ground truths. The sinograms calculated from the ground truths by GATE were used as the measurements for our tests. All the tests were performed using MATLAB on a PC with an Intel i7-3770 3.40 GHz and 8 GB RAM. The system matrix G used for matching the dimensions of the simulated sinogram was generated by the Fessler tomography toolbox [44].
Both the EM and our model-based TV method were applied to reconstruct the activity maps from the Monte Carlo simulation measurements. The ability of two frameworks (the EM and model-based TV method) to reconstruct a PET activity map was compared in this experiment. This ability is quantified using the relative errors bias, variance, and contrast recovery coefficient (CRC) of the region of interest (ROI) calculated as where u i is the ith pixel of ground truth u andû i is the ith pixel of the reconstructed images. S is the mean activity of the ROI and B is the mean activity of the white matter region (background) in the reconstructed image. The bias and variance are used to evaluate the accuracy of the reconstruction and CRC is used to indicate the relative contrast of the ROI in the reconstructed images. In our method, the TV factor (the value of weighting parameter μ) is an important parameter, which can be optimized to smoothen the results. To study the impact of the initial condition of the TV factor in our method, we plotted the bias curve and variance curve versus the TV factor, which are shown in Fig 2. When the TV factor increases, the bias and variance curves first decrease and then increase. Therefore, we can choose the nadir value of the curve as the optimized TV factor, which is 0.0025 in Poisson-TV and 4 in LS-TV for all experiments at all count levels.
The activity maps of Hoffman brain phantom and Zubal phantom reconstructed by the EM and the proposed TV methods are presented in Figs 3 and 4. The sinograms used in Figs 3 and 4 constitute high count level data; the count is approximately 10 7 . This clearly indicates that the quality of recovery of both Poisson-TV and LS-TV is much higher than that of the EM method. In addition, the reconstructed images also indicate that the edge of the reconstruction produced by the proposed methods is sharper. However, the LS-TV method yields sharper edge information than does Poisson-TV. To evaluate the accuracy of the reconstruction, we marked the specific region of the images (Figs 3(b) and 4(b)) for quantitative analysis. In addition, we present the profiles of the reconstructed images by all methods in the Fig 5. It is clear that the results of Poisson-TV give the closest fit to the ground truth.    To evaluate the performance of the Poisson-TV and LS-TV methods at different count (dose) levels, five counting level values (total number of photon counts in the reconstruction plane: 5 Ã 10 5 , 1 Ã 10 6 , 3 Ã 10 6 , 6 Ã 10 6 and 9 Ã 10 6 ) were simulated in this experiment.
The activity maps of the Hoffman brain phantom (Fig 6(a)) and Zubal phantom (Fig 6(b)) reconstructed by EM, Poisson-TV, and LS-TV from the data with different numbers of counts are given in Fig 6. Tables 1 and 2 list the bias and variance of EM, Poisson-TV, and LS-TV for the different counting level data. This demonstrates the statistical analysis of the reconstructed results, which confirms that our model yields better estimates in terms of bias when the counting level changes. The robustness of the proposed method to noise was also evaluated by sinogram data with different counting levels (a lower count means more serious noise). As indicated in Tables 1 and 2, the variance of the Poisson-TV and LS-TV methods is 22% to 74% less than that of the EM method at all counting levels. In other words, the proposed method is more robust to noise than EM. Moreover, Table 1 also shows that the bias and variance for LS-TV and Poisson-TV decrease faster than those of EM when the counting level increases. For a full comparison, the bias, variance, and CRC curves of EM, Poisson-TV, and LS-TV in ROI 1, 2, and 3 for the brain and Zubal phantoms are given in Fig 7. Fig 7(a) and 7(b) show the bias curves. The bias curve for the EM was not improved significantly when the counting level increased, whereas the bias curve for our method decreased . Fig 7(c) and 7(d) show the variance curves, which exhibit the same trend as Fig 7(a) and 7(b). Fig 7(e) and 7(f) show the CRC curves. All the CRC curves increase when the counting level increases, whereas the CRC curves for the proposed method increase faster and higher, and LS-TV gives the best CRC in all methods. Furthermore, we plotted the relative error of reconstruction after each iteration of the algorithms to evaluate the convergence speed of the Poisson-TV, LS-TV, EM, and IST methods. Here, the IST method was used to evaluate the effectiveness of the ADM to solve the TV-problem. The results are shown in Fig 8. The y-axis represents relative error and the x-axis represents the iteration times of the methods. The black curve represents Poisson-TV, the green curve LS-TV, the yellow curve EM, and the red curve the IST method. For all methods, the relative error approaches close to 0 within 20 iterations (since the stopping threshold is 10 −3 for all experiments). With the help of the ADM, the relative error of the Poisson-TV and LS-TV decreases to approximately 0.1 with only two iterations, while the relative error of IST decreases to approximately 0.22 and the relative error of EM is 0.65. This shows that the curves of both Poisson-TV and LS-TV decrease faster than those of IST and EM. All the results demonstrate that the convergence speed of the ADM is faster than that of the IST method and much faster than that of the EM method.  Fig 9. To evaluate the detailed information of the reconstruction, the marked areas (marked by a red rectangle) are given in Fig 9(b) and zoomed for comparison. It is clearly indicated that the quality of the recovery, in particular the detailed information yielded by the Poisson-TV and LS-TV methods, is better than that of the EM method.

Clinical case
In this section, we validate the proposed method using clinical patient data. The clinical patient data comprised a PET scan acquired from a volunteer at a local hospital. The PET system used was a Hamamatsu SHR-22000 whole body PET scanner. It has 32 crystal rings and can be operated in 2D or 3D mode. The trans-axial resolution of the central FOV is about 3.7 mm. The patient lay on the bed and the main target of the scanning was the thorax. The total scanning time was 40 min. Random, normalization, dead time, scatter, and attenuation correction were applied to the measurement data using the programs provided by the scanner prior to reconstruction. The measurement data were stored in the sinogram model, the size of the reconstructed images was 128 Ã 128.
To evaluate the effectiveness of the proposed method, the Poisson-TV, LS-TV, and EM algorithms were used in this experiment. Their results are shown in Fig 10. In Fig 10(a), the difference between these three methods is not sufficiently clear, however, all the algorithms can provide a clear reconstruction. However, when we zoom in on the highlighted region (marked by a red rectangle) in the reconstructed images, it is clear that, as compared to EM, both the Poisson and LS TV algorithm can provide better edge information and better visual quality, in particular for the small highlighted point marked by the red arrow. In the results of the EM algorithm, this point almost disappears in the noise background. In contrast, we can easily locate this point in the results of Poisson-TV and LS-TV. In the comparison of the LS-TV and Poisson-TV, since the results of the LS-TV method suffer from some undesirable artifacts (such as edge over-smoothness and staircase effects), Poisson-TV yields a much better visual quality than LS-TV and recovers more detail of textures in the reconstructed image.

Conclusion
This paper presented a TV-constraint reconstruction algorithm for low dose PET image reconstruction. According to the property of the PET measured data, two types of TV-constraint models, LS-TV and Poisson-TV, were established. To efficiently obtain high quality reconstructed images, the ADM is used. To evaluate the effectiveness of the proposed LS-TV and Poisson-TV methods, PET data measured from a Monte Carlo simulation and real scanning were used. As compared with the traditional EM method, the proposed TV-constraint models are able to reconstruct the images with higher accuracy and much clearer structures. The ADM also provides a faster convergence rate as compared to the IST method for solving the TV problem. Both Poisson-TV and LS-TV can provide a high visual quality and are robust to noise corruptions at a low dose level; Poisson-TV gives a more accurate estimation with lower bias and variances, whereas the CRC of the LS-TV is better. This means that, although both the proposed TV methods can provide a highly accurate and robust reconstruction, Poisson-TV is more suitable for kinetic parameter estimation, which requires a highly accurate estimation with low bias and variance for PET images. In contrast, LS-TV is more suitable for tumor diagnosis, which requires a high contrast for PET images in order to locate the tumor.