Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2895983

Formats

Article sections

- Abstract
- Introduction
- II. Physical Model Formulations
- III. Component Sinogram Restoration
- IV. Regularization Design
- V. Simulation Studies
- VI. Conclusion and Discussion
- REFERENCES

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 November 1.

Published in final edited form as:

Published online 2009 March 24. doi: 10.1109/TMI.2009.2018283

PMCID: PMC2895983

NIHMSID: NIHMS207164

Joonki Noh, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA (Email: ude.hcimu@knoojhon).

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

Dual-energy (DE) X-ray computed tomography (CT) has been found useful in various applications. In medical imaging, one promising application is using low-dose DECT for attenuation correction in positron emission tomography (PET). Existing approaches to sinogram material decomposition ignore noise characteristics and are based on logarithmic transforms, producing noisy component sinogram estimates for low-dose DECT. In this paper, we propose two novel sinogram restoration methods based on statistical models: penalized weighted least square (PWLS) and penalized likelihood (PL), yielding less noisy component sinogram estimates for low-dose DECT than classical methods. The proposed methods consequently provide more precise attenuation correction of the PET emission images than do previous methods for sinogram material decomposition with DECT. We report simulations that compare the proposed techniques and existing approaches.

THE combination of positron emission tomography (PET) and X-ray computed tomography (CT) in a single scanner has provided a variety of significant advantages in nuclear medicine [1]. First, PET/CT provides reasonably accurate alignment of functional and anatomical information. In oncology imaging, for example, PET/CT improves the identification and localization of lesions. Second, CT transmission images can be used for attenuation correction of PET emission images. CT-based attenuation correction (CTAC) for PET provides several benefits over the conventional attenuation correction by PET transmission scans [1]–[3]. For instance, low noise attenuation correction factors (ACFs) are provided by CTAC without lengthy transmission scans, and post-injection biases are avoided.

For CTAC, since X-ray source spectra for CT transmission scans typically have a broad range of energies (30 keV ~ 140 keV), one must transform CT values to the linear attenuation coefficients (LACs) evaluated at the PET energy (511 keV). Various approaches to this transformation have been suggested in the literature and these can be roughly categorized into two groups [4]: segmentation based methods and scaling methods. Segmentation based approaches first separate the CT image into regions associated with different material types such as soft tissue and bone and then replace the segmented areas with proper LACs evaluated at 511 keV based on the material types. It is difficult to make a clear segmentation of some material types [5], hampering the acceptance of the segmentation based methods.

In linear scaling, CT values are multiplied by the ratio of the LACs of water: at the CT energies and at the PET energy. However, linear scaling provides poor estimates for the LACs of bone minerals at 511 keV. Bilinear scaling resolves this problem by using two different scaling factors for different ranges of CT values [6]. One scaling factor considers water–air mixtures whereas the other considers water–bone mixtures. For objects containing materials with high atomic numbers such as iodine contrast agents, bilinear scaling can introduce quantitative errors in ACFs and these errors propagate into the reconstructed PET images [3], [7].

The classical CTAC approaches reviewed above use a single X-ray source spectrum.^{1} As alternatives, dual-energy (DE) CT-based methods, also known as dual-kVp methods, that use two different X-ray spectra have drawn attention in the literature. DECT exploits the energy dependence of LACs for the basis material characterization by collecting two sets of transmission scans [8]. Requiring no segmentation or scaling, DECT can eliminate one potential source of errors in CTAC unlike SECT. For attenuation correction of single photon emission computed tomography (SPECT), similar ideas were suggested in [9]–[11]. In DECT based methods, estimates of separate component images associated with two basis materials are reconstructed first from sinogram measurements and combined then to form ACFs at 511 keV. We now review previous methods for DECT imaging and also the use of DECT for attenuation correction in PET.

Early methods for exploiting two different energies in X-ray CT decomposed the energy dependence of LACs into two components corresponding to two types of interactions of photons: photoelectric absorption and Compton scattering [8], [12]–[14]. For the same decomposition, by singular value decomposition (SVD), [15] showed that complete energy dependent information can be achieved by collecting two sets of measurements using two different incident source spectra when the scanned object contains no *k*-edge materials near the effective energies of the source spectra.

Prior to the 1990s, reconstruction methods based on filtered backward projection (FBP) were predominant in DECT. In the early 1990s, a few algebraic iterative algorithms such as [16]–[18] were proposed for DECT. However, these methods did not account for noise statistics. Since additional X-ray scans in DECT introduce higher radiation doses than a single transmission scan, it is desirable to reduce the radiation dose as much as possible for clinical purposes. Statistical image reconstruction methods built on appropriate physical and statistical models can suppress noise, enabling the use of lower radiation doses for DECT.

Statistical approaches to reconstructing DECT images have been explored for monochromatic measurement models. For instance, [19]–[22] proposed penalized weighted least square (PWLS) methods in DECT to reconstruct soft tissue and bone mineral images with constraints in the object domain. These monochromatic methods did not fully exploit the energy dependence of LACs and required additional beam hardening corrections. Some maximum likelihood (ML) algorithms for image reconstruction in DECT were derived from measurement models considering the polychromatic nature of X-ray spectra in [23]–[25]. For attenuation correction in PET, iterative reconstruction algorithms based on PWLS approaches from polychromatic measurement models were developed recently [3], [7]. CTAC based on statistical image reconstruction in DECT provides more accurate attenuation corrections for PET than CTAC by bilinear scaling in SECT and by the FBP image reconstruction in DECT [7].

Previous DECT based methods for CTAC first reconstruct component CT images and then estimate ACFs for the PET emission images by synthesizing the obtained CT images. However, if the primary purpose of DECT is PET attenuation correction, then component images are not necessary. A synthesized sinogram at 511 keV suffices. In addition, since PWLS and PL methods are derived from proper statistical models, they provide less noisy component sinogram estimates than the classical sinogram decomposition, producing more precise ACFs in low-dose DECT [26]. Our PWLS method in Section III is straightforward to develop but is based on a simple model of the statistical properties of measurements in the projection domain, providing a suboptimal solution in terms of noise reduction for a given low radiation dose. Thus as an alternative, we also propose a PL method to estimate component sinograms. Our proposed methods generalize the previous sinogram restoration approaches developed in [27] for SECT to methods for DECT.

The remainder of this paper is organized as follows. We first introduce physical model formulations for polychromatic measurements and the object being scanned in Section II. We then review conventional approaches to decomposing component sinograms and propose two statistically motivated methods: PWLS and PL to estimate component sinograms from multiple-kVp measurements in Section III. Section IV discusses the design of regularization penalties for achieving approximately uniform spatial resolution and for matching the resolutions of component sinogram estimates. Simulations to compare the proposed DECT based methods and existing approaches are provided in Section V. Finally conclusions and discussions are presented in Section VI. Mathematical details are given in Appendix I and II.

We consider a general measurement model where multiple sets of polychromatic measurements are collected for *M*_{0} different incident spectra and forward projections (line integrals) are recorded for *N _{d}* radius-angle pairs for each incident spectrum, forming

For *m* = 1,..., *M*_{0} and *i* = 1,..., *N _{d}* let

$$\mathsf{E}\left[{y}_{mi}\right]={\stackrel{\u2012}{y}}_{mi}\triangleq \int {I}_{mi}\left(\mathcal{E}\right)\phantom{\rule{thickmathspace}{0ex}}\text{exp}\left(-{\int}_{{\mathcal{L}}_{i}}\mu (\overrightarrow{x},\mathcal{E})\mathrm{d}\ell \right)\mathrm{d}\mathcal{E}+{r}_{mi}$$

(1)

where ${I}_{mi}\left(\mathcal{E}\right)$ denotes the product of the *m*th incident source spectrum and the detector gain for the *i*th ray, and ${\int}_{{\mathcal{L}}_{i}}\cdot \mathrm{d}\ell $ is the line integral along the *i*th ray. $\mu (\overrightarrow{x},\mathcal{E})$ denotes the LAC of the object being scanned at the spatial location $\overrightarrow{x}$ and the photon energy $\mathcal{E}.\phantom{\rule{thickmathspace}{0ex}}{r}_{mi}$ denotes additive background contributions, for example room background, dark current, and scatter. We also can use *r _{mi}* to model electronic noise in a shifted Poisson approach. We treat ${I}_{mi}\left(\mathcal{E}\right)$ and

Since the number of measurements is finite whereas the LAC of the scanned object is a continuous function of $\overrightarrow{x}$ and $\mathcal{E}$, we parameterize $\mu (\overrightarrow{x},\mathcal{E})$ using a basis material decomposition. We model the LAC using a set of basis functions that are separable in the space and energy domains [19], [22], [32] as follows:

$$\mu (\overrightarrow{x},\mathcal{E})=\sum _{l=1}^{{L}_{0}}{\beta}_{l}\left(\mathcal{E}\right){\rho}_{l}\left(\overrightarrow{x}\right)$$

(2)

where ${\beta}_{l}\left(\mathcal{E}\right)$ denotes mass attenuation coefficient (MAC) and ${\rho}_{l}\left(\overrightarrow{x}\right)$ is the unknown density map of the *l*th material type. *L _{0}* is the number of material types comprising the object being scanned. In DECT, we usually have

Combining the measurement model in (1) and the object model in (2) yields the following simplified expression for the ensemble mean of measurements:

$${\stackrel{\u2012}{y}}_{mi}={I}_{mi}{e}^{-{f}_{mi}\left({\mathit{s}}_{i}\right)}+{r}_{mi}$$

(3)

for *m* = 1,...,*M*_{0} and *i* = 1,...,*N _{d}*, where

$${f}_{mi}\left({\mathit{s}}_{i}\right)\triangleq -\text{log}\phantom{\rule{thinmathspace}{0ex}}{v}_{mi}\left({\mathit{s}}_{i}\right),$$

(4)

$${v}_{mi}\left({\mathit{s}}_{i}\right)\triangleq \int {p}_{mi}\left(\mathcal{E}\right){e}^{-\sum _{l=1}^{{L}_{0}}{\beta}_{l}\left(\mathcal{E}\right){\mathit{s}}_{li}}\mathrm{d}\mathcal{E}.$$

(5)

The nonlinear function *f _{mi}*(

$${I}_{mi}\triangleq \int {I}_{mi}\left(\mathcal{E}\right)\mathrm{d}\mathcal{E},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{p}_{mi}\left(\mathcal{E}\right)\triangleq \frac{{I}_{mi}\left(\mathcal{E}\right)}{{I}_{mi}},$$

(6)

$${\mathit{s}}_{i}\triangleq {({s}_{1i},\dots ,{s}_{{L}_{0}i})}^{T},\phantom{\rule{1em}{0ex}}{s}_{li}\triangleq {\int}_{{\mathcal{L}}_{i}}{\rho}_{l}\left(\overrightarrow{x}\right)\mathrm{d}\ell .$$

(7)

The nonlinear function *f _{mi}*(

Usually in DECT, we reconstruct the component density images, and ${\rho}_{1}\left(\overrightarrow{x}\right)$ and ${\rho}_{2}\left(\overrightarrow{x}\right)$ from the DE sinograms. For the purpose of PET attenuation correction, however, it is sufficient to have sinogram-domain estimates of the component material integrals ${\left\{{\mathit{s}}_{i}\right\}}_{i=1}^{{N}_{d}}$. Therefore we focus on estimating ${\left\{{\mathit{s}}_{i}\right\}}_{i=1}^{{N}_{d}}$ hereafter.

We discuss methods to recover component sinograms, i.e., ${\left\{{s}_{li}\right\}}_{i=1}^{{N}_{d}}$ for *L*_{0} material types from noisy measurements, i.e., ${\left\{{y}_{mi}\right\}}_{i=1}^{{N}_{d}}$ for *M*_{0} incident source spectra in this section. We first review classical material sinogram decomposition, and then propose two statistically principled approaches, PWLS and PL, for estimating component sinograms with improved accuracy and/or precision.

Given the noisy measurement *y _{mi}*, the conventional method for estimating values of the nonlinear function

$${\widehat{f}}_{mi}\triangleq -\text{log}\left\{\text{smooth}\left(\frac{{y}_{mi}-{r}_{mi}}{{I}_{mi}}\right)\right\}$$

(8)

where smoothing in the radial (detector) direction is often applied to reduce noise [35]. Equating (4) and (8) yields a system of *M*_{0} nonlinear equations and *L*_{0} unknowns for the *i*th ray, where the *l*th unknown variable is *s _{li}*. Since we usually have the same number of source spectra and material types, that is

$${\widehat{\mathit{s}}}_{i}\triangleq {\mathit{f}}_{i}^{-1}\left({\widehat{\mathit{f}}}_{i}\right),\phantom{\rule{1em}{0ex}}i=1,\dots ,{N}_{d}$$

(9)

where ${\mathit{f}}_{i}\triangleq {({f}_{1i},\dots ,{f}_{{M}_{0}i})}^{T}$. This is called the conventional sinogram decomposition or sinogram preprocessing approach in DECT. Note that the measurement noise was ignored in (8) and (9), yielding noisy estimates of component sinograms and hampering their acceptance for low-dose DECT.

To reduce the noise in the classical sinogram decomposition, one could estimate ${\left\{{\mathit{s}}_{i}\right\}}_{i=1}^{{N}_{d}}$ by minimizing the following function:

$$\Lambda \left(\mathbf{s}\right)\triangleq \sum _{i=1}^{{N}_{d}}{({\widehat{\mathit{s}}}_{i}-{\mathit{s}}_{i})}^{T}{\mathit{Q}}_{i}({\widehat{\mathit{s}}}_{i}-{\mathit{s}}_{i})+R\left(\mathbf{s}\right)$$

(10)

where a possible choice of the weight matrix for the *i*th ray is ${\mathit{Q}}_{i}={\text{cov}}^{-1}\left({\widehat{\mathit{s}}}_{i}\right)$. Its approximation can be found by the method proposed in [36]. The *L*_{0} × *N _{d}* sinogram matrix

Instead of solving the *M*_{0} nonlinear equations, we propose to estimate component sinograms from the ${\widehat{f}}_{mi}$ values in (8) by minimizing a PWLS cost function. Subject to the nonnegativity constraint on the entries of *s*_{i}, we have

$${\widehat{\mathbf{s}}}_{\mathrm{PWLS}}=\underset{\mathbf{s}\in {\mathbf{R}}^{{L}_{0}\times {N}_{d}}}{\text{arg min}}\Phi \left(\mathbf{s}\right)$$

(11)

and

$$\Phi \left(\mathbf{s}\right)\triangleq \sum _{i=1}^{{N}_{d}}\frac{1}{2}{\left({\widehat{\mathit{f}}}_{i}-{\mathit{f}}_{i}\left({\mathit{s}}_{i}\right)\right)}^{T}{\mathit{D}}_{i}\left({\widehat{\mathit{f}}}_{i}-{\mathit{f}}_{i}\left({\mathit{s}}_{i}\right)\right)+R\left(\mathbf{s}\right)$$

(12)

where $\mathbf{s}\triangleq {\left\{{\mathit{s}}_{i}\right\}}_{i=1}^{{N}_{d}}$ denotes the *L*_{0} × *N _{d}* sinogram matrix. For the

$$\text{var}\left({\widehat{f}}_{mi}\right)\approx \frac{\text{var}\left({y}_{mi}\right)}{{({\stackrel{\u2012}{y}}_{mi}-{r}_{mi})}^{2}}\approx \frac{1}{{y}_{mi}}.$$

(13)

For the roughness penalty function *R*(**s**) in (12), we use

$$R\left(\mathbf{s}\right)\triangleq \sum _{l=1}^{{L}_{0}}{\gamma}_{l}\sum _{k=1}^{K}\frac{1}{2}{\left({\left[\mathbf{C}{\stackrel{\u02c7}{\mathit{s}}}_{l}\right]}_{k}\right)}^{2}$$

(14)

where *γl* denotes a regularization parameter controlling the tradeoff between data fidelity and roughness penalty. **C** is a second-order difference matrix and the column vector ${\stackrel{\u02c7}{\mathit{s}}}_{l}$ denotes the *l*th component sinogram, i.e., ${\stackrel{\u02c7}{\mathit{s}}}_{l}\triangleq {({s}_{l1},\dots ,{s}_{l{N}_{d}})}^{T}$. We choose **C** to regularize component sinograms only in the radial direction [38]; thus *K* ≈ *N _{d}*.

We use the optimization transfer principle (OTP) [39] to minimize Φ(**s**). In the framework of OTP, we design a sequence of separable quadratic surrogates (SQSs) satisfying the surrogate conditions.^{3} Nonnegativity constraints on the sinogram matrix **s** are easily imposed since the surrogates for PWLS data fitting term and the roughness penalty term are additively separable. The SQSs allow a simultaneous update of *s _{li}* for

$${s}_{li}^{(n+1)}={\left[{s}_{li}^{\left(n\right)}-\frac{1}{{h}_{li}}\frac{\partial \Phi \left({s}_{li}^{\left(n\right)}\right)}{\partial {s}_{li}}\right]}_{+}$$

(15)

where ${\left[a\right]}_{+}(\triangleq \text{max}(a,0))$ enforces the nonnegativity constraint on *s _{li}*, and ${s}_{li}^{\left(n\right)}$ denotes the estimate of

Although the PWLS approach described in the previous section is statistically motivated, it requires the logarithmic transformation in (8) to obtain ${\widehat{f}}_{mi}$. Thus the solution provided by (12) is based on an incomplete model of the data statistics and can be suboptimal in terms of noise reduction for a given radiation dose [26]. To further improve the estimated component sinograms, we now propose a PL method that uses the raw measurements, i.e., ${\left\{{y}_{mi}\right\}}_{m=1,i=1}^{{M}_{0},{N}_{d}}$.

For simplicity, we assume that the measurement *y _{mi}* obeys a Poisson distribution with the ensemble mean in (1), i.e.,

$${y}_{mi}~\text{Poisson}\phantom{\rule{thickmathspace}{0ex}}\left({\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)\right).$$

(16)

One can easily generalize this to include additive electronic noise via the shifted Poisson approach [28], [42]. Based on (16), we define the PL cost function as follows:

$$\Psi \left(\mathbf{s}\right)\triangleq -L\left(\mathbf{s}\right)+R\left(\mathbf{s}\right)=\sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{L}_{mi}\left({\mathbf{s}}_{i}\right)+R\left(\mathbf{s}\right)$$

(17)

where denotes the sinogram matrix. The negative log-likelihood is given by

$${L}_{mi}\left({\mathit{s}}_{i}\right)\triangleq {\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)-{y}_{mi}\phantom{\rule{thinmathspace}{0ex}}\text{log}\phantom{\rule{thinmathspace}{0ex}}{\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)$$

(18)

where constants independent of **s**_{i} are ignored and ${g}_{mi}\left(x\right)\triangleq x-{y}_{mi}\phantom{\rule{thinmathspace}{0ex}}\text{log}\phantom{\rule{thinmathspace}{0ex}}x$ is a convex function with respect to *x*.

With the same form of the roughness penalty function, *R*(**s**), as (14), we estimate component sinograms by performing the following PL minimization:

$${\widehat{\mathbf{s}}}_{\mathrm{PL}}=\underset{\mathbf{s}\in {\mathbf{R}}^{{L}_{0}\times {N}_{d}}}{\text{arg min}}\Psi \left(\mathbf{s}\right)$$

(19)

subject to the nonnegativity of the entries of *s*_{i}. We exploit OTP to achieve a sequence of SQSs for ${\widehat{\mathbf{s}}}_{PL}$. Since the PL data fidelity term *L*(**s**) and penalty term are additively separable, SQSs can be derived and the nonnegativity constraint is easily imposed in the framework of OTP. After applying OTP, we arrive at the following equation for updating an estimate of component sino-grams in the *n*th step. For any *l* and *i*, we have

$${s}_{li}^{(n+1)}={\left[{s}_{li}^{\left(n\right)}-\frac{1}{{d}_{li}}\frac{\partial \Psi \left({s}_{li}^{\left(n\right)}\right)}{\partial {s}_{li}}\right]}_{+}.$$

(20)

The curvature *d _{li}* is derived in (62) of Appendix I-B

There is no straightforward way to show that the updates provided by (20) monotonically decrease the PL cost function Ψ(**s**) at every iteration because of the approximations used to make the PL curvature *d _{li}* precomputable. However, if the approximations used are reasonably accurate or if

After estimating all component sinograms for *L*_{0} material types, i.e., ${\left\{{\widehat{s}}_{li}\right\}}_{i=1}^{{N}_{d}}$ for *l* = 1,..., if desired, one can use the sinogram estimates to reconstruct *L*_{0} component images, e.g., soft tissue and bone mineral in DECT. A straightforward approach is to apply FBP to each component sinogram for estimating the corresponding basis material image ${\rho}_{l}\left(\overrightarrow{x}\right)$. This approach usually gives less noisy estimates of component images than the classical sinogram decomposition combined with FBP [26]. In addition, it is less computationally expensive than fully iterative methods for reconstructing component images, e.g., [3], [7], [23]. However, in this paper, we focus on using the estimated component sinograms ${\left\{{\mathit{s}}_{i}\right\}}_{i=1}^{{N}_{d}}$ to compute ACFs, enhancing the quality of the PET emission images.

Matching the resolution of ACFs to that of the PET images is important to avoid artifacts [43] and requires appropriate regularization parameters in (14). We analyze the local impulse response (LIR) of the proposed PL component sinogram estimates, showing that the component sinograms do not have spatially uniform resolution, and do not have matched spatial resolutions. Therefore, we design modified regularizing penalty functions that provide approximate resolution uniformity and match, extending the ideas in [44] to DECT transmission tomography.

The LIR measures the changes in the estimated component sinograms induced by the perturbation of a particular element of a component sinogram. Extending [44] to our PL sinogram restoration problem for DECT, we define the following LIR, focusing on *L*_{0} = *M*_{0} = 2:

(21)

where ${\ell}_{l}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$ describes the responses that appear on all component sinograms when we place an impulse at the *j*th element of the *l*th component sinogram and is a lexicographically ordered 2*N _{d}* × 1 vector. $\stackrel{\u2012}{y}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$ denotes a 2

From (14) and (17), it can be shown that (21) reduces to the following LIR in DECT, when an impulse is placed on the first component sinogram:

$${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})={[\mathbf{F}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})+\mathbf{R}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})]}^{-1}\phantom{\rule{thickmathspace}{0ex}}\mathbf{F}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2}){e}_{1}^{j}$$

(22)

where ${e}_{1}^{j}$ denotes a 2*N _{d}* × 1 unit vector containing 1 at the position corresponding to the

$$\mathbf{F}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\triangleq \left[\begin{array}{cc}\hfill {\mathbf{D}}_{1}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\hfill & \hfill {\mathbf{D}}_{12}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\hfill \\ \hfill {\mathbf{D}}_{12}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\hfill & \hfill {\mathbf{D}}_{2}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\hfill \end{array}\right]$$

(23)

$$\mathbf{R}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\triangleq \left[\begin{array}{cc}\hfill {\gamma}_{1}{\mathbf{C}}^{T}\mathbf{C}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\gamma}_{2}{\mathbf{C}}^{T}\mathbf{C}\hfill \end{array}\right]$$

(24)

where **C** is from (14) and

$$\begin{array}{cc}\hfill {\mathbf{D}}_{1}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})& ={\text{diag}}_{i=1}^{{N}_{d}}\left\{\sum _{m=1}^{2}\frac{{({\stackrel{\u2012}{y}}_{mi}-{r}_{mi})}^{2}}{{\stackrel{\u2012}{y}}_{mi}}{\left(\frac{{w}_{mi,1}}{{v}_{mi}}\right)}^{2}\right\}\hfill \\ \hfill {\mathbf{D}}_{12}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})& ={\text{diag}}_{i=1}^{{N}_{d}}\left\{\sum _{m=1}^{2}\frac{{({\stackrel{\u2012}{y}}_{mi}-{r}_{mi})}^{2}}{{\stackrel{\u2012}{y}}_{mi}}\cdot \frac{{w}_{mi,1}{w}_{mi,2}}{{v}_{mi}^{2}}\right\}\hfill \\ \hfill {\mathbf{D}}_{2}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})& ={\text{diag}}_{i=1}^{{N}_{d}}\left\{\sum _{m=1}^{2}\frac{{({\stackrel{\u2012}{y}}_{mi}-{r}_{mi})}^{2}}{{\stackrel{\u2012}{y}}_{mi}}{\left(\frac{{w}_{mi,2}}{{v}_{mi}}\right)}^{2}\right\}\hfill \end{array}$$

and for *l* = 1,2,

$${w}_{mi,l}\left({\mathit{s}}_{i}\right)\triangleq \int {p}_{mi}\left(\mathcal{E}\right){\beta}_{l}\left(\mathcal{E}\right){e}^{-\beta {\left(\mathcal{E}\right)}^{T}{\mathit{s}}_{i}}\mathrm{d}\mathcal{E}.$$

For given $\mathcal{E},\phantom{\rule{thickmathspace}{0ex}}\beta \left(\mathcal{E}\right)\triangleq {({\beta}_{1}\left(\mathcal{E}\right),\dots ,{\beta}_{{L}_{0}}\left(\mathcal{E}\right))}^{T}$ denotes a MAC vector. Appendix II provides the detailed derivation of (22) for the PL method.

From (22), we conclude that the LIR in (22) is shift variant since the FIM contains the block matrices ${\mathbf{D}}_{1}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$ and ${\mathbf{D}}_{2}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$ whose entries vary along the diagonal, i.e., as the index for rays changes. Therefore, the PL method with conventional quadratic penalty functions in (14) yields nonuniform spatial resolution. The LIR in (22) also reveals coupling of the two component sinogram estimates induced by the terms, ${\mathbf{D}}_{12}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$ in the FIM. Thus, a perturbation of the soft tissue sinogram affects both the soft tissue and bone estimates. These interactions are due to the coupling of the two component sinograms through ${f}_{mi}({s}_{1i},{s}_{2i})$. The conventional sinogram decomposition also has similar cross-coupling. We discuss next a method to mitigate the spatial nonuniformity of the component sinogram estimates by designing spatially variant penalty functions.

We now modify the penalty matrix **R** in (24) into a new spatially variant penalty matrix to make the LIR of the PL method have approximately uniform resolution and matched spatial resolution of the estimated component sinograms, extending [44]. This also simplifies the choice of the regularization parameters, γ_{1} and γ_{2}.

Assuming that the block matrices in the FIM vary slowly along the diagonal, we approximate the FIM near the *j*th element as follows:

$$\mathbf{F}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\approx \left[\begin{array}{cc}\hfill {\left[{\mathbf{D}}_{1}\right]}_{j}{\mathbf{I}}_{{N}_{d}}\hfill & \hfill {\left[{\mathbf{D}}_{12}\right]}_{j}{\mathbf{I}}_{{N}_{d}}\hfill \\ \hfill {\left[{\mathbf{D}}_{12}\right]}_{j}{\mathbf{I}}_{{N}_{d}}\hfill & \hfill {\left[{\mathbf{D}}_{2}\right]}_{j}{\mathbf{I}}_{{N}_{d}}\hfill \end{array}\right]\triangleq {\mathbf{S}}_{j}\otimes {\mathbf{I}}_{{N}_{d}}$$

(25)

where, for example, ${\left[{\mathbf{D}}_{1}\right]}_{j}$ denotes the *j*th diagonal entry of ${\mathbf{D}}_{1}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2}),\phantom{\rule{thickmathspace}{0ex}}{\mathbf{I}}_{{N}_{d}}$ is an *N _{d}* ×

$${\mathbf{S}}_{j}\triangleq \left[\begin{array}{cc}\hfill {\left[{\mathbf{D}}_{1}\right]}_{j}\hfill & \hfill {\left[{\mathbf{D}}_{12}\right]}_{j}\hfill \\ \hfill {\left[{\mathbf{D}}_{12}\right]}_{j}\hfill & \hfill {\left[{\mathbf{D}}_{2}\right]}_{j}\hfill \end{array}\right].$$

(26)

Substituting (25) into (22) and simplifying yield

$${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\approx \left({\mathbf{S}}_{j}^{-{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right){\left[{\mathbf{I}}_{2{N}_{d}}+\left({\mathbf{S}}_{j}^{-{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right)\mathbf{R}\times \left({\mathbf{S}}_{j}^{-{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right)\right]}^{-1}\left({\mathbf{S}}_{j}^{{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right){e}_{1}^{j}$$

(27)

where ${\mathbf{S}}_{j}^{{\scriptstyle \frac{1}{2}}}$ denotes a square root factorization of the matrix **S**_{j} defined in (26).

Instead of a more standard penalty matrix whose nonzero block matrices are Toeplitz in (24), we choose a new penalty matrix given by

$$\stackrel{~}{\mathbf{R}}\triangleq \left[\begin{array}{cc}\hfill {\gamma}_{1}{\stackrel{~}{\mathbf{C}}}_{1}^{T}{\stackrel{~}{\mathbf{C}}}_{1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\gamma}_{2}{\stackrel{~}{\mathbf{C}}}_{2}^{T}{\stackrel{~}{\mathbf{C}}}_{2}\hfill \end{array}\right]$$

(28)

where ${\stackrel{~}{\mathbf{C}}}_{1}$ and ${\stackrel{~}{\mathbf{C}}}_{2}$ are spatially variant second-order difference matrices for the first material type (soft tissues) and the second material type (bone minerals), respectively. To provide resolution uniformity of component sinogram estimates, we define ${\stackrel{~}{\mathbf{C}}}_{1}$ and ${\stackrel{~}{\mathbf{C}}}_{2}$ as follows:

$${\stackrel{~}{\mathbf{C}}}_{1}\triangleq \mathbf{C}{\stackrel{~}{\mathbf{D}}}_{1},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\stackrel{~}{\mathbf{D}}}_{1}\triangleq {\text{diag}}_{i=1}^{{N}_{d}}\left\{\sqrt{{\left[{\mathbf{D}}_{1}\right]}_{i}}\right\}$$

(29)

$${\stackrel{~}{\mathbf{C}}}_{2}\triangleq \mathbf{C}{\stackrel{~}{\mathbf{D}}}_{2},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{\stackrel{~}{\mathbf{D}}}_{2}\triangleq {\text{diag}}_{i=1}^{{N}_{d}}\left\{\sqrt{{\left[{\mathbf{D}}_{2}\right]}_{i}}\right\}.$$

(30)

By replacing the conventional penalty matrix **R** with the new penalty matrix $\stackrel{~}{\mathbf{R}}$, we have a useful approximation^{4} to obtain ${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$

$$\left({\mathbf{S}}_{j}^{-{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right)\stackrel{~}{\mathbf{R}}\left({\mathbf{S}}_{j}^{-{\scriptstyle \frac{1}{2}}}\otimes {\mathbf{I}}_{{N}_{d}}\right)\approx \Gamma \otimes {\mathbf{C}}^{T}\mathbf{C}$$

(31)

where the regularizing parameter matrix is defined as

$$\Gamma \triangleq \left[\begin{array}{cc}\hfill {\gamma}_{1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\gamma}_{2}\hfill \end{array}\right].$$

(32)

For small and γ_{1} and γ_{2}, substituting (31) into (27) yields the following approximation to the LIR:

$${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})\approx {[{\mathbf{I}}_{2{N}_{d}}+\Gamma \otimes {\mathbf{C}}^{T}\mathbf{C}]}^{-1}{e}_{1}^{j}$$

(33)

indicating approximately uniform resolution of component sinogram estimates since **C**^{T }**C** is Toeplitz. Using (33), we are able to tabulate the relationship between the spatial resolution of a component sinogram estimate, usually quantified by full-width half-maximum (FWHM), and the regularizing parameter γ_{l}. For a given target FWHM, we can then determine the corresponding regularizing parameter. The assumption that γ_{1} and γ_{2} are small values is reasonable since we will use a hybrid approach combining the penalizing methods with small regularizing parameters and postsmoothings in Section V as recommended in [45].

We performed several simulations to evaluate the proposed PWLS and PL component sinogram estimates for attenuation correction of the PET emission images. First, we compared the classical DECT sinogram decomposition and the statistically motivated PWLS and PL using a NCAT phantom consisting of soft tissues and bone minerals [46]. Second, we compared the classic bilinear scaling with a single-kVp spectrum, and the PWLS and PL methods with dual-kVp spectra using the same NCAT phantom but containing iodine contrast agents. In all the results below, we used the modified regularizer in (28) for the PWLS and PL methods.

Fig. 1 shows two true component densities: soft tissues and bone minerals of the NCAT phantom used in simulations. This phantom contains 512 × 512 pixels and the pixel size is 0.1 × 0.1 cm^{2}. Fig. 2 shows two source spectra that are incident on the NCAT phantom having 80 kVp and 140 kVp, where two dashed vertical lines at ${\stackrel{\u2012}{\mathcal{E}}}_{1}=57\phantom{\rule{thickmathspace}{0ex}}\mathrm{keV}$ and ${\stackrel{\u2012}{\mathcal{E}}}_{2}=72\phantom{\rule{thickmathspace}{0ex}}\mathrm{keV}$ denote their effective energies in (61), respectively. For simplicity, we used the same incident spectra for each ray, i.e., *I _{mi} = I_{m}*, which ignores the effects of bow tie filters. Since clinical SECT scans have in the order of 10

Two component densities of the NCAT phantom used in simulations: (a) the density map of soft tissues and (b) the density map of bone minerals.

Two incident spectra ${I}_{m}\left(\mathcal{E}\right)$ against energy $\mathcal{E}$ [keV] for *m* = 1,2: 80 kVp (top) and 140 kVp (bottom). The dashed vertical lines indicate the effective energies, ${\stackrel{\u2012}{\mathcal{E}}}_{1}$ and ${\stackrel{\u2012}{\mathcal{E}}}_{2}$, respectively.

To provide approximately uniform spatial resolution of the component sinogram estimates, we used the proposed penalty matrix $\stackrel{~}{\mathbf{R}}$ in (28) with small amount of regularization by choosing γ_{1} = γ_{2} = 2^{−8} for the PWLS and PL methods. We then applied a postsmoothing filter to the restored component sinograms so that the three DECT based methods have matched spatial resolution [45]. The cost functions of the PWLS and PL methods were checked at every iteration, verifying that the algorithms in (15) and (20) monotonically decreased the corresponding cost functions, respectively. After estimating the component sinograms, we computed the PET ACFs as follows. For the *i*th ray

$${\mathrm{ACF}}_{i}\triangleq \text{exp}{\phantom{\mid}\left(\sum _{l=1}^{{L}_{0}}{\beta}_{l}\left(\mathcal{E}\right){\widehat{s}}_{li}\right)\mid}_{\mathcal{E}=511\phantom{\rule{thickmathspace}{0ex}}\mathrm{keV}}.$$

(34)

We applied these ACFs to the PET sinogram and applied FBP to reconstruct the emission images as shown in Fig. 3. The normalized root mean squared error (NRMSE) of the PET image with ACF by conventional DECT decomposition was 12%, whereas the DECT-PWLS and DECT-PL methods yielded a lower NRMSE of 7.4%. These are global values over the whole PET image. Fig. 4 shows the component CT sinograms (soft tissue and bone) restored by DECT-PL and the corresponding FBP reconstructed component CT images. Further comparisons of the component CT sinograms and reconstructed images are reported in [26].

True PET emission image and reconstructed PET images by three DECT based attenuation corrections: (a) true PET emission image, (b) reconstructed PET image with ACF by the conventional DECT sinogram decomposition, (c) reconstructed PET image with ACF by **...**

To compare the classical bilinear scaling with a single-kVp spectrum and our PWLS and PL methods with dual-kVp spectra, we placed a small amount of iodine contrast agents in three areas of the NCAT phantom. Two of them were added into regions in the heart and one was placed near the border of the left lung. Fig. 5(a) shows the resulting true CT density image. We assumed that the contrast was diluted and one of them, in the anterior part of the heart, is indistinguishable from the soft tissues surrounding it in the true PET emission image shown in Fig. 5(b). To quantify the effects of errors in the estimated ACFs on the reconstructed PET emission images, we again used the noiseless PET image. The same regularizing parameters and spatially variant penalty matrices were applied and the same postsmoothing filters were used as in the previous section.

True CT density (left) and PET image (right) containing iodine contrast agents in the three regions. Two iodine contrast agents (red arrows) in the center of the NCAT phantom correspond to heart and the other one (green arrow) in the left side is associated **...**

In SECT based methods combined with bilinear scaling, sinograms were first restored by three different approaches: 1) the standard sinogram preprocessing, 2) PWLS method, and 3) PL method. CT images were then reconstructed by FBP before performing bilinear scaling. We call these three approaches SECT-BS, SECT-PWLS-BS, SECT-PL-BS, respectively. For the SECT based methods, we used a spectrum with 140 kVp whose shape is the same as ${I}_{2}\left(\mathcal{E}\right)$. To ensure that we would not bias the results in favor of the DECT based approaches, we set the number of incident photons per ray to be 5 × 10^{5} for generating the SECT data, more than twice the total photons in the DECT scans.

After compensating for attenuation in the PET images by six competing methods, three of which are based on DECT and the remaining three are based on SECT and bilinear scaling, FBPs produce the reconstructed PET emission images shown in Fig. 6 with their NRMSEs given in Table I. The DECT-PWLS and DECT-PL methods yielded lower NRMSE values than the classical DECT sinogram decomposition in the attenuation correction of the PET images when iodine contrast agents are present. Similarly, the statistically motivated SECT based approaches, SECT-PWLS-BS and SECT-PL-BS, provide better attenuation corrections for the PET images than the standard sinogram preprocessing in SECT. Table I also shows that the DECT methods have lower NRMSE values than their SECT counterparts.

Reconstructed PET emission images by six competing attenuation correction methods in the presence of iodine contrast agents: (a) conventional DECT sinogram decomposition, (b) DECT-PWLS method, (c) DECT-PL method, (d) SECT-BS method, (e) SECT-PWLS-BS method, **...**

We also compared the local NRMSEs of four selected local regions in the reconstructed PET emission images. Fig. 7 shows these four local areas chosen from the true PET emission image. Table II shows that the DECT-PWLS and DECT-PL provided the lowest NRMSEs. Note that Region 2 contains the area where diluted iodine contrast agent was indistinguishable from local soft tissues in the true PET emission image and Region 4 contains soft tissues only.

Four selected regions for local NRMSE analysis marked by white boxes. The region numbers were counted clockwise from Region 1 to 4.

We conducted additional simulations to investigate the bias and variance of the reconstructed PET images corrected by the DECT based methods. We synthesized 50 Poisson distributed realizations of two sets of CT measurements in the projection domain from the NCAT phantom having iodine contrast. After correcting attenuation by CTACs and reconstructing the PET emission images by FBPs for all realizations, we evaluated two quantities: normalized root mean squared bias (NRMSB) and normalized root mean variance (NRMV) of the reconstructed PET emission images, summarized in Table III. Recall that the noiseless PET image was used for the simulations, producing small NRMVs. Without increasing variances, the PWLS and PL methods based on DECT yielded PET reconstructions having lower biases than did the other competing methods: DECT sino-gram material decomposition and SECT based approaches.

Biases and Standard Deviations of the Reconstructed PET Images When Iodine Contrast Agents are Present

Since the true CT density in Fig. 5(a) and PET image in Fig. 5(b) contain three different material types; soft tissues, bone minerals, and iodine whereas the object model for the DECT based CTACs in (2) assumed two basis materials, there may exist bias caused by model mismatch. To scrutinize this effect induced by iodine contrast, we increased the sizes of the three iodinated regions and repeated the simulations without Poisson noise. Fig. 8 shows horizontal profiles obtained from DECT-PL and SECT-PL-BS at two different vertical locations, one of which corresponds to the iodine contrast agent placed in the anterior part of phantom's heart [red arrow in Fig. 8(a)] and the other is associated with the posterior one of the heart [red arrow in Fig. 8(b)]. The green arrow in Fig. 8(b) marks the iodine contrast agent in the left lung. Other two DECT based CTACs and two SECT based CTACs had very similar results to those of DECT-PL and SECT-PL-BS, respectively, so are not shown. Fig. 8(a) and (b) suggest that, although the DECT based CTACs cause some biases in the reconstructed PET emission images, they are more robust to model mismatch than their SECT based counterparts. We will explore image domain statistical approaches for DECT based CTACs in the future to mitigate these biases caused by model mismatch around iodinated contrast agents.

Errors in X-ray CT-based attenuation correction of PET emission data will propagate into the reconstructed PET emission images. Such errors can arise from several sources, one of which is the inability of a SECT scan to reliably distinguish materials of differing densities and effective atomic numbers, e.g., bone and iodinated contrast agents. Although these materials can have the same reconstructed values in a CT image in terms of Hounsfield units (HU), they will have significantly different LACs at the annihilation photon energies of 511 keV used for PET. It has been recognized for some time [8] that two CT scans acquired with different spectral distributions can be used to estimate spatial density patterns of two basis functions or material components, albeit at a significant amplification of image noise. This increase is not surprising given the large overlap of two spectra. For attenuation correction of PET emission data, however, the increase in noise in some ways is not as critical due to the generally larger noise levels in PET relative to CT. In addition, there is noise reduction for CTAC due to two different forms of signal averaging. First, we do not need the separate component sinograms (or images), but rather just the sum of the two components after scaling to 511 keV. Second, while CT images are typically reconstructed on a 512 × 512 grid, PET images use a 128 × 128 grid for the same field of view (FOV).

With the considerations listed above, it becomes feasible to use DECT for attenuation correction of PET emission data although noise amplification is still an issue. In addition, an important consideration in PET/CT imaging is radiation dose to the patient. Thus methods that further reduce noise in the DECT scan components are essential for a low radiation dose.

We proposed novel PWLS and PL methods for statistical sinogram restoration in DECT used for attenuation correction in PET. The goal is to reduce bias when materials with higher atomic numbers such as iodine or metallic objects are present in the patient, which are scaled incorrectly by the standard SECT bilinear scaling approach. We also designed spatially variant penalty functions that generate component sinograms with approximately uniform spatial resolution. These methods produced more accurate ACFs than conventional approaches, reducing the overall NRMSE compared to CTAC estimates using conventional SECT and DECT methods. In addition, the statistically motivated DECT methods reduced bias when iodine is present without unduly increasing noise in the final PET image (Tables II and III). From the simulations, there did not seem to be a significant difference between the results for the DECT-PWLS and DECT-PL methods. This may be due to the level of noise and the relatively simple model used in the simulations.

In the simulations of Section V, after downsampling, the two sets of CT rays synthesized in a parallel-beam geometry matched those of the PET sinogram. To cope with more practical CT scans, e.g., axial or helical CT scans, the proposed DECT based CTACs can be modified as follows. First, the component sinograms are decomposed in the CT projection domain by the PWLS or PL method, and then the estimated sinograms are combined to synthesize a monochromatic CT sinogram at 511 keV. Second, this sinogram is reconstructed using conventional axial or helical FBP to form an attenuation map at 511 keV. Third, this attenuation map is reprojected to produce ACFs that match any PET geometry. Such backprojection/reprojection steps are routinely used in CTAC for PET.

Several important considerations are not addressed in this study. One that is well known is the potential for patient motion between two sets of CT scans, which could lead to significant artifacts. Using simultaneous or near-simultaneous acquisition of the two CT data sets, e.g., by fast kVp switching [49], can mitigate these artifacts. Furthermore, CTAC for PET requires a lower spatial resolution than diagnostic CT images, so some effects of small motions might be reduced by downsampling the CT sinograms to PET resolution. However, such downsampling might not completely suppress these motion-related artifacts because of the nonlinearities of polychromatic CT, and it is possible that effects akin to the exponential edge-gradient [50] might persist. In addition, patient motion between the CT scans and the PET scan is well known to cause other artifacts in the reconstructed PET emission images, but this problem occurs even for SECT based CTACs and is beyond the scope of this work. Another consideration is that knowledge of the X-ray spectrum for each ray or measurements of the function ${f}_{mi}\left({\mathit{s}}_{i}\right)$ in (4) is needed, which may be challenging to determine. Finally, the cases where improved accuracy in PET/CT imaging is necessary have to be delineated. There are clinical scenarios where accurate estimation of tracer uptake is not needed for PET imaging. Cases where improved estimation of PET tracer up-take by CTAC with DECT will most likely improve patient outcomes are with evaluation of responses to therapy where bone is involved and/or contrast agent is used and reduction of artifacts from prostheses and other objects.

The authors would like to thank anonymous reviewers for valuable comments on this paper.

This work was supported by the National Institute of Health under Grant 1R01CA115870.

This appendix derives (15) and (20). To obtain surrogates for the PWLS cost function $\Phi \left(\mathbf{s}\right)$ and PL cost functions $\Psi \left(\mathbf{s}\right)$, we consider the data fidelity terms and penalty terms separately. We first derive the algorithm for PWLS and then move to the algorithm for PL.

First, we express the data fidelity term of Φ in (12) as

$${\Phi}_{d}\left(\mathbf{s}\right)\triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{\varphi}_{mi}\left({\mathit{s}}_{i}\right){y}_{mi}$$

(35)

where ${\varphi}_{mi}\left({\mathit{s}}_{i}\right)\triangleq \frac{1}{2}{({\widehat{f}}_{mi}-{f}_{mi}\left({\mathit{s}}_{i}\right))}^{2}$ is the LS term for the *m*th spectrum and *i*th ray. By a Taylor series expansion of ϕ_{mi}(**s**_{i} about ${\mathit{s}}_{i}^{\left(n\right)}$, we have an inequality

(36)

where **W**_{mi} is a *L*_{0} × *L*_{0} positive definite matrix satisfying the condition that for any **s**_{i} in the matrix sense^{5} and ${\nabla}^{2}$ is a Hessian matrix. $\text{Let}\phantom{\rule{thickmathspace}{0ex}}{\varphi}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ denote the right-hand side of (36). Since ${\varphi}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ then majorizes ϕ_{mi}(**s**_{i}, we define a surrogate function for Φ(**S**) in the th step, given by

$${\Phi}_{d}^{\left(n\right)}\left(\mathbf{s}\right)\triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{\varphi}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right){y}_{mi}$$

(37)

where the exact form of **W**_{mi} is determined below.

Second, we consider the additively separable penalty term in (14). By applying De Pierro's additive convexity trick [51], we define a surrogate function for the penalty term as

$$\begin{array}{cc}\hfill {R}^{\left(n\right)}\left(\mathbf{s}\right)& \triangleq \sum _{l=1}^{{L}_{0}}\sum _{i=1}^{{N}_{d}}{\gamma}_{l}{r}_{li}^{\left(n\right)}\left({s}_{li}\right)\hfill \\ \hfill {r}_{li}^{\left(n\right)}\left({s}_{li}\right)& \triangleq \sum _{k=1}^{K}\frac{{a}_{ki}}{2}{\left(\frac{{c}_{ki}}{{a}_{ki}}\left({s}_{li}-{s}_{li}^{\left(n\right)}\right)+{\left[\mathbf{C}{\stackrel{\u02c7}{\mathit{s}}}_{l}^{\left(n\right)}\right]}_{k}\right)}^{2}\hfill \end{array}$$

(38)

where denotes the element of **C** at (*k, i*) and *a _{ki}* is given by ${a}_{ki}\triangleq \mid {c}_{ki}\mid \u2215{\sum}_{{i}^{\prime}=1}^{{N}_{d}}\mid {c}_{k{i}^{\prime}}\mid $ if ${c}_{ki}\ne 0;{a}_{ki}\triangleq 0$ if

By adding two surrogate functions ${\Phi}_{d}^{\left(n\right)}\left(\mathbf{s}\right)$ and *R*^{(n)}(**s**), we have a SQS for the PWLS cost function defined as follows:

$${\Phi}^{\left(n\right)}\left(\mathbf{s}\right)\triangleq {\Phi}_{d}^{\left(n\right)}\left(\mathbf{s}\right)+{R}^{\left(n\right)}\left(\mathbf{s}\right)$$

(39)

satisfying the surrogate conditions [41]. It can be checked that is Φ^{(n)}(**s**) is additively separable with respect to the index *l* and *i*. Therefore, Φ^{(n)} majorizes Φ(**s**). By differentiating Φ^{(n)} with respect to *s _{li}*, equating it to zero, and simplifying, one arrives at the PWLS algorithm in (15).

Finally, we need to determine **W**_{mi}. From the definition of ϕ_{mi}(**s**_{i}, its Hessian is given by

(40)

Since *f _{mi}(s_{i})* is nonnegative, monotonically increasing, and concave, we have an inequality

(41)

We define and further simplify the above inequality into

(42)

where ${\mathit{s}}_{i}^{\text{max}}\triangleq {\text{arg max}}_{{\mathit{s}}_{i}}\Vert {\nabla}^{2}{f}_{mi}\left({\mathit{s}}_{i}\right)\Vert $. For our spectra in Fig. 2, we found that ${\mathit{s}}_{i}^{\text{max}}=0$ for any *i*, but we could not prove that this will always hold. By defining **W**_{mi} from the right side of the inequality above, we have the first piece of the curvature *h _{li}* as follows:

$${w}_{mi}\triangleq \left({\Vert {\mathit{b}}_{mi}\Vert}^{2}+{\left[{\widehat{f}}_{mi}\right]}_{+}\Vert -{\nabla}^{2}{f}_{mi}\left(0\right)\Vert \right){y}_{mi}$$

(43)

producing

$${h}_{li}\triangleq \sum _{m=1}^{{M}_{0}}{w}_{mi}+{\gamma}_{l}\cdot {\stackrel{~}{c}}_{i}$$

(44)

where ${\stackrel{~}{c}}_{i}\triangleq {\sum}_{k=1}^{K}\mid {c}_{ki}\mid \cdot {\sum}_{{i}^{\prime}=1}^{{N}_{d}}\mid {c}_{k{i}^{\prime}}\mid $ and *c _{ki}* denotes the element of

First, the data fidelity term of Ψ(**s**) is expressed as

$$-L\left(\mathbf{s}\right)=\sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{L}_{mi}\left({\mathit{s}}_{i}\right)\triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{g}_{mi}\left({\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)\right)$$

(45)

where *g _{mi}(x) = x — y_{mi}* log

$$\begin{array}{cc}\hfill {\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)& =\int {I}_{mi}\left(\mathcal{E}\right){e}^{-{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}}\mathrm{d}\mathcal{E}+{r}_{mi}\hfill \\ \hfill & =\int \left[\frac{{I}_{mi}\left(\mathcal{E}\right)}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)}\right]{t}_{mi}({\mathit{s}}_{i},\mathcal{E}){b}_{mi}\times \left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)\mathrm{d}\mathcal{E}\hfill \end{array}$$

(46)

where

$${t}_{mi}({\mathit{s}}_{i},\mathcal{E})\triangleq {e}^{-{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}}+\frac{{r}_{mi}}{{I}_{mi}\left(\mathcal{E}\right)},\phantom{\rule{1em}{0ex}}{b}_{mi}({\mathit{s}}_{i},\mathcal{E})\triangleq \frac{{\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)}{{t}_{mi}\left({\mathit{s}}_{i}\mathcal{E}\right)}.$$

Note that the term inside the square bracket behaves as a probability density function (PDF) since it is nonnegative and integrates into unity over $\mathcal{E}$. By Jensen's inequality and the convexity of *g _{mi}(x)*, we have an inequality

$${g}_{mi}\left({\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)\right)\le \int \frac{{I}_{mi}\left(\mathcal{E}\right){g}_{mi}\left({t}_{mi}({\mathit{s}}_{i},\mathcal{E}){b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)\right)}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)}\mathrm{d}\mathcal{E}$$

(47)

on the basis of De Pierro's multiplicative convexity trick [52]. $\text{Let}\phantom{\rule{thickmathspace}{0ex}}{l}_{mi,1}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ denote the right-hand side of (47). Since ${l}_{mi,1}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ majorizes *g _{mi}*(

$${L}_{1}^{\left(n\right)}\left(\mathbf{s}\right)\triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{l}_{mi,1}^{\left(n\right)}\left({\mathit{s}}_{i}\right).$$

(48)

However, it is difficult to minimize ${L}_{1}^{\left(n\right)}\left(\mathbf{s}\right)$ directly.

Second, we want to seek a quadratic surrogate for by ${L}_{1}^{\left(n\right)}\left(\mathbf{s}\right)$ by applying the same idea as (36). To do that, we focus on

$${h}_{mi}(a,b,\mathcal{E})\triangleq {g}_{mi}\left({t}_{mi}({\mathit{s}}_{i},\mathcal{E}){b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)\right)$$

(49)

where two variables are defined as $a\triangleq {\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}$ and $b\triangleq {\mathit{s}}_{i}^{\left(n\right)}$. By a Taylor series expansion of ${h}_{mi}(a,b,\mathcal{E})$ with respect to *a*, we have an inequality,

$${h}_{mi}(a,b,\mathcal{E})\le {h}_{mi}\left({\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)},b,\mathcal{E}\right)+{\stackrel{.}{h}}_{mi}\left({\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)},b,\mathcal{E}\right)\times \left(a-{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)}\right)+\frac{1}{2}{\stackrel{\u02c7}{c}}_{mi}{\left(a-{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)}\right)}^{2}$$

(50)

where ${\stackrel{\u02c7}{c}}_{mi}$ is a positive constant satisfying the condition that ${\stackrel{\u02c7}{c}}_{mi}\ge {\ddot{h}}_{mi}({\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i},b,\mathcal{E})$ for any ** s_{i}**. Here ${\stackrel{.}{h}}_{mi}(a,b,\mathcal{E})$ and ${\ddot{h}}_{mi}(a,b,\mathcal{E})$ denote the first- and second-order derivatives of ${h}_{mi}(a,b,\mathcal{E})$ with respect to

$$\begin{array}{cc}\hfill {L}_{2}^{\left(n\right)}& \triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{l}_{mi,2}^{\left(n\right)}\left({\mathit{s}}_{i}\right),\hfill \\ \hfill {l}_{mi,2}^{\left(n\right)}& \triangleq \int \frac{{I}_{mi}\left(\mathcal{E}\right)}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)}{h}_{mi}^{\left(n\right)}\left({\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i},{\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)\mathrm{d}\mathcal{E}\hfill \end{array}$$

(51)

where the exact form of ${\stackrel{\u02c7}{c}}_{mi}$ is discussed below. Since ${L}_{2}^{\left(n\right)}\left(\mathbf{s}\right)$ is quadratic with respect to *s*_{i}, we have many algorithms to minimize it, for example, the coordinate descent algorithm. An alternative choice is finding a SQS function for ${L}_{2}^{\left(n\right)}\left(\mathbf{s}\right)$ to further simplify the optimization.

Thirdly, we seek a SQS for ${L}_{2}^{\left(n\right)}\left(\mathbf{s}\right)$ by applying De Pierro's additive convexity trick [51]. In this case, ${h}_{mi}^{\left(n\right)}(a,b,\mathcal{E})$ plays a key role to derive the SQS function. We express ${\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}$ in the following way:

$${\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}=\sum _{l=1}^{{L}_{0}}{\pi}_{l}\left(\mathcal{E}\right)\left(\frac{{\beta}_{l}\left(\mathcal{E}\right){s}_{li}-{\beta}_{l}\left(\mathcal{E}\right){s}_{li}^{\left(n\right)}}{{\pi}_{l}\left(\mathcal{E}\right)}+{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)}\right)$$

(52)

where the coefficient ${\pi}_{l}\left(\mathcal{E}\right)$ can be defined as

$${\pi}_{l}\left(\mathcal{E}\right)\triangleq \frac{{\beta}_{l}\left(\mathcal{E}\right)}{\sum _{{l}^{\prime}=1}^{{L}_{0}}{\beta}_{{l}^{\prime}}\left(\mathcal{E}\right)}\ge 0,\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\sum _{l=1}^{{L}_{0}}{\pi}_{l}\left(\mathcal{E}\right)=1.$$

(53)

Since ${h}_{mi}^{\left(n\right)}(a,b,\mathcal{E})$ has a quadratic form with respect to *a* and ${\pi}_{l}\left(\mathcal{E}\right)$ acts as a probability mass function (PMF), by Jensen's inequality, we have

$${h}_{mi}^{\left(n\right)}\left({\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i},{\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)\le \sum _{l=1}^{{L}_{0}}{\pi}_{l}\left(\mathcal{E}\right)\cdot {h}_{mi}^{\left(n\right)}\left(\frac{{\beta}_{l}\left(\mathcal{E}\right){s}_{li}-{\beta}_{l}\left(\mathcal{E}\right){s}_{li}^{\left(n\right)}}{{\pi}_{l}\left(\mathcal{E}\right)}+{\beta}^{T}\left(\mathcal{E}\right){\mathit{s}}_{i}^{\left(n\right)},{\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right).$$

(54)

$\text{Let}\phantom{\rule{thickmathspace}{0ex}}{k}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ denote the right-hand side of (54). Since ${k}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right)$ majorizes ${h}_{mi}^{\left(n\right)}(a,b,\mathcal{E})$, we can define a surrogate function for ${L}_{2}^{\left(n\right)}\left(\mathbf{s}\right)$ in the th step, given by

$$\begin{array}{cc}\hfill {L}_{3}^{\left(n\right)}\left(\mathbf{s}\right)& \triangleq \sum _{m=1}^{{M}_{0}}\sum _{i=1}^{{N}_{d}}{l}_{mi,3}^{\left(n\right)}\left({\mathit{s}}_{i}\right)\hfill \\ \hfill {l}_{mi,3}^{\left(n\right)}\left({\mathit{s}}_{i}\right)& \triangleq \int \frac{{I}_{mi}\left(\mathcal{E}\right)}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)}{k}_{mi}^{\left(n\right)}\left({\mathit{s}}_{i}\right)\mathrm{d}\mathcal{E}.\hfill \end{array}$$

(55)

Therefore, we have

$$-L\left(\mathbf{s}\right)\le {L}_{1}^{\left(n\right)}\left(\mathbf{s}\right)\le {L}_{2}^{\left(n\right)}\left(\mathbf{s}\right)\le {L}_{3}^{\left(n\right)}\left(\mathbf{s}\right)$$

(56)

and ${L}_{3}^{\left(n\right)}\left(\mathbf{s}\right)$ is a SQS for –*L*(**s**) by the construction.

Fourthly, we combine the surrogate function for the data fidelity term, ${L}_{3}^{\left(n\right)}\left(\mathbf{s}\right)$ and the surrogate function for the penalty term,*R ^{(n)}*(

$${\Psi}^{\left(n\right)}\left(\mathbf{s}\right)\triangleq {L}_{3}^{\left(n\right)}\left(\mathbf{s}\right)+{R}^{\left(n\right)}\left(\mathbf{s}\right)$$

(57)

satisfying the surrogate conditions. It can be checked that Ψ^{(n)}(**s**) is additively separable with respect to the index *l* and *i*. By differentiating Ψ^{(n)}(**s**) with respect to *s _{li}*, equating it to zero, and simplifying, one arrives at the same form as the PL algorithm in (20) but having

$${d}_{li}=\sum _{m=1}^{{M}_{0}}\int \frac{{I}_{mi}\left(\mathcal{E}\right){\beta}_{l}^{2}\left(\mathcal{E}\right){\stackrel{\u02c7}{c}}_{mi}}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right){\pi}_{l}\left(\mathcal{E}\right)}\mathrm{d}\mathcal{E}+{\gamma}_{l}\cdot {\stackrel{~}{c}}_{i}.$$

(58)

Finally, we now discuss how to determine the curvature ${\stackrel{\u02c7}{c}}_{mi}$. When measurements are Poisson distributed, the optimal ${\stackrel{\u02c7}{c}}_{mi}$ can be found in [53], yielding the fastest convergence. Instead of the optimal curvature, however, we pursue a precomputable curvature to reduce computational costs per iteration in this paper. Assuming that the curvature of *L _{mi}(s_{i})* varies slowly around minimizer, we have an approximation for the curvature as follows:

$${\stackrel{\u02c7}{c}}_{mi}\approx {y}_{mi}$$

(59)

where a similar idea can be found in [54].

By substituting (59) into (58), the curvature *d _{li}* is

$${d}_{li}=\sum _{m=1}^{{M}_{0}}\int {\beta}_{l}\left(\mathcal{E}\right)\sum _{j=1}^{{L}_{0}}{\beta}_{j}\left(\mathcal{E}\right)\frac{{I}_{mi}\left(\mathcal{E}\right){y}_{mi}}{{b}_{mi}\left({\mathit{s}}_{i}^{\left(n\right)},\mathcal{E}\right)}\mathrm{d}\mathcal{E}+{\gamma}_{l}\cdot {\stackrel{~}{c}}_{i}.$$

(60)

We replace $\mathcal{E}$ with the effective energy of the *m*th incident spectrum at the *i*th ray, defined as an expectation of energy $\mathcal{E}$

$${\stackrel{\u2012}{\mathcal{E}}}_{mi}\triangleq \frac{\int \mathcal{E}{I}_{mi}\left(\mathcal{E}\right)\mathrm{d}\mathcal{E}}{\int {I}_{mi}\left(\mathcal{E}\right)\mathrm{d}\mathcal{E}}.$$

(61)

This produces an approximate to the curvature *d _{li}* as follows:

$${d}_{li}\approx \sum _{m=1}^{{M}_{0}}{\beta}_{l}\left({\stackrel{\u2012}{\mathcal{E}}}_{mi}\right)\sum _{j=1}^{{L}_{0}}{\beta}_{l}\left({\stackrel{\u2012}{\mathcal{E}}}_{mi}\right){y}_{mi}+{\gamma}_{l}\cdot {\stackrel{~}{c}}_{i}.$$

(62)

Now we have the PL algorithm in (20).

If we place an impulsive perturbation in the first sinogram, the LIR from (21) can be expressed as

(63)

First from (3), it can be checked that a first-order partial derivative of * _{mi}*(

$$\frac{\partial {\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{li}}=-({\stackrel{\u2012}{y}}_{mi}\left({\mathit{s}}_{i}\right)-{r}_{mi})\frac{{w}_{mi,l}\left({\mathit{s}}_{i}\right)}{{v}_{mi}\left({\mathit{s}}_{i}\right)}$$

(64)

for any *m, l*, and *i*( = *j*). Given *l* = 1, we have a 2*N _{d}* × 1 vector and by stacking these up for

$$\frac{\partial \stackrel{\u2012}{y}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})}{\partial {s}_{1j}}=-\left[\begin{array}{c}\hfill {\text{diag}}_{i=1}^{{N}_{d}}\left\{({\stackrel{\u2012}{y}}_{1i}\left({\mathit{s}}_{i}\right)-{r}_{1i})\frac{{w}_{1i,1}\left({\mathit{s}}_{i}\right)}{{v}_{1i}\left({\mathit{s}}_{i}\right)}\right\}{e}^{j}\hfill \\ \hfill {\text{diag}}_{i=1}^{{N}_{d}}\left\{({\stackrel{\u2012}{y}}_{2i}\left({\mathit{s}}_{i}\right)-{r}_{2i})\frac{{w}_{2i,1}\left({\mathit{s}}_{i}\right)}{{v}_{2i}\left({\mathit{s}}_{i}\right)}\right\}{e}^{j}\hfill \end{array}\right]$$

(65)

where *e ^{j}* denotes a

Second, we consider two *N _{d}* × 2

(66)

where $\Psi ({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2},y)$ denotes the PL cost function expressed in terms of three variables. The Hessians are defined in terms of these three variables. For example, ${\nabla}^{200}\Psi $ is the second-order derivative of Ψ with respect to ${\stackrel{\u02c7}{\mathit{s}}}_{1}$ and ${\nabla}^{020}\Psi $ is the second-order derivative of Ψ with respect to ${\stackrel{\u02c7}{\mathit{s}}}_{2}$. ${\nabla}^{110}\Psi $ is the Hessian of Ψ in terms of ${\stackrel{\u02c7}{\mathit{s}}}_{1}$ and ${\stackrel{\u02c7}{\mathit{s}}}_{2}$. By expressing Ψ in terms of the negative Poisson log-likelihood function and penalty function, and evaluating (66) at $\stackrel{\u2012}{y}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$, we now have to necessary components for the LIR.

Finally, combining these two pieces in (65) and (66) yields the following expression for ${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})$:

$${\ell}_{1}^{j}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})=-{\left[\begin{array}{cc}\hfill {\nabla}^{200}L+{\gamma}_{1}{\mathbf{C}}^{\mathbf{T}}\mathbf{C}\hfill & \hfill {\nabla}^{110}L\hfill \\ \hfill {\nabla}^{110}L\hfill & \hfill {\nabla}^{020}L+{\gamma}_{2}{\mathbf{C}}^{\mathbf{T}}\mathbf{C}\hfill \end{array}\right]}^{-1}\times \left[\begin{array}{c}\hfill {\nabla}^{101}L\frac{\partial \stackrel{\u2012}{y}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})}{\partial {s}_{1j}}\hfill \\ \hfill {\nabla}^{011}L\frac{\partial \stackrel{\u2012}{y}({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2})}{\partial {s}_{1j}}\hfill \end{array}\right]$$

where *L* denotes the negative Poisson log-likelihood function expressed in terms of $({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2},y)$ and all Hessians of *L* are evaluated at ${\stackrel{\u02c7}{\mathit{s}}}_{1}={\widehat{\stackrel{\u02c7}{\mathit{s}}}}_{1}\left(\stackrel{\u2012}{y}\right),{\stackrel{\u02c7}{\mathit{s}}}_{2}={\widehat{\stackrel{\u02c7}{\mathit{s}}}}_{2}\left(\stackrel{\u2012}{y}\right)$, and *y = y̵*. Evaluating all necessary Hessians of *L* by chain rules allows the LIR in (22). The following three partial derivatives are useful to obtain this final expression for the LIR:

$$\frac{\partial L({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2},y)}{\partial {s}_{li}}=\sum _{m=1}^{2}\left(\frac{{y}_{mi}}{{\stackrel{\u2012}{y}}_{mi}}-1\right)\frac{\partial {\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{li}}$$

(67)

$$\frac{{\partial}^{2}L({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2},y)}{\partial {s}_{li}\partial {y}_{mi}}=\frac{1}{{\stackrel{\u2012}{y}}_{mi}}\frac{\partial {\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{li}}$$

(68)

$$\frac{{\partial}^{2}L({\stackrel{\u02c7}{\mathit{s}}}_{1},{\stackrel{\u02c7}{\mathit{s}}}_{2},y)}{\partial {s}_{li}\partial {s}_{{l}^{\prime}i}}=-\sum _{m=1}^{2}\frac{{y}_{mi}}{{\stackrel{\u2012}{y}}_{mi}^{2}}\frac{\partial {\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{li}}\frac{\partial {\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{{l}^{\prime}i}}+\sum _{m=1}^{2}\left(\frac{{y}_{mi}}{{\stackrel{\u2012}{y}}_{mi}}-1\right)\frac{{\partial}^{2}{\stackrel{\u2012}{y}}_{mi}({s}_{1i},{s}_{2i})}{\partial {s}_{li}\partial {s}_{{l}^{\prime}i}}$$

(69)

where *l,ĺ* = 1,2, *i* = 1,...,*N _{d}*, and

The material in this paper was presented in part at SPIE Medical Imaging 2008: Physics of Medical Imaging.

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

^{1}We call the conventional CTACs single-energy computed tomography (SECT) based methods to distinguish them from DECT based methods.

^{2}For simplicity, we focus on the 2-D static cases in this paper. The developed techniques can be extended to the helical and cone-beam cases.

^{3}These are also known as majorization conditions [40], [41].

^{4}The off-diagonal blocks are neglected in the approximation used in (31).

^{5}For simplicity, we denote a nonnegative definite matrix **A** as **A**
0.

Joonki Noh, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA (Email: ude.hcimu@knoojhon).

Jeffrey A. Fessler, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA (Email: ude.hcimu.scee@relssef).

Paul E. Kinahan, Department of Radiology, University of Washington, Seattle, WA 98195 USA (Email: ude.notgnihsaw.u@nahanik).

1. Beyer T, Townsend DW, Brun T, Kinahan PE, Charron M, Roddy R, Jerin J, Young J, Byars L, Nutt R. A combined PET/CT scanner for clinical oncology. J. Nuc. Med. 2000 Aug.41(8):1369–1379. [PubMed]

2. Kinahan PE, Townsend DW, Beyer T, Sashin D. Attenuation correction for a combined 3d PET/CT scannter. Med. Phys. 1998 Oct.25(10):2046–2053. [PubMed]

3. Kinahan P, Fessler JA, Alessio A, Lewellen TK. Quantitative attenuation correction for PET/CT using iterative reconstruction of low-dose dual-energy CT. Proc. IEEE Nuc. Sci. Symp. Med. Im. Conf. 2004;5:3285–3289.

4. Kinahan PE, Hasegawa BH, Beyer T. X-ray based attenuation correction for PET/CT scanners. Seminars in Nuclear Medicine. 2003;33(3):166–179. [PubMed]

5. Robinson PJ, Kreel L. Pulmonary tissue attenuation with computed tomography: Comparison of inspiration and expiration scans. J. Comp. Assist. Tomogr. 1979;3:740–748. [PubMed]

6. Blankespoor SC, Wu X, Kalki K, Brown JK, Cann CE, Hasegawa BH. Attenuation correction of SPECT using X-ray CT on an emission-transmission CT systems: Myocardial perfusion assessment. IEEE Trans. Nucl. Sci. 1996 Aug.43(4):2263–2274.

7. Kinahan PE, Alessio AM, Fessler JA. Dual energy CT attenuation correction methods for quantitative assessment of response to cancer therapy with PET/CT imaging. Technol. Cancer Res. Treatment. 2006 Aug.5(4):319–328. [PubMed]

8. Alvarez RE, Macovski A. Energy-selective reconstructions in X-ray computed tomography. Phys. Med. Biol. 1976 Sep.21(5):733–744. [PubMed]

9. Gingold EL, Hasegawa BH. Systematic bias in basis material decomposition applied to quantitative dual-energy x-ray imaging. Med. Phys. 1992 Jan.19(1):25–33. [PubMed]

10. Hasegawa BH, Lang TF, Brown JK, Gingold EL, Reilly SM, Blankespoor SC, Liew SC, Tsui BMW, Ramanathan C. Object-specific attenuation correction of SPECT with correlated dual-energy X-ray CT. IEEE Trans. Nucl. Sci. 1993 Aug.40(4):1242–1252.

11. Guy MJ, Castellano-Smith IA, Flower MA, Flux GD, Ott RJ, Visvikis D. DETECT-dual energy transmission estimation CT-for improved attenuation correction in SPECT and PET. IEEE Trans. Nucl. Sci. 1998 Jun.45(3):1261–1267.

12. Macovski A, Alvarez RE, Chan J, Stonestrom JP, Zatz LM. Energy dependent reconstruction in X-ray computerized tomography. Comput. Biol. Med. 1976 Oct.6(4):325–336. [PubMed]

13. Marshall WH, Alvarez RE, Macovski A. Initial results with prereconstruction dual-energy computed tomography (PREDECT) Radiology. 1981 Aug.140(2):421–430. [PubMed]

14. Stonestrom JP, Alvarez RE, Macovski A. A framework for spectral artifact corrections in X-ray CT. IEEE Trans. Biomed. Eng. 1981 Feb.28(2):128–141. [PubMed]

15. Lehmann LA, Alvarez RE. Energy-selective radiography: A review. In: Kareiakes COJ, Thomas S, editors. Digital Radiography: Selected Topics. Plenum; New York: 1986. pp. 145–188.

16. Kotzki PO, Mariano-Goulart D, Rossi M. Prototype of dual energy X-ray tomodensimeter for lumbar spine bone mineral density measurements; Choice of the reconstruction algorithm and first experimental results. Phys. Med. Biol. 1992 Dec.37(12):2253–2265. [PubMed]

17. Michael GJ. Tissue analysis using dual energy CT. Australasian Phys. Eng. Sci. Med. 1992 Mar.15(1):25–37. [PubMed]

18. Markham C, Fryar J. Element specific imaging in computerised tomography using a tube source of X-rays and a low energy-resolution detector system. Nucl. Instrum. Meth. Phys. Res. A. 1993 Jan.324(1–2):383–388.

19. Clinthorne NH. A constrained dual-energy reconstruction method for material-selective transmission tomography. Nucl. Instrum. Meth. Phys. Res. A. 1994 Dec.353(1):347–348.

20. Sukovic P, Clinthorne NH. Data weighted vs. non-data weighted dual energy reconstructions for X-ray tomography. Proc. IEEE Nucl. Sci. Symp. Med. Im. Conf. 1998;3:1481–1483.

21. Sukovic P, Clinthorne NH. Design of an experimental system for dual energy X-ray CT. Proc. IEEE Nucl. Sci. Symp. Med. Im. Conf. 1999;2:1021–1022.

22. Sukovic P, Clinthorne NH. Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography. IEEE Trans. Med. Imag. 2000 Nov.19(11):1075–1081. [PubMed]

23. Fessler JA, Elbakri I, Sukovic P, Clinthorne NH. Maximum-likelihood dual-energy tomographic image reconstruction. Proc. SPIE Med. Imag. Image Proc. 2002;1:38–49.

24. Xu J, Frey EC, Taguchi K, Tsui BMW. A Poisson likelihood iterative reconstruction algorithm for material decomposition in CT. Proc. SPIE Med. Imag. Phys. Med. Imag. 2007:65101Z.

25. O'sullivan JA, Benac J. Alternating minimization algorithms for transmission tomography. IEEE Trans. Med. Imag. 2007 Mar.26(3):283–297. [PubMed]

26. Noh J, Fessler JA, Kinahan PE. Low-dose dual-energy computed tomography for PET attenuation correction with statistical sinogram restoration. Proc. SPIE Med. Imag. Phys. Med. Imag. 2008:691312.

27. La Riviere PJ, Bian J, Vargas PA. Penalized-likelihood sino-gram restoration for computed tomography. IEEE Trans. Med. Imag. 2006 Aug.25(8):1022–1036. [PubMed]

28. Snyder DL, Helstrom CW, Lanterman AD, Faisal M, White RL. Compensation for readout noise in CCD images. J. Opt. Soc. Am. A. 1995 Feb.12(2):272–283.

29. Ruth C, Joseph PM. Estimation of a photon energy spectrum for a computed tomography scanner. Med. Phys., vol. 1997 May;24(5):695–702. [PubMed]

30. Colijn AP, Beekman FJ. Accelerated simulation of cone beam X-ray scatter projections. IEEE Trans. Med. Imag. 2004 May;23(5):584–590. [PubMed]

31. Elbakri IA, Fessler JA. Statistical image reconstruction for polyenergetic X-ray computed tomography. IEEE Trans. Med. Imag. 2002 Feb.21(2):89–99. [PubMed]

32. Sukovic P, Clinthorne NH. Basis material decomposition using triple-energy X-ray computed tomography. IEEE Instrum. Measurement Technol. Conf., Venice, Italy. 1999;3:1615–1618.

33. Cardinal HN, Fenster A. An accurate method for direct dual-energy calibration and decomposition. Med. Phys. 1990 May;17(3):327–341. [PubMed]

34. Cardinal HN, Fenster A. Analytic approximation of the log-signal and log-variance functions of x-ray imaging systems, with application to dual-energy imaging. Med. Phys. 1991 Sep.18(5):867–887. [PubMed]

35. Hsieh J. Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise. Med. Phys. 1998 Nov.25(11):2139–2147. [PubMed]

36. Fessler JA. Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography. IEEE Trans. Imag. Process. 1996 Mar.5(3):493–506. [PubMed]

37. Fessler JA. Hybrid Poisson/polynomial objective functions for tomographic image reconstruction from transmission scans. IEEE Trans. Imag. Process. 1995 Oct.4(10):1439–1450. [PubMed]

38. Daube-Witherspoon ME, Carson RE. Investigation of angular smoothing of PET data. IEEE Trans. Nucl. Sci. 1997 Dec.44(6-2):2494–2499.

39. Jacobson MW, Fessler JA. An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms. IEEE Trans. Imag. Procress. 2007 Oct.16(10):2411–2422. [PMC free article] [PubMed]

40. Lange K, Hunter DR, Yang I. Optimization transfer using surrogate objective functions. J. Computat. Graph. Stat. 2000 Mar.9(1):1–20.

41. Hunter DR, Lange K. A tutorial on MM algorithms. Amer. Stat. 2004 Feb.58(1):30–37.

42. Snyder DL, Hammoud AM, White RL. Image recovery from data acquired with a charge-coupled-device camera. J. Opt. Soc. Am. A. 1993 May;10(5):1014–1023. [PubMed]

43. Chatziioannou A, Dahlbom M. Detailed investigation of transmission and emission data smoothing protocols and their effects on emission images. IEEE Trans. Nucl. Sci. 1996 Feb.43(1):290–294.

44. Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs. IEEE Trans. Imag. Process. 1996 Sep.5(9):1346–1358. [PubMed]

45. Stayman JW, Fessler JA. Compensation for nonuniform resolution using penalized-likelihood reconstruction in space-variant imaging systems. IEEE Trans. Med. Imag. 2004 Mar.23(3):269–284. [PubMed]

46. Segars WP, Tsui BMW. Study of the efficacy of respiratory gating in myocardial SPECT using the new 4-D NCAT phantom. IEEE Trans. Nucl. Sci. 2002 Jun.49(3):675–679.

47. Whiting BR, Massoumzadeh P, Earl OA, O'sullivan JA, Snyder DL, Williamson JF. Properties of preprocessed sinogram data in x-ray computed tomography. Med. Phys. 2006 Sep.33(9):3290–3303. [PubMed]

48. Lasio GM, Whiting BR, Williamson JF. Statistical reconstruction for x-ray computed tomography using energy-integrating detectors. Phys. Med. Biol. 2007 Apr.52(8):2247–2266. [PubMed]

49. Kalender WA, Perman WH, Vetter JR, Klotz E. Evaluation of a prototype dual-energy computed tomographic apparatus. 1. Phantom studies. Med. Phys. 1986 May;13(3):334–339. [PubMed]

50. Joseph PM, Spital RD. The exponential edge-gradient effect in x-ray computed tomography. Phys. Med. Biol. 1981 May;26(3):473–487. [PubMed]

51. De Pierro AR. A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography. IEEE Trans. Med. Imag. 1995 Mar.14(1):132–137. [PubMed]

52. De Pierro AR. On the relation between the ISRA and the EM algorithm for positron emission tomography. IEEE Trans. Med. Imag. 1993 Jun.12(2):328–333. [PubMed]

53. Fessler JA. Statistical image reconstruction methods for transmission tomography. In: Sonka M, Fitzpatrick JM, editors. Handbook of Medical Imaging, Volume 2. Medical Image Processing and Analysis. SPIE; Bellingham: 2000. pp. 1–70.

54. Erdoğgan H, Fessler JA. Ordered subsets algorithms for transmission tomography. Phys. Med. Biol. 1999 Nov.44(11):2835–2851. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |