Open Access
Add to BibSonomyAbstract
A multi-source inverse-geometry CT (MS-IGCT) system consists of a small 2D detector array and multiple x-ray sources. During data acquisition, each source is activated sequentially, and may have random source intensity fluctuations relative to their respective nominal intensity. While a conventional 3rd generation CT system uses a reference channel to monitor the source intensity fluctuation, the MS-IGCT system source illuminates a small portion of the entire field-of-view (FOV). Therefore, it is difficult for all sources to illuminate the reference channel and the projection data computed by standard normalization using flat field data of each source contains error and can cause significant artifacts. In this work, we present a raw data normalization algorithm to reduce the image artifacts caused by source intensity fluctuation. The proposed method was tested using computer simulations with a uniform water phantom and a Shepp-Logan phantom, and experimental data of an ice-filled PMMA phantom and a rabbit. The effect on image resolution and robustness of the noise were tested using MTF and standard deviation of the reconstructed noise image. With the intensity fluctuation and no correction, reconstructed images from simulation and experimental data show high frequency artifacts and ring artifacts which are removed effectively using the proposed method. It is also observed that the proposed method does not degrade the image resolution and is very robust to the presence of noise.
© 2015 Optical Society of America
1. Introduction
The multi-source inverse-geometry CT (MS-IGCT) system consists of a small 2D detector array and multiple x-ray sources. A series of sources are aligned transaxially and each source acquires projection data sequentially for a small portion of the entire field-of-view (FOV). A second series of sources are aligned at a longitudinal offset relative to the first series, and increase the volumetric coverage and provide additional sampling of the object, which reduces cone-beam artifacts. Using a small 2D detector array can reduce the detector cost and the amount of scatter, and by modulating the x-ray flux of each source independently, a concept called virtual bowtie, the effective dose to the patient can be minimized [1–5].
During a scan, X-ray source intensity fluctuations often occur and these must be accounted for in order to avoid artifacts in the reconstructed image. Conventional 3rd generation CT systems have a single X-ray source and one or more reference channels at the edges of the detector array (outside the FOV) as shown in Fig. 1(a). The reference channels monitor the source intensity and allow any fluctuations to be properly accounted for by normalizing the raw data during image reconstruction. However, in the MS-IGCT system, multiple sources cannot directly illuminate a common reference channel in the detector array because each source covers only a small portion of the entire FOV, as shown in Fig. 1(b). The sources at the edge of the FOV are the only ones that can directly illuminate detector-edge reference channels. Therefore, projection data from all other sources will include source-fluctuation error terms that cannot be properly normalized in a conventional manner. As discussed in detail below, the projection data of the MS-IGCT system are rebinned into a full fan-beam data set for image reconstruction and data from all sources (and their fluctuation errors) are combined together in the rebinning process. Figure 2(a) shows a profile of the rebinned projection from a uniform water phantom with and without source intensity fluctuation.
Fig. 2 (a) Profile of rebinned projection with and without source intensity fluctuation, (b) profile of filtered projection, and reconstructed water phantom, (c) with and (d) without source intensity fluctuation. The display window of the uniform water phantom is [0.18 0.22] cm−1.
The projection data of each source were generated with random ±10% source intensity fluctuation that introduces discontinuities along the rebinned projection profile. The discontinuities are amplified by the ramp filter as shown in Fig. 2(b). As a result, the reconstructed image has significant high frequency artifacts as shown in Fig. 2(c). For the comparison, a reference image is shown in Fig. 2(d).
In this paper, we present a raw-data normalization method that eliminates image artifacts caused by unintended source-intensity fluctuations in MS-IGCT. This is an amplification of preliminary work presented in an SPIE conference abstract [6]. We first model the source intensity fluctuation by modifying Beer’s law and analyze the sampling pattern of the MS-IGCT system in 2D Radon space. The proposed method is validated using noiseless and noisy simulation data and its effect on image quality is analyzed using MTF and standard deviation of the reconstructed noise image. Projection data were acquired using an experimental MS-IGCT system for an ice-filled PMMA phantom and a post-mortem rabbit. Reconstructed images with and without the proposed method are presented.
2. Materials and methods
2.1. System overview
The MS-IGCT system includes 32 sources arranged in a 16 by 2 array with 16 sources in the x (transaxial) direction and 2 sources in the z direction [7–10]. The array is made of four sub-arrays (modules) that are tiled in the x direction, each with 4 by 2 sources [7–10]. In the composite array, focal spots are separated in z by 100 mm and in x by 25 mm. Each row of 16 transaxial sources produces a 22 cm diameter FOV. The detector uses modules from the newer GE CT750HD detector because of its speed (i.e., 7.5 kHz frame rate), detection efficiency (using GE’s Gemstone material), and improved low-signal performance. Each module consists of 64 by 16 cells of about 1.1 mm × 1.0 mm each. 12 modules are tiled in z producing a total sensitive area of about 70 mm in the x direction and 192 mm in the z direction [4,5]. Figure 3 shows a picture of the integrated MS-IGCT gantry.
2.2. Raw data normalization
2.2.1. Modelling source intensity fluctuation
We modeled x-ray transmission with source intensity fluctuation as:
where Nin is the expected number of incident photons, Ni,j,k, and pi,j,k are the detected number of photons and projection data for the ith source in the jth view at the kth detector channel, and αi,j is the corresponding intensity fluctuation parameter which varies randomly for different view and sources. Note that all detector channels see the same source intensity fluctuation.After taking the negative logarithm in Equation (1), the projection data pi,j,k are expressed as
where εi,j = log(αi,j) is the error caused during jth view by ith source intensity fluctuation. A goal of the data-normalization method is to accurately estimate the error terms.2.2.2. Normalization method
Referring to Fig. 1(b), it is assumed that source 1 directly illuminates reference channels at the edge of the MS-IGCT detector so that the intensity fluctuations (and the associated error term) of source 1 can be directly monitored. As shown in Fig. 4, the projection data of adjacent sources in the MS-IGCT system share an overlap region in 2D Radon space, and therefore the projection data of source 1 can be used as a reference to estimate the error terms of source 2. For example, the projection data of sources 1 and 2 within the overlap region are expressed as:
where p1,j,k, p1,j,k+1, p1,j+1,k and p1,j+1,k+1 are the nearest 4 points of p2,l,m within the overlap region as shown in Fig. 4. Note that the radial distance and view angle of each point in 2D Radon space are pre-calculated using the known system geometry.Fig. 4 Sampling of 2D Radon space for the first three sources: source 1 (S1), source 2 (S2), and source 3 (S3).
To calculate an error estimate, , we use the following steps:
- Step 1. Estimate from p1,j,k, p1,j,k+1, p1,j+1,k, p1,j+1,k+1 within the overlap region using bi-linear interpolation.
- Step 2. Calculate detector channel dependent error, , by subtracting from p2,l,m.
- Step 3. Perform steps 1 and 2 for different detector channels, m, of source 2, and then average the over the detector channels to form .
By performing steps 1–3 for the projection data of source 2 for all views, error terms for all source 2 views can be estimated. Then, raw data of source 2 are normalized correctly by subtracting the estimated error terms from the projection data of source 2. In the same manner, raw data of the next adjacent sources are normalized using corrected sources as a reference. Note that using corrected source as a reference for the next uncorrected source will introduce an accumulating error over all corrected sources. The effect of this error propagation on image quality is examined as discussed below.
2.3. Computer simulations
The proposed method was validated using simulated noisy and noiseless projection data. An analytical Shepp Logan phantom [11] and a uniform water phantom (μ = 0.2 cm−1, diameter = 20 cm) were simulated, and projection data were produced by calculating the line integrals along each ray through the phantom. The finite focal spot and detector cells were modeled by averaging projection data from discrete sub-sources and sub-detectors spanning the finite apertures. Note that the simulation assumed a hypothetical mono-energetic beam and an ideal photon counting detector. Source intensity fluctuations were modeled from a uniform distribution with maximum values of ±10% of the actual source intensity. We also simulated the case where the sources have random flux differences from each other, but the errors are consistent over all views. For noisy projections, photon arrivals were modeled from a Poisson distribution with mean of 5000 photons across all detector channels. The typical SNR of clinical CT head scans is around 100, which is comparable to our noise simulation setting. We only modeled Poisson-distributed photon noise; this should also be relevant to other random noise in the imaging system. For image reconstruction, corrected projection data were rebinned into a full fan beam geometry [5], and then a direct fan beam reconstruction algorithm was used. As a reference, images were reconstructed from projection data without source intensity fluctuation. The simulation and reconstruction parameters are summarized in Table 1. Note that the simulation parameters match the real system geometry.
Table 1. Simulation and reconstruction parameters
The proposed fluctuation-correction method incorporates a spatial interpolation of the estimated errors, which introduces error propagation to other sources and can also introduce blurring. Accumulation of statistical errors, if significant, should be visible in images of the uniform phantom. To examine the effect of the proposed method on image resolution, projection data of a 0.5 mm diameter sphere centered at (1 mm, 1 mm) and (70 mm, 1 mm) were generated with and without source fluctuation errors. The proposed method was applied to the projections that included errors, and the MTF was calculated from the reconstructed image. For reference, the MTF was also calculated from a reconstructed image formed using projection data without errors. For the quantitative comparison, mean square error (MSE) between the reference MTF and the MTF with correction was calculated. To test the performance of the proposed method with noisy data, noise-only projection data (i.e., no object) with error terms were generated using uniform noise (i.e., 5000 photons per ray), and then reconstructed after using the proposed method. A reference noise image was also reconstructed, and then the standard deviation (STD) across the FOV was compared. The STD was calculated from 100 noise realizations using a 10 by 10 pixel ROI, and the MSE of STD from the reference noise image and the noise image with correction was calculated.
2.4. Experimental setup and phantom scan
Experimental data were acquired for an ice-filled PMMA cylinder and for a post-mortem rabbit in a water-filled PMMA cylinder that was frozen to keep the rabbit stationary relative to the rotating PMMA cylinder. For data acquisition, the objects were rotated while each source illuminates the objects sequentially. The objects were rotated because gantry rotation was not available for these experiments. Each source acquired 270 views over 360 degrees of the object rotation, the sources operated at 80 kVp and 100–300 mA, with each source activated individually for 10–100 ms, and the total exposure time was 8 seconds. In addition, flat field data (air shots) for each source were acquired, and beam hardening correction was performed by using a first order polynomial fit, where the fitting parameters were estimated by comparing the analytical projection data of a uniform cylinder phantom with the scan data of a uniform PMMA cylinder phantom [12]. Then, projection data with and without the proposed correction were reconstructed using reconstruction parameters in Table 1.
3. Results
Figure 5(a) shows simulated reconstructed images and corresponding central profiles for a uniform water phantom without noise (infinite flux). The figure shows reference results with no source fluctuation and results with source fluctuation both with and without fluctuation correction. The images are displayed with a window level of [0.18 0.22] cm−1. Figure 5(b) shows a similar set of results for a Shepp Logan phantom. These images are displayed with a window level of [1 1.05] cm−1. The images from data containing errors and without correction show high frequency artifacts. However, with the proposed method, the artifacts were removed effectively. Figures 6(a) and 6(b) show the reconstructed images and corresponding central profiles for noisy data, demonstrating the robustness of the proposed method.
Fig. 5 (a) Reconstructed images and central profiles of a uniform water phantom and a Shepp Logan phantom without correction. Noiseless simulation of (a) a uniform water phantom and (b) a Shepp Logan phantom. The display windows are [0.18 0.22] cm−1 and [1 1.05] cm−1 for the uniform water phantom and the Shepp Logan phantom, respectively.
Fig. 6 (a) Reconstructed images and central profiles of a uniform water phantom and a Shepp Logan phantom with and without correction. Noiseless simulation of (a) a uniform water phantom and (b) a Shepp Logan phantom. The display windows are [0.18 0.22] cm−1 and [1 1.05] cm−1 for the uniform water phantom and the Shepp Logan phantom, respectively.
We also examined the artifact pattern when the source intensity varied randomly from source to source but were consistent over all views of each source. As shown in Fig. 7, the artifacts generated were rings and bands, as expected when errors are consistent across views. This is also the case with detector channel errors in conventional systems. The proposed method effectively removed the artifacts.
Fig. 7 Reconstructed noise images of a Shepp Logan phantom with and without correction. The display window is [11.05] cm−1.
Figure 8(a) shows the reference point-spread function (PSF) and the PSF with source intensity fluctuation and correction, both near isocenter and off-center. The image size of the PSF is 29 × 29 mm2 (100 × 100). Figure 8(b) compares the MTF near isocenter and off-center from the reference image and the image with correction. It is observed that the proposed method does not degrade the image resolution since it performs simple subtraction of the error terms from projection data. It is only the error term that is interpolated. Note that the MSE between the reference MTF and the MTF with correction is 2.5 × 10−4 near isocenter and 2.0 × 10−4 at off-center, respectively. While the source intensity error terms can be estimated reliably with noiseless data, the estimation performance of the error terms might be degraded when noise is present.
To examine the effect of noise, the STD was calculated for several points across the FOV for both the reference noise image and the noise image with source intensity fluctuation and correction. Figure 8(c) shows the reference noise image, noise image with correction, and STD across the FOV. The MSE of STD between the reference noise image and noise image with correction is 1.4 × 10−6. It is observed that the STD from both noise images matched well, which demonstrates that any estimation errors caused by the proposed method are negligible compared to CT image noise. Note that the increasing STD for the larger radial distance is an inherent property of the direct fan beam reconstruction algorithm [13–17].
Figures 9(a) and 9(b) show experimental results for the ice-filled PMMA phantom and for the frozen rabbit. Each figure shows the reconstruction without correction, the reconstruction with correction, and central profiles with and without correction. Before the correction, both images show ring type artifacts, which degrades the uniformity and low contrast detectability in the reconstructed images. The proposed method removed the ring type artifacts effectively and improved the uniformity and low contrast detectability of the reconstructed images.
Fig. 8 (a) Reference PSF and PSF with correction near isocenter and off-center. (b) Comparison of corresponding MTF near isocenter and off-center. (c) Reference noise image, noise image with correction, and comparison of STD across the FOV from both noise images.
Fig. 9 Reconstructed images and central profiles of (a) an ice-filled PMMA phantom and (b) a rabbit without and with correction. The display windows is [0.15 0.25] cm−1 for the ice-filled PMMA phantom and a rabbit.
4. Conclusions
In this work, we investigated the effect of source intensity fluctuations on image quality in MS-IGCT and presented a raw data normalization method. In simulations, the source intensity fluctuations were generated from a uniform distribution and added to the noiseless and noisy simulation data. Without correction, the reconstructed images showed high-frequency image artifacts and rings or bands, depending on whether the errors varied or were consistent across views. These artifacts were effectively removed by the proposed method.
To test the effect of the proposed method on image resolution, MTF near isocenter and off-center were compared between a reference (no source fluctuation) image and a corrected image. As expected, the reference-image and corrected-image MTFs matched well. We also calculated the STD across the FOV from the reference noise image and corrected noise and observed that the STD of both noise images showed good agreement. The MTF study and noise simulation demonstrate that the proposed method does not significantly degrade the image resolution or noise performance.
The proposed method was also tested using experimental data of an ice-filled PMMA phantom and a frozen rabbit. While our simulation assumed that the maximum intensity fluctuation of each source is the same (i.e., 10%), the real data (i.e., air scan data of each source) showed that the peak variation of intensity for each source varies from 2% to 6% and are sometimes consistent over views, and therefore the reconstructed images showed both ring and high frequency artifacts. The proposed method removed these artifacts effectively and improved the uniformity and low contrast detectability in the reconstructed images.
Acknowledgments
The authors of this paper are grateful for the editing help provided by Denise Largen. Research reported in this publication was supported in part by Grant Number R01 EB006837 from the National Institute of Biomedical Imaging and Bioengineering (NIH-NIBIB), in part by GE Global Research, and in part by GE Healthcare. This research was also supported by the MSIP (Ministry of Science, ICT and Future Planning), Korea, under the IT Consilience Creative Program ( NIPA-2014-H0201-14-1002) supervised by the NIPA (National IT Industry Promotion Agency). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of Biomedical Imaging and Bioengineering or the National Institutes of Health.
References and links
1. B. De Man, S. Basu, D. Beque, B. Claus, P. Edic, M. Iatrou, J. LeBlanc, B. Senzig, R. Thompson, M. Vermilyea, C. Wilson, Z. Yin, and N. J. Pelc, “Multi-source inverse geometry CT: a new system concept for x-ray computed tomography,” Proc. SPIE 6510, 65100H (2007). [CrossRef]
2. B. De Man, S. Basu, P. Fitzgerald, D. Harrison, M. Iatrou, K. Khare, J. LeBlanc, B. Senzig, C. Wilson, Z. Yin, and N. J. Pelc, “Inverse geometry CT: the next generation CT architecture?” in Proceedings of IEEE NSS-MIC (IEEE, 2007), pp. 2715–2716.
3. T. G. Schmidt, R. Fahrig, and N. J. Pelc, “An inverse-geometry volumetric CT system with a large-area scanned source: a feasibility study,” Med. Phys. 9(31), 2623–2627 (2004). [CrossRef]
4. B. De Man, A. Caiafa, Y. Cao, K. Frutschy, D. Harrison, L. Inzinna, R. Longtin, B. Neculaes, J. Reynolds, J. Roy, J. Short, J. Uribe, W. Waters, X. Zhang, Y. Zou, B. Senzig, and N. J. Pelc, “Multi-source inverse-geometry CT: from system concept to research prototype,” in Proceedings of IEEE NSS-MIC (IEEE, 2009), pp. 3531–3533.
5. J. Baek, B. De Man, J. Uribe, R. Longtin, D. Harrison, J. Reynolds, B. Beculaes, K. Frutschy, L. Inzinna, A. Caiafa, B. Senzig, and N. J. Pelc, “A multiple-source inverse-geometry CT system: initial results with an 8 spot x-ray source array,” Phys. Med. Biol. 59(5), 1189–1202 (2014). [CrossRef] [PubMed]
6. J. Baek and N. J. Pelc, “Data normalization method for a multisource inverse geometry CT system,” Proc. SPIE 8313, 83132B (2012). [CrossRef]
7. K. Frutschy, B. De Man, P. Edic, B. Lounsberry, B. Neculaes, Y. Zou, L. Inzinna, K. Conway, X. Zhang, Y. Cao, S. Gunturi, A. Caiafa, and J. Reynolds, “X-ray multisource for medical imaging,” Proc. SPIE 7258, 725822 (2009). [CrossRef]
8. K. Frutschy, B. Neculaus, L. Inzinna, A. Caiafa, J. Reynolds, Y. Zou, X. Zhang, S. Gunturi, Y. Cao, B. De Man, D. McDevitt, R. Roffers, B. Lounsberry, and N. J. Pelc, “High power distributed x-ray source,” Proc. SPIE 7622, 76221H (2010). [CrossRef]
9. B. Neculaes, L. Inzinna, A. Caiafa, J. Reynolds, K. Frutschy, B. De Man, J. Uribe, R. Longtin, D. Harrison, and N. J. Pelc, “A high power distributed x-ray source for multisource inverse geometry CT,” Proc. SPIE 7622, 76221H (2010)
10. J. Uribe, J. Reynolds, L. Inzinna, R. Longtin, D. Harrison, B. De Man, B. Neculaes, A. Caiafa, W. Waters, K. Frutschy, B. Senzig, J. Baek, and N. J. Pelc, “Multi-source inverse-geometry CT: prototype system integration,” in Proceedings of IEEE NSS-MIC (IEEE, 2010), pp. 2578–2581.
11. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, 1988).
12. J. Hsieh, Computed Tomography Principles, Design, Artifacts, and Recent Advances (SPIE Publications, 2003).
13. K. E. Bennet and R. L. Byer, “Fan-beam-tomography noise theory,” J. Opt. Soc. Am. A 3(5), 624–633 (1986). [CrossRef]
14. K. E. Bennet and R. L. Byer, “Optical tomography: experimental verification of noise theory,” Opt. Lett. 9(7), 270–272 (1984). [CrossRef]
15. K. E. Bennet, G. W. Faris, and R. L. Byer, “Experimental optical fan beam tomography,” Appl. Opt. 23(16), 2678–2685 (1984). [CrossRef]
16. G. L. Zeng, “Nonuniform noise propagation by using the ramp filter in fan-beam computed tomography,” IEEE Trans. Med. Imag. 23(6), 690–695 (2004). [CrossRef]
17. L. Zhu and J. Starlack, “A practical reconstruction algorithm for CT noise variance maps using FBP reconstruction,” Proc. SPIE 6510, 651023 (2007). [CrossRef]
