PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 November 1.
Published in final edited form as:
PMCID: PMC2895983
NIHMSID: NIHMS207164

Statistical Sinogram Restoration in Dual-Energy CT for PET Attenuation Correction

Joonki Noh, Jeffrey A. Fessler, Fellow, IEEE, and Paul E. Kinahan, Senior Member, IEEE

Abstract

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.

Keywords: Index Terms—Attenuation correction, dual-energy (DE) X-ray computed tomography (CT), low radiation dose, positron emission tomography/computed tomography (PET/CT), sinogram restoration

Introduction

A. Background

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.

B. Literature Review

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.

II. Physical Model Formulations

A. General Measurement Model

We consider a general measurement model where multiple sets of polychromatic measurements are collected for M0 different incident spectra and forward projections (line integrals) are recorded for Nd radius-angle pairs for each incident spectrum, forming M0 sinograms.2

For m = 1,..., M0 and i = 1,..., Nd let ymi denote the measurement for the mth incident spectrum and the ith ray. We assume that ymi is a random variable whose ensemble mean mi is defined by the following underlying physics:

E[ymi]=ymiImi(E)exp(Liμ(x,E)d)dE+rmi
(1)

where Imi(E) denotes the product of the mth incident source spectrum and the detector gain for the ith ray, and Lid is the line integral along the ith ray. μ(x,E) denotes the LAC of the object being scanned at the spatial location x and the photon energy E.rmi denotes additive background contributions, for example room background, dark current, and scatter. We also can use rmi to model electronic noise in a shifted Poisson approach. We treat Imi(E) and rmi as known (or separately calibrated) nonnegative quantities by methods proposed in [28]–[30]. By modeling the polychromatic source spectra in (1), the sinogram restoration methods in Section III can correct for beam hardening artifacts. Therefore separate beam hardening correction steps are not needed [23], [31]. We refer to {ymi}i=1Nd as the sinogram measurements for the mth incident source spectrum. The LAC μ(x,E) is the property we want to estimate in a CT scan, but that for PET imaging it is a confounding aspect (albeit monochromatically at E=511keV) that we want to remove [4].

B. Basis Material Decomposition

Since the number of measurements is finite whereas the LAC of the scanned object is a continuous function of x and E, we parameterize μ(x,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:

μ(x,E)=l=1L0βl(E)ρl(x)
(2)

where βl(E) denotes mass attenuation coefficient (MAC) and ρl(x) is the unknown density map of the lth material type. L0 is the number of material types comprising the object being scanned. In DECT, we usually have L0 = 2, e.g., soft tissues and bone minerals. Other material decompositions are possible, e.g., [8], and could be used in our approaches. The goal in DECT is to estimate {ρl(x)}l=12 from {ymi}i=1Nd for M0 = 2 incident spectra.

C. Measurement Model Reformulation

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

ymi=Imiefmi(si)+rmi
(3)

for m = 1,...,M0 and i = 1,...,Nd, where

fmi(si)logvmi(si),
(4)
vmi(si)pmi(E)el=1L0βl(E)slidE.
(5)

The nonlinear function fmi(si) characterizes the beam hardening caused by polychromatic source spectra and it can be measured using calibration phantoms [33], [34]. We define the total intensity Imi for the mth incident spectrum and the ith ray, and the sinogram vector si as follows:

ImiImi(E)dE,pmi(E)Imi(E)Imi,
(6)
si(s1i,,sL0i)T,sliLiρl(x)d.
(7)

The nonlinear function fmi(si) is monotonically increasing and concave. These two properties play key roles in the development of our sinogram restoration algorithms as shown in Appendix I.

III. Component Sinogram Restoration

Usually in DECT, we reconstruct the component density images, and ρ1(x) and ρ2(x) 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 {si}i=1Nd. Therefore we focus on estimating {si}i=1Nd hereafter.

We discuss methods to recover component sinograms, i.e., {sli}i=1Nd for L0 material types from noisy measurements, i.e., {ymi}i=1Nd for M0 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.

A. Conventional Sinogram Decomposition

Given the noisy measurement ymi, the conventional method for estimating values of the nonlinear function fmi is to invert (3) as follows:

f^milog{smooth(ymirmiImi)}
(8)

where smoothing in the radial (detector) direction is often applied to reduce noise [35]. Equating (4) and (8) yields a system of M0 nonlinear equations and L0 unknowns for the ith ray, where the lth unknown variable is sli. Since we usually have the same number of source spectra and material types, that is M0 = L0, solving this nonlinear system of equations produces the following estimate of component sinograms:

s^ifi1(f^i),i=1,,Nd
(9)

where fi(f1i,,fM0i)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 {si}i=1Nd by minimizing the following function:

Λ(s)i=1Nd(s^isi)TQi(s^isi)+R(s)
(10)

where a possible choice of the weight matrix for the ith ray is Qi=cov1(s^i). Its approximation can be found by the method proposed in [36]. The L0 × Nd sinogram matrix s is obtained by concatenating si long the ray index and R(s) is a roughness penalty function for the component sinograms. Although (10) could reduce the noise somewhat, the main source of noise amplification is the inverse performed in (9). We present next two alternative approaches that consider noise characteristics and avoid the inverse in (9), thus providing better estimates of component sinograms.

B. Penalized Weighted Least Square

Instead of solving the M0 nonlinear equations, we propose to estimate component sinograms from the f^mi values in (8) by minimizing a PWLS cost function. Subject to the nonnegativity constraint on the entries of si, we have

s^PWLS=arg minsRL0×NdΦ(s)
(11)

and

Φ(s)i=1Nd12(f^ifi(si))TDi(f^ifi(si))+R(s)
(12)

where s{si}i=1Nd denotes the L0 × Nd sinogram matrix. For the ith ray, Didiagm=1M0{ymi} is the M0 × M0 weight matrix. If the measurements approximately follow Poisson distribution and rmi is small, then the above Di is a reasonable choice since an approximate variance of f^mi is [36], [37]

var(f^mi)var(ymi)(ymirmi)21ymi.
(13)

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

R(s)l=1L0γlk=1K12([Csˇl]k)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 sˇl denotes the lth component sinogram, i.e., sˇl(sl1,,slNd)T. We choose C to regularize component sinograms only in the radial direction [38]; thus KNd.

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 sli for l = 1,...,L0 and i = 1,...,L0. After applying OTP, we arrive at the following equation for updating an estimate of component sino-grams in the nth step:

sli(n+1)=[sli(n)1hliΦ(sli(n))sli]+
(15)

where [a]+(max(a,0)) enforces the nonnegativity constraint on sli, and sli(n) denotes the estimate of sli in the nth iteration. The curvature hli is defined in (44) of Appendix I-A. We precompute it before iterating. It can be shown that the associated surrogate function satisfies the surrogate conditions. Thus the update provided by (15) decreases the PWLS cost function every iteration. Appendix I-A gives the detailed derivations.

C. Penalized Likelihood

Although the PWLS approach described in the previous section is statistically motivated, it requires the logarithmic transformation in (8) to obtain 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., {ymi}m=1,i=1M0,Nd.

For simplicity, we assume that the measurement ymi obeys a Poisson distribution with the ensemble mean in (1), i.e.,

ymi~Poisson(ymi(si)).
(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:

Ψ(s)L(s)+R(s)=m=1M0i=1NdLmi(si)+R(s)
(17)

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

Lmi(si)ymi(si)ymilogymi(si)
(18)

where constants independent of si are ignored and gmi(x)xymilogx 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:

s^PL=arg minsRL0×NdΨ(s)
(19)

subject to the nonnegativity of the entries of si. We exploit OTP to achieve a sequence of SQSs for 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 nth step. For any l and i, we have

sli(n+1)=[sli(n)1dliΨ(sli(n))sli]+.
(20)

The curvature dli 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 dli precomputable. However, if the approximations used are reasonably accurate or if dli provides a sufficiently large curvature that ensures the surrogate conditions, we can expect monotonicity of Ψ(s) when we implement (20). This can be empirically checked by evaluating Ψ(s) in each iteration for given data. In Section V, the simulations corresponding to the PL method exhibited empirical monotonicity with the iteration in (20).

After estimating all component sinograms for L0 material types, i.e., {s^li}i=1Nd for l = 1,..., if desired, one can use the sinogram estimates to reconstruct L0 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 ρl(x). 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 {si}i=1Nd to compute ACFs, enhancing the quality of the PET emission images.

IV. Regularization Design

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.

A. Local Impulse Response of Component Sinograms

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 L0 = M0 = 2:

equation image
(21)

where lj(sˇ1,sˇ2) describes the responses that appear on all component sinograms when we place an impulse at the jth element of the lth component sinogram and is a lexicographically ordered 2Nd × 1 vector. y(sˇ1,sˇ2) denotes a 2Nd × 1 vector containing lexicographically ordered ymi(si) for m = 1,2 and i = 1,...,Nd.y is similarly defined with ymi. An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0002.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0003.jpg denote column gradient operator and row gradient operator with respect to y, respectively. Here, sˇl denotes the th component sinogram and si is the th column of the sinogram matrix defined as s[s1,,sNd]

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:

1j(sˇ1,sˇ2)=[F(sˇ1,sˇ2)+R(sˇ1,sˇ2)]1F(sˇ1,sˇ2)e1j
(22)

where e1j denotes a 2Nd × 1 unit vector containing 1 at the position corresponding to the jth element of the first component sinogram where 1 ≤ jNd. By replacing e1j with e2j that is similarly defined, we can also obtain 2j(sˇ1,sˇ2). The Fisher information matrix (FIM) in the sinogram domain and the penalty matrix in (22) have the following forms:

F(sˇ1,sˇ2)[D1(sˇ1,sˇ2)D12(sˇ1,sˇ2)D12(sˇ1,sˇ2)D2(sˇ1,sˇ2)]
(23)
R(sˇ1,sˇ2)[γ1CTC00γ2CTC]
(24)

where C is from (14) and

D1(sˇ1,sˇ2)=diagi=1Nd{m=12(ymirmi)2ymi(wmi,1vmi)2}D12(sˇ1,sˇ2)=diagi=1Nd{m=12(ymirmi)2ymiwmi,1wmi,2vmi2}D2(sˇ1,sˇ2)=diagi=1Nd{m=12(ymirmi)2ymi(wmi,2vmi)2}

and for l = 1,2,

wmi,l(si)pmi(E)βl(E)eβ(E)TsidE.

For given E,β(E)(β1(E),,βL0(E))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 D1(sˇ1,sˇ2) and D2(sˇ1,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, D12(sˇ1,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 fmi(s1i,s2i). 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.

B. Spatially Variant Penalty Design

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 jth element as follows:

F(sˇ1,sˇ2)[[D1]jINd[D12]jINd[D12]jINd[D2]jINd]SjINd
(25)

where, for example, [D1]j denotes the jth diagonal entry of D1(sˇ1,sˇ2),INd is an Nd × Nd identity matrix, [multiply sign in circle] denotes a Kronecker product, and the 2 × 2 matrix Sj is defined as

Sj[[D1]j[D12]j[D12]j[D2]j].
(26)

Substituting (25) into (22) and simplifying yield

1j(sˇ1,sˇ2)(Sj12INd)[I2Nd+(Sj12INd)R×(Sj12INd)]1(Sj12INd)e1j
(27)

where Sj12 denotes a square root factorization of the matrix Sj 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

R~[γ1C~1TC~100γ2C~2TC~2]
(28)

where C~1 and 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 C~1 and C~2 as follows:

C~1CD~1,D~1diagi=1Nd{[D1]i}
(29)
C~2CD~2,D~2diagi=1Nd{[D2]i}.
(30)

By replacing the conventional penalty matrix R with the new penalty matrix R~, we have a useful approximation4 to obtain 1j(sˇ1,sˇ2)

(Sj12INd)R~(Sj12INd)ΓCTC
(31)

where the regularizing parameter matrix is defined as

Γ[γ100γ2].
(32)

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

1j(sˇ1,sˇ2)[I2Nd+ΓCTC]1e1j
(33)

indicating approximately uniform resolution of component sinogram estimates since CT 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].

V. Simulation Studies

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.

A. DECT Based Attenuation Corrections

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 cm2. 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 E1=57keV and E2=72keV denote their effective energies in (61), respectively. For simplicity, we used the same incident spectra for each ray, i.e., Imi = Im, which ignores the effects of bow tie filters. Since clinical SECT scans have in the order of 106 photons per ray [47], [48], to simulate DECT with low radiation dose, we set the number of incident photons per ray for I1(E) and I2(E) to be 2.8 × 104 and 2 × 105, respectively. This very low value for I1(E) helps contain the extra dose of the lower energy scan. Using the true component densities in Fig. 1, we synthesized DECT measurements using Poisson random variables whose ensemble means follow the physical measurement model in (1) in a parallel-beam geometry. After generating measurements with a high spatial resolution in the projection domain, we downsampled them to make a typical resolution in PET/CT sinogram, producing 256 (radial direction) × 200 (angular direction) samples with 0.2 cm radial spacing.

Fig. 1
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.
Fig. 2
Two incident spectra Im(E) against energy E [keV] for m = 1,2: 80 kVp (top) and 140 kVp (bottom). The dashed vertical lines indicate the effective energies, E1 and E2, respectively.

To provide approximately uniform spatial resolution of the component sinogram estimates, we used the proposed penalty matrix 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 ith ray

ACFiexp(l=1L0βl(E)s^li)E=511keV.
(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].

Fig. 3
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 ...
Fig. 4
Restored component CT sinograms (top) and the corresponding component CT images (bottom) by DECT-PL: (a) restored soft tissue CT sinogram, (b) restored bone CT sinogram, (c) reconstructed soft tissue CT image by FBP, and (d) reconstructed bone CT image ...

B. Comparison With Bilinear Scaling

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.

Fig. 5
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 I2(E). 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 × 105 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.

Fig. 6
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, ...
TABLE I
Global NRMSEs of the Reconstructed PET Images by Six CTACs When Iodine Contrast Agents are Present

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.

Fig. 7
Four selected regions for local NRMSE analysis marked by white boxes. The region numbers were counted clockwise from Region 1 to 4.
TABLE II
Local NRMSEs of the Reconstructed PET Images When Iodine Contrast Agents are Present

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.

TABLE III
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.

Fig. 8
Horizontal profiles of the reconstructed PET emission images corrected by DECT-PL and SECT-PL-BS at two different vertical locations: (a) the anterior part of phantom's heart and (b) the posterior part of phantom's heart. The colored arrows correspond ...

VI. Conclusion and Discussion

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 fmi(si) 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.

Acknowledgment

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.

Appendix I

Algorithm Derivations

This appendix derives (15) and (20). To obtain surrogates for the PWLS cost function Φ(s) and PL cost functions Ψ(s), 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.

A. PWLS Algorithm

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

Φd(s)m=1M0i=1Ndϕmi(si)ymi
(35)

where ϕmi(si)12(f^mifmi(si))2 is the LS term for the mth spectrum and ith ray. By a Taylor series expansion of ϕmi(si about si(n), we have an inequality

equation image
(36)

where Wmi is a L0 × L0 positive definite matrix satisfying the condition that An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0012.jpg for any si in the matrix sense5 and 2 is a Hessian matrix. Letϕmi(n)(si) denote the right-hand side of (36). Since ϕmi(n)(si) then majorizes ϕmi(si, we define a surrogate function for Φ(S) in the th step, given by

Φd(n)(s)m=1M0i=1Ndϕmi(n)(si)ymi
(37)

where the exact form of Wmi 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

R(n)(s)l=1L0i=1Ndγlrli(n)(sli)rli(n)(sli)k=1Kaki2(ckiaki(slisli(n))+[Csˇl(n)]k)2
(38)

where denotes the element of C at (k, i) and aki is given by akickii=1Ndcki if cki0;aki0 if cki = 0. Then we have R(n)(s) that majorizes R(s).

By adding two surrogate functions Φd(n)(s) and R(n)(s), we have a SQS for the PWLS cost function defined as follows:

Φ(n)(s)Φd(n)(s)+R(n)(s)
(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 sli, equating it to zero, and simplifying, one arrives at the PWLS algorithm in (15).

Finally, we need to determine Wmi. From the definition of ϕmi(si, its Hessian is given by

equation image
(40)

Since fmi(si) is nonnegative, monotonically increasing, and concave, we have an inequality

equation image
(41)

We define An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0017.jpg and further simplify the above inequality into

equation image
(42)

where simaxarg maxsi2fmi(si). For our spectra in Fig. 2, we found that simax=0 for any i, but we could not prove that this will always hold. By defining Wmi from the right side of the inequality above, we have the first piece of the curvature hli as follows:

wmi(bmi2+[f^mi]+2fmi(0))ymi
(43)

producing

hlim=1M0wmi+γlc~i
(44)

where c~ik=1Kckii=1Ndcki and cki denotes the element of C at (k, i). Now we complete the derivation of the PWLS algorithm in (15).

B. PL Algorithm

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

L(s)=m=1M0i=1NdLmi(si)m=1M0i=1Ndgmi(ymi(si))
(45)

where gmi(x) = x — ymi log x is a convex function with respect to x. By combining (1) and (2), we have an expression for the ensemble mean of measurements

ymi(si)=Imi(E)eβT(E)sidE+rmi=[Imi(E)bmi(si(n),E)]tmi(si,E)bmi×(si(n),E)dE
(46)

where

tmi(si,E)eβT(E)si+rmiImi(E),bmi(si,E)ymi(si)tmi(siE).

Note that the term inside the square bracket behaves as a probability density function (PDF) since it is nonnegative and integrates into unity over E. By Jensen's inequality and the convexity of gmi(x), we have an inequality

gmi(ymi(si))Imi(E)gmi(tmi(si,E)bmi(si(n),E))bmi(si(n),E)dE
(47)

on the basis of De Pierro's multiplicative convexity trick [52]. Letlmi,1(n)(si) denote the right-hand side of (47). Since lmi,1(n)(si) majorizes gmi(mi(si)), we can define a surrogate function for –L(s) in the nth step, given by

L1(n)(s)m=1M0i=1Ndlmi,1(n)(si).
(48)

However, it is difficult to minimize L1(n)(s) directly.

Second, we want to seek a quadratic surrogate for by L1(n)(s) by applying the same idea as (36). To do that, we focus on

hmi(a,b,E)gmi(tmi(si,E)bmi(si(n),E))
(49)

where two variables are defined as aβT(E)si and bsi(n). By a Taylor series expansion of hmi(a,b,E) with respect to a, we have an inequality,

hmi(a,b,E)hmi(βT(E)si(n),b,E)+h.mi(βT(E)si(n),b,E)×(aβT(E)si(n))+12cˇmi(aβT(E)si(n))2
(50)

where cˇmi is a positive constant satisfying the condition that cˇmih¨mi(βT(E)si,b,E) for any si. Here h.mi(a,b,E) and h¨mi(a,b,E) denote the first- and second-order derivatives of hmi(a,b,E) with respect to a, respectively. Lethmi(n)(a,b,E) denote the right-hand side of (50). hmi(n)(a,b,E) then majorizes hmi(a,b,E). Thus we define a surrogate function for L1(n)(s) in the nth step, given by

L2(n)m=1M0i=1Ndlmi,2(n)(si),lmi,2(n)Imi(E)bmi(si(n),E)hmi(n)(βT(E)si,si(n),E)dE
(51)

where the exact form of cˇmi is discussed below. Since L2(n)(s) is quadratic with respect to si, we have many algorithms to minimize it, for example, the coordinate descent algorithm. An alternative choice is finding a SQS function for L2(n)(s) to further simplify the optimization.

Thirdly, we seek a SQS for L2(n)(s) by applying De Pierro's additive convexity trick [51]. In this case, hmi(n)(a,b,E) plays a key role to derive the SQS function. We express βT(E)si in the following way:

βT(E)si=l=1L0πl(E)(βl(E)sliβl(E)sli(n)πl(E)+βT(E)si(n))
(52)

where the coefficient πl(E) can be defined as

πl(E)βl(E)l=1L0βl(E)0,l=1L0πl(E)=1.
(53)

Since hmi(n)(a,b,E) has a quadratic form with respect to a and πl(E) acts as a probability mass function (PMF), by Jensen's inequality, we have

hmi(n)(βT(E)si,si(n),E)l=1L0πl(E)hmi(n)(βl(E)sliβl(E)sli(n)πl(E)+βT(E)si(n),si(n),E).
(54)

Letkmi(n)(si) denote the right-hand side of (54). Since kmi(n)(si) majorizes hmi(n)(a,b,E), we can define a surrogate function for L2(n)(s) in the th step, given by

L3(n)(s)m=1M0i=1Ndlmi,3(n)(si)lmi,3(n)(si)Imi(E)bmi(si(n),E)kmi(n)(si)dE.
(55)

Therefore, we have

L(s)L1(n)(s)L2(n)(s)L3(n)(s)
(56)

and L3(n)(s) is a SQS for –L(s) by the construction.

Fourthly, we combine the surrogate function for the data fidelity term, L3(n)(s) and the surrogate function for the penalty term,R(n)(s) in (38). Then we have a SQS for the PL cost function defined as follows:

Ψ(n)(s)L3(n)(s)+R(n)(s)
(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 sli, equating it to zero, and simplifying, one arrives at the same form as the PL algorithm in (20) but having

dli=m=1M0Imi(E)βl2(E)cˇmibmi(si(n),E)πl(E)dE+γlc~i.
(58)

Finally, we now discuss how to determine the curvature cˇmi. When measurements are Poisson distributed, the optimal 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 Lmi(si) varies slowly around minimizer, we have an approximation for the curvature as follows:

cˇmiymi
(59)

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

By substituting (59) into (58), the curvature dli is

dli=m=1M0βl(E)j=1L0βj(E)Imi(E)ymibmi(si(n),E)dE+γlc~i.
(60)

We replace E with the effective energy of the mth incident spectrum at the ith ray, defined as an expectation of energy E

EmiEImi(E)dEImi(E)dE.
(61)

This produces an approximate to the curvature dli as follows:

dlim=1M0βl(Emi)j=1L0βl(Emi)ymi+γlc~i.
(62)

Now we have the PL algorithm in (20).

Appendix II

Derivation of Local Impulse Response

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

equation image
(63)

First from (3), it can be checked that a first-order partial derivative of ymi(si) with respect to sli is given by

ymi(s1i,s2i)sli=(ymi(si)rmi)wmi,l(si)vmi(si)
(64)

for any m, l, and i( = j). Given l = 1, we have a 2Nd × 1 vector and by stacking these up for m = 1,2 and i = 1,..., Nd

y(sˇ1,sˇ2)s1j=[diagi=1Nd{(y1i(si)r1i)w1i,1(si)v1i(si)}ejdiagi=1Nd{(y2i(si)r2i)w2i,1(si)v2i(si)}ej]
(65)

where ej denotes a Nd × 1 unit vector whose jth entry is 1.

Second, we consider two Nd × 2Nd matrices An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0020.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0021.jpg in an augmented matrix. Extending the trick used in [44, Sec. III], it can be shown that the Nd × Nd augmented matrix is expressed as follows. For any measurements y(≥ 0)

equation image
(66)

where Ψ(sˇ1,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, 200Ψ is the second-order derivative of Ψ with respect to sˇ1 and 020Ψ is the second-order derivative of Ψ with respect to sˇ2. 110Ψ is the Hessian of Ψ in terms of sˇ1 and sˇ2. By expressing Ψ in terms of the negative Poisson log-likelihood function and penalty function, and evaluating (66) at y(sˇ1,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 1j(sˇ1,sˇ2):

1j(sˇ1,sˇ2)=[200L+γ1CTC110L110L020L+γ2CTC]1×[101Ly(sˇ1,sˇ2)s1j011Ly(sˇ1,sˇ2)s1j]

where L denotes the negative Poisson log-likelihood function expressed in terms of (sˇ1,sˇ2,y) and all Hessians of L are evaluated at sˇ1=sˇ^1(y),sˇ2=sˇ^2(y), 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:

L(sˇ1,sˇ2,y)sli=m=12(ymiymi1)ymi(s1i,s2i)sli
(67)
2L(sˇ1,sˇ2,y)sliymi=1ymiymi(s1i,s2i)sli
(68)
2L(sˇ1,sˇ2,y)slisli=m=12ymiymi2ymi(s1i,s2i)sliymi(s1i,s2i)sli+m=12(ymiymi1)2ymi(s1i,s2i)slisli
(69)

where l,ĺ = 1,2, i = 1,...,Nd, and m = 1,2. To obtain 2j(sˇ1,sˇ2), one can follow the same procedure as above.

Footnotes

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.

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

2For 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.

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

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

5For simplicity, we denote a nonnegative definite matrix A as A An external file that holds a picture, illustration, etc.
Object name is nihms-207164-ig0013.jpg 0.

Contributor Information

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

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

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

REFERENCES

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]