Abstract
Soft tissue images from portable cone beam computed tomography (CBCT) scanners can be used for diagnosis and detection of tumor, cancer, intracerebral hemorrhage, and so forth. Due to large field of view, X-ray scattering which is the main cause of artifacts degrades image quality, such as cupping artifacts, CT number inaccuracy, and low contrast, especially on soft tissue images. In this work, we propose the X-ray scatter correction method for improving soft tissue images. The X-ray scatter correction scheme to estimate X-ray scatter signals is based on the deconvolution technique using the maximum likelihood estimation maximization (MLEM) method. The scatter kernels are obtained by simulating the PMMA sheet on the Monte Carlo simulation (MCS) software. In the experiment, we used the QRM phantom to quantitatively compare with fan-beam CT (FBCT) data in terms of CT number values, contrast to noise ratio, cupping artifacts, and low contrast detectability. Moreover, the PH3 angiography phantom was also used to mimic human soft tissues in the brain. The reconstructed images with our proposed scatter correction show significant improvement on image quality. Thus the proposed scatter correction technique has high potential to detect soft tissues in the brain.
1. Introduction
Cone beam computed tomography (CBCT) scanners have been initially developed for dental applications. Dental CBCT images usually represent only bone structures due to limited contrast resolution of X-ray detectors. In recent years, flat panel detectors (FPDs) were tremendously improved in their capability of detecting small differences in attenuation of X-ray beam, thus providing soft tissue detectability, such as tumor, muscle, and intracerebral hemorrhage (ICH) [1–4]. This particular FPD model has been introduced in portable CBCT scanners with large field of view (FOV) to cover the head region. The portable CBCT scanner can be freely moved where rapid diagnosis is needed such as in the emergency and operation rooms [2–4]. The scanner can perform patient screening before treatment, during operation, and after operation. For having a large FPD, soft tissue images from CBCT scanner are usually degraded as more artifacts occur such as lag, glare [5], motion artifact [4, 5], beam hardening effects (BHE) [4–6], and X-ray scattering effects [5–7]. However, it is well known that the most critical artifact is caused by the X-ray scattering effect which reduces image quality. X-ray scatter signals directly affect contrast and CT numbers of soft tissue images. The number of X-ray scatter signals is increased as an increase in FOV and object thickness [7].
There are many methods to reduce X-ray scatter signals in FPDs, such as antiscatter grid plate [8], beam stop array plate [9], primary modulation methods [10], Monte Carlo simulation (MCS) software [11, 12], and X-ray scattering models [13–18]. The antiscatter grid plate can reduce X-ray scatter signals to increase image quality; however, it is not sufficient for improving image quality on soft tissue images [8], and the apparent dose increase is always undesired. A beam stop array plate [9] can be used for measuring X-ray scatter signals; however, this method is not practical in clinics. The primary modulation method was to place a high frequency attenuation plate between an X-ray source and an object to obtain a modulated projection image, and then this modulated projection image was filtered to eliminate low-frequency components due to X-ray scatter signals by a high-pass filter [10]. Monte Carlo simulation software based on Geant4 libraries [11, 12] has been used widely to design the simulated system for estimating high accuracy X-ray scatter signals. Generally, the MCS software usually utilized expensive computation time. The technique proposed by Star-Lack et al. [14] described efficient X-ray scattering correction on the large ellipse phantom using asymmetric kernels. Since our work aims to reduce the X-ray scatter signals in the patient’s head only, symmetric kernels can be assumed. Several iterative techniques were proposed for scatter correction including subtraction from measured data [14, 15] and deconvolution by the statistical method [13, 16, 17]. The problem of the subtraction method is that the corrected value can become negative as discussed in [16].
In this work, we propose the X-ray scattering correction method for improving soft tissue images on the large flat panel detector of portable CBCT. We apply the deconvolution technique based on the maximum likelihood expectation maximization (MLEM) method [13, 16, 17] for estimating the primary signal. First, the MCS software is used to form the shape and the amplitude of the X-ray scatter signals as a point spread function (PSF) and a scatter fraction (SF) function according to actual parameters of the CBCT system. Second, we match equivalent thickness in projection images to obtain subprojection images. Third, those subprojection images are convolved with the prepared kernels or PSFs and combined to obtain a convolved projection image. Finally, the convolved projection image is used in the MLEM method to estimate the primary signal, and the process is iterated until convergence to FBCT data. Once the primary signal is estimated, the reconstruction algorithm based on filtered backprojection (FBP) technique [19] is employed. In the experiment, our proposed technique will be tested with the QRM-ConeBeam phantom by QRM GmbH, Germany [20], and the PH3 angiographic CT head phantom ACS by Kyoto Kagaku Co., Ltd., Japan [21]. The experimental results are compared with FBCT data.
2. Materials and Methods
2.1. Materials
Our prototype portable CBCT scanner prototype consists of the flat panel detector (Varian PaxScan 4030CB) and the high frequency X-ray source which shoots X-ray pulses by synchronizing with exposure time of the flat panel detector. The distance from source to detector (DSD) and the distance from source to an object (DSO) are 500 mm and 786 mm, respectively. The dynamic gain was operated for the pixel size of 0.388 mm. This system used 90 kVp, 9 mA, and the total filtration of 0.5 mm Aluminum. A full rotation of scanning was achieved. In this work, we used PSF and SF from MCS software to simulate X-ray scatter signals by using PMMA sheets with a pencil beam as shown in Figure 1.
The MCS software computed the X-ray scatter signals by simulating a pencil beam with the PMMA sheets; one set of the PMMA sheets was 38 mm thick and the total of 10 sets was used. The measured X-ray scatter signals at each thickness are interpolated to achieve X-ray scatter signals at thickness of 1 mm. We used 2 phantoms in the experiments: the QRM-ConeBeam phantom by QRM GmbH, Germany [20], is a cylindrical tissue equivalent phantom at 120 kVp with the diameter of 160 mm and the height of 160 mm and the PH3 angiographic CT head phantom ACS by Kyoto Kagaku Co., Ltd., Japan [21], is a life-size adult phantom made of Urethane base resin (SZ-50) and Epoxy base resin as shown in Figure 2.
(a)
(b)
2.2. Method
According to the fundamental law of physics, while the X-ray beam travels to penetrate each layer inside an object, there exist three phenomena of interactions: photoelectric effect, Compton scattering, and Rayleigh Scattering of X-ray photons. The Compton scattering effect is essential for cross section images reconstruction. The X-ray beam that penetrates an object and goes directly to the flat panel detector is desired; however, scattered X-rays from other directions are usually combined at the same sensor position of the flat panel detector. This degrades image quality, such as cupping artifacts, low contrast, and inaccurate CT numbers. We can write the measured X-ray signal model as follows: where is the measured signal, is the primary signal which is the expected signal to be used in cross section image reconstruction, is the scatter signal, and denotes the pixel coordinate in the projection image. The scatter signal can be written in the form of the primary signal convolved with the kernel function, , as follows: where denotes a 2D convolution operator. The statistical method is based on Bayes’ rule as discussed in [13]. The average value of the statistical method is considered by the X-ray photons behind an object as Poisson distribution [13]. The maximum likelihood is used to estimate the primary signal, , and can be written in the form of the MLEM [13, 16, 17] algorithm as follows: where is the estimated primary signal at the th iteration. The kernel measurement at different thickness, , is obtained according to the thickness of the PMMA sheet, so is divided according to thickness, . We can rewrite the above equation as follows:The scatter correction method starts by creating a database of kernels according to thickness of PMMA sheets using the MCS software and initializing the primary signal, . The next step is to match each pixel of the primary signal with the estimated PMMA equivalent thickness; that is, the projection image data are divided into different groups or subprojection image data according to thickness. The subprojection image data sets of the primary signal are convolved with the kernels and summed up to obtain the convolved projection image. Then the convolved projection image is used to perform deconvolution by the MLEM algorithm. Finally, we check the condition for convergence. If it does not converge, we will return to the thickness mapping process again. In summary, all steps of the scatter correction can be illustrated as in Figure 3.
2.2.1. Kernel Measurements
Kernel measurements in this work are obtained by measuring different thickness of the object. Generally, kernels can be measured from experiments or simulated by the MCS software. For accuracy of measuring the X-ray scatter signal, we used the MCS software to simulate the scatter signals according to the actual CBCT scanner system: 90 kVp Voltage as a pencil beam and the total flat filtration of 5.5 mm Al. The pencil beam penetrates through PMMA sheets. We used the total of 10 sets of PMMA sheets with the thickness of 38 mm each set. The measured scatter signals are interpolated to obtain the scatter signal at each 1 mm thickness, so the total of 380 scatter signals is obtained. The example of scatter signals using this work is shown in Figure 4. They are normalized to obtain the point spread function (PSF) according to thickness, . The kernel at each thickness can be described as follows: where and denote the spatial coordinate, is PMMA thickness, and is the compensating value from the experiments, which depends on thickness. The values of used in this work are shown in Table 1. Moreover, the scatter fraction values are measured as an average at the center region of each kernel for a suitable size of region of interest (ROI). In the experiments, we employ linear curve fitting to the measured average value of scatter fraction. The scatter fraction, , can be described as follows: where and are coefficients of the linear function. The SF value depends on object thickness and highly affects the amount of scatter correction.
2.2.2. Thickness Map Measurements
We match the projection image data to obtain the subprojection image data set according to PMMA thickness. First, we create the log signal function from pure PMMA plates for thickness mapping. We divide thickness of the projection image, , by using Beer’s law [11, 12, 19] as follows: where is the linear attenuation coefficient at the effective energy, is the measured primary signal without attenuation, and is the primary signal. We divide the projection image data of the object into the subprojection image data set as the log signal function of PMMA sheet thickness. The measured primary signal at different PMMA thickness from the MCS software can create the log signal function as the relationship between actual PMMA thickness and the log signal value. This function is used for transferring the log signal value to the equivalent value of PMMA thickness.
2.2.3. Evaluation
Evaluation in the experimental results starts as checking the performance of the MLEM method whether it converges according to the log likelihood function, [13]. For simplicity, we ignore the insignificant terms as follows:To verify our proposed algorithm, we will compare both projection images and cross section images with the FBCT data. The CT number in the reconstructed images is calibrated by using the water average value in the CT number section of the QRM-ConeBeam phantom as shown in Figure 5(a). Note that the water average value is measured by using FBCT instead of CBCT to avoid scatter effects. We can normalize the Hounsfield unit (HU) as follows:where is the average value of any material in the cross section images and is the water average value from FBCT. To measure the low contrast detectability, we use Sections A, B, and C of the QRM-ConeBeam phantom. Section A has higher contrast between inserts and background than Sections B and C, while Section C has the lowest contrast. Here we measure the density value at numbers 1 to 4 as shown in Figure 5(b). Different gray level values represent different density values.
(a)
(b)
After estimating scatter correction, the contrast is increased; however, the noise signal value in the projection images is increased as well. In this study, we measure the contrast value between two different inserts and its contrast to noise ratio (CNR) as follows:where is the average value in the background region and and are the standard deviation (STD) value of any insert material and the background, respectively. However, we also measure percentage of cupping in the cross section images in order to evaluate the proposed method and the remaining beam hardening effect as follows:where is an average CT number value from four peripheral ROIs in the uniform areas and is the average CT number value at the center.
3. Results
All parameters used in this study are shown in Tables 1 and 2. Table 2 shows the implementation parameters including the coefficient values of the SF function: , ; number of groups: 5; and the thickness of each group: 40 mm. From the experiment, we tried dividing into different groups and found that five groups with 40 mm thickness in each group were sufficient for acceptable estimation of the primary signal. Moreover, according to the experiment, the proposed method did not introduce any additional artifacts.
3.1. Scatter Correction Results in the Projection Image
The scatter correction results in the projection images were compared with the FBCT data. FBCT was obtained from the narrow-collimated scan of the CBCT system. The collimation was made of 3 mm thick leaded blades with 3 mm opening in the vertical direction. Due to narrow collimation, the FBCT profile data contain less scatter signal than the profile data obtained from CBCT. In the experiment, we used the PMMA sheet (thickness of 60 mm) that is attached with the lead sheet (thickness of 3 mm) for measuring the scatter signal value and comparing its profile with FBCT as shown in Figure 6.
(a)
(b)
(c)
Ideally, the profile data obtained from FBCT should be almost identical to the actual primary signal; therefore, the FBCT data are used as the benchmark in comparison with the estimated data from CBCT. However, FBCT acquisition is limited to only a small strip around the center of the X-ray beam. Figure 7 shows comparison among uncorrected CBCT and corrected CBCT and FBCT profiles of the projection data acquired from the QRM-ConeBeam phantom. The log likelihood values calculated using (8) are plotted versus the number of iterations as shown in Figure 8. The log likelihood values seem to converge as the number of iterations increases. One of the stopping criteria can be monitored from the convergence of the log likelihood values. After the convergence of the log likelihood value, the reconstructed cross section images of the proposed scatter correction method should be close to the cross section image of the FBCT.
(a)
(b)
(c)
3.2. Measurements in the CT Number Section
Corrected projection images are reconstructed by the filtered backprojection method [19] using the Shepp Logan filter with the cutoff at 0.6 and the voxel size of 0.3 mm. The cross section images in the CT number section are measured and compared with FBCT data. Figure 9 shows the cross section images with and without scatter correction and their profiles comparison. All reconstructed images in Figure 9 are displayed with the window width and level () of 1500 and 1000, respectively.
(a)
(b)
(c)
(d)
The profile data comparison is plotted across the bone insert and the air hole in the CT number section as shown in Figure 9(d). The proposed profile is almost identical to the profile from FBCT as the cupping effect is reduced. The CT number values are calculated by using the average value of the water insert according to the QRM-ConeBeam phantom’s specifications [20]. Since the X-ray spectrum we used cannot clearly discriminate the density value of water and soft tissue (background), we used the soft tissue value instead. In the experiment, we measured the CT numbers of three inserts up to 10 iterations as shown in Table 3.
The accuracy of the CT number value at each insert increases as the number of iterations increases; however, once convergence is reached the accuracy of CT number value stays the same even when we perform more number of iterations. These CT number values are considered to be valid for diagnosis after convergence is achieved. However, noise is also increased after scatter correction. We evaluated the contrast value between bone and soft tissue inserts and the influence of noise by measuring the contrast to noise ratio as shown in Table 4. In Table 5, to calculate percentage of cupping, we compare the center average value with the average value of four peripheral areas using (11).
From Table 5, the cupping artifacts are significantly reduced, comparing with FBCT data; however, the small effect from beam hardening is still present even with FBCT. As the log likelihood value is decreased, the CT number values of bone, air, and soft tissue as well as the contrast values are significantly improved; however, noise is also increased. In this study, we did not perform noise suppression after scatter correction, so this would affect the CNR.
3.3. Measurement Results in the Low Contrast Sections and the PH3 Phantom
There are three sections of low contrast detectability in the QRM-ConeBeam phantom, namely, Sections A, B, and C. They have the same pattern as shown in Figure 5(b), but the density values of the inserts within each section are different. The results of three low contrast sections after 5 iterations of scatter correction are reconstructed by the FBP method using the Hamming filter with the cutoff frequency of 0.6 and the voxel size of 0.6 mm and then compared with FBCT as shown in Figure 10. We measured the CT numbers in each section according to numbers 1 to 4 as shown in Figure 5(b) and the errors were calculated by using the absolute different value of the corrected results with FBCT.
(a) Section A, FBCT
(b) Section A, without correction
(c) Section A, with correction
(d) Section B, FBCT
(e) Section B, without correction
(f) Section B, with correction
(g) Section C, FBCT
(h) Section C, without correction
(i) Section C, with correction
The measured HU values in three sections of low contrast are shown in Tables 6–8. Note that one of the reasons that the error after scatter correction in Section A seems to be a little higher than Sections B and C might be because the position of Section A is further away from the central ray than Sections B and C.
Results in three sections of low contrast indicate better image quality. The low contrast values in three sections can be discriminated; thus the contrast is significantly improved. Moreover, some inserts appear after correction. In addition to the QRM-ConeBeam phantom, we applied the proposed method to another phantom using the same parameters and the results are shown in Figures 11 and 12. The reconstructed images in these two figures are displayed with the window width and level of 900 and −100, respectively. Figure 11 shows the cross section images of the PH3 phantom along with profile data comparison. The ventricles region in the PH3 angiographic CT head phantom is significantly improved as shown in Figure 11(b). The corrected profile in Figure 11(d) is very close to the FBCT profile. Figure 12 shows the reconstructed images with and without scatter correction in the coronal and sagittal planes. Although image quality is improved, the beam hardening effect still presents in these planes.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4. Discussion
Although the log likelihood function seems to converge as the number of iterations increases, image quality after scatter correction may not be approaching FBCT exactly according to the log likelihood function. The small change in the log likelihood value after several iterations may result in overcorrection in the reconstructed images as indicated by the percentage of cupping in Table 5. For example, the percentage of cupping at 10 iterations is smaller than that of FBCT. Therefore, to ensure convergence and avoid overcorrection, only 5 iterations seem to be sufficient (Figures 10–12). In this work, we attempt to improve image quality of low contrast by applying the proposed method to the three sections of low contrast in the QRM-ConeBeam phantom. The results of three sections in Figure 10 are significantly improved; that is, visibility of some inserts in Section C can be observed. Moreover, in Table 5, the percentage of cupping artifacts in FBCT is about 10% which means that other causes of cupping artifacts besides scattering are present. One of them could be beam hardening effects, which is ignored in this study.
5. Conclusions
In this paper, we propose the scatter correction method to improve image quality of soft tissue images acquired from portable CBCT. Our proposed technique is based on estimation of X-ray scatter signals using the MLEM method and kernel modeling with Monte Carlo simulation. By benchmarking with FBCT, the scatter correction results with CBCT show significant improvement on image quality. With the QRM-ConeBeam phantom, the cupping artifacts are reduced, CT numbers of inserts are approaching FBCT, contrast is increased, and low contrast detectability becomes more apparent. Moreover, scatter correction in the reconstructed images of the PH3 angiographic CT head phantom can bring out soft tissue structures prominently. Therefore, our proposed scatter correction method has a high possibility to detect soft tissue images using portable CBCT. For future work, we will test our proposed technique on real patient data.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to thank the Faculty of Medicine, Prince of Songkla University, for research funding and feedbacks of the portable CBCT prototype.