## 1 Introduction
Note although the notations are unified for different reconstruction problems in this paper, the meaning of these notations, such as \(\mathbf{x}\) and \(\mathbf{y}\), can be different for different imaging modalities. Vectors \(\mathbf{y}\) and \(\mathbf{x}\) are related through a system matrix \(A\); see Equation (4) below for some examples. For tomographic reconstruction problems, matrix \(A\) is usually assumed known so its estimation is not covered by this paper. Rather, we focus on how to estimate \(\mathbf{x}\) from the observed \(\mathbf{y}\) and the known system matrix \(A\). We denote the estimate of \(\mathbf{x}\) by \(\hat{\mathbf{x}}\). Statistical reconstruction \(\hat{\mathbf{x}}\) obtained by maximum penalized likelihood (MPL) (also known as maximum a posteriori (MAP)) is defined by  
\[\hat{\mathbf{x}}=\arg\max_{\mathbf{x}\geq 0}\Psi(\mathbf{x}) \tag{1}\]where \(\Psi(\mathbf{x})\) is an objective function derived from the probability distribution for measurements and the penalty function. When the \(y_{i}\)'s are assumed independent (given \(\mathbf{x}\)), the penalized likelihood objective function is  
\[\Psi(\mathbf{x})=l(\mathbf{x})-hJ(\mathbf{x}) \tag{2}\]  
where \(l(\mathbf{x})\) is the log-likelihood function given by  
\[l(\mathbf{x})=\sum_{i=1}^{n}l_{i}(\mu_{i}(\mathbf{x});y_{i}) \tag{3}\]  
Here \(h>0\) is the smoothing parameter and \(J(\mathbf{x})\) is the penalty function used to smooth \(\hat{\mathbf{x}}\). In Equation (3), \(l_{i}\) denotes the log-density function for measurement \(y_{i}\), and \(\mu_{i}\) is a function of \(\mathbf{x}\in\mathds{R}_{+}^{p}\) (here \(\mathds{R}_{+}^{p}\) denotes the non-negative orthant of \(\mathds{R}^{p}\)) representing the mean measurement of camera bin \(i\). Examples of \(\mu_{i}\) include  
\[\mu_{i}(\mathbf{x})=\left\{ \begin{array}{ll}\eta_{i}(\mathbf{x})+r_{i}&\text{emission} \\ b_{i}e^{-\eta_{i}(\mathbf{x})}+r_{i}&\text{transmission}\end{array} \right. \tag{4}\]  
where \(\eta_{i}(\mathbf{x})=A_{i}\mathbf{x}\) with \(A_{i}\) being the \(i\)th row of matrix \(A\), \(b_{i}\) is the known blank scan counts of the \(i\)th detector and \(r_{i}\) the known mean background counts. Another example is polyenergetic transmission scans (such as X-ray CT) where  
\[\mu_{i}(\mathbf{x})=\sum_{m=1}^{M}b_{im}e^{-A_{i}\mathbf{x}_{m}}+r_{i} \tag{5}\]  
and here \(\mathbf{x}_{m}=(x_{1m},\ldots,x_{pm})^{T}\) denotes the attenuation map corresponding to the \(m\)-th energy spectrum, \(\mathbf{x}\) is a vector formed by the \(\mathbf{x}_{m}\)'s and \(b_{im}\) is the blank scan count from energy spectrum \(m\). In Equation (3) the notation \(l_{i}(\mu_{i};y_{i})\) is used to emphasize that \(l_{i}\) is a function of \(\mu_{i}\) and it also involves measurement \(y_{i}\). We can also write this function as \(l_{i}(\eta_{i})\) or \(l_{i}(\mathbf{x})\) in different contexts when there is no ambiguity. However, the functional properties of \(l_{i}\) may change with respect to its different arguments. For example, if assuming \(y_{i}\) follows a Poisson distribution for either emission or transmission scans, then  
\[l_{i}(\mu_{i})=-\mu_{i}+y_{i}\log\mu_{i} \tag{6}\]  
This is clearly a concave function of \(\mu_{i}\) for both emission and transmission cases. However, for \(l_{i}(\mathbf{x})\) (treated as a function of \(\mathbf{x}\)), it may be no longer concave for transmission but still concave for emission scans. Concavity is an important property exploited by the optimization transfer algorithms. Let \(\mathbf{\mu}\) be an \(n\)-vector of all \(\mu_{i}\). The first term of Equation (2), _i.e._, \(l(\mathbf{x})\), measures similarity between \(\mathbf{y}\) and \(\mathbf{\mu}\). Different probability distributions have been used to model \(y_{i}\) even under the same imaging modality.