Raw data normalization for a multi source inverse geometry CT system.

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.


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][2][3][4][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 3 rd 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. 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.

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][8][9][10]. The array is made of four sub-arrays (modules) that are tiled in the x direction, each with 4 by 2 sources [7][8][9][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.

Modelling source intensity fluctuation
We modeled x-ray transmission with source intensity fluctuation as: where N in is the expected number of incident photons, N i, j,k , and p i, j,k are the detected number of photons and projection data for the i th source in the j th view at the k th 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 Eq. (1), the projection data p i, j,k are expressed as where ε i, j = log(α i, j ) is the error caused during j th view by i th source intensity fluctuation. A goal of the data-normalization method is to accurately estimate the error terms.

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 p 1, j,k , p 1, j,k+1 , p 1, j+1,k and p 1, j+1,k+1 are the nearest 4 points of p 2,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.
Step 3. Perform steps 1 and 2 for different detector channels, m, of source 2, and then average theε 2,l,m over the detector channels to formε 2,l .
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.

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

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.    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. 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. 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 mm 2 (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.

Results
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][14][15][16][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.

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 offcenter 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