Abstract
In a typical positron emission tomography/computed tomography (PET/CT) system, the attenuation correction is necessary for PET image reconstruction, which involves a transformation from the CT Hounsfield units (HU) to its linear attenuation coefficient (LAC) at 511 keV. This transformation is usually aided by an empirical bilinear function, followed by the forward projection of the transformed attenuation image. In this paper, we propose a direct method that calculates attenuation factors from CT projections, without using a reconstructed CT image. In this method, the human body is considered as a mixture of three distinct components: air, water and bone. Then, we estimate the proportions of these three components along each x-ray path and restore the attenuation factor at 511 keV with the known water and bone LACs. Our numerical results show that the proposed method produces as accurate estimation as the conventional HU mapping method.
Introduction
In quantitative positron emission tomography (PET), the attenuation correction is important to reconstruct a radioactive tracer distribution accurately [11]. For simple or small objects, object geometry and tissue approximation can be utilized to estimate linear attenuation coefficients (LACs). For example, in small animal PET, the object is assumed to be filled with water. In brain imaging, Siemens’ SMART Neuro AC algorithm generates an attenuation map (the so-called μ map) based on an uncorrected PET scan [27]. However, for large animal or human imaging, an accurate attenuation correction is necessary. In a PET/CT system, CT scanning provides the anatomical background for image registration and μ map for the attenuation correction. Compared to a nuclear transmission scan with an isotope at or near the energy of the nuclear tracer, the CT scan is faster and more accurate [13, 17].
Note that CT projections and PET projections are obtained under different energy levels (40–140 keV and 511 keV), the LAC obtained from CT projections should be transformed to that at 511 keV for PET image reconstruction. There are mainly three methods to implement this transformation: linear scaling, segmentation, and multi-linear scaling [13, 14]. In the first method, the entire CT image is multiplied by the ratio of the attenuation coefficients of water at 511 keV and the average CT energy. This method provides poor estimation for bone materials. The second method segments the CT image into different tissue types and replaces them with theoretical LACs at 511 keV. This segmentation-based method is not widely used since some material types are difficult to be distinguished [19]. The most widely used method is the multi-linear scaling method that utilizes thresholds to break the CT image into different tissue types, and applies different scaling factors. It has been shown that this method achieves reasonable results for low-Z materials. However, it suffers from significant biases for high-Z materials, like bone and contrast agents [1 , 19].
Dual-energy CT (DECT) provides an alternative approach to deal with this problem. DECT enables biological structures to be decomposed into two basis ‘materials’, either physical effects such as photoelectric and Compton effects, which are two different mechanisms in the attenuation process, or two basis materials such as water and bone [3]. Examples using DECT in the PET attenuation correction were reported for density/atomic number decomposition [9], Compton/photoelectric decomposition [24], and water-bone/water-iodine decomposition [14]. These decomposition methods lead to a reduction of contrast agent artifacts in PET images. DECT significantly decreases the severity of beam hardening or streak artifacts as compared to conventional CT. DECT can be implemented in several ways, such as dual-layer detection, dual-source scanning, or fast kVp-switching [12]. However, the dual-energy attenuation correction method generally introduces extra radiation exposure compared to the single-kVp scanning and with an additional instrumentation cost. Even though a recent paper [25] reported that DECT can achieve good PET image quality with the same radiation dose as a single CT scan, a good DECT noise suppression technique is required and the x-ray spectra should be optimized carefully.
In the literature, most CT-based attenuation correction methods focus on the transformation in the image domain except for [18] which decomposes CT projections into basis material projections with DECT. In this paper, we focus on the direct transformation in the projection domain with only single-kVp CT. In this direct method, the human body is considered as a mixture of three distinct components: air, water and bone. Their effective proportions are estimated along x-ray paths for restoration of the corresponding 511keV-LAC projections. In section 2, we review the currently used bilinear mapping transformation for PET/CT, and describe our method based on different mixture models. In the third section, numerical tests are reported. In the last section, we discuss the relevant issues and conclude the paper.
Theory and method
HU value transformation
In commercial PET/CT scanners, a single CT scan is performed to obtain a CT image and then the CT image is transformed to the attenuation map at 511 keV. Since the x-ray beam contains photons of different energies, the reconstructed CT image is an average of LACs over the energy spectrum of the x-ray beam, and represented in Hounsfield units (HU):
We propose a mixture model of three distinct materials: water, air and bone, to support the bilinear mapping function from HU to LAC511 (Equation 2). Materials that have lower LACs than water are considered as combinations of water and air, while those that have higher LACs are assumed to be mixtures of water and bone so that the LAC for a given material type is the linear interpolation between either water and air or water and bone. This is a physical interpretation of the bilinear function HU-LAC511.
Let μ
A
(E) , μ
w
(E) and μ
B
(E) denote the LACs of air, water and bone respectively at energy E. Based on our assumption, the LAC at position x is a linear function of μ
A
(E) , μ
w
(E) and μ
B
(E); that is,
Let S (E) and I
0 denote the x-ray spectrum and intensity respectively. They can be measured using a standard method such as [20]. Then, the detected x-ray signal at a certain detector position is expressed as
According to Equation (4),
If a region in a patient only contains materials whose attenuation coefficients are less than the water value, then Equation (5) becomes
Similarly, if the region only contains materials denser than water; e.g., the brain, we can discard the air term and solve the following optimization problem
Since Equations (9) and (10) are monotonic with respect to p
w
in each associated term, the solution can be easily found using a numerical method such as bisection search. In our study, we utilized the Matlab (The MathWorks, Natick, MA) function fminbnd() to solve these problems. Finally, the estimated weighting factors are used to construct the sinogram at 511keV:
Because the line integral paths in a CT scan are different from the LORs in a PET scan in most cases, an interpolation step in the sinogram domain is required to find attenuation correction factors associated with LORs in the PET scan.
In some parts of patient body, there exist materials that have LACs either higher or lower than the water value. Theoretically, a DECT scan instead of a conventional CT scan is required to solve Equation 3accurately. Since DECT is not always available, here we use a segmentation and interpolation method to estimate the weighting factors in Equation 3.
We would like to segment the sinogram data into two categories: regions with and without bone structures respectively. Unfortunately, bone segmentation in projection domain is not a commonly used technique in CT application. However, the lung and rib segmentation from projection data has been well developed in computer-aided diagnosis (CAD) [2 , 21]. This makes the bone segmentation in chest scanning possible.
First, bone and rib cage segmentation can be performed in the 2D projection for each angle, which is similar to what is done for CAD. Then, we can map the segmentation result to the sinogram domain and separate the sinogram into the areas with and without bone structures. The bone-free area could be treated as a mixture of air and water (note that in this case the weight for water could be greater than 1;that is, the extrapolation is allowed), and we can solve the problem for the weighting factors p
A
and p
W
according to Equation (9). After that, the leftover area that combines three ingredients: air, water and bone, can be resolved based on the knowledge obtained from the bone-free area. Specifically, the weighting factors of air can be interpolated from p
A
values in the bone-free region, and p
W
and p
B
are estimated by solving the following optimization problem
To test our algorithm, we performed several simulation studies with or without noise. The experimental settings are described as follows.
The x-ray spectrum (shown in Fig. 1) was generated with an open source software package called Spectrum GUI [28]. Using this software, we simulated the GE Maxiray 125 tube with 1 mm air filter in all the experimental settings. The highest x-ray energy was 140 keV. The x-ray spectrum was divided into channels with a step size 1 keV. Parallel-beam projections were synthesized from 180 equi-angular orientations from 0° to 180°, and 400 equally distributed detector cells per view with a total length 40 cm. In each energy channel, we simulated corresponding CT measurements. The detected photons at different energies were summed up to generate the final data. And we added Poisson noise into thedataset.
Three digital phantoms were designed to test the three mixture models. The materials used for the phantoms are listed in Table 1. To calculate the LACs, we used Jaroslaw Tuszynski’s open source software [29] which is based on NIST’s tables [30]. The LAC of air is 0. Moreover, we tested the water-bone model with human abdomen CT images.
Numerical results
HU-LAC511 transform
Before we perform tests with noisy data, we established a bilinear transform HU-LAC511 with noise-free data. Data used for image reconstruction were first corrected by a cubic function similar to that in [10] for suppression of beam hardening artifacts.
The fitted HU-LAC511 transform is expressed as
We first tested the water-bone model using the digital phantom #1 as shown in Fig. 2(a), which contains only materials of LACs higher than the water value. Table 2 indicates the exact contents of Phantom #1. Poisson noise was added to projection data. Fig. 3 shows the estimated sinogram at 511 keV via the proposed method based on the water-bone model in comparison with the ground truth which was obtained by projecting the true μ map, which is calculated from the phantom components, along the LORs. The maximum relative error between these two sinograms is 1.97% , and the PSNR (peak signal-to-noise ratio) across the sinogram region is 51.98 dB. Similarly, the digital phantom #2 with materials of LACs lower than the water value was designed to test the water-air model, as shown in Fig. 2(b) and Table 3. Figure 4(a) shows the estimated sinogram obtained via the water-air model. The maximum relative error is 1.71% , and the PSNR is 52.15 dB.
Air-water-bone model
To test the air-water-bone model, we used a lung phantom shown in Fig. 2(c). It consists of water, LN-450 (lung tissue substitute), and cortical bone. Assuming that the projection data involving bone were already known from a bone-segmentation procedure, we first calculated the air components in the bone-free areas. Then, we implemented a bilinear interpolation method to estimate the weighting factors of air in the bone areas. For each pixel in the bone area, we searched in horizontal and vertical directions for the nearest pixels in the bone-free domain. Then, p A at this pixel was a weighted sum of the neighboring p A values.
For the conventional HU to μ map transform, the image was reconstructed via filtered back-projection, and transformed to the μ map according to Equations (1) and (14). Then, the sinogram was obtained by projecting the μ map with the same projector. With our proposed method, we followed the same steps used in the aforementioned studies with water-air and water-bone models. Also, we reconstructed the μ map from the estimated sinogram. Figure 5 shows the estimated weighting factors of air, water and bone.
For comparison, we followed the traditional reconstruction and transformation method using Equation (14). Figure 6 displays the two sinograms and their differences from the truth. The PSNRs for the proposed method is 37.14 dB, while that for the reconstruction and transformation method is 23.19 dB. Figure 7 are the μ maps from the HU-LAC511 transformation and the estimated sinogram by proposed method. In Fig. 8, two profiles in the sinograms are compared respectively along the vertical and horizontal lines in Fig. 6.
Abdomen phantom
Furthermore, spectral CT datasets of the human abdomen were used to evaluate our method. We first calculated the water and bone basis material coefficients from a spectral CT datasets obtained with a GE Discovery CT750 CT scanner and used them to construct monochromatic images with available theoretical water and bone LACs. Then, these images were re-projected by the same projector as in section 2.5 with Poisson noise added.
The overall structure of one slice CT image can be seen in Fig. 10(c). Since the abdomen contains only small air regions, we chose the water-bone model for simplicity instead of the air-water-bone model, which definitely is more accurate but requires more work. Figure 9 displays the estimated sinogram with the water-bone model and the one obtained from the HU-LAC511 transform. The PSNRs of the two sinograms are 34.12 dB and 24.74 dB respectively. Figure 10 shows the PET attenuation map, which demonstrated a good matching up with the true μ map. Moreover, Fig. 11 shows two profile comparisons in the sinogram domain. The results demonstrated that the proposed method performs better than the traditional method.
Discussions and conclusion
In the first two experiments, we demonstrated the Water-Air and Water-Bone models with digital phantoms. The results showed a good matching between the proposed method and the ground truth. The small errors are mainly from the bilinear model, which cannot exactly represent the LACs of all the materials. The effect of this inaccuracy on the PET images will require evaluation of real components of the human body. Because the bilinear relationship between HU value and LAC511 has been successfully used in most commercial PET/CT systems, we believe the material decomposition model for our method is acceptable and would not have any significant adverse effect on real PET/CT imaging. Yet, real data experiments need to be performed as a follow-up study.
For the three-mixture model, theoretically, we cannot solve the weighting factors based on a single energy scan. However, it is possible to estimate air components in the bone areas from bone-free areas by linear interpolation. Though it’s not accurate, our simulation result showed that the errors were still acceptable (the PSNR was better than the conventional reconstruction and HU-LAC511 transformation method).
In experiments with the lung phantom and the abdomen spectral CT data, we compared our method with the current multi-linear scaling approach in commercial PET/CT system. Both methods gave good estimation on sinograms at 511 keV (see Figs. 6 and 9). It should be noted that in real PET/CT system the intersect paths in a CT scan and the lines of response in a PET scan are described in different coordinates. Hence, the μ map obtained from the HU-LAC511 transform needs to be re-projected through the LORs. For the estimated sinogram of our method, such kind of LOR resampling is also required before it is used for PET image reconstruction. However, one advantage of our method is that it doesn’t need to reconstruct the image and thus can suppress the propagation of noise in the reconstruction, transformation and re-projection steps. More importantly, our proposed direct method enables the LAC511 estimation from truncated CT projections, which makes interior PET/CT possible for dose-saving and other benefits. Carefully inspecting profiles in Figs. 8 and 11, we would see that the one obtained by the proposed method follows the ground truth much better than the HU-LAC511 transformation.
In the abdomen experiment, bone segmentation seems impossible and no reports have been published for such a technique. However, the water-bone model has generated a very good estimation despite the existence of some small air regions. This indicates the robustness of the proposed method. Generally speaking, the selection of different model should be based on the anatomy. For instance, we believe water-bone model should work well for head and abdomen scanning. Even though some parts of them contain air, e.g. noise airway, they can be ignored as demonstrated by this experiment. And for thoracic cavity scanning, the airways cannot be ignored and we must choose the air-water-bone model.
A potential application of the proposed approach is the attenuation correction in interior PET/CT. The concept of interior tomography has been extended from CT reconstruction to SPECT, MRI, and phase-contrast tomography [22, 23]. However, the idea of interior PET/CT is missing. In such a case, the CT data are truncated and we cannot get a full CT image of the patient (interior CT has been proved in [26], but the accuracy of the area outside of the ROI is not guaranteed. Thus, the projection calculated from the reconstructed image is not accurate.). The proposed direct recovering method can handle this problem.
It should be pointed out that our proposed method differs from that in [18]. Although both methods target the sinogram restoration in the projection domain, a DECT is required for material decomposition in [18], while our method utilizes a surface scan to provide a whole body contour, without any additional x-ray dose involved. Furthermore, an angular discrepancy between dual-energy projections could degrade the accuracy of the basis material decomposition in the projection domain.
In conclusion, we’ve developed a direct approach to restore the sinogram at 511 keV in PET/CT system based on three possible models: air-water, water-bone and water-air-bone mixture model. This method avoids the reconstruction of CT image and will benefit interior PET/CT. We demonstrated its accuracy and stability through numerical phantom studies. Further studies should include interpolation methods to suit the sinogram to PET LOR reconstruction and the effects of the mixture model’s inaccuracy on the PET images.
