Robust Framework for PET Image Reconstruction Incorporating System and Measurement Uncertainties

In Positron Emission Tomography (PET), an optimal estimate of the radioactivity concentration is obtained from the measured emission data under certain criteria. So far, all the well-known statistical reconstruction algorithms require exactly known system probability matrix a priori, and the quality of such system model largely determines the quality of the reconstructed images. In this paper, we propose an algorithm for PET image reconstruction for the real world case where the PET system model is subject to uncertainties. The method counts PET reconstruction as a regularization problem and the image estimation is achieved by means of an uncertainty weighted least squares framework. The performance of our work is evaluated with the Shepp-Logan simulated and real phantom data, which demonstrates significant improvements in image quality over the least squares reconstruction efforts.


Introduction
Positron Emission Tomography (PET) is one of the most important medical imaging modality which provides in vivo functional information of biological organs. It utilizes the idea of injecting chemical compounds tagged with positron emitting isotopes into a body to acquire complete coincidence data, which records the concentration information of the isotope distributions at specific locations within the body. The radioactivity images are then reconstructed based on these photon counting measurements.
The reconstruction of isotope concentration distribution is an ill-posed inverse problem. Most of current approaches to tackle this problem can be classified into two general categories, namely the analytical methods, which rely on the inversion of Radon transform, and the iterative approaches, which are based on a statistical description of the physical problem [1][2][3][4][5][6]. Because of the random nature of the radioactive disintegration, the tomographic data are noisy, and therefore it is straight-forward to regard PET reconstruction as a statistical estimation problem. Such approaches, when reconstructing PET images, need to introduce modeling of the data statistics and to make use of some prior information about the PET imaging system (often referred to as system probability matrix, which represents the probability that an emission event will be detected). For instance, Poisson/Gaussian assumptions on photon counting measurement data may be employed to deal with measurement uncertainties, thus constraining the solution space of reconstruction problem in maximum likelihood/least squares based frameworks [7][8][9][10][11].
However, so far, the common feature of all statistical-based methods for PET image reconstruction is that the system response model is assumed to be exactly known a priori. In real situations, however, it is almost impossible to have the exact system model information because real imaging systems are subject to a number of complicated physical effects (such as positron range, photon emission angle, detector sensitivity normalization factors, intercrystal scatter et al.) [12][13][14]. On the other hand, it has been acknowledged that the quality of the system model largely determines the quality of the final reconstructed images [15][16][17], and the importance of incorporating PET system uncertainties into the reconstruction framework is well recognized yet seldom addressed [14,14,[18][19][20][21].
In this paper, we investigate the application of the uncertainty weighted least squares principle to PET image reconstruction. Our algorithm, which is based on a min-max formulation, allows the simultaneous incorporation of system model and measurement statistical uncertainties, thus providing a more robust and accurate solution.

PET System and Measurement Modeling
In the PET measurement, initially, when a positron-emitting nuclide decays in the body, the nucleus rids of itself of excess positive charge by emitting a positron, which almost immediately loses its energy by collisions in the surrounding tissues and then combines with an electron and annihilates. Two back-to-back gamma rays of equal energy are then generated. These photon pairs can be detected externally by two opposite detectors using a coincidence technique, forming a coincidence event. These acquired coincidence data record the concentration information of the isotope distributions at specific locations within the body.
However, in reality, the coincidence events may also include those that the two gamma rays originating from two unrelated position annihilations are detected within the coincidence resolving time and those that the annihilation photons have an interaction of Compton scattering and lose their directional and energy information before they arrive to the detector system. The former are referred as random coincidences and the latter are called scattered coincidences. To recover a reasonable image, random and scattered photon pairs should be subtracted from coincidence events.
In practice, the emission sinogram data y, collected in 2-D mode, is a 2-D projection matrix by scanning all detector bins at each angle. The 2-D projection matrix can be transformed into a 1-D vector in the lexicographic ordering. So y is a m by 1 vector with y~fy j Dj~1,:::,mg and m is the total number of bins of the system. x~fx i Di~1,:::,ng is a n by 1 vector with n the total number of image voxels, which represents the unknown radioactivity of emission object in voxel i. The relationship between the projection data and emission object is given through an affine transform: where D is the m by n system matrix that gives the probability of a photon emitted from voxel i being detected in projection bin j. The value of detection probability matrix D depends on various factors: such as the geometry of the detection system, detector efficiency, attenuation effects, dead time correction factors, and the scattering of photons from one crystal to another. It is typically difficult or even impossible to obtain an ideal known D matrix for any real world situation. The errors of system model are broadly classified into two groups -deterministic and statistical [12]. The deterministic errors arise because of the wellknown ambiguities (e.g. the geometric ambiguity, attenuation, positron range) presented in the formulation of the model. The statistical errors are results of the random nature of photon detecting, i.e. statistical variations in detector-pair sensitivity, temporo-spatially variant response of the detector caused by the combined effects of intercrystal penetration, cross-talk and the statistical fluctuations in the photomulitplier tube [14,15,21]. Here, the goal of our work is to recover the unknown activity distribution x based on the noisy measurements y and the system model D with random errors.

Objective Function Formulation
A way to solve the reconstruction problem is to use the statistical principles, where the objective function is: with L(x,y) being the log-likelihood, R(x) denoting a regularizing penalty term and b is a hyperparameter that controls the resolution of the reconstructed image. Then the focus of PET imaging becomes to estimate the isotope concentration x from the noisy measurement data y such thatx x~argminW(x). Please note here, in the statistical image reconstruction algorithms, the system matrix D links a tomographic image with the measurements. While the ideal D is almost impossible to obtain, the performance of estimators designed without considering these uncertainties can be severely degraded and sometimes even unacceptable(such as for small animal PET imaging). For example, if the actual system matrices were DzdD, until now, all aforementioned works are based on D alone, without taking the existence of dD into account. And this inexactness may seriously affect the accuracy and reliability of the estimation results. Here, we introduce an uncertainty weighted least square framework which considers the statistical variations of the system model D and the measurement data y.
In order to handle the uncertainty issues of the system and the measurements, a min-max cost function formulation can be adopted to achieve robust solution [22] for (2): where the notation ExE 2 Q is defined as the square of the weighted Directly solving (4) will take too much storage space and computational time. Similar to the penalized weighted method proposed by Fessler [9], we have adopted an iterative algorithm to get the convergent solution of the uncertainty penalized weighted least squares (UPWLS, or called robust least squares, RLS) framework based on state space description of PET imaging. In the following section, we will derive the UPWLS (RLS) formulation with Table 1 giving the notation of related symbols and abbreviations.

UPWLS Framework for PET Imaging
Deterministic Interpretation of PET Imaging. Here, the stationary PET inverse problem is considered. In static imaging case the concentration x is assumed to be nonvarying, which means Discretizing it, we have Please note here, x(t) means activity distribution after t th time (updated) step and u(t) represents the uncertainties of the state updating process. Together with PET observations, the PET imaging can thus be interpreted in state space description as: where v models the measurement noise. Here we treat u(t) and v(t) as random variables with mean and covariance matrix as Uncertain State Space Model for PET imaging. Now, let us consider the case of uncertainty in matrix D, the state space equations in the subsection above could be rewritten with the uncertainty item dD as: Here we introduce the norm-bounded uncertainty model [23] for the system uncertainties as where EDEƒ1 is an arbitrary contraction, fM,E d g are known real constant matrices of proper dimensions that specify how the uncertain parameters in D enter the matrix D. The case of accurate models could be obtained by setting E d~0 . For above uncertain state space description, once having a priori estimatex x(tDt) together with variance P(tDt) for the state x(t) (x x(tDt) is the estimation of x(t) at step t, with corresponding variances P(tDt).), and giving measurements y, updating the estimate of state variable fromx x(tDt) tox x(tDtz1) could be realized by solving The estimation criterion is to minimize the worst possible effects of the disturbances on the signal estimation errors, which ensures that if the disturbances are small, the estimation errors will be as small as possible. This characteristic makes the method to be appropriate for some practical problems that disturbanfces exist in both system and measurements.

UPWLS Solution for PET Image Reconstruction
To solve the objective function (15), we model uncertainties in system and measurement with a norm-bounded structure [23] as where EDEƒ1 is an arbitrary contraction, fH,E d ,E y g are known real constant matrices of proper dimensions that specify how the uncertain parameters in D enter the matrices D and y. Consider about the min-max cost function (15), by setting a corresponding relationship as it is easy to get the objective function: Which is just the equivalent form of (4) under the uncertainty model (16). Here Hy is defined as According to the result published in 2001 [22], we can obtain the following unique solution for (18) by solving: where parameterl l could be calculated by minimizing G(l) defined in Appendix S1 over the interval l[(l l ,?) (further detail about how to define G(l) and how to decidel l will be given in Appendix S1) with In practice, we can choose a reasonable approximation forl l, which is to set it equal to a multiple of its low bound l l as [22] l l~( where aw0 is designed by the user, that could be chosen to be time variant as well. For any determinedl l with the corresponding relationship in equation (17), if letR R~Ŵ W {1~R {l l {1 MM T , equation (20) could be written as By settingP and letx we can get Further more, we can obtain a form for iteration of P(tz1Dtz1) if let that is According to the conclusions above, an iterative process of UPWLS estimation for the state space based PET reconstruction would be summarized as(Please see Appendix S2 for a brief derivation): Given the uncertain model: with known D, M, E d and covariance matrix R, as defined in (12)- (14). Initialization: with given initial state estimationx x(0) ( the initial activity distributions are zero or set to the results from filtered back projection method for faster convergence), covariance matrix P(0) and measurement y(0), here, P(0) is initialized with the inversion of penalty term Q, the initial penalty Q we used is a simple quadratic smoothness penalty. Set a to get l l~(1za)EMR {1 ME and calculateR R~R{l l {1 MM T , then make initialization according to (32) Update step from a priori estimationx x(tDt) and P(tDt): with Correction step with given measurement y:

Digital Phantom Data
The well-known and widely used Shepp-Logan synthetic emission phantom ( Fig. 1) with known radioactivity concentrations has been used to evaluate our algorithm. To generate realistic data, we simulate the emission coincidence events during prompt windows and delayed  windows respectively. The prompt data has to be modified to subtract the effects of the random events. Taking these effects into account, the measured sinogram y is created based on the equation: y prompt~P oissonfy true z60%:y true z20%:y true g y delay~P oissonf60%:y true g y~y prompt {y delay Here, we model the random and scatter events to be uniform field of 60 percents and 20 percents respectively. y prompt is the number of coincident photon pairs collected in the prompt windows, while y delay is the number of coincident photon pairs collected in the delay windows. The total photon counts are set to be 100K. Then we generate 50 realizations of pseudorandom emission measurements independently. In order to investigate how the noise in the system matrix affects the reconstructed image, we generate a noisy system matrix D through changing the factors related to attenuation effects with the mean relative error in the range of 0%-15%, where the error is defined by with N p is the total number of pixels. Finally, each noisy sinogram is reconstructed with three different algorithms: the popular EM and PWLS+CG (where the penalty term is a commonly used standard quadratic smoothing penalty) [9], and the proposed UPWLS method. These processes are executed iteratively until it meets the convergence criterion, which is defined by two consecutive normalized errors x(tz1) and x(t) through Ex(tz1){x(t)Evz with z being a small constant, and x defines the normalized error between the estimated and the exact value with where J t i is the estimated value, J r i is the corresponding true value, and i indicates the pixel.
The images of the mean pixel values obtained by the three algorithms (EM [7], PWLS+CG [9,10], UPWLS) based on noisy system matrix are shown in Fig. 2, while the horizontal profiles of the 18 th row through the sample mean are plotted in Fig. 3. The figures show obviously that the PWLS+CG results have large bias, while the UPWLS framework seems free of such a bias.
A detailed statistical analysis on the estimation results against the ground truth phantom map is performed. Let N p be the total number of pixels andx x i be the final reconstruction result of pixel i respectively, and x m be the mean value of the ground truth through all the pixels, then we have the following error definitions: The bias and variance of errors are averaged over the 50 reconstructions to give the estimates bias and variance at different noise levels which are summarized in Table 2. EM and PWLS perform well in noise-free case, but degradation of the image quality is observed in the noisy system model. The UPWLS framework results demonstrate the bias and standard variation remain more stable over the changing system matrix, which can be observed more clearly by percentage. For example, the bias is improved by 25.21% and 35.48% in average over EM and PWLS algorithms, and the standard deviation is improved by 26.63% and 23.56% in average over EM and PWLS algorithms, respectively, for the case with model error of 6%. Overall, these figures and results illustrate that it is possible that small noise errors lead to large estimation errors for traditional methods. On the other hand, very stable results are obtained with UPWLS framework, showing its desired robustness for real-world problems.

Real Phantom Data
The real data set used in this study was acquired on Hamamatsu SHR-22000 scanner [24] using a 6-spheres phantom, which is   usually employed to measure the recovery coefficients. The SHR-22000 scanner is designed as a whole body imaging system. It has a 838mm detector ring diameter with a patient aperture of 600mm, an axial field of view (FOV) of 224mm, operates in 2D or 3D mode. For the phantom, there are six circular regions of different diameters. These sphere objects have diameters of 37mm, 28mm, 22mm, 17mm, 13mm, 10mm and are inserted in a circular cylinder with diameter of 200mm corresponding to a volume of 9300ml, as shown in Fig. 4. The phantom filled with pure water was located at the center of both transaxial and axial FOV in the scanner using the patient bed. Transmission data were acquired in 2-D and prompt-delayed coincidence model using rotating 68Ge rod sources with total activity around 148MBq. A long blank scan was first acquired for 60 minutes. We injected F-18 concentration with initial activity of 107.92 kBq/ml into the six spheres. A 120minutes scan was then performed. Fig. 4 shows the sinograms obtained from the emission scan. Here, system matrix D ij is computed by with a single ray approximation model, which can be approximately viewed as the length of intersection between the jth pixel and the ith projection ray, i.e. D ij~lij . Longer l ij indicates larger detection sensitivity.The random events have been removed by utilizing delayed window coincidence technique. Conventional EM methods, PWLS and our algorithm as described in the previous section have been applied to recover images from the noisy data as shown in Fig. 5. Along with the vertical profiles (34# col), as shown in Fig. 6, it is evident that UPWLS results are faster, more stable, and can achieve more robust convergence of the estimates. Since PWLS converges more slowly, a longer computation time is thus expected. In a further quantitative analysis of the algorithms, the mean concentration values with standard derivation of the estimates are summarized in Table 3. It can seen that the variance performance for the UPWLS is better than EM, while the mean performance is slightly worse. Two reasons lead to such result: (1) The parameters used in UPWLS are not yet optimal; (2) The UPWLS does not incorporate the Poisson model of the measurement data while EM does.

Discussion
There are many different approaches for PET image reconstruction from projections, most of which are based on exactly known system model or system matrix. Several attempts were made to tackle this problem based on a specific type of modelling error, such as positron range, non-colinearity of the photon pair, and depth of interaction effect et al. Unlike these existing efforts that are limited to single type of modelling error, we propose an UPWLS method to handle statistical uncertainties of the system.
The P (0), E d , and R matrices that are used in the UPWLS framework are chosen according to the confidence measurements on state x, system model and measurement noise respectively. For example, if we know that the noise in the system model is smaller than the measurement noise, we should make the E d matrix smaller than the R matrix, which de-emphasizes the importance of the uncertainties of system model relative to the measurement noise. Further, we believe that any prior knowledge of the system model and measurement data should enable us to achieve higher estimation efficiency and more robust results.
In our current implementation, the parameters are set to some empirically fixed values in all experiments: H~1,E dd iagf0:5|10 {4 g: Ideally, these parameters could be adaptively updated during the estimation process. The simulation experiments are designed to show the robustness and accuracy of the proposed method, and the physical phantom experiment is used to show its efficiency and accuracy for the real world problem. Overall, our experiment results reveal that it is possible that small noise errors in system matrix may lead to large estimation errors for exactly known model-based schemes, while the UPWLS framework produces consistent results even with highly noisy system matrix, which promises robustness for real-world problems.
On the other hand, current methods requires huge storage and expensive computation as today's PET scanners have a large set of detector pairs and the inversion of the system matrix. With the method mentioned in the Appendix S3, we are able to get the inversions faster and more accurate. Furthermore, another fast algorithm, called fast state space filter has been developed. The most interesting thing is that our framework naturally allows the combination of the reconstruction process with the data acquisition task.
As a continuation of this work, based on the tracer kinematics equation coming from compartment model, we can recover dynamic changes of tracer density in a continuous time domain for dynamic PET with an uncertain system model.
Detailed investigations on these issues are underway.

Conclusion
In this paper, we have developed an uncertainty weighted least squares framework for the estimation of activity distribution from PET measured data. The approach enables us to incorporate the system uncertainties into PET image reconstruction, thus providing more robust and stable results. Analytical and experimental results with Shepp-Logan simulation phantom data and real PET measurements demonstrate the power of the proposed method.

Supporting Information
Appendix S1 Derivations for UPWLS solution.