Reconstruction for 3D PET Based on Total Variation Constrained Direct Fourier Method

This paper presents a total variation (TV) regularized reconstruction algorithm for 3D positron emission tomography (PET). The proposed method first employs the Fourier rebinning algorithm (FORE), rebinning the 3D data into a stack of ordinary 2D data sets as sinogram data. Then, the resulted 2D sinogram are ready to be reconstructed by conventional 2D reconstruction algorithms. Given the locally piece-wise constant nature of PET images, we introduce the total variation (TV) based reconstruction schemes. More specifically, we formulate the 2D PET reconstruction problem as an optimization problem, whose objective function consists of TV norm of the reconstructed image and the data fidelity term measuring the consistency between the reconstructed image and sinogram. To solve the resulting minimization problem, we apply an efficient methods called the Bregman operator splitting algorithm with variable step size (BOSVS). Experiments based on Monte Carlo simulated data and real data are conducted as validations. The experiment results show that the proposed method produces higher accuracy than conventional direct Fourier (DF) (bias in BOSVS is 70% of ones in DF, variance of BOSVS is 80% of ones in DF).


Introduction
Positron emission tomography (PET) is a nuclear medicine imaging method based on coincidence detection of photon pairs emitted from positron annihilation events [1]. Apart from operating in two-dimensional (2D) data acquisition, PET are often used in three-dimensional (3D) mode, where the interplane septa are removed. The 3D data, in contrast with the 2D data acquired, approximate line integrals of the radioactive tracer distribution along all possible lines of response (LOR's) which are not restricted to transaxial planes [2][3][4]. The 2D acquisition to 3D acquisition leads to a significant improvement of the scanner sensitivity, due to the increased number of measured line of response [5][6][7].
In 3D PET, the amount of data collected by scanner is extremely large, therefore, the focus in reconstruction has been largely on the reduction of the computation cost. To date, significant improvements have been achieved by using the two following approaches. There have been considerable efforts aimed at speeding iterative reconstruction methods with a combination of computer acceleration techniques (e.g. graphics processing unit) [8,9], and intelligent coding (e.g. precomputing factors or implementing [3,10] the calculations of the system matrix in the reconstruction process [11][12][13][14]). A viable alternative to the above method is a class of approaches known as direct Fourier (DF) strategy, which has long been studied by many researchers [15]. The primary emphasis has been on rebinning the 3D data in plane integrals or in a 2D data set [16][17][18], and the 3D projection algorithm (3DRP) [19], and the direct Fourier reconstruction with Fourier reprojection (3D-FRP) [20], and the Fourier-wavelet restoration technique [21].
While direct Fourier method is fast and easy to implement [22][23][24][25], the reconstruction accuracy suffers from performance limitations due to the fact that the Fourier basis cannot adequately represent spatially inhomogeneous data, like that typically found in the images [26,27]. For example, tumors or cancerous tissue, because they absorb most of the fluorine radioisotope, the intensity values of these regions are larger than surrounding materials [28,29]. In the case, the tumors in the images that are reconstructed by direct Fourier method tends to be blurred and illegible [30,31]. On the other hand, for the organisms imaged by PET, including tumors, their intensity values are homogeneous within that region, i.e., locally piece-wise constant [32][33][34][35]. This induces us to incorporate the TV regularization into DF framework to improve the reconstruction accuracy. Total variation was first proposed by Rudin, Osher and Fatemi (ROF) [36]. Unlike other regularization methods [37][38][39], it avoids over-smoothing the inhomogeneous region by minimizing the gradient field of the image. To be precise, TV preserves the locally piece smooth region but also encourages sharp variations near the edges [40][41][42][43][44].
Thus, the problem of reconstruction can be formulated to be an optimization problem, whose objective function consists of TV norm of the reconstructed image and the data fidelity term measuring the consistency between the reconstructed image and sinogram. This paper applies BOSVS [45] which called the Bregman operator splitting algorithm with variable step size to solve the optimization problem [46,47]. The algorithm uses the variable Barzilai-Borwein (BB) step instead of the fixed BB step used in original Bregman Operator Splitting (BOS) algorithm, the stepsize rule starts with a inial value, and increases the nominal step until termination condition is satisfied. By combining a variable splitting and the alternating direction method of multipliers (ADMM) with a Barzilai-Borwein approximation to the Hessian, the convergence of reconstruction turns out to be faster.
We evaluate the quality of the proposed method by using Monte Carlo simulated data and real patient data. After rebinning the 3D data, we compare performance in terms of contrast recovery, noise performance, performance on detecting a small target or regions, initialization issues and robustness.

Methods
The whole framework includes two parts: 3D data rebinning and 2D image reconstruction. This section is organized as follows: First of all, the detail about 3D data in PET scanner is given, and then we adopt the Fourier rebinning algorithm (FORE) to sort the 3D data into a stack of ordinary 2D data sets as sinogram data. Then, we introduce the total variation (TV) constrained direct Fourier reconstruction schemes. At last, we apply BOSVS to solve the resulting problem.

3D PET Imaging Model
3D Data for PET Scanner. Our study uses the same geometry and data parametrization as those in [5] for the approximate Fourier rebinning algorithm FORE. As shown in Fig 1, the scanner is a truncated cylinder with radius R(R = r f ) and length L, with its axis along the z axis. Line integrals of the tracer distribution f(x, y, z) are measured for all LOR's connecting detector A and detector B on the surface of the cylinder. After normalization by a factor 1= these data are given by where we have simplified notations by replacing the limits AE ffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 f À s 2 q of the integral over l by ±1, using the fact that f = 0 outside the cylinder of radius r f . The variable z is the axial coordinate of the point mid-way between detector A and detector B and the axial angular variable is δ = tanθ where θ is the angle between the LOR's and the transaxial plane [5,48].
Rebinning 3D Data into 2D. The principle of a rebinning algorithm is summarized in Fig  2. For a fixed pair (z, δ), the set of data is referred to as an oblique sinogram. Each oblique sinogram is parameterized with the standard variables s and ϕ: s is the signed distance between the z axis and the LOR's, ϕ is the angle between the x axis and the projection of the LOR's onto a transaxial plane. The ranges of these variables are −r f s r f and 0 ϕ 2π.
When a scanner is operated in 2D mode, the measured LOR's are confined to lie in a transaxial plane, so that a 2D data set can be also described by Eq (1) but with δ = 0 pðs; ; zÞ 2D ¼ pðs; ; z; 0Þ ð2Þ These 2D data contain, for each slice z, one ordinary (or direct) sinogram which can be independently reconstructed by 2D reconstruction methods. Not surprisingly, this slice by slice reconstruction of a three-parameter data set is considerably faster than the reconstruction of the four-parameter 3D data set with the 3D direct reconstruction algorithm such as 3DRP.
This observation suggests an alternative approach to 3D reconstruction, in which the 3D data are not directly reconstructed, but serve to estimate 2D data from which the image can then be recovered using any 2D reconstruction algorithm such as DF. Therefore, we use the FORE(Fourier rebinning algorithm) as a method to estimate p 2D (s, ϕ, z) from p(s, ϕ, z, 0) [5]. 2D PET Imaging Model DF methods follow directly from the ideal procedure suggested by the central slice theorem applied to the discrete case. And in a real case, only a finite number of projection are taken, giving a finite number of F(u, v) samples located on a uniformly distributed polar grid for the central slice theorem.
Thus, in order to obtain an approximation of f(x, y) by 2D inverse Fourier transform of F(u, v), we have first to interpolate in Fourier domain from the radial points to the points on a Cartesian grid. Basically, direct Fourier reconstruction methods are three steps methods: First, 1D discrete Fourier transform (through FFT algorithm) of the parallel projections taken at angles θ 1 , θ 2 , Á Á Á, θ M . Then, after the interpolation from Polar to Cartesian grid, the process is completed by applying 2D inverse Fourier transform (using FFT).
Unfortunately, the result of such a straight method suffer from artifacts due to interpolation in Fourier plane and to aliasing (since function f(x, y) is not band limited). Nevertheless, due to its lower computational complexity O(N 2 logN), compared to that of the widespread Filtered Back Projections O(N 3 ), the method has been object of research and various techniques has been proposed in order to improve its performance [15]. Likewise, our proposed model will be based on the three steps and a generalized flowchart can be found in

Problem Formulation
Total variation regularization is a recently image processing techniques, which has been shown to be very successful in many image processing applications such as PET and MRI [49]. And recent research has confirmed that for image restoration the use of total variation (TV) regularization instead of the l 1 term makes the recovered image quality sharper by preserving the edges or boundaries more accurately, which is essential to characterize images. The advantages of TV minimization stem from the property that it can recover not only sparse signals or images, but also dense staircase signals or piecewise constant images.
As was commonly known, the fundamental equation for Fourier transform can be shown as: where F represents the Fourier transform, f is the spatial frequency information matrix, u is the reconstructed image. Thus the problem can be formulated as: Here α is the weighting factor of TV, and following the standard treatment we will vectorize an two-dimensional image u into one-dimensional column vector, i.e. u 2 C N , where N is the number of pixels of the image u. Then, the (isotropic) TV norm is defined by where for each i = 1, ‥, N, D i 2 R 2 × N has two nonzero entries in each row corresponding to finite difference approximations to partial derivatives of u at the i-th pixel along the coordinate axes.

Solution
We introduce BOSVS to numerically solve the total variation based image reconstruction problem (4), which uses variable splitting method to reduce computational cost and the Barzilai-Borwein step size selection method is adopted for much faster convergence. [45] Firstly, we introduce auxiliary variables w i to transform D i out of the non-differentiable norms which is clearly equivalent to the original problem (4) as they share the same solutions u. Then, we define G(u) as GðuÞ ¼ 1 2a kFu À f k 2 2 , the augmented Lagrangian can be defined as: where λ 2 C 2N is a Lagrange multiplier, and β is the penalty factor. The method of multiplier iterates the minimizations of Lagrangian L in Eq (7) with respect to (w, u) and the updates of the multipliers λ Thus, the problem divides into several subproblems, which can be efficiently solved.
For w subproblem, it can be get from Problem (11) has a closed form solution, 1-D Shrinkage can be applied, which for given vector b and scalar μ > 0, has solution shrinkage operator For u subproblem, Using Bregman operator splitting to G(u) [45]: where δ is a constant step size. Then we drop constants in updating u k+1 and complete square, as a result u can be written as below: Take the Barzilai-Borwein step size method in the selection of step size: Using Hessian information in the selection of the step size δ by letting 1 The step size δ isn't a constant any more, which can reduce computational cost and accelerate the convergence. The BB step size can be written as below: It converges if d 1 L G (L G : Lipschitz constant for rG) as proved in [47]. In our experiment, the initial value of δ is chosen to be 0.8 and the iteration is stopped when δ < 10 −5 .
Finally, u can be calculated as below: The computation time for a single slice with a size of 192 by 192 pixels and iterating 100 times turns out to be 1.297345 seconds.

Results
In order to demonstrate the performance of our proposed method, we design two groups of experiments using Monte Carlo simulated data and real data respectively. All the codes in the following experiments are implemented in Matlab R2011a (MathWorks Corporation, USA) and run in a desktop computer with i5 Intel Core CPU and 6 GB memory.

Monte Carlo Simulation
PET Simulation Setup. Monte Carlo simulations can provide a relatively accurate reference for the development and assessment of new image reconstruction algorithms. Simulations in our study are performed using toolbox GATE. The simulated PET scanner is Hamamatsu SHR22000 and the phantoms chosen here for simulations are Zubal thorax phantom [50] and Hoffman brain phantom. The simulation outputs are stored in sinogram. Before reconstruction, random correction, normalization correction, attenuation correction and scatter correction are performed properly. The sinogram of brain phantom are reconstructed as images of size 64 Ã 64 and sinogram of zubal phantom are reconstructed as images of size 128 Ã 128.
Reconstruction Results Comparison of Different Methods. Since these experiments are based on Monte Carlo simulations, we can get the true activity distributions at any time exactly. In order to analyze the reconstruction results quantitatively, we define the measurements as follows: variance where u i ,û i and u n represent the estimated activity value of pixel i, the true activity value of pixel i, and the average estimated value of all the n pixels in one ROI respectively. Furthermore, we compute the contrast recovery coefficient (CRC), which is calculated by where S is the mean activity of the region of interest and B is the mean activity of the white matter region (background) in the reconstructed image. Fig 4 gives the true activity distributions and reconstructed images of the center slice by DF and BOSVS. Table 1 calculates the bias and variance between true activity distributions and reconstruction results by DF and BOSVS respectively. Noise Toleration. We provide four Monte Carlo simulated sinogram for DF and BOSVS reconstruction in different kinds of counting rates and the results are illustrated as Fig 8. Obviously we can conclude that BOSVS results in a higher accuracy than DF. BOSVS can also have a good reconstructed result even in low counting rate.
The calculated bias and variance between reconstruction results and true activity distributions through all four kinds of counting levels are plotted in Fig 9, and some results are  Table 2 for better understanding. From the quantitative analysis in the experiment, the bias and variance fall with the counting rate increasing. Moreover, BOSVS has a better performance than DF whatever the counting rate is.
Comparison of Different ROIs. Different regions in phantom may have different physiological and physical property, and reconstruction of a small target in the presence of background activity can be a challenge in PET imaging, as small targets may either be visually obstructed by noise in a regularized reconstruction or be over-smoothed by regularization. In this section we compare two reconstruction methods aimed at different regions. Take ROI1 and ROI2 for example, Fig 10(A) gives the image of bias and variance versus counting rate in different ROIs of phantom, and more details are shown in Table 2. Fig 10(B) compares the contrast recovery coefficient (CRC) versus counting rate by different methods for two regions of  interest (ROI). The plots show that the proposed method achieves a better performance than DF in different regions, even sinogram is in low counting level. Effect of Parameters. From our proposed model, α is the weighting factor of TV, which can be optimized to smooth results. Fig 11(A) shows the bias and variance when TV factor α is changed. By increasing the TV weight factor, we can see that the bias and variance of reconstructed results decrease first and then after one optimal point, bias and variance will increase. As we know, TV can smooth the image, to get the balance of data fidelity and TV, we choose α = 1.5 for brain phantom and α = 1.1 for zubal phantom to get the best reconstructed results.
The other important parameter in the proposed model is the initialization of step size δ, Fig  11(B) gives the relative error between reconstruction results and true activity distributions versus iterations in BOSVS with different initialization of step size δ. We can see that the step size nearly has nothing to do with the convergence of the reconstruction algorithm. It just matters how fast the algorithm converge at the beginning, and at last the iterations will stop at the same time with all kinds of step size initialization. So in this experiment, we just choose the step size initialization δ as 0.8 both for brain phantom and Zubal phantom.

Real Data
Real data experiments are performed on CTI ECAT data, and 47 slices real patient data are obtained. The scanner used is CTI ECAT PET system, which contains three rings of 48 block detectors each, for an axial field of view of 16.2 cm. The raw sinogram, along with the attenuation and efficiency factors, are radial samples by 192 angular samples (over 180 degrees), and will be reconstructed to images of 128 Ã 128. The values of sinogram range from 30.76 to 98.28, and sum is 1.75756 × 10 6 . Sinogram have been pre-corrected for delayed coincidences, and four slices of reconstructed images are presented for visual judgment. Note that all these data are publicly available in Fessler's website. Patient records/information has been anonymized and de-identified prior to analysis.  We can obviously see that BOSVS perform far better than traditional direct Fourier. Fig 13 gives the details of our selected area in the image reconstructed by two methods, we mark out our interest area of all four slices with the red rectangles in the reconstruction result by BOSVS, then compare the details of the same area of reconstruction results by two methods, we can conclude that BOSVS shows better recovery of activity distributions than DF even in small regions.

Discussion
This paper introduces total variation to traditional PET direct Fourier reconstruction methods to preserve the edges or boundaries more accurately and keep images more smooth. The main advantage of the proposed method is its capability of generating a high-resolution reconstruction results. We discuss the performance of the proposed method with DF in term of benefit of total variation regularization, performance of a small target or regions, initialization issues and robustness in the following context.

Benefit of Total Variation Regularization
The experiments on Monte Carlo simulation data tested the reconstruction performance of traditional DF reconstruction method and the proposed method. Both in the brain phantom and zubal phantom, DF gives a relatively less bias and variance. In contrast, the proposed method exhibits high accuracy and stability, as indicated in Figs 4, 6, 7 and Table 1. Meanwhile, DF does not perform well when our data is in low counting rate, DF gives a poor bias and variance when the counting rate of data are far less than 1 × 10 6 , and DF even could not give an acceptable reconstruction results in low counting level. Whereas the bias and variance of our proposed method are far better as shown in the Table 2, total variation shows strong noise toleration and reconstruction performance.
The effectiveness of our method with total variation in estimating reconstruction performance is also extensively examined in the real data experiments. The reconstructed images indicate that the accuracy of results are greatly improved using our proposed method. Four slices real patient data, obtained from CTI ECAT PET system, with high level noises and kinds of physical factors, can be recovered with higher quality than DF reconstruction methods, as shown in Figs 12 and 13. It is clear that the results achieved through BOSVS have less background noise and present a much clearer contour than that of DF. However, it can be noticed that there are some rough color lumps which is very likely to appear when we add TV regularization into our reconstruction algorithm since TV based algorithms tend to over-smooth image details and textures. [51]

Performance of a Small Target or Regions
We compare the reconstructed images by the DF and proposed methods in Fig 7 and we apply CRC as our measurement for reconstruction of regions. The image reconstructed by DF has a CRC of 0.55 for brain phantom and 0.58 for zubal phantom. In comparison, the BOSVS preserves the contrast (CRC: brain 0.61 zubal 0.66). And BOSVS demonstrates a more robust performance of CRC versus counting rate than DF. These results indicate that the proposed method can work robustly for reconstruction for both large and small targets.

Initialization Issues
The number of iterations, the measurement, and the initialization of the external noise and all parameters are all important in our proposed methods. The measurements represent information about the examined object, but our method does not explicitly require the type of measurements. On the other hand, in order to acquire a strain image of high quality, we simply perform 100 iterations for the scheme. The experiments on Monte Carlo simulation data and real patient data indicate that a convenient guideline for the initialization of parameter in proposed method can be easily concluded as follows: • Since real patient data are originally noisy and we make kinds of experiments in different counting rates, the external noise can be initialized to zero.
• Most of the parameters of BOSVS such as u, w, λ always change with iterations, and BOSVS is almost an automated algorithm. The initialization of these parameters will not affect the reconstruction performance. So in our proposed methods, we initialize these parameters as 0, and the values of parameters will be adjusted by the algorithm with iterations. • TV weighting factor α and the initialization of step size δ has been discussed in Fig 11. The initialization should be adjusted according to each case, and should be carefully set by the designer to obtain the desired estimation performance.

Robustness of the Proposed Method
The proposed method can better deal with internal and external disturbances/noise, and only require a TV weighting factor and a initialization of step size for the iterative recovery process.
It makes no assumptions about external noise statistics, which makes it more appropriate for certain practical problems where the disturbances/uncertainties are unknown. The robustness of the proposed method to noise was evaluated by employing different counting rate simulated sinogram data (For brain phantom, bias and variance of BOSVS are almost 60% of ones of DF), and we also examine the sensitivity of step size initialization for BOSVS(No matter what the initial value is, the iteration will stop at no more than 50th loop).

Conclusion
This paper presents a total variation constrained direct Fourier reconstruction algorithm for 3D PET. By introducing TV regularization to traditional PET imaging model, we can enable smoother results and better denoising performance. Monte Carlo simulation and real patient data show the proposed method can implement the reconstruction in higher accuracy compared with traditional direct Fourier method, and it requires few iterations and little computational cost.