## 1 Introduction
For example, for emission tomography, if assuming the Poisson model for \(y_{i}\) (_i.e._, \(y_{i}\sim\text{Poisson}(\mu_{i})\)) then \(l_{i}\) is given by Equation (6), or if considering the weighted least squares then  
\[l_{i}=-(y_{i}-\mu_{i})^{2}/w_{i} \tag{7}\]  
where \(w_{i}\) is the weight. When \(w_{i}=\mu_{i}\) we have the weighted least squares model as suggested in [11]. Another example in emission (or transmission) tomography is the randoms-precorrected PET scan (assume no scattering to simplify). In this context, the observed measurements are \(y_{i}=y_{i}^{\text{prompt}}-y_{i}^{\text{Delay}}\)where \(y_{i}^{\text{Prompt}}\) and \(y_{i}^{\text{Delay}}\) (both unavailable directly) denotes the number of coincidences of the prompt and delayed windows respectively. Although we can assume \(y_{i}^{\text{Prompt}}\sim\text{Poisson}(A_{i}\mathbf{x}+r_{i})\) and \(y_{i}^{\text{Delay}}\sim\text{Poisson}(r_{i})\) and that they are independent, the exact distribution of \(y_{i}\) cannot be derived directly (e.g., [9]). An approximate probability model suggested in [9] is the shifted Poisson distribution, namely \(y_{i}+2r_{i}\sim\text{Poisson}(A_{i}\mathbf{x}+2r_{i})\), which gives  
\[l_{i}=-(A_{i}\mathbf{x}+2r_{i})+(y_{i}+2r_{i})\log(A_{i}\mathbf{x}+2r_{i}) \tag{8}\]  
or the weighted least squares given by  
\[l_{i}=-(y_{i}-A_{i}\mathbf{x}_{i})^{2}/(A_{i}\mathbf{x}+2r_{i}) \tag{9}\]  
Note that the shifted Poisson approximation matches the first two moments with the true probability model for \(y_{i}+2r_{i}\) when both the prompt and delayed measurements are assumed independent and follow Poisson distributions. In this paper, we present and discuss several important non-negatively constrained penalized likelihood reconstruction algorithms. When designing a reconstruction algorithm in tomographic imaging, one considers the following important issues: (i) the algorithm is computationally efficient, and ideally it involves only forward-projection (e.g., \(A\mathbf{x}\)) and back-projection (e.g., \(A^{T}\mathbf{y}\)) operations; (ii) the algorithm can be easily applied to different measurement probability models and imaging modalities; (iii) the algorithm can impose the non-negativity constraint; (iv) the algorithm converges fast. Our discussions on the algorithms in this paper will mainly focus on these points. In tomographic imaging, it is important to produce smoothed reconstructions as severe noise in a reconstruction can cause false diagnoses. Smoothing can generally be achieved by one of the following five practices: (i) early termination of the iterations (e.g., [12]); (ii) MPL reconstructions with an appropriate smoothing parameter (e.g., [13]); (iii) functional representation of the unknown image by a set of smooth basis functions (e.g., [14]); (iv) post smoothing of the reconstruction within each iteration (e.g., [15]) or after all iterations ([16]); and (v) pre-smoothing of the camera data (_i.e._, sinogram) followed by filter backprojection (FBP) (e.g., [17,18]). We focus on the penalized likelihood approach to smoothing in this paper. In Equation (2), the smoothing parameter \(h\) balances two conflicting targets: fidelity of the \(\mu_{i}\)s to the \(y_{i}\)s and smoothness of \(\mathbf{x}\). Although an appropriate choice of \(h\) is important for achieving a reconstruction with _balanced_ fidelity and smoothness, we will not consider how to estimate \(h\) in this paper. A penalty function \(J(\mathbf{x})\) is used to smooth or regulate the estimate \(\hat{\mathbf{x}}\). Usually, \(J(\mathbf{x})\) takes the form of  
\[J(\mathbf{x})=\sum_{j=1}^{p}\rho(C_{j}\mathbf{x}) \tag{10}\]  
where \(C_{j}\mathbf{x}\) represents a neighborhood operation (such as the first or second order difference) on pixel \(j\), and function \(\rho(\cdot)\) measures the magnitude of \(C_{j}\mathbf{x}\). A common choice of \(\rho\) is the quadratic function: \(\rho(v)=\frac{1}{2}v^{2}\). Generally, a quadratic penalty tends to produce images with over-smoothed edges.