|Home | About | Journals | Submit | Contact Us | Français|
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 . 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 –. 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 : 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 , 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 . 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 , .
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 . 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 –. 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 , –. For the same decomposition, by singular value decomposition (SVD),  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 – 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, – 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 –. For attenuation correction in PET, iterative reconstruction algorithms based on PWLS approaches from polychromatic measurement models were developed recently , . 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 .
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 . 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  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 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 y̵mi is defined by the following underlying physics:
where denotes the product of the mth incident source spectrum and the detector gain for the ith ray, and is the line integral along the ith ray. denotes the LAC of the object being scanned at the spatial location and the photon energy 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 and rmi as known (or separately calibrated) nonnegative quantities by methods proposed in –. 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 , . We refer to as the sinogram measurements for the mth incident source spectrum. The LAC is the property we want to estimate in a CT scan, but that for PET imaging it is a confounding aspect (albeit monochromatically at ) that we want to remove .
Since the number of measurements is finite whereas the LAC of the scanned object is a continuous function of and , we parameterize using a basis material decomposition. We model the LAC using a set of basis functions that are separable in the space and energy domains , ,  as follows:
where denotes mass attenuation coefficient (MAC) and 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., , and could be used in our approaches. The goal in DECT is to estimate from for M0 = 2 incident spectra.
for m = 1,...,M0 and i = 1,...,Nd, where
The nonlinear function fmi(si) characterizes the beam hardening caused by polychromatic source spectra and it can be measured using calibration phantoms , . We define the total intensity Imi for the mth incident spectrum and the ith ray, and the sinogram vector si as follows:
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.
Usually in DECT, we reconstruct the component density images, and and 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 . Therefore we focus on estimating hereafter.
We discuss methods to recover component sinograms, i.e., for L0 material types from noisy measurements, i.e., 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.
Given the noisy measurement ymi, the conventional method for estimating values of the nonlinear function fmi is to invert (3) as follows:
where smoothing in the radial (detector) direction is often applied to reduce noise . 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:
where . 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 by minimizing the following function:
where a possible choice of the weight matrix for the ith ray is . Its approximation can be found by the method proposed in . 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.
Instead of solving the M0 nonlinear equations, we propose to estimate component sinograms from the values in (8) by minimizing a PWLS cost function. Subject to the nonnegativity constraint on the entries of si, we have
where denotes the L0 × Nd sinogram matrix. For the ith ray, 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 is , 
For the roughness penalty function R(s) in (12), we use
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 denotes the lth component sinogram, i.e., . We choose C to regularize component sinograms only in the radial direction ; thus K ≈ Nd.
We use the optimization transfer principle (OTP)  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:
where enforces the nonnegativity constraint on sli, and 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.
Although the PWLS approach described in the previous section is statistically motivated, it requires the logarithmic transformation in (8) to obtain . 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 . To further improve the estimated component sinograms, we now propose a PL method that uses the raw measurements, i.e., .
For simplicity, we assume that the measurement ymi obeys a Poisson distribution with the ensemble mean in (1), i.e.,
where denotes the sinogram matrix. The negative log-likelihood is given by
where constants independent of si are ignored and 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:
subject to the nonnegativity of the entries of si. We exploit OTP to achieve a sequence of SQSs for . 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
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., 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 . This approach usually gives less noisy estimates of component images than the classical sinogram decomposition combined with FBP . In addition, it is less computationally expensive than fully iterative methods for reconstructing component images, e.g., , , . However, in this paper, we focus on using the estimated component sinograms 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  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  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  to our PL sinogram restoration problem for DECT, we define the following LIR, focusing on L0 = M0 = 2:
where 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. denotes a 2Nd × 1 vector containing lexicographically ordered for m = 1,2 and i = 1,...,Nd.y is similarly defined with ymi. and denote column gradient operator and row gradient operator with respect to y, respectively. Here, denotes the th component sinogram and si is the th column of the sinogram matrix defined as
where denotes a 2Nd × 1 unit vector containing 1 at the position corresponding to the jth element of the first component sinogram where 1 ≤ j ≤ Nd. By replacing with that is similarly defined, we can also obtain . The Fisher information matrix (FIM) in the sinogram domain and the penalty matrix in (22) have the following forms:
where C is from (14) and
and for l = 1,2,
From (22), we conclude that the LIR in (22) is shift variant since the FIM contains the block matrices and 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, 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 . 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 . 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:
where, for example, denotes the jth diagonal entry of is an Nd × Nd identity matrix, denotes a Kronecker product, and the 2 × 2 matrix Sj is defined as
where 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
where and 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 and as follows:
By replacing the conventional penalty matrix R with the new penalty matrix , we have a useful approximation4 to obtain
where the regularizing parameter matrix is defined as
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 .
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 . 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 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 and 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 , , to simulate DECT with low radiation dose, we set the number of incident photons per ray for and to be 2.8 × 104 and 2 × 105, respectively. This very low value for 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.
To provide approximately uniform spatial resolution of the component sinogram estimates, we used the proposed penalty matrix 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 . 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
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 .
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.
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 . 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.
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.
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.
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  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 , 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  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 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 and PL cost functions , 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
where is the LS term for the mth spectrum and ith ray. By a Taylor series expansion of ϕmi(si about , we have an inequality
where Wmi is a L0 × L0 positive definite matrix satisfying the condition that for any si in the matrix sense5 and is a Hessian matrix. denote the right-hand side of (36). Since then majorizes ϕmi(si, we define a surrogate function for Φ(S) in the th step, given by
where the exact form of Wmi is determined below.
where denotes the element of C at (k, i) and aki is given by if if cki = 0. Then we have R(n)(s) that majorizes R(s).
By adding two surrogate functions and R(n)(s), we have a SQS for the PWLS cost function defined as follows:
satisfying the surrogate conditions . 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
Since fmi(si) is nonnegative, monotonically increasing, and concave, we have an inequality
We define and further simplify the above inequality into
where . For our spectra in Fig. 2, we found that 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:
where and cki denotes the element of C at (k, i). Now we complete the derivation of the PWLS algorithm in (15).
First, the data fidelity term of Ψ(s) is expressed as
Note that the term inside the square bracket behaves as a probability density function (PDF) since it is nonnegative and integrates into unity over . By Jensen's inequality and the convexity of gmi(x), we have an inequality
on the basis of De Pierro's multiplicative convexity trick . denote the right-hand side of (47). Since majorizes gmi(y̵mi(si)), we can define a surrogate function for –L(s) in the nth step, given by
However, it is difficult to minimize directly.
Second, we want to seek a quadratic surrogate for by by applying the same idea as (36). To do that, we focus on
where two variables are defined as and . By a Taylor series expansion of with respect to a, we have an inequality,
where is a positive constant satisfying the condition that for any si. Here and denote the first- and second-order derivatives of with respect to a, respectively. denote the right-hand side of (50). then majorizes . Thus we define a surrogate function for in the nth step, given by
where the exact form of is discussed below. Since 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 to further simplify the optimization.
Thirdly, we seek a SQS for by applying De Pierro's additive convexity trick . In this case, plays a key role to derive the SQS function. We express in the following way:
where the coefficient can be defined as
Since has a quadratic form with respect to a and acts as a probability mass function (PMF), by Jensen's inequality, we have
denote the right-hand side of (54). Since majorizes , we can define a surrogate function for in the th step, given by
Therefore, we have
and is a SQS for –L(s) by the construction.
Fourthly, we combine the surrogate function for the data fidelity term, 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:
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
Finally, we now discuss how to determine the curvature . When measurements are Poisson distributed, the optimal can be found in , 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:
where a similar idea can be found in .
If we place an impulsive perturbation in the first sinogram, the LIR from (21) can be expressed as
First from (3), it can be checked that a first-order partial derivative of mi(si) with respect to sli is given by
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
where ej denotes a Nd × 1 unit vector whose jth entry is 1.
Second, we consider two Nd × 2Nd matrices and 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)
where denotes the PL cost function expressed in terms of three variables. The Hessians are defined in terms of these three variables. For example, is the second-order derivative of Ψ with respect to and is the second-order derivative of Ψ with respect to . is the Hessian of Ψ in terms of and . By expressing Ψ in terms of the negative Poisson log-likelihood function and penalty function, and evaluating (66) at , we now have to necessary components for the LIR.
where L denotes the negative Poisson log-likelihood function expressed in terms of and all Hessians of L are evaluated at , 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:
where l,ĺ = 1,2, i = 1,...,Nd, and m = 1,2. To obtain , one can follow the same procedure as above.
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.
4The off-diagonal blocks are neglected in the approximation used in (31).
5For 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).