Patient-specific scatter correction in clinical cone beam computed tomography imaging made possible by the combination of Monte Carlo simulations and a ray tracing algorithm

Abstract

Purpose. Cone beam computed tomography (CBCT) image quality is limited by scattered photons. Monte Carlo (MC) simulations provide the ability of predicting the patient-specific scatter contamination in clinical CBCT imaging. Lengthy simulations prevent MC-based scatter correction from being fully implemented in a clinical setting. This study investigates the combination of using fast MC simulations to predict scatter distributions with a ray tracing algorithm to allow calibration between simulated and clinical CBCT images. Material and methods. An EGSnrc-based user code (egs_cbct), was used to perform MC simulations of an Elekta XVI CBCT imaging system. A 60keV x-ray source was used, and air kerma scored at the detector plane. Several variance reduction techniques (VRTs) were used to increase the scatter calculation efficiency. Three patient phantoms based on CT scans were simulated, namely a brain, a thorax and a pelvis scan. A ray tracing algorithm was used to calculate the detector signal due to primary photons. A total of 288 projections were simulated, one for each thread on the computer cluster used for the investigation. Results. Scatter distributions for the brain, thorax and pelvis scan were simulated within 2% statistical uncertainty in two hours per scan. Within the same time, the ray tracing algorithm provided the primary signal for each of the projections. Thus, all the data needed for MC-based scatter correction in clinical CBCT imaging was obtained within two hours per patient, using a full simulation of the clinical CBCT geometry. Conclusions. This study shows that use of MC-based scatter corrections in CBCT imaging has a great potential to improve CBCT image quality. By use of powerful VRTs to predict scatter distributions and a ray tracing algorithm to calculate the primary signal, it is possible to obtain the necessary data for patient specific MC scatter correction within two hours per patient.

Cone beam computed tomography (CBCT) image quality is limited by a large scatter contamination in the projection images. The scattered photons reduce contrast and introduce the so-called cupping artefact in the reconstructed CBCT images [1].

Several different approaches have been taken to compensate for the scattered photons in CBCT imaging, as presented in an extensive review by Rührnschopf and Klingenbeck [2,3]. The scatter compensation techniques either concern the prevention of scattered photons from being detected on the CBCT scanner, or try to compensate for the scattered photons in the projection images through data processing after the acquisition of the CBCT projection images.

Different ways of performing post-acquisition scatter compensation have been investigated over the recent years. The simplest approaches assume the scatter contamination to constitute a constant fraction of the total signal in each projection image, and this technique is currently used to compensate for scattered photons on the Elekta XVI CBCT imaging system [4]. More sophisticated scatter models are continuously developed in the pursuit of better image quality in CBCT imaging. One of the promising techniques is the use of Monte Carlo (MC) simulations to predict the scatter distributions in CBCT projection images based on the physical properties of the imaging object itself.

The use of MC-based scatter correction in CBCT imaging is limited by the extensive simulation time required. In one of the first MC studies of scattered photons in CBCT imaging, Jarry et al. reported a simulation time of 430 hours [5]. To reduce this simulation time, they suggested a reduction of the number of simulated projection images, as well as an increase in the size of the phantom voxels and/or simulated detector pixels. This suggestion was adopted in a more recent study by Poludniowski et al. [6], who used a small number of fixed detector points to reduce the required simulation time. To recover the full scatter distributions, interpolation between the fixed detector points was used. This approach was justified by the assumption that spatial frequencies in scatter distributions are low – an assumption which may or may not be valid based on the specific imaging geometry as pointed out by Mainegra-Hing and Kawrakow [7].

Recent advancements in computing power as well as MC algorithms have made full scale MC simulations of scattered photons in CBCT imaging possible in a time frame not prohibitive for clinical applications [7,8]. This work presents a MC model to predict the patient-specific scatter distributions in CBCT imaging. No assumption is made regarding the spatial frequency of the scatter distributions, nor is the phantom voxel size increased. To allow easy comparison and calibration against clinical CBCT projection images, the CBCT detector signal due to primary photons is calculated as well as the detector signal due to scattered photons. For efficiency, the primary signal is calculated using a ray tracing algorithm which has no statistical uncertainty on the calculation result.

Material and methods

CBCT projection images based on primary photons were calculated using a ray tracing algorithm, while MC simulations were used to calculate scatter images. Projection images were calculated in angular steps of 1.25° covering a full revolution, and CBCT reconstructions were made using an FDK type reconstruction algorithm [9]. Calculation phantoms were made from planning CT images of three patients – one brain, one thorax and one pelvis scan – who had previously received radiotherapy (RT) in our clinic.

Calculation geometry

The calculation geometry was specified to resemble an Elekta XVI CBCT imaging system operated with the S20 and F0 filters [4]. This setting operates without a bowtie filter, and with the kV detector panel on the central axis of the kV x-ray tube. The S20 filter allows a cylinder of 27 cm diameter to be reconstructed, and uses an x-ray cone and fan beam angle of 15.75°. The detector panel operates with 512 × 512 pixels of 0.8 × 0.8 mm², totalling an active detector size of 40.96 × 40.96 cm². The x-ray field is thus slightly larger than the active detector region. The source to axis distance of the XVI unit is 100 cm, and the source to detector distance is 153.6 cm.

The x-ray source used to calculate the primary and scatter signal at the CBCT detector position was a 60 keV monoenergetic and uniform x-ray source. This energy was chosen according to Spezi et al. who reported the mean energy of the XVI unit to be 61.2 keV when using the F0 and S20 collimators in combination [10]. Using a monoenergetic x-ray source further simplified the calculations using the ray tracing algorithm.

MC simulations

Scattered photons in the CBCT projection images were simulated using the recent EGSnrc user code egs_cbct [7,8]. The accuracy of EGSnrc-based simulation of scattered photons in CBCT imaging has been shown in previous studies comparing MC simulations with direct measurements of scattered photons [5,11]. Scatter signals were normalised against a simulated projection image with no object between the source and detector (referred to as a blank scan).

Since electrons deposit their energy locally and do not contribute to the detector signal in CBCT imaging, MC calculations were performed simulating only photon transport to improve the simulation efficiency. Photon cross sections from the XCOM database were used, as provided by the EGSnrc software. Bound Compton scattering was simulated without rejections at run time (the norej option in egs_cbct), and Rayleigh scattering was turned on. Variance reduction techniques (VRTs) were used to improve the scatter calculation efficiency in the MC simulations [7]. In particular, particle splitting using fixed splitting numbers was used in combination with an edge preserving smoothing algorithm, path length biasing, and delta transport for photons not aimed at the detector.

The detector signal was simulated as air kerma, using forced detection for efficient recording of the detector signal. Scatter distributions were calculated to within 2% statistical uncertainty for all 288 simulated projection images and all three patient phantoms.

Ray tracing algorithm

The used ray tracing algorithm has been developed specifically to be able to read the same voxelised phantoms as any C++-based EGSnrc user code. The ray tracing algorithm provided the attenuation between the source and each detector pixel, and could be converted to the detector kerma through an exponential transform and multiplication by the blank scan from the MC simulations. No sampling was involved in the ray tracing calculation, and hence no statistical uncertainties needed to be considered in relation to the primary signal.

CBCT reconstruction

The calculated primary and scatter signal in the CBCT projection images were added to produce the total calculated CBCT projection images. To investigate the influence of scattered photons on reconstructed CBCT image quality, three different calculated CBCT images were reconstructed for each of the patient phantoms. The first CBCT reconstruction was based on the total signal without any scatter compensation. The second CBCT reconstruction used a constant scatter subtraction, where 20% of the total signal in the central detector pixel was subtracted from all pixels in each projection image prior to reconstruction. This subtraction resembles the default scatter correction in the Elekta XVI reconstruction software, as currently used in our clinic. The last CBCT reconstruction was based on the projection images with only primary photons recorded. These images show the potential image quality of CBCT images if an ideal scatter correction method were available.

Hardware

MC simulations were performed on a medium sized computer cluster made from 24 Intel Xeon X5650 2.66 GHz CPUs. Each CPU allows 12 threads to be processed in parallel and hence 288 projection images could be simulated in parallel. The ray tracing algorithm and FDK reconstruction algorithm were run on the cluster front end PC equipped with 2 AMD Opteron 254 2.8 GHz processors. The 2011 purchase price of the cluster and front end PC was €31 500.

Results

Calculated projection images

Figure 1 shows an example of the calculated CBCT projection images for the three patient phantoms. The calculated scatter distributions vary slowly over the projection images. To quantify this variation, pixel values in each scatter image were binned in histograms (Figure 1G–I). The ratio between the minimum and maximum value of the 95% most central pixel values in the histograms was used as a measure of the scatter variation. This ratio was found to be 1.4 for the brain CBCT projection image (Figure 1d), 1.6 for the lung CBCT projection image (Figure 1e), and 1.5 for the pelvis CBCT projection image (Figure 1f).

Figure 1. Example CBCT projection images of a brain (A, D), thorax (B, E) and pelvis (C, F) scan. Primary images (A–C) were calculated using the ray tracing algorithm, and have been log transformed. Scatter images (D–E) were calculated by MC simulations, and the images show air kerma calculated at the detector plane. Histograms (G–I) were calculated from both the primary and scatter signal, and normalised to unit area in the displayed air kerma range. For the brain scan, the primary photons in the air region surrounding the head gave rise to a signal in the histogram at higher air kerma values than shown in (G).

Calculation time

The required time to simulate the scatter distributions in the CBCT projection images was 23 minutes for the brain scan, 50 minutes for the thorax scan and 120 minutes for the pelvis scan (Table I). The calculation time for the ray tracing algorithm only depends on the number of voxels in the calculation phantom and the detector resolution. The present phantom calculations required 40 minutes to calculate the primary signal in all 288 projections for any one of the three calculation phantoms. FDK reconstructions were completed in 13 minutes per scan.

Table I. Scatter simulation time.

	Brain CBCT	Thorax CBCT	Pelvis CBCT
Nhist	5 × 10^4	5 × 10^4	10^5
σmax	1.3%	1.5%	2.0%
TCPU	23 min	50 min	120 min

Simulated number of histories (N_hist), highest statistical uncertainty (σ_max) and CPU time (T_CPU) for the MC simulations of CBCT scatter distributions.

Reconstructed CBCT image quality

Figure 2 shows the reconstructed CBCT images for all three calculation phantoms. For each CBCT scan, three images with different scatter contamination are shown. To allow comparison between the reconstructed images, level and windowing were calculated as the CBCT slice mean and three times the standard deviation, respectively.

Figure 2. CBCT images reconstructed from the calculated brain (A–C), thorax (D–F) and pelvis (G–I) scans. Reconstructions are based on projection images calculated in steps of 1.25°, and reconstructed using an FDK type algorithm. For each scan, three images are reconstructed. The first reconstruction is based on the total signal from both primary and scattered photons (A, D, G). The second is based on a constant scatter subtraction of 20% in each projection image prior to reconstruction (B, E, H), and the third is based on the primary photons only (C, F, I). The displayed CBCT slices are examples, and not chosen centrally from the datasets.

The calculated brain CBCT images (Figure 2a–c) show little effect of the scatter compensation. Gradients around the bony structures are steeper in the scatter-free CBCT image compared to the CBCT images based on the total signal with and without the constant scatter subtraction. No significant cupping artefact was observed even in the CBCT image reconstructed based on the total signal.

In the thorax CBCT images (Figure 2d–f) steeper gradients are found around the bony structures and between different tissues in the mediastinum when scatter correction is applied. The CBCT image based on the constant scatter subtraction shows some improvement compared to the reconstructed image without scatter compensation, but the largest difference is observed when the scattered photons are removed completely prior to reconstruction. Gradients around major vessels in the lung tissue are steeper in the scatter-free image, compared to the CBCT images with and without the constant scatter subtraction applied.

The pelvis CBCT images reconstructed from the total signal with or without the constant scatter subtraction (Figure 2G–H) show a very pronounced cupping artefact. This cupping is removed when only the primary photons are used for the reconstruction (Figure 2i), and soft tissue details can be discriminated when they are not obscured by cupping.

Discussion

With a required calculation time of no more than two hours per patient, this study showed that clinical use of MC-based scatter correction is feasible. Previous limitations of prohibitively long simulation times [5] or use of sparse views and heavy interpolation on the detector panel [6] were circumvented by use of the latest MC software with powerful VRTs for scatter calculations in combination with a ray tracing algorithm for calculations of the primary signal. With the combined primary and scatter signal it is possible to calibrate the simulated detector signal against a set of clinical projection images. This calibration is what allows the simulated scatter distribution to be subtracted from clinical projection images to improve clinical CBCT image quality.

Calculated projection images

The simulated detector signal due to scattered photons was found to vary with the position on the detector panel in each projection image. This variation must be considered before deciding what statistical uncertainty is acceptable for the MC simulations of the scatter distributions, since smaller statistical uncertainties come at the cost of longer simulation times. In the present investigation, the uncertainty was more than one order of magnitude smaller than the variation between the max and min values of the scatter distributions. This was deemed sufficient to resolve whatever structure was present in the scatter distributions. Furthermore, the number of simulated projection images should reflect the smooth structure, and the employed 288 projections are similar to the number of views used in previous studies [5,6]. The ultimate test of whether the MC simulations are sufficiently accurate should be based on the image quality of scatter corrected clinical CBCT images of known phantoms.

Calculation time

The different simulation times required to reach the 2% statistical uncertainty in the scatter distributions are explained by different x-ray attenuation of the three calculation phantoms. The brain scan has the shortest path length, and hence the shortest required simulation time. For the thorax scan, the increased volume of the body required a longer simulation to reach the desired uncertainty. The pelvis scan required the longest time, caused by an increased density of the anatomy compared to the thorax scan. The fact that the pelvis scan required a longer simulation time compared to the brain and thorax scan showed that considerable self-absorption of scatter in the patient anatomy is present. This was also reflected in the increased number of primary histories used for the pelvis scan compared to the thorax and brain scans (Table I).

Ray tracing calculations of the primary signal required 40 minutes to complete all 288 projection images in the present investigation. This time could be reduced by optimising the code for faster calculations, or by taking advantage of the computer cluster to calculate one projection image on each thread available. The FDK reconstruction time was 13 minutes, and it is noted that both ray tracing and FDK reconstruction can be performed in less than 10 seconds on a GPU [12].

Reconstructed CBCT image quality

While it is beyond the scope of this investigation to fully account for the differences in image quality between the three reconstructed CBCT images with different scatter subtractions, a few comments are in place. When no scattered photons were included in the projection images prior to reconstruction, higher contrast was found between regions of varying density. This was mostly noted for the pelvis CBCT images, where the cupping artefact obscured all low-contrast objects in the reconstructed CBCT images with scattered photons in the projection images. In the lung CBCT image, the same trend was observed where removal of scattered photons yielded better homogeneity and steeper gradients between tissues of different densities. For the brain CBCT images little effect of the scatter correction was observed. The severity of scatter contamination has been reported to increase with increasing object size [1,11], which is in agreement with the present investigation where the scatter artefacts are most severe in the larger regions of the body. Furthermore, the largest magnitude of scatter in the brain scan was found to be outside of the patient head, and thus with no impact on the region of interest in the reconstructed image.

Clinical considerations

In order to use the calculated scatter distributions to correct clinical CBCT images, it is necessary to normalise the calculated signal against the measured signal of the clinical CBCT projection images. This is where the total images are needed, and the reason why the MC simulations should be combined with a ray tracing algorithm to provide all the necessary data for clinical CBCT scatter correction efficiently. The isocentre for RT should be used to position the isocentre of the calculation geometry to ensure that calculations are performed in the same reference system as the clinical images are acquired.

Due to the slowly varying structure in the scatter distributions, it is our belief that only one calculation of the scatter signal is required to correct all clinical CBCT images acquired for the same patient during fractionated RT. A proper alignment of the calculated projection images against the clinical projection images is paramount, but the anatomical changes which occur during the RT treatment are believed to cause only minor changes in the scatter distributions. This hypothesis should be tested before clinical implementation. In cases where large anatomical changes are observed during fractionated RT, a resimulation of the treatment plan based on a new CT image is likely to be performed. This CT image could be used for a new calculation to account for changes in the scatter distributions.

Before clinical scatter correction in CBCT imaging based on the proposed combination of a ray tracing algorithm and MC simulations can be realised, the model must be thoroughly validated through phantom studies. This is required to ensure that the MC-based scatter corrected CBCT images represent the imaging objects better than CBCT images made with the currently available scatter correction methods, and will thus include a thorough investigation of the CBCT image quality. During our work, we have observed that the CBCT detector panel on our linear accelerator is saturated whenever there is a direct line of sight between the source and detector. This saturation and the potential non-linear detector response need to be included properly in the scatter correction model to perform the correct normalisation between the calculated and the clinical projection images. Other clinical factors such as detector scatter glare, beam hardening and the inclusion of a bow tie filter for some CBCT scans should also be considered in future work to ensure a successful application of the calculated data to clinical CBCT images. However, these effects are small compared to that of efficient scatter compensation, and none of the effects are expected to cause a major increase in the required simulation time. Beam hardening is readily included in the simulations by exchange of the monoenergetic source by a realistic x-ray spectrum. Scatter glare can be deconvolved as shown by Poludniowski et al. [13]. Bow tie filtration can be included in a source simulation which is a one-off simulation that does not extend the calculation time needed for MC-based scatter correction.

With the use of a medium sized computer cluster and the latest advancements in VRTs, it is possible to use MC simulations to predict patient-specific scatter distributions in CBCT imaging within a time frame not prohibitive for clinical applications. No assumptions were made on the spatial distribution of the scattered photons, nor was the patient geometry downsampled to reduce the required simulation time in the present investigation. The detector signal from primary photons was calculated using a ray tracing algorithm. In combination, the MC simulations and ray tracing algorithm provided the total detector signal required for calibrations against a clinical CBCT imaging system within two hours on our computer cluster.

Declaration of interest: This work is supported by Elekta and CIRRO – The Lundbeck Foundation Center for Interventional Research in Radiation Oncology and The Danish Council for Strategic Research. RST acknowledges financial support from Oticonfonden. The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.