TIGRE-VarianCBCT for on-board cone-beam computed tomography, an open-source toolkit for imaging, dosimetry and clinical research

We presented TIGRE-VarianCBCT, an open-source toolkit Matlab-GPU for Varian on-board cone-beam CT with particular emphasis to address challenges in raw data preprocessing, artifacts correction, tomographic reconstruction and image post-processing. The aim of this project is to provide not only a tool to bridge the gap be-tween clinical usage of CBCT scan data and research algorithms but also a framework that breaks down the imaging chain into individual processes so that research effort can be focused on a specific part. The entire imaging chain, module-based architecture, data flow and techniques used in the creation of the toolkit are presented. Raw scan data are first decoded to extract X-ray fluoro image series and set up the imaging geometry. Data conditioning operations including scatter correction, normalization, beam-hardening correction, ring removal are performed sequentially. Reconstruction is supported by TIGRE with FDK as well as a variety of iterative algorithms. Pixel-to-HU mapping is calibrated by a Catphan TM 504 phantom. Imaging dose in CTDIw is calculated in an empirical formula. The performance was validated on real patient scans with good agreement with respect to vendor-designed program. Case studies in scan protocol optimization, low dose imaging and iterative algorithm comparison demonstrated its substantial potential in performing scan data based clinical studies. The toolkit is released under the BSD license, imposing minimal restrictions on its use and distribution. The toolkit is accessible as a module at https://github.com/CERN/TIGRE.


Introduction
The clinical application of on-board cone-beam computed tomography (CBCT) in image-guided radiotherapy (IGRT) has significantly increased the accuracy of tumour localization, bringing substantial improvement to clinical outcome [1,2]. Currently, on-board CBCT has become a critical part of IGRT treatment, and for each patient scans are recurrently performed over the treatment course as per institutional routine [3]. These intra-fraction CBCT images as well as inter-fraction series offer significant online and longitudinal information towards precise and adaptive radiotherapy (ART) [4][5][6].
Even though on-board CBCT has become the golden standard for patient setup verification [7], its image quality is relatively coarse and far below the reference level in diagnostic and simulation CT. This adversely limits its clinical role and hinders to unleash its potential, especially in online ART which requires high accuracy and precision in online dose re-calculation [4,8].
The poor imaging quality of on-board CBCT can be partially attributed to proprietary imaging programs embedded on linac systems. Although these programs have been task-specifically optimized to pipeline the imaging workflow, the overall imaging chain is generally obsolete and lags behind the pace of technology development. For example, while it is well known that iterative reconstruction outperforms FDK especially in low dose scenarios [9][10][11], FDK is still the predominant algorithm in on-board CBCT.
Nevertheless, most of these open-source tools are focused on the core process of image reconstruction, leaving other processes in imaging chain rarely addressed. For on-board CBCT, the embedded imaging programs are fully closed, and scan data are implicitly described. The result is a lack of tools available for users to revisit scan data to explore new opportunities in on-board CBCT imaging. It has been a dilemma of new researchers and students that even with scan data at hand, the gaps from raw scan data of X-ray exposures to clinical images in HU numbers still remain unfilled yet.
This work is to present an open-source and versatile toolkit TIGRE-VarianCBCT that allows to perform off-line CBCT imaging using scan data from Varian on-board CBCT. The particular emphasis as well as the key originality of the proposed toolkit is to fully address challenges in raw data preprocessing, artifacts correction, tomographic reconstruction and image post-processing.
To make it user friendly, cross-platform, fast to run and communitysupported, the toolkit was developed in Matlab as an independent and task-specific module to TIGRE. The latest version is open to access via https://github.com/CERN/TIGRE with documentation, installation and demos.
As an in-depth introduction to the open-source toolkit, this paper serves as a step-by-step roadmap to walk through the imaging chain from raw CBCT scan data to volumetric images in HU numbers. The rest of the paper is organized as follows: the geometry and imaging chain of on-board CBCT are briefly described, and then the toolkit structure and key algorithms are detailed; to validate and illustrate the streamlined workflow, a demo project is shown with a code snippet template; to preview the easy access this toolkit provides to low-dose CBCT research, three proof-of-concept case studies are presented; the scalability, limitations and future vision are discussed in the last section.

Geometry of On-board CBCT
The overview of on-board CBCT on a typical Varian linac is shown in Fig. 1. The kilo-voltage (kV) X-ray source and flat panel detector (FPD) are mounted on two positioning unit (PU) arms attached to the linac gantry. The relative vertical and lateral positions between the FPD between X-ray source are adjustable by PU arms to allow for different fieldof-view sizes and patient clearance. During CBCT scan, a patient lying on a flat-top couch is delivered to the imaging centre. Then, the X-ray source and FPD are fully extended to planned positions, and the cone-shaped beam irradiated from the X-ray source traverses the patient following the Beer-Lamber law and is captured by the FPD. At the same process, the X-ray source and FPD rotate simultaneously around the isocentric axis in a circular trajectory to acquire X-ray radiographic images at different sampling points over the scan range.

Imaging chain and key algorithms
According to the vendor-provided technical guide [21], the overall imaging chain of on-board CBCT from raw X-ray fluoro image series to volumetric CBCT slices can be summarized into four major modules as illustrated in Fig. 2.

Imager calibration
The first module, imager calibration, is part of routine quality control (QC) of the imaging system rather than CBCT scan, but it has a substantial impact on image qualities. The key processes include FPD dynamic gain or response uniformity correction, dark noise or offset subtraction, defect pixel removal, etc. During scan data exporting, the calibration files acquired in the latest valid calibration are exported along in the same process.
Since imager calibration is automatically performed as default in CBCT scan, there was not any imager calibration specific function in the current toolkit. Despite this, to mitigate steaking/ring artifacts induced by out-of-date calibration -which is not rare in clinic, denoising in both projection and image domains were devised in the following modules.

Data conditioning
Data conditioning is to preprocess raw scan data and prepare linear integrals of attenuation for image reconstruction. The preprocessing includes corrections for five physical phenomena in data acquisition: Xray scatter, angle-dependent geometrical shifts in gantry rotation, beam shaper induced field non-uniformity, pulse-to-pulse tube output fluctuation and poly-energetic spectrum induced beam-hardening. These operations are essential to enhance overall image qualities, especially in suppressing shading/cupping artifacts and improving HU accuracy.
Scatter. To mitigate scattered X-rays, a 1D anti-scatter grid (ASG) is mounted in front of the FPD as shown in Fig. 1(c). In addition to this hardware-based technique, a software-based scatter correction scheme is also applied. Briefly, it is an adaptive scatter kernel superposition based iterative algorithm [22] that takes into consideration of object thickness variations, object edges and ASG effects in Eq.(1) and the point-spreading effect in FPD in Eq.(2) as.
where γ is the proportionality coefficient, τ is the water-equivalent thickness, I p is the primary signal, i is the thickness group index, R i is the mask for thickness group i, g i is the form function of scatter kernels for thickness group i as shown in Fig. 3(b), c e i is the edge-modified amplitude factor for thickness group i, and k is the anti-scatter grid response function; ̅̅̅̅̅̅̅̅ ̅ where n c is the empirical normalization factor that is equal to 0.04, and (a 1 , a 2 , a 3 , a 4 , a 5 ) are the fitting parameters. Note that all the parameter values are saved in the accompanying calibration file (SC/Calibration. xml). As implementation of this algorithm, the scatter correction code in this toolkit first parsed the calibration file and automatically imported all the parameter values. The iteration number for scatter correction was 8, which was adequate for convergence. A representative X-ray image and the estimated scatter-to-primary map were shown in Fig. 3(a)&(c). Besides, since scatters are mostly low frequency signal, input images were first downsampled by 12 as default to achieve commensurate benefit in running time. Also, to further speed up the procedure, GPUbased fast Fourier transform was also activated by default.
Geometric shifts. As shown in Fig. 1(c), a Titanium foil filter as thin as 0.89 mm is mounted behind the tube window to absorb low-energy photons. Along the foil filter, a beam shaper or bowtie filter is placed along the beam path to shape the beam profiles. The filter is made of aluminium and its profile thickness varies from 2 mm to 28 mm. The use of bowtie brings several benefits to on-board CBCT, such as reducing patient skin dose and lateral X-ray scatter. Varian on-board CBCT provides two bowties, one for scans with detector-at-centre (Full-Fan mode) that applies to protocols of Thorax/Pelvis/PelvisLarge and the other one for scans with detector-off-centre (Half-Fan mode) that applies to protocols of ImageGently/Head/Spotlight. In this toolkit, in order to compensate for the beam field nonuniformity induced by filters, object-attenuated X-ray intensities were normalized by airscan intensities acquired in routine QC as: where I represents the X-ray intensity in patient scan, I air represents the X-ray intensity in airscan, and F c represents a factor that accounts for Xray tube output fluctuation. During data acquisition, beam pulse output is monitored by an air chamber close to the tube window as in Fig. 1(c). Herein the pulse-topulse fluctuations were compensated by a scaling factor as: where R and R air represent the air chamber readings in patient and air scans respectively. During gantry rotation for CBCT data acquisition, the heavy load and Illustration of scatter estimation in this toolkit: (a) raw X-ray intensity map extracted from a real patient CBCT scan, where X-ray field non-homogeneity is caused by the bowtie; (b) kernel used for fast superposition-based scatter estimation, where the anisotropic effect is caused by the anti-scatter grid; (c) estimated scatter-to-primary ratio (SPR) map of (a), where severe scattering is induced by denture implants and the skull.
gravity effect induce subtle variations in relative imaging geometry. This angle-dependent geometrical instability will lead to slight shifts in nonuniform beam profiles, and if left uncorrected a shading pattern known as crescent artifacts are present. To improve the accuracy of the signal intensities measured and normalize the projections, a sector-based rotating air norm scheme was used as.
where β represents the projection angle, (β 1 , β 2 ) the control points in airscan calibration, and (ω, 1-ω) the bilinear interpolation weights. The vendor recommends to perform monthly air-scan calibration. During air-scan calibration, the gantry rotates in clock-wise (CW) and then counter-clockwise (CCW) directions respectively, and X-ray intensity maps are acquired at intervals of 36 degrees, i.e., 10 equal-spaced control points. The latest air-scan maps are exported alongside with every CBCT scan data.
Beam-hardening modelling & correction. Incident X-ray beams from the tube are poly-energetic over a spectrum rather in single energy photons. For X-ray beams used in kV on-board CBCT, low-energy photons are more easily absorbed, which is known as beam hardening. If this effect is left uncorrected, HU value accuracy will be compromised with cupping profiles, the most common artifacts, in thick objects.
In this toolkit, a classic spectrum modelling algorithm [23] was implemented and the beam-hardening correction shown in Fig. 4(a) generally fell into 5 steps: • Step 1): the factory-defined calibration file (ASC/Calibration.xml) was first parsed to load necessary data. The data included the normalized X-ray spectrum, coil filter material and thickness, bowtie material and shape parameters, attenuation coefficients of typical Step 2): the X-ray spectrum filtered by the coil filter was modelled by the Beer-Lamber law as where μ f is the attenuation coefficient of the coil filter, L f is coil filter thickness, S 0 is the X-ray tube spectrum, and S 1 is the filtered spectrum. • Step 3): to take into account of the cone beam geometry and bowtie effect, a pixel-by-pixel look-up table (LUT) of the bowtie-transvers path length was built as in Eq. (8), and a group of bowtie-filtered X-ray spectra were computed as in Eq. (9): where l is bowtie transverse path length, (u,v) is detector unit coordinate, S 2 is bowtie transvers path length dependent spectrum.
• Step 4): a logarithmic attenuation LUT w for water over the thickness of 0 to 600 mm were computed in Eq.(10), and therefore an analytical LUT p for logarithmic projection signal intensity were derived. The LUT p were dependent on and water thickness and the bowtie transverse path length which was tabulated in LUT b . The ideal projection signal intensity over water thickness can be precomputed in Eq. (11) as where t is water thickness. The beam hardening effect of concave projection signal trend is illustrated in Fig. 4(d). • Step 5): a group of 1st order models fitted by the logarithmic projection signals over the first 6-cm thick water were built, and the concave non-linear projection signals were remapped to its linear response values as in Eq. (12), where P is the measured projection signal, func is the linear remapping function, and P' is the corrected projection signal. A representative profile of the CTP486 test module in Catphan TM 504 was shown in Fig. 4 (e). The beam hardening correction detailed above was a pure software operation, and in current version only water was accounted as the calibration material. Due to the computing intensive nature, the beam hardening correction process was relatively slow and took about 10 min to complete for a typical dataset with 500 frames of X-ray acquisitions. Since beam hardening model data are applicable to the same scan protocol, the process can be accelerated by loading precomputed models.
Auxiliary corrections. In addition to the abovementioned procedures, the toolkit proposed two more corrections towards artifacts suppression. One was to remove redundant frames that were over-sampled due to gantry acceleration and deceleration in starting and finishing phases of data acquisition. The defined threshold in angular interval was 95 %, that is, if the angular position between two neighbouring frames was smaller than 95 % of the nominal sampling interval, the next frame would be removed from the scan data sequence. The other correction was to apply a lateral trimmed mean filter or median filter along the uaxis in projection domain. While the spatial resolution was reduced to some extent, the effect of noise reduction and ring artifacts removal was significant and therefore was recommended to perform.

Image reconstruction
The third module is image reconstruction, the core process calculating tomographic slices from preprocessed projections. In this toolkit, projection and geometry data were structured in the TIGRE style to fully utilize the provided various algorithms. As the reconstruction backbone, TIGRE provided not only classic FBP/FDK but also 13 iterative algorithms. These iterative algorithms generally fell into 5 classes: gradient descent-based algebraic algorithms (SART, OS-SART-SIRT), statistical algorithm (ML-EM), Krylov-subspace based algorithm (CGLS), total variation regularization -proximal based algorithms (FISTA, SART-TV) and total variation regularization -POCS based algorithm (ASD-POCS, OS-ASD-POCS, B-ASD-POCS-β, PCSD, AwPCSD, Aw-ASD-POCS). These algorithms were out of the box with well-encapsulated input and output interfaces.

Post-processing and auxiliary functions
The last module is to post-process the reconstructed slices in image domain. It consists of three major features: image denoising, field-ofview (FOV) masking, and pixel-to-HU mapping.
In this toolkit, slice-by-slice median and mean filtering were devised, which was adequate to reduce random noise in most cases. For sparseview scans, 3D TV denoising was recommended to further suppress streaking artifacts. For FOV masking, both geometric FOV and selfdefined masks were available to crop the 3D image matrix. An image domain ring artifacts removal method was also provided as an alternative to lateral filtering in projection domain. The method was based on polar-cartesian coordinate transformation, and rings could be easily identified and filtered out in polar coordinate system. For pixel-to-HU mapping, we adopted the same method as in the technical guide by scanning a vendor-provided Catphan TM 504 (The Phantom Laboratory) calibration. As shown in Fig. 5, there were 10 cylindrical inserts made of 9 representative materials, and a ROI was drawn within each insert to calculate the mean pixel value. Then, a linear mapping model was fitted from the pixel-HU pairs. This feature was introduced mainly to facilitate extensive HU-based application, such as cross-modality registration and radiotherapy dose calculation.
In addition, two important auxiliary features were provided. One was to export volumetric images as compact nrrd data. This facilitated interactions with medical image analysis tools, such as 3D-Slicer and MITK. The other was to calculate imaging dose index. For Varian onboard CBCT, imaging dose of each scan was defined in CT Dose Indexweighted (CTDIw) as: • N proj (13) where DF norm represents the normalized dose conversion factor per 100 mAs and per projection, N proj represents the projection number [21]. Note that the value of DF norm , which could be calibrated with a 16 cm/ 32 cm CTDI phantom, was stored in the associated scan logfile (scan. xml). Note that TIGRE itself also provided some smart tools that were off-the-shelf and very helpful, such as slice-by-slice plotting, noise simulator and quality metrics.

Streamlined workflow
The streamlined workflow from data importing to image analysis is illustrated in Fig. 6. The end-to-end workflow generally fell into three steps.
The first step was to load raw scan data and perform data conditioning. All the processes were encapsulated in a single VarianData-Loader function. Within VarianDataLoader, raw scan data including Xray exposure image series, geometry look-up tables and calibration files were first decoded by specific functions, and then fed into the data conditioning pipeline. It is worth noting that all the above-mentioned processes were encapsulated as individual functions. This provided substantial convenience for users to develop task-specific correction algorithm. Also, if a CUDA-compatible GPU was found, most of the computing intensive operations would automatically switch to run on GPU to improve efficiency. After this step was finished, the projection data and geometry parameters were ready for TIGRE reconstruction. The second step was to select a TIGRE-provided reconstruction algorithm, pass the data as input and run. Note that while all the iterative algorithms in TIGRE provided default hyperparameter values, users should keep in mind that these values were for reference only, and all the hyperparameters should be well fine-tuned to maximize the performance. The third step was to use tools and metrics for image postprocessing and analysis. This was a free end allowing users to explore.
In addition, this toolbox contained a demo to illustrate how to use and help pages for all the functions. In help pages, technical reference for original algorithms, detailed explanation on every parameter were included for advanced options as well as the contact of the source code contributor for any feedback.

Implementation examples
For Varian TrueBeam-series linacs, the default source-axis-distance (SAD) and source-imager-distance (SID), vertical distances to source and FPD, are 100 cm and 150 cm respectively. The SAD and SID values can be defined in either treatment plans or control panel to allow for suitable field-of-view (FOV) and patient clearance. The pixel area of the FPD (PaxScan 4030CB, Varex Imaging) is about 40 cm in width (u-axis) and 30 cm in height (v-axis). When the central beam projects exactly to the detector centre, the scan FOV is about 26 cm in diameter, which is usually enough for paediatric patients as well as head and neck (HN) sites. In this case, a short scan range of 200 degrees known as Full-Fan Mode is adequate to acquire sufficient projections, and therefore is recommended as the default trajectory. For thorax, abdomen and pelvis sites, the FPD can be shifted along the u-axis by an offset as large as 16 cm to maximize the FOV diameter to about 46 cm. This scan mode is known as Half-Fan Mode, and a full scan range of 360 degrees is required to acquire adequate projections.
Depending on the combination of source voltage, detector fan type and angular scan range, seven imaging protocols are pre-configured by the vendor as standard protocols for clinical application -Head, Image Gently, Short Thorax, Spotlight, Thorax, Pelvis and Pelvis Large. In addition, 4D CBCT is a licensed feature.
As for reconstruction, the vendor-embedded algorithm in on-board CBCT is the classic FDK, and iterative reconstruction (iCBCT) [11] is a licensed feature. To validate the imaging chain and exhibit the versatility of the toolkit, four case studies were present below by means of offline imaging. The raw data were from scans in various protocols for either clinical treatment or quality assurance at our institution (XXX, XXX), and the usage had been approved by IRB.

Demo study: template and validation Motivation
To demonstrate how to use the toolkit and validate its performance, a typical code snippet was given below illustrating the easiness from loading CBCT scan data to obtain clinical HU images. The sample snippet could be utilized as a template for new users to get started with their own tasks.

Data, method and code snippet
We collected raw on-board CBCT scan data of three patients treated at our institution. The patients were scanned on the same TrueBeam linac in protocols of Head, Thorax and Pelvis Large respectively. The entire workflow of raw data loading and conditioning, FDK-reconstruction, image denoising and HU mapping were performed using the following code snippet.

Results
The typical reconstructed images of the Head, Thorax and Pelvis Large scans are shown in Fig. 7. The overall qualities of the heads displayed in a wide 2200-HU display window in Fig. 7(a)&(g) were comparable, the tissue contrast was similar and the temporal and occipital bony structures were clearly reconstructed with sharp boundaries. While the difference in details and HU inconsistency could be identified in a narrow 500-HU display window, the noise level in both images were acceptable with negligible artifacts. The profiles in Fig. 7(m) overlapped in most part except in air, indicating that the vendor-designed program applied a HU cut-off at − 1000.
Compared with Head scan images, the reconstructed images of Thorax in Fig. 7 Fig. 7

(c)&(f)&(i)&
(l) exhibited severe streaking artifacts, smashed bony boundaries and compromised tissue contrast. This is mainly attributed to patient respiratory movement during data acquisition, which took at least 60 s to complete a full-arc scan. It also indicates the clinical significance in effective patient respiratory management for lung and pelvis tumour treatment. Besides, HU inconsistency and profile discrepancies in Fig. 7 (n)&(o) can be identified, especially in bone regions. This can be attributed to deviations in image chain and implementation between the vendor programme and this toolkit, especially in the processes of detector response calibration, signal denoising, polychromatic X-ray spectrum modelling and beam hardening correction, reconstruction, et al. The technical issues and pertinent hyperparameter values are not fully addressed in detail, which may somehow lead to HU shifts.

Case study A: image quality comparison of different Head scan protocols Motivation
There were three vendor-provided protocols to perform head scans: standard Head, Image Gently, and Spotlight. Despite different settings in source voltage, these protocols were all designed in Full-Fan mode (200degree short-arc scan range) rather than Half-Fan (360-degree full-arc scan range). It comes as common sense that given the same sampling rate, the larger scan range is used, the better image quality is achieved. Therefore, it is not rare that some institutions have a self-designed fullfan head scan protocol as an in-house option to acquire high-quality CBCT images. This protocol was usually derived from the default head protocol (short-arc) and typically entitled as High Quality to provide high-quality guidance images for SRS/SRT treatment.
While several protocols were available for therapists to select for head scans, the optimal option in terms of quality-to-dose balance was unclear yet. In this case study, we compared image qualities and imaging doses of the four abovementioned protocols to indicate the optimal trade-off between quality and dose.

Data and method
The head section of an Alderson Radiation Therapy phantom (Radiology Support Devices Inc.) was scanned on a Varian Edge linac at our institution with protocols of standard Head (100 kVp, short-arc), High Quality (100 kVp, full-arc), Image Gently (80 kVp, short-arc), and Spotlight (125 kVp, short-arc). Tomographic images were reconstructed by FDK, and then denoised by intra-slice median filtering.
Image quality was evaluated in an ROI-based approach-one ROI was contoured within maxillary sinus to represent air and the other within brainstem to represent soft-tissue, and then signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) were calculated as following: where (μ, σ) represent the mean and standard deviation over the ROI, and (T, A) denote the tissue and air ROIs respectively. In addition, the relative SNR and CNR values normalized by protocol-specific CTDI were calculated as nSNR and nCNR to represent image quality to imaging dose ratio. Fig. 8 shows the representative slices of the head phantom that were scanned in different head protocols. Compared with the reference image

e)&(h)&(k) and an obese rib mets patient in (c)&(f)&(i)&(l).
The scan protocols were standard Head, Thorax, and Pelvis Large respectively. The first two-row images of (a)-(f) were reconstructed by the vendor-provided program, and the middle two-row of (g)- of standard Head in Fig. 8(a), the images of High Quality in Fig. 8(b) and Spotlight in Fig. 8(c) exhibited reduced noise and about 20 % improvement in SNR and CNR properties. However, this 20 % improvement in SNR and CNR introduced 80 % increase in CTDI to High Quality, and about 300 % increase in CTDI to Spotlight. In addition, while the profile in Fig. 8(e) indicates that Image Gently presented improved bone to soft-tissue contrast, the subjective perception of Fig. 8 (d) as well as the SNR and CNR values in Fig. 8(f) were corrupted with too much noise.
It is interesting to note, while the SNR and CNR values of High Quality and Spotlight images were equivalent, the CTDI of High Quality was only about 45 % of Spotlight. Whereas Image Gently image presented the lowest SNR and CNR values, its nSNR and nCNR values were the largest, which indicates the highest radiation efficacy.

Case study B: Impact of projection reduction on image quality Motivation
As shown in the case study above, the self-designed High Quality protocol could deliver equivalent image quality with Spotlight by using less imaging dose. As a result, it was intuitive to raise the idea that maybe scan protocols could be optimized to achieve uncompromised image quality but with reduced dose. A straightward strategy was to reduce the total number of projections, which would proportionally decrease incident radiation dose.
In this case study, our aim was to clinically investigate the trend of image quality degradation over projection reduction, providing reference to re-tune relevant imaging parameters.

Data and method
The patient scan data reconstructed as in Fig. 7(a) was used in this study. The scan was acquired in the standard Head protocol, where the original projection number was 500. As a retrospective study, the fullstack projection series was downsampled to simulate virtual scans that were performed in increased angular intervals. The downsampling factors used herein were 2, 4 and 8, and correspondingly the imaging dose were reduced to 1/2, 1/4 and 1/8 of the original. These dose reduction protocols were denoted as Low Dose (LD), Ultra-Low Dose (ULD) and Under-Dosed (UD) hereafter. The projection numbers of the simulated LD, ULD and UD scans were 250, 125 and 63, which were equivalent to 0.8-degree, 1.6-degree and 3.2-degree in angular sampling intervals. The LD, ULD and UD scans as well as the original full-stack scan were all reconstructed by FDK and then denoised slice by slice with a 3x3 median filter.
For image quality evaluation, SNR and CNR in the dose-reduction images as well as the original full-stack image (benchmark) were calculated in a similar ROI-based approach. In addition, the peak-signalto-noise ratio (PSNR) and structure similarity (SSIM) values relative to the benchmark were also calculated as.
where I ref denotes the benchmark image reconstructed from the original full-stack scan, and MSE represents mean-squared error. (2μ x μ x + C 1 ) where (μ, μ y ), (σ x , σ y ), and σ xy are the local means, standard deviations, and cross-covariance for images ×,y; (C 1 , C 2 ) are constants, and C 1 = (0.01*L) 2 , C 2 = (0.03*L) 2 , L = 3000. Fig. 9 shows the representative slices of the same patient virtually scanned in different protocols. As the projections used for reconstruction became sparse, image noise and radial streaking artifacts significantly increased. Likewise in Table 1, while the CTDI in UD was only 1/8 of the benchmark, the SNR and CNR decreased by 72 % and 67 % respectively. Table 1 also indicates a general trend that when acquired projections were halved off, SNR and CNR values in reconstructed images declined by about 1/3, and SSIM values decreased about 10 %. The finding revealed in this patient-based study [24] is consistent with the phantom based study that CBCT dose can be significantly reduced by protocol optimization.

Case study C: IR based image quality enhancement of ULD-Image Gently for paediatric patients Motivation
Section 3.3 indicates that when views of data acquisition were reduced, tomographic image quality would dramatically deteriorate. However, it is important to note the prerequisite that all the tomographic images were reconstructed by the same algorithm (FDK in Section 3.3). When different algorithms were used involved, image quality could vary to some extent. Fortunately, the proposed toolkit provides seamless access to a large group of iterative reconstruction (IR) algorithms in TIGRE. The following study demonstrated the image quality enhancement of IR algorithms in a typical low-dose scenario.

Data and method
The scan data was from a 9-year-old non-Hodgkin lymphoma (NHL) paediatric patient treated at our institution. The CBCT scan was performed using the standard Image Gently protocol (kVp = 80, scan range = 200 degrees, N proj = 500). Herein, we used the same approach of 4x downsampling to generate ultra-low dose (ULD) scan data (125 views in total). As benchmark, the original scan data was reconstructed by FDK. As comparison, the ULD scan data was reconstructed by FDK as well as four typical IR algorithms. The IR algorithms were OS-SART, MLEM, CGLS and OS-ASD-POCS, which were accelerated by GPU (NVidia Quadro P5000). For image quality evaluation, we calculated SNR and CNR values in an ROI-based approach similar as above-one ROI in nasal cavity and the other within CTV. The PSNR values were also calculated between test images and the benchmark.

Results
Fig . 10 shows the representative slices of the paediatric patient virtually scanned and reconstructed in different approaches. Compared with the benchmark slice in Fig. 10(a) that was from the standard Image Gently scan, the slice from the ULD scan in Fig. 10(b) was rather noisy with severe steaking artifacts. This is consistent with the finding in Session 3.3.3 that when FDK was used for reconstruction, image quality would significantly decrease over projection reduction.
Compared with the slice reconstructed by FDK, the slices by IR algorithms exhibited substantial noise and artifacts reduction. Also, we can see that that the quality enhancement performance among IR algorithms were not the same. As shown in Table 2, the SNR and CNR in the MLEM slice were 140 % and 196 % higher than those in the ULD-FDK slice, and 28 % and 49 % higher than those in the benchmark slice. Likewise, the SNR and CNR values in the ASD-POCS slice were 419 % and 429 % higher than those in the ULD-FDK slice, and 172 % and 184 % higher than those in the benchmark slice. Among the low-dose algorithms, MLEM and ASD-POCS algorithms delivered the highest structure similarity and the SSIM values were 12.7 % and 15.6 % higher than the low-dose FDK, while the CGLS exhibited the lowest SSIM value with 0.710 (6.07 % lower than the low-dose FDK).
We compared the time efficiency of different algorithms in TIGRE on the same computational platform. Table 2 evidently shows that despite GPU acceleration the running time of IR algorithms were still much longer than FDK.

Discussion
The initiative of this project was to provide not only a tool to bridge the gap between clinical usage of CBCT scan data and research algorithms but also a platform that breaks down the imaging chain into individual processes so that research effort can be focused on a specific part. Even though some vendors provide similar software for research purpose such as iTools software [25,26] by Varian Medical System, to our knowledge these tools are strictly authorized to limited users. On the contrary, the proposed toolkit was open source under the BSD license, imposing minimal restrictions on its use and distribution. Also, the  toolkit was designed in a module-oriented paradigm with minimal data formatting. This essentially facilitates task-specific secondary development and helps in keeping the project convenient to maintain. From the beginning of this work, we focused on reproducing the vendor-designed imaging chain as much as possible. Therefore, a lot of efforts had been directed to reading technical whitepapers, decoding data files, parsing parameters, discussing with service engineers, et al. This brings another key feature that all the parameters in data acquisition and most hyperparameters in data conditioning and reconstruction are intrinsically initialized by pertinent values from scan data itself. This initialization strategy paves the way for new users to pick up, and in the same process, provides options for experienced users to fine-tine the imaging chain.
It is important to note two major limitations of this work. First, as the name indicates, this toolkit was for on-board CBCT on Varian linacs only. This is attributed to the limitation of our institution that we are a Varian reference site and we do not have any linac from other vendors. The project might be renamed in future, say IGRT instead, if community contributors added new features for on-board CBCT from other vendors. Second, despite the efforts in reproducing the vendor-designed imaging chain, the imaging chain we built was still not and would be impossibly identical to the vendor version, as indicated in the demo case study. This limitation comes from two aspects. On the one hand, while this opensource toolkit was developed by us, the official technical references were not open to us. Some published studies were lacking in detail, some manuals were out-of-date, some documents were inaccessible, some details were inconsistent across different references, and some parameters were hard-coded or encrypted somehow and somewhere. To make up missing pieces, we implemented relevant processes based either on our own research experience-oriented interpretation or pertinent code in other tools. For example, the beam-hardening correction function was mainly referred to the counterpart in MIRT [18]. On the other hand, the embedded program for on-board CBCT control is fully closed, making it impossible to access any intermediate outputs. This set a large barrier to perform cross-validation, which adversely hindered us fine-tuning hyperparameters. This became a non-negligible issue in some artifacts correction algorithms that required hyperparameters to optimize correction strength.
Also, there are two minor limitations of the current version TIGRE-VarianCBCT. First, as several technical references put, most artifacts correction algorithms designed by the vendor were plain and simple.  Since this alpha version was mainly to reproduce the vendor imaging chain, the artifacts corrections in this work are classic as well. New algorithms will be introduced in future versions, probably with community support. Second, the current toolkit is coded in Matlab and dependent on the TIGRE-Matlab package in image reconstruction. This is a limitation as Matlab itself is not free. Since the TIGRE-python package has been officially released, part of future efforts will be directed to developing a python version to further broaden its accessibility. As for its potential application, three typical areas can be easily listed: 1) development of new CBCT imaging techniques, 2) teaching demonstration of CT imaging as well as IGRT, 3) data manipulation for virtual clinical trials and data mining.
First, this toolkit was primarily designed to facilitate research activities in on-board CBCT imaging and would be helpful for users to explore novel solutions in many issues, including scan protocol optimization, artifacts suppression, iterative algorithms, et al. For instance, the study [27] clearly shows that accurate modelling of bowtie placement is crucial in avoiding crescent artifacts, and it is indicated [28] that a phantom-based bowtie normalization strategy prevails over the abovementioned airscan method in mitigating the effects of bowtie misplacement and beam-hardening. We believe that validation and comparison studies can be facilitated by the proposed toolkit, and maybe a hybrid method can be developed for improved performance.
Second, compared with other tools coded in C/C++, the toolkit was in Matlab, making it very user-friendly. Instructors in medical imaging/ physics could use it as a convenient tool to demonstrate an entire imaging chain in CT and image guidance in IGRT. For instance, the study [29] evaluated the performance of a commercial iterative algorithm in log file-based dose verification. We believe the proposed toolkit could provide more algorithms to perform a comparative study.
Third, patient scan data are growing rapidly from daily treatments, and these data are invaluable assets for clinical studies. The case studies had demonstrated an unparalleled advantage in simulating virtual patient scans by data re-sampling. Compared with either physical or computational phantom scans, these virtual scans are from real patients and acquired without any extra X-ray exposure. This could allow users to perform clinical trials in a retrospective approach, such as imaging parameter fine-tuning and personalized imaging. Likewise, large datasets can also be formed for machine learning.

Conclusions
We presented TIGRE-VarianCBCT: an open-source Matlab-GPU toolkit for Varian on-board cone-beam CT. The toolkit streamlined the entire process to generate volumetric CBCT images from raw scan data. The performance was validated on patient scan data with good agreement with respect to vendor-designed program. Case studies in scan protocol optimization, low dose imaging and iterative algorithm comparison demonstrated its substantial potential in performing scan data based clinical studies.

Ethical statement
The usage of all the patient scan data were approved by the IRB at Peking University Cancer Hospital.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.