|Home | About | Journals | Submit | Contact Us | Français|
Photon-counting computed tomography (PCCT) is an emerging imaging technique that enables multi-energy imaging with only a single scan acquisition. To enable multi-energy imaging, the detected photons corresponding to the full x-ray spectrum are divided into several subgroups of bin data that correspond to narrower energy windows. Consequently, noise in each energy bin increases compared to the full-spectrum data. This work proposes an iterative reconstruction algorithm for noise suppression in the narrower energy bins used in PCCT imaging. The algorithm is based on the framework of prior image constrained compressed sensing (PICCS) and is called spectral PICCS; it uses the full-spectrum image reconstructed using conventional filtered back-projection as the prior image. The spectral PICCS algorithm is implemented using a constrained optimization scheme with adaptive iterative step sizes such that only two tuning parameters are required in most cases. The algorithm was first evaluated using computer simulations, and then validated by both physical phantoms and in-vivo swine studies using a research PCCT system. Results from both computer-simulation and experimental studies showed substantial image noise reduction in narrow energy bins (43~73%) without sacrificing CT number accuracy or spatial resolution.
Dual-energy CT has become a valuable tool in diagnostic imaging (McCollough et al. 2015); it is used for virtual mono-energetic imaging, automated bone/plaque removal in CT angiography, virtual non-contrast-enhanced/non-calcium imaging, kidney stone characterization, and discrimination of gout from pseudo gout. Currently, dual-energy CT uses energy-integrating detectors (EIDs) and provides two sets of projection data associated with two different x-ray spectra. Dual-energy CT has been used for two basis-material decomposition (Alvarez et al. 1976; Heismann et al. 2003; Alvarez 2011) and three basis-material decomposition (Liu et al. 2009; Yu et al. 2009; Long et al. 2014; Mendonça et al. 2014).
Photon-counting CT (PCCT) is conceptually advantageous when compared to dual-energy CT (Atak et al. 2015; Faby et al. 2015). Such a technique uses a photon-counting detector (PCD) that is capable of resolving the energy information of each incident x-ray photon. Consequently, multiple sets of projection data with different x-ray spectra can be obtained from a single scan, which enables multi-basis material decomposition (Roessl et al. 2007; Schlomka et al. 2008; Kappler et al. 2013; Faby et al. 2015; Li et al. 2015). PCCT is not currently used in clinical practice due to limitations in PCD materials at high x-ray flux (Taguchi et al. 2013). Several research PCCT systems have been developed by various vendors and institutes (Schlomka et al. 2008; Shikhaliev 2008; Iwanczyk et al. 2009; Bennett et al. 2014), some of which are capable of handling high x-ray flux (Köhler 2004; Xu et al. 2013; Kappler et al. 2014; Gutjahr et al. 2016; Yu et al. 2016).
In particular, PCCT is able to produce narrow-energy-bin images with improved energy resolution. These improved-energy-resolution images are further processed to generate clinically useful information. For example, they may be used to generate mono-energetic images, iodine maps, or virtual non-contrast images, via multi-basis material decomposition; they may also be used to generate images of enhanced soft tissue contrast via various combination schemes. However, because the number of photons used in a narrow energy bin decreases, the associated image noise increases (Roessl et al. 2007; Leng et al. 2011; Alvarez 2013). The increased image noise induces substantial difficulty in subsequent processing. For example, increased image noise would result in a more ill-posed inversion process of a calibration matrix in image-based material decomposition. Therefore, reducing noise in narrow energy bins is an important task for PCCT imaging. To this end, traditional denoising methods may be applied either in the image domain (Vogel et al. 1996; Tomasi et al. 1998; Buades et al. 2005; Borsdorf et al. 2008; Li et al. 2014), or in the projection domain (Manduca et al. 2009; Tang et al. 2012). However, these image processing methods only work in one domain, and the amount of noise reduction may be limited if image details are to be preserved; oftentimes, it is necessary to make a tradeoff between noise reduction and spatial resolution.
Noise reduction may be better achieved via iterative image reconstruction methods, such as compressed-sensing based methods (Song et al. 2007; Chen et al. 2008; Sidky et al. 2008; Sidky et al. 2009; Defrise et al. 2011; Loris et al. 2011; Sidky et al. 2012; Xu et al. 2012), and statistical model-based methods (Sauer et al. 1993; Bouman et al. 1996; Elbakri et al. 2002; Ramani et al. 2012). These methods incorporate data consistency and/or statistical information into the noise reduction process. However, they consider each set of projection data individually and do not utilize data redundancy between different energy bins. Such redundancy is abundant in PCCT, as all sets of PCCT projection data are from a single scan.
To exploit data redundancy in the energy domain for the task of noise reduction, Leng et al. developed a noise reduction method called local highly constrained back-projection reconstruction (HYPR-LR, (Leng et al. 2011)). This method first creates a composite image using images from all energies so as to achieve the minimum noise level. The composite image is then used to reduce noise in individual energy bin images, while preserving CT number accuracy and spatial resolution in each energy bin image. HYPR-LR is straightforward, efficient and effective, but the noise in each energy bin image is limited by the noise level of the composite image.
This work, which is a continuation of the HYPR-LR method, develops an iterative image reconstruction algorithm for noise reduction in narrow energy bins by making use of images from full-spectrum projection data. The algorithm is based on the PICCS (prior image constrained compressed sensing) framework (Chen et al. 2008; Leng et al. 2008; Szczykutowicz et al. 2010; Lubner et al. 2011; Chen et al. 2012; Lauzier et al. 2012). To distinguish it from PICCS and to emphasize the use of data redundancy in the energy domain, we call this method spectral PICCS. For conciseness, spectral PICCS sometimes is also abbreviated as SPICCS in figure legends.
In this paper, we describe the spectral PICCS algorithm and demonstrate its performance qualitatively and quantitatively via computer simulations and experiments using a research whole-body PCCT system (Gutjahr et al. 2016; Yu et al. 2016; Yu et al. 2016). Both computer-simulation and experimental studies showed substantial noise reduction in narrow energy bin images with preserved CT number accuracy and spatial resolution. In some cases, the noise performance in narrow energy bin images was better than that in the corresponding full-spectrum images, which is a fundamental advantage compared to HYPR-LR.
In this section, the spectral PICCS algorithm is explained and relevant notation is introduced. Then details of the computer simulations are provided, and the figures-of-merit (FOM) used for quantitative assessments are defined. Finally, a brief description of the research PCCT system is provided, along with details of the physical phantom and in-vivo swine studies. The computer simulations used realistic x-ray interaction models but ideal detector response, whereas the phantom and swine studies were performed on a research PCCT system.
The PICCS reconstruction framework (Chen et al. 2008) is described as below. Let x̲ and b̲ be the column vectors denoting the object of interest and its corresponding projection data, respectively. Additionally, let A be the matrix representing the forward projection and x̲p be the prior image. Let
Here c is a weighting factor between (0,1] and TV stands for total variance. The PICCS reconstruction result, denoted as x̲*, is a solution to the constrained minimization problem
Equation (2) can be solved by an iterative method, as presented in (Chen et al. 2008), with each iteration consisting of two steps. First, data consistency Ax̲ = b̲ is enforced using the algebraic reconstruction technique (ART) (Kak et al. 2001). Second, TV summation expressed in Equation (1) is minimized using a gradient descent algorithm. These two steps alternate until the difference between the images reconstructed from two neighboring iterations is below a predefined threshold.
The spectral PICCS algorithm follows the PICCS framework; see the diagram in Figure 1. We consider two sets of projection data, one being bin data that belongs to a relatively narrow energy window, and the other being full spectrum data that contains all detected x-ray photons. The task is to perform image reconstruction for the bin data.
To start, two images were reconstructed using the conventional filtered-back-projection (FBP) algorithm (Kak et al. 2001), one for the bin data, denoted as , and the other for the full spectrum data, denoted as . The FBP images of the bin data and full spectrum data were used as the seed image and the prior image of the spectral PICCS algorithm, respectively. Next, the seed bin image was updated by an outer iterative process that contained three steps, namely the simultaneous ART (SART) step, the projection-onto-convex (POC) step, and the TV minimization step. Note that each SART or TV minimization step was also iterative, which we call inner iterative processes so as to distinguish them from the outer iterative process.
The SART step updated the bin image every projection view, and traversed all projection views in a random fashion to accelerate convergence (Köhler 2004). The forward projection of the SART step was implemented using the Joseph method (Joseph 1982). The relaxation parameter used for the inner SART iteration (equivalent to the β in (Sidky et al. 2008)) was 1/k, where k was the iteration number of the outer iterative process. The rationale for choosing 1/k was based on the rules for the relaxation parameter presented in Section 2.1 in (Bertsekas 2015), which would guarantee convergence if ART alone was considered. The POC step ensured that the bin image (in terms of linear attenuation coefficient) was non-negative (Sidky et al. 2008).
The TV minimization step was implemented using the scheme proposed by (Sidky et al. 2009) except that the step size of each TV iteration was calculated using a backtracking line search algorithm [Algorithm 3.1] (Nocedal et al. 2006). For each outer iterative process, the number of TV iterations was fixed, e.g., 50.
The spectral PICCS algorithm stopped when one of the following two criteria was met: 1) when the maximum outer iteration number, denoted as Nstop, was reached, or 2) when the difference image between two consecutive outer iterations became smaller than a predefined threshold. For this 2nd condition, let IFBP be the bin image reconstructed using the FBP algorithm, and let and be the bin images obtained from the SART step in the (k-1)th and kth outer iterations, respectively. We define
for 1<k<Nstop, where ‖·‖ stands for the Frobenius norm (L2 norm). We call rSART normalized image update, and use it for quantification of image update. Note that rSART is a relative value that is normalized by the seed image, hence the choice of values for rSART is relatively stable across different objects and scanning techniques.
Overall, the spectral PICCS algorithm was mainly controlled by three parameters, i.e., the iteration number for the TV minimization step, the threshold (normalized image update) for the 2nd stopping criterion, as well as the weighting factor c.
Two digital phantoms were designed for our computer simulations, a characterization phantom and a water phantom. The characterization phantom was composed of seven materials, including water, calcium of 300mg/ml concentration, iodine of 20mg/ml concentration, C1, C2, C3, and C4 (see the central plane of the characterization phantom in Figure 2). C1, C2, and C3 were defined according to (Woodard et al. 1986); they were representative of muscle, grey matter, and white matter, respectively. C4 was defined according to (Kramer et al. 1982), and was representative of bone. The calcium and iodine portions were spherical balls of diameter 40 mm, whereas C1, C2, and C3 were spherical balls of diameter 20 mm. The centers of all these five objects were located on a circle of diameter 120 mm. C4 was a spherical shell of diameter 200 mm and radial width of 10 mm; it was concentric to the 120 mm circle. The free space inside C4 was filled with water.
The water phantom consisted of a cylindrical wall of C4 filled with water. The C4 wall was of longitudinal length 20 mm and radius 200 mm with radial width 10 mm. Therefore, the cross section of the water phantom was exactly the same as the central cross section of the characterization phantom (Figure 2), except that the water phantom did not contain C1, C2, C3, the calcium disk, or the iodine disk.
Projection data of the two digital phantoms were obtained by a simulation tool called DRASIM (Stierstorfer 2000). For pre-patient filters, a 2 mm thick aluminum filter mimicking tube housing was included, and a 0.9 mm thick titanium filter and a 1.5 mm thick aluminum filter were added. Three x-ray interactions, including photoelectric, Compton, and Rayleigh, were taken into account. The detector utilized CdTe based sensors, and was assumed to have an ideal response. Consequently, physical effects such as charge sharing, pulse pile-up, and K-escape were not modeled.
For data acquisition, a full circular trajectory of radius 595 mm was used. This trajectory consisted of 1152 equiangular x-ray source positions in one rotation. This geometry is similar to that used on the research PCCT system. The characterization phantom was placed such that its central cross section was located on the trajectory plane, and the center of the C4 shell was at the center of the circular trajectory. The water phantom was placed such that its middle plane was on the trajectory plane, and its axis coincided with the rotation axis of the circular trajectory. The detector was one circular row composed of 480 channels with quarter-detector pixel shift, which corresponded to 0.5 mm slice thickness and 0.7 mm in-plane pixel pitch at the iso-center. This configuration provided a field-of-view (FOV) of diameter 330 mm at the iso-center.
All projection data were acquired using 1s rotation, 550 mA, and 140 kV. Four energy thresholds were used: 20, 54, 64, and 84 keV. We reconstructed images for four energy bins: [20, 54) keV, [54, 64) keV, [64, 84) keV, and [84, 140] keV. We considered the energy threshold data within [20, 140] keV as the full spectrum data, and used the corresponding FBP images as the prior images for the spectral PICCS image reconstruction of the four energy bins.
Both seed (bin) and prior (full spectrum) images were reconstructed by the conventional FBP algorithm with the ramp filter apodized by a Hanning window within the Nyquist frequency band associated with the image pixel size. For spectral PICCS, the number of TV iterations was set to 50, whereas the maximum outer iteration number was set to 100 (Nstop).
For both digital phantoms, a series of weighting factors c from 0.1 to 1 with a step size of 0.1 was used. All c values were used to investigate the impact of c on spectral PICCS image quality, whereas the c value of 0.5 was used for other quantitative evaluation. The final spectral PICCS images were obtained once the normalized image update, i.e., rSART in Equation (3), dropped below 0.05%. The image pixel size was 0.5 mm × 0.5 mm, whereas the image thickness was 0.5 mm. Prior to image reconstruction, water beam hardening correction was performed for all simulated projection data using the method presented in (Hsieh 2009). After reconstruction, the linear attenuation coefficient images were converted to Hounsfield Unit (HU) using the standard conversion formula and the linear attenuation coefficient of water at the effective energy corresponding to the respective energy bin or threshold.
Spectral PICCS was quantitatively evaluated in terms of convergence, CT number accuracy, spatial resolution, noise, as well as the impact of the weighting factor c.
First, the convergence was analyzed using the normalized image update rSART (Equation (3)). Second, the CT numbers of water, calcium and iodine were measured using disk-like regions of interest (ROIs) that were of radius 25 mm (Figure 3A). The measured CT numbers were compared between the FBP and spectral PICCS images for evaluation of CT number bias. Third, spatial resolution was measured using modulation transfer functions (MTF). To generate an MTF curve, an edge spread function corresponding to the calcium edge was calculated, which was differentiated to generate a line spread function. The line spread function was then averaged radially and Fourier transformed. The final MTF curve was obtained by normalizing the Fourier transformed curve to the magnitude of its zero-frequency component.
Fourth, for a fair noise comparison, spatial resolution was matched between FBP and spectral PICCS images. As will be reported later, the spectral PICCS images had better spatial resolution than the FBP images. To this end, Gaussian low-pass filters were applied to the spectral PICCS images such that the resultant MTF curves corresponding to the calcium disk best matched those measured from the FBP images in terms of root mean square errors. Noise was then measured in the water ROI (Figure 3A) and compared between pairs of FBP and spectral PICCS images that had the matched MTF curves.
Finally, the impact of c on spectral PICCS image quality was investigated in terms of 1) noise power spectrum (NPS) using the water phantom, 2) MTF using the calcium edge of the characterization phantom, 3) noise using the water ROI of the characterization phantom, and 4) CT numbers using the water, calcium and iodine ROIs (Figure 3A). To measure NPS, 40 Poisson realizations were inserted to the same set of projection data, and 40 sets of noisy images were reconstructed. These 40 sets of images were further used to generate 20 sets of difference images. In each difference image, 24 square ROIs of size 10 mm × 10 mm located on a radius of 60 mm were drawn (Figure 3B), and an NPS curve was calculated for each ROI using the method presented in (Siewerdsen et al. 2002). The final NPS curve was obtained by averaging 24×40 NPS curves. The measured NPS curves, MTF curves and noise were studied for different c values.
The research PCCT system used the same platform as the 2nd generation dual-source dual-energy (DSDE) CT system (SOMATOM Definition Flash, Siemens Healthcare GmbH, Forchheim, Germany). It consisted of two x-ray sources that were about 90 degrees apart. One source was coupled to an EID and the other to a PCD; see Figure 4A. We refer to the former source-detector assembly as the EID subsystem, and the latter as the PCD subsystem.
The EID consisted of 64 rows of width 0.6 mm (at iso-center), providing a longitudinal coverage of 38.4 mm. The PCD consisted of 32 rows of width 0.5 mm (at iso-center), providing a longitudinal coverage of 16 mm. The FOVs of the EID and PCD subsystems were 500 mm and 275 mm in diameter, respectively. Unlike the 2nd generation DSDE CT system, the EID and PCD subsystems were not simultaneously energized. However, due to the limited FOV size of the PCD, in-plane data truncation occurred when the PCD subsystem imaged a large object. To avoid such artifacts when imaging a large object, a separate low-dose data-completion scan (DCS) was performed using the EID subsystem. The truncated PCD projection data were estimated and completed using the DCS data. A more detailed description about the DCS method can be found in (Yu et al. 2016).
The research PCCT system provided two data acquisition modes, namely the macro mode and the chess mode. The macro mode allowed two energy thresholds, whereas the chess mode allowed four energy thresholds. In this work, we will only explain how the chess mode works, as all our experimental data were acquired using this mode. For details of the macro mode, we refer to (Yu et al. 2016).
A macro pixel was the smallest readout unit of the PCD, and its size defined the detector row width. Each macro pixel was composed of 4×4 subpixels (Figure 4B), and each subpixel registered two energy thresholds. In the chess mode, the 16 subpixels were interlaced into two groups, with one group registering two energy thresholds, and the other registering two different energy thresholds. As a result, the chess mode allowed 4 energy thresholds. With such a configuration, the chess mode only uses half dose for each threshold (Yu et al. 2016), and presents a unique noise correlation between different energy bins (Schmidt et al. 2015). Let E1, E2, E3 and E4 be the 4 thresholds in an ascending order, and let Ekv be the tube potential. The chess mode generated 4 sets of bin data with energy ranges of [E1, E2), [E2, E3), [E3, E4) and [E4, Ekv], as well as 3 sets of threshold data with energy ranges of [E1, Ekv], [E2, Ekv], [E3, Ekv]. In this work, the goal was to apply the spectral PICCS algorithm to the bin data by using FBP images of the [E1, Ekv] threshold data as prior images.
The research PCCT system provided two types of scan protocols, namely the head and body scan protocols. The head scan protocols assumed a 275 mm scan FOV, and thus did not require a DCS, whereas the body scan protocols required a DCS. In this work, the physical phantom studies used the head scan protocols, whereas the swine studies used the body scan protocols.
The head portion of an anthropomorphic torso phantom was scanned on the research PCCT system using the head scan protocols. The head was placed such that it was entirely within the PCD FOV. The head was scanned in the chess mode using an axial acquisition, head protocol, 140 kV and 280 effective mAs (CTDIvol = 71 mGy). The four thresholds used in the chess mode were 25, 48, 65, and 84 keV. For each set of bin or threshold data, 2304 fan-beam projections were generated. These projection data were preprocessed by a workstation provided by the vendor, which included water beam hardening correction, scatter correction and the logarithm operation.
Spectral PICCS image reconstruction was performed for the central plane using the two central rows of projection data, which resulted in a slice thickness of 1 mm. All images were reconstructed using a pixel size of 0.5 mm × 0.5 mm. The seed and prior images were generated by the FBP algorithm using a Hanning apodization window. For the spectral PICCS algorithm, the iteration number of the TV minimization step and the maximum outer iterations were 50 and 100 respectively, the same as those used in the simulation studies. The weighting factor c = 0.5 was used for the spectral PICCS algorithm. For comparison, TV images were also reconstructed by setting c = 1.0. All spectral PICCS (and TV) reconstructions were stopped when the normalized image update (rSART) became smaller than 0.05%.
With approval from our institutional animal care and use committee, a 40 kg female pig was scanned on the research PCCT system. The neck and thorax regions were scanned in the chess mode with an axial acquisition, body protocol. As the entire pig was placed on the patient table, in-plane data truncation occurred in all scans, and DCSs were performed for data truncation corrections.
One neck scan and two thoracic scans were performed at 140 kV. The neck scan used energy thresholds 25, 45, 65 and 85 keV, as well as 440 effective mAs (CTDIvol = 111 mGy). The two thoracic scans were performed right after a contrast enhanced scan, hence there were iodine residuals in the heart. Both thoracic scans were performed with breathing suspended using the same energy thresholds 20, 25, 57 and 77 keV, except that the first thoracic scan used 550 effective mAs (CTDIvol = 60.7 mGy), whereas the second used 276 effective mAs (CTDIvol = 30.5 mGy).
Each scan generated 2304 projections, which were further preprocessed by a workstation provided by the vendor. Besides those preprocessing steps applied to the head phantom data, truncation correction using the DCS method was also applied to the pig data. Note that the output projection data of the DCS method were in parallel-beam geometry, indicating that a rebinning process was involved in the preprocessing steps. Unlike in the simulation and head phantom studies that images were reconstructed using fan-beam data, the pig images were reconstructed using parallel-beam data. However, the reconstruction parameters and kernels (Hanning apodization) used for the pig studies were the same as those used for the simulation and head phantom studies.
For comparison, HYPR-LR images of the swine neck were also generated from the FBP images. For each energy bin, both bin and prior images were low-pass filtered and the ratio of the low-pass filtered bin images over the low-pass filtered prior images were calculated as the weighting images. The final HYPR-LR images were obtained as the multiplication of the prior images and the weighting images (Leng et al. 2011). The HYPR-LR images were then compared to the FBP and spectral PICCS images both qualitatively and quantitatively.
Reconstruction speed of spectral PICCS was evaluated on a personal computer that was equipped with dual core CPUs (Intel(R) Xeon(R) CPU E5520 @ 2.27 GHZ) and 6 GB RAM. Graphic processing units (GPU) were not used. The operation system was 64-bit Windows 7 Professional, and the algorithm was implemented in C. Average reconstruction time of 10 outer iterations was recorded and calculated for each object. The size of the projection data and image is summarized in Table 4 in the Appendix.
The spectral PICCS images of the characterization phantom converged for all energy bins for the 1st 100 iterations; see the convergence curve of Bin-1 images in Figure 5 as an example. The final spectral PICCS images for Bin-1 to Bin-4 were generated using 22, 19, 19 and 24 iterations, respectively; see Figure 6. The average time taken by one iteration was 33 seconds; see Table 4 in the Appendix. The spectral PICCS images of different energy bins demonstrated different features, indicating that the spectral PICCS preserved energy information. Take the beam hardening artifacts indicated by the dashed rectangle in Figure 6 as an example; they are dominant in Bin-1 images but almost vanish in the Bin-4 image. Compared to the FBP images, it is obvious that the spectral PICCS images have less noise. These observations are consistent with our quantitative assessment below.
In terms of CT number accuracy, the differences between the FBP and spectral PICCS images for water, calcium and iodine are shown in Figure 7. All differences converged to a number smaller than 2 HU. Around 20 iterations, which corresponded to the generation of the images presented in Figure 6, the CT number differences between FBP and spectral PICCS were less than 3 HU for all materials and energy bins. Therefore, we conclude that the spectral PICCS algorithm provided satisfactory CT number accuracy.
In terms of spatial resolution, the MTF curves were measured for Bin 2, 3, and 4. The MTF curve of the Bin-1 image was excluded as the calcium edge was contaminated too much by beam hardening artifacts. The widths of the MTF curves of the FBP images at its 50% peak values were 0.258 mm−1 for all three bins, whereas those of the spectral PICCS images were 0.398, 0.375, and 0.367 mm−1 for Bin 2, Bin 3, and Bin 4, respectively. As an example, the MTF curves of Bin 4 are shown in Figure 8. As we can see, substantial enhancement in spatial resolution was achieved in the spectral PICCS images relative to the FBP images that were reconstructed with a Hanning apodization window.
For a fair noise comparison between FBP and spectral PICCS, spatial resolution was matched. The MTF curves corresponded to the low-passed spectral PICCS images fit the MTF curve of the FBP image almost perfectly; see the MTF curves of Bin 4 in Figure 8 as an example. With matched spatial resolution, image noise was measured in the water ROI (Figure 3A) for both the FBP and spectral PICCS images (Table 1). Due to beam hardening artifacts in the Bin-1 images, neither spatial resolution nor noise measurement was performed for the Bin-1 images. Compared to the FBP images, a noise reduction between 55% and 60 % was achieved in the spectral PICCS images. Furthermore, the spectral PICCS images had even less noise than the prior FBP image, which cannot be achieved by the HYPR-LR method (Leng et al. 2011).
Finally we report on the impact of the weighting factor c on spectral PICCS image quality. The NPS curves corresponding to different c values in Bin 3 are plotted in Figure 9A. The NPS peaks of both the spectral PICCS and FBP images appeared at the same frequency location (0.185 lp/mm). In addition, the areas under the NPS curves of the spectral PICCS images were substantially smaller than that of the FBP image. It was clear that the area under the NPS curve achieved a minimum value when c was 0.5, indicating maximum noise reduction.
The results from the characterization phantom were consistent with the NPS results from the water phantom. The MTF curves corresponding to the calcium edge in Bin 4 were plotted in Figure 9B. The spectral PICCS MTF curves did not change much across difference c values, and were much wider than the FBP MTF curve. Also, we note that the measured noise in the water ROI of the characterization phantom was minimum when c was 0.5 (without spatial resolution matching).
In addition, we compared the CT numbers between the FBP images and the spectral PICCS images using different c values. It turned out that the maximum differences in CT numbers of water, calcium, and iodine were 3.4, 3.2 and 3.5 HU, respectively.
In summary, according to the results from both the water and characterization phantoms, a c value of 0.5 was optimal in terms of noise performance of the spectral PICCS algorithm, which is consistent to what was suggested in (Lauzier et al. 2012). Thus, all c values in the other spectral PICCS reconstructions of this work were set to 0.5.
The spectral PICCS images of the head phantom were obtained at the 39th, 36th, 39th, and 36th iterations for Bin 1 – 4, respectively. The average time taken by one iteration was 98 seconds; see Table 4 in the Appendix. All images converged for the 1st 100 iterations.
The reconstruction results of the head phantom are shown in Figure 10. The spectral PICCS images preserved spectral information. For example, compared to the Bin-1 images, substantial reductions in beam hardening and calcium blooming were achieved in Bin 4; see the regions indicated by arrows in Figure 10. Further, compared to the FBP images, considerable noise reduction was achieved in the spectral PICCS images. The amount of noise reduction can be better appreciated in a zoomed in region; see Figure 11. The region used in this figure is indicated by the dashed rectangle in the bottom right of Figure 10. Noise measurements are presented in Table 2; these results were measured in the region indicated by the solid rectangle as shown in Figure 10. Noise reductions ranging from 43% to 53% were achieved. Note that these results were obtained without matching spatial resolution, and thus are conservative. It is also important to point out that this amount of noise reduction was achieved without introducing patchy artifacts, which are evident in the TV reconstructions (Figure 10 and Figure 11, right columns).
The spectral PICCS images of the pig neck were obtained at the 55th, 50th, 50th, and 51th iterations for Bin 1, 2, 3, and 4, respectively. The average time taken by one iteration was 220 seconds; see Table 4 in the Appendix. Compared to the FBP seed images, the spectral PICCS images achieved considerable noise reduction without introducing patchy artifacts (Figure 12B and E). In particular, spectral PICCS images recovered some subtle structures that were less visible in the FBP images; see the structures in the dashed red circles in Figure 12. In contrast, the same structures were noticeably contaminated by patch artifacts in the TV image (Figure 12D). The HYPR-LR images also achieved considerable noise reduction and had recovered subtle structures, however, the amount of noise reduction in HYPR-LR images were observed less than in the spectral PICCS images. These visual inspections were confirmed by our quantitative analysis below.
For quantitative analysis, we measured CT numbers and noise in the pig neck images. For CT numbers, we plotted the profiles along the line depicted in Figure 12A. The small disks above the pig neck were calcium (left four) and iodine solutions (right four) with different concentrations. For each bin, the spectral PICCS, FBP, and HYPR-LR profiles matched each other except that the spectral PICCS profile appeared smoother indicating lower image noise; see Figure 13 as an example. In addition, the slope of these profiles also overlap well among the three algorithms, which indicates that the spectral PICCS, FBP, and HYPR-LR profiles had similar spatial resolution.
The noise was measured in a region of interest (ROI) indicated by the white circle (Figure 12A), and the results are shown in Table 3. For Bin 2, 3 and 4, the spectral PICCS images achieved 49–73% noise reduction relative to the FBP images. In addition, the spectral PICCS images in Bin 2, 3 and 4 had less noise than that of the prior images. The HYPR-LR images also achieved substantial noise reduction, but their noise in bin images was higher than the prior images as expected, except for Bin 1. In Bin 1, the HYPR-LR image showed lowest noise among the three algorithms and was even lower than that in the prior image, which is considered the theoretical limit. The exact reason was unclear, however we suspected this was at least partially related to the fact that Bin 1 image, with the lowest energy range, had most contamination from any non-ideal physical effects, such as charge-sharing, k-escape and beam hardening. Note that bone beam hardening correction was not applied to the FBP images.
First, we report reconstruction results from Bin-1, starting with the 550 mAs scans. Because of the narrow energy window, the FBP image demonstrated excessive noise (Figure 14A) such that most details of the soft tissue were lost. In contrast, the spectral PICCS image (Figure 14B) recovered most of the soft tissue details. This Bin 1 image was obtained using 31 iterations, with each iteration took about 185 seconds on average; see Table 4 in the Appendix. The correctness of the recovered soft tissue was validated by the prior image (Figure 14C). We further subtracted the FBP Bin-1 image from the spectral PICCS Bin-1 image. The mean CT number measured inside the black ROI shown in Figure 14D was −0.7 HU, which indicates that the prior image did not induce bias in the spectral PICCS image.
Several other observations were also made in Figure 14. First, although the energy range of Bin 1 was located in the lower end of the full spectrum that was used for the prior image, the contrast of myocardium was not noticeably enhanced neither in the FBP Bin-1 image nor in the spectral PICCS image. Second, although soft tissue details were recovered, the noise texture in the spectral PICCS image was not as natural when compared to the prior image. Finally, it should be pointed out that when the dose was decreased by half (276 mAs), both FBP and spectral PICCS failed.
Furthermore, we compared the spectral PICCS images of the 276 mAs scans to the FBP images of the 550 mAs scans; see Figure 15. The spectral PICCS reconstruction results of Bin 2, 3 and 4 were obtained using 27, 27 and 26 iterations, respectively. We measured noise in both FBP and spectral PICCS images using the white circular ROIs as indicated in the bottom row of Figure 15. For FBP, the noise was 25.0, 30.5 and 16.0 HU in Bin-2, −3 and −4 images, respectively. For spectral PICCS, the noise was 9.3, 14.2 and 7.7 HU in Bin-2, −3 and −4 images, respectively. As the appearance of the spectral PICCS half-dose images was comparable to that of the FBP full-dose images, we conclude that use of spectral PICCS could allow at least 50% dose reduction when compared to FBP. Note that the myocardium in the FBP images were different from that in the spectral PICCS images in terms of both shape and contrast, because the projection data used for the two reconstructions were acquired in-vivo during different time periods.
Spectral PICCS used three main tuning parameters: the iteration number for the TV minimization step, the threshold for the normalized-image-update stopping criterion, and the weighting factor c. In this work, these parameters were set to be identical for the simulation and experimental studies, and satisfactory spectral PICCS images were obtained for each case, which indicated that these parameters were reasonable for a range of object sizes and scan techniques. This is important for translation to clinical practice, as the parameters did not require substantial tuning efforts for different imaging scenarios. Regarding the weighting factor c, it was shown in Section 3.1 that it had negligible impact on the shape of NPS and MTF curves for high contrast objects, but substantially changed the noise magnitude, which reached a minimum for c = 0.5. These findings, together with those reported in (Lauzier et al. 2012), that c has a noticeable impact on noise texture and that c around 0.5 results in a noise texture similar to that of FBP images, suggest fixing c at 0.5 for the spectral PICCS algorithm. Consequently, two tuning parameters should be satisfactory in most scenarios for the spectral PICCS algorithm.
It is worth mentioning that the x-ray spectrum was non-overlapping between different energy bins in our computer simulation studies, but were overlapping in the real data studies due to detrimental physical effects such as charge sharing and K-escape. Nevertheless, the spectral PICCS produced satisfactory reconstruction results. We also showed that spectral PICCS worked well in a case of an extremely narrow energy window (relative to the range of x-ray spectrum for medical usage). For the 550 mAs pig-thorax scans, spectral PICCS was able to recover almost all the soft tissue details that were lost in the FBP Bin-1 images ([20, 25] keV), without introducing bias. Note that the Bin-1 projection data from the 550 mAs scans were contaminated by various detrimental physical effects such as charge sharing and K-escape. Consequently, Bin-1 data may have used more x-ray photons than what would be expected in an ideal scenario, which could help allay photon starvation. On the other hand, the detrimental physical effects may have increased the effective energy of Bin-1, which is likely the reason that the myocardium contrast in Bin 1 was not noticeably better than that in the prior image.
In the simulation studies of the characterization phantom (Section 2.2.4), the spatial resolution of SPICCS was higher than the corresponding FBP images. To match the spatial resolution between these two algorithms, a Gaussian filter was applied to the SPICCS images so that the spatial resolution of the calcium disk was matched between the two algorithms. Unlike FBP algorithms, spectral PICCS is non-linear, and its images may have different spatial resolution for edges of different contrast levels (Li et al. 2014; Yu et al. 2015). Consequently, although the spatial resolution corresponding to the calcium disk of the characterization phantom was matched between the spectral PICCS and FBP images, it might not be matched for low contrast objects such as C2. The choice of the calcium disk used in this work for spatial resolution matching was mainly based on the fact that its edge was long and of high contrast such that MTF curves could be accurately measured. Also, the spatial resolution improvement shown in this case was measured with piece-wise constant objects, which were suitable for TV-based algorithms. The amount of spatial resolution improvement may be different in cadaveric and animal studies where complex objects and textures exist. Results in our cadaveric and animal studies demonstrated that at least SPICCS didn’t degrade spatial resolution in these scenarios (Figure 13).
Recently a Poly-United-Iterative Reconstruction algorithm involving a concept of prior image (Poly-UIR, (Xi et al. 2015) was proposed for spectral CT. The presented work is different from the Poly-UIR algorithm. First, although the mathematical formulation is similar, the optimization method was different. The Poly-UIR algorithm uses an unconstrained optimization method called split Bregmen, whereas the spectral PICCS uses a constrained optimization method. Although the convergence rate was not reported in (Xi et al. 2015), an unconstrained optimization method in general requires considerably more iterations than a constrained optimization method. In addition, compared to the three tuning parameters (including c) in spectral PICCS, Poly-UIR requires five tuning parameters. Second, projection data used for algorithm demonstration are different as well. For Poly-UIR, projection data were first obtained from a dual-energy CT system capable of fast kV switching, which were further forward projected via computer simulation programs to simulate a PCCT system with a down-sampling rate (60 views per 360 degrees), whereas in the presented work, both simulation and real data from a research PCCT system were used with a typical CT sampling rate (1000~2000 views per 360 degrees). Third, image quality was assessed differently. The work of (Xi et al. 2015) focused on noise reduction with preserved similarity between the Poly-UIR images and the corresponding reference images, whereas our goal in this work was to demonstrate noise reduction with preserved CT number and spatial resolution.
Finally, we would like to point out a couple of potential limitations of the spectral PICCS algorithm in its current form. First, its forward/back projectors, as well as the TV process are all in 2D, and thus the algorithm is for 2D images. The performance of the algorithm for 3D images is currently under development. Second, current spectral PICCS does not involve statistical noise modeling, which could be important in cases where the number of x-ray photons is limited. To this end, many statistical image reconstruction algorithms developed for dual-energy CT, such as those mentioned in the introduction (Sauer et al. 1993; Bouman et al. 1996; Elbakri et al. 2002; Ramani et al. 2012) as well as the dPIRPLE algorithm developed at the Johns Hopkins University (Dang et al. 2014), may be incorporated into the spectral PICCS algorithm.
Our intention of the spectral PICCS algorithm was to reduce noise in narrow energy bins for improved performance in further spectral imaging applications such as multi-material decomposition and K-edge imaging. In future work, we are interested in combining the spectral PICCS algorithm with our general volume constrained material decomposition algorithm (Li et al. 2015). Our ultimate goal is, however, is to incorporate material decomposition into the spectral PICCS framework such that noise reduction and material decomposition can be simultaneously achieved in one step, similar to the idea proposed in (Barber et al. 2015). To further improve the performance, the statistical material decomposition methods that are developed for dual-energy CT, such as those with direct noise modeling for individual set of spectral data (Long et al. 2014; Niu et al. 2014) as well as those involving noise correlation modeling between different sets of spectral data (Kalender et al. 1988; Zhang et al. 2014; Brown et al. 2015; Liu et al. 2015), are inspiring and worth investigation in future work.
A full-spectrum image-constrained reconstruction algorithm called spectral PICCS has been developed. For both ideal spectra separation and overlapping spectra, the algorithm converged and was able to accurately preserve spectral information: the CT number difference between the spectral PICCS images and the FBP images was no larger than 4 HU in water, calcium and iodine for all cases. Furthermore, the spectral PICCS algorithm was able to achieve substantial noise reduction without sacrificing spatial resolution. The amount of noise reduction relative to the FBP images were 1) 55 – 60% in our computer simulation studies with spatial resolution matching and 2) 43 – 73% in our experimental studies without spatial resolution matching. In some of the computer-simulation and swine studies, bin images generated from the spectral PICCS had less noise than the prior image, indicating a fundamental advantage when compared to the HYPR-LR method. We also have shown that the spectral PICCS algorithm could achieve at least 50% dose reduction when used in conjunction with a research PCCT system.
The authors would like to thank Drs. Steffen Kappler, Ahmed Halaweish, Andre Henning for their onsite support of the research photon-counting CT system. We also would like to thank Dr. Erik Ritman, Dr. Ahmed Halaweish, Steven Jorgensen and Jill Anderson for their contributions in the swine studies. The first author would like to thank Dr. Frederic Noo for pointing out the relaxation parameter of the ART algorithm. The project described was supported by Grant numbers R01 EB016966 and C06 RR018898 from the National Institute of Health, in collaboration with Siemens Healthcare. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of Health. The equipment and concepts described in this work are based on a research scanner and are not commercially available.
|Objects||Image Size||Number of|