VOXSI: A voxelized single- and dual-energy CT scenario generator for quantitative imaging.

BACKGROUND AND
PURPOSE: Dedicated CT simulation models have the potential to investigate several acquisition, reconstruction, or post-processing parameters without giving any radiation dose to patients. A software program was developed for the simulation and the analysis of single-energy and dual-energy CT images. Simulation and analysis functionalities of the software are described.
MATERIALS AND METHODS: In the software, named VOXSI (VOXelized CT SImulator), the X-ray source, user specified simulation geometry, CT setup and the detector energy response can be varied. CT image reconstructions can be performed with an implementation of the ASTRA toolbox. In the DECT post processing toolkit, GUI tools are provided to calculate effective atomic number, relative electron density, pseudo-monoenergetic images, and material map images. Quantitative CT number validation, based on a RMI 467 tissue characterization phantom model, was performed between experimental and simulated CT scans at three different X-ray tube potentials (80, 120, and 140 kVp) with a third generation CT scanner.
RESULTS: Overall, a good agreement was found for the mean CT numbers of the RMI 467 inserts. For all energies, the maximum difference in CT numbers between experimental and simulated data was below 17 HU for the soft tissues and below 48 HU for the osseous tissues.
CONCLUSION: The software's simulation algorithm showed a good agreement between the CT measurements and CT simulations of the RMI 467 phantom at different energies. The capabilities of the software are demonstrated by an elaborated dual-energy CT research example.

Introduction

X-ray imaging is widely used to acquire images in many fields, such as oncology and radiotherapy, where imaging plays a crucial role in diagnosing and treating a wide range of cancerous diseases. An image modality such as single-energy (SE) computed tomography (CT) is frequently used for radiology and radiotherapy treatment planning [1]. Besides conventional SECT imaging, interest in dual-energy CT (DECT) for radiology and radiotherapy is increasing [2]. A growing body of literature suggests that DECT is beneficial for subjective and objective image quality, soft tissue characterization, projection or image based metal artifact reduction, brachytherapy dose calculations, particle therapy dose calculations, and small animal radiotherapy dose calculations [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17].

Multiple aspects of CT imaging should be investigated to establish an optimal scanning protocol for each specific treatment site and patient size. Unfortunately, changing hardware components, or varying imaging parameters is not always possible in a clinical workflow. A dedicated CT simulation model has the potential to investigate scanning protocols, reconstruction algorithms, and beam hardening effects without giving radiation dose to patients or requiring the development of prototypes and acquisition of phantom scans.

In the literature, simulation tools for different imaging modalities were published. For example, tools were presented to simulate planar X-ray images [18], [19], [20]. Ay and Zaidi developed a fan-beam CT and cone-beam CT X-ray CT simulator of mathematical volumes based on Monte Carlo (MC) simulations [21]. Other researchers applied a ray-tracing approach algorithm to simulate fan-beam CT images [22], [23]. Segars et al. developed a parallel-beam, fan-beam and cone-beam CT projector to simulate realistic humanoid anatomies from mesh volumes, without an analysis toolkit that provides tools to process the simulated data [24], [25]. Another simulation tool, ImaSim simulates planar kilo- and megavoltage, fan-beam CT, and cone-beam CT (without scatter) images of simulation geometries formed by combining only a few mathematical volumes [26].

Current CT simulation models use either MC simulations or a ray-tracing approach. MC simulates individual particle interactions and requires much more processing power compared to the ray-tracing alternative which ignores individual photon interactions to speed up the calculation time.

In this paper, a new fan-beam CT scenario generator accessible via a graphical user interface (GUI) is described. The simulation software has the flexibility to use a variety of voxelized geometries, and to modify X-ray spectra, geometrical setup, detector responses, and CT image reconstruction parameters. The analysis software toolkit includes DECT analysis capabilities to examine relevant parameters such as the relative electron density (RED), effective atomic number (Zeff), pseudo-monoenergetic energy images (PMIs), and material maps. In material maps, unique material flags are assigned to the voxels to obtain a spatial material distribution. Compared with other available CT simulation packages, our software allows the user to perform a fast multi-threaded CT simulation, from voxelized user defined geometries, while many simulation and reconstruction parameters can be varied. No other simulation packages are currently available that allow the user to perform different DECT analyses. For validation of our simulation model, the model was compared with measurements at three different X-ray tube potentials.

Materials and methods

A schematic overview of the different tasks and toolkits embedded in the software are depicted in Fig. 1. The simulation software, named VOXSI (VOXelized CT SImulator), was developed in MATLAB® (The MathWorks Inc, Natick, MA) and runs on 64-bit Windows and Linux computers.

Fig. 1

The flowchart of VOXSI. The X-ray source, simulation geometry, CT setup, and the detector energy response (top row) must be defined before a simulation can be started. After the simulation, the projections are exportable, or can be reconstructed using the implemented ASTRA toolbox and MATLAB’s FBP algorithm. As an alternative, the open-source image reconstruction toolkit (RTK) can be used externally (not implemented in the GUI) to reconstruct the projections (dotted line). The SECT and DECT image analysis toolkit provides tools to analyze the reconstructed CT images. Snapshots of VOXSI’s interface are added in Supplementary material 1.

X-ray source

As X-ray spectrum generator, SpekCalc was implemented in VOXSI’s GUI [27]. SpekCalc is a validated simulation software package which calculates filtered polychromatic X-ray spectra (40–300 kVp) for X-ray tubes with a tungsten anode via a deterministic/MC model. The model calculates the bremsstrahlung spectrum in combination with a MC pre-calculated energy and depth distribution of electrons in a tungsten target [28], [29]. Additionally, measured spectra, spectra simulated in an external package, monoenergetic photons, and characteristic gamma rays emitted by radionuclides can be used as input of VOXSI’s CT simulation algorithm.

Voxelized geometry

VOXSI was designed such that only voxelized geometries are used as input for the CT simulation algorithm, e.g. DICOM CT scans or user specified voxelized geometries. Delineation tools, provided in the GUI, such as a brush (draw/erase) and intensity based region growing help the user to assign materials in the voxels of a DICOM CT scan. The user specified geometries should have unique integer values assigned to every voxel. In VOXSI, these unique integer values can be linked to a material in the material database.

The 3D humanoid XCAT phantom, which also has 4D capabilities, serves as a realistic and flexible voxelized phantom [25]. The activity mode of the XCAT phantom produces a humanoid phantom with unique organ labels, defined in the XCAT parameter file, which can be loaded in VOXSI’s geometry editor. This methodology is also applicable on the 4D mouse and rat whole body phantoms: MOBY and ROBY [30].

CT setup and detector response

VOXSI simulates a non-helical fan-beam CT setup, although clinical CT scanners often scan in helical mode. A non-helical approach is assumed because a voxelized geometry consists of discrete slices. Different fan-beam setup parameters can be varied in the software; source/detector to origin distance, fan-beam angle, curved or planar detector, number of detectors, and number of projections.

An essential input of VOXSI is the detector energy response, defined as the ratio of the absorbed and incident photon energy. The detector response curve can often be obtained from the manufacturers by request. As an alternative, the detector response curve can be obtained by performing MC simulations in external software using the detector specifications [30], [31]. Whichever way it was obtained, the requested, pre-calculated, or even measured detector energy response curve can be loaded in VOXSI. As standard, an ideal detector response with a 100% absorption for all incident photon energies, is provided in the software.

Simulation and CT image reconstruction

A multithreaded ray-tracing approach was implemented in VOXSI to reduce calculation time. E.g. 45 s of computing time was needed to simulate a 512 × 512 × 1 geometry, a 120 kVp X-ray spectrum, 1024 detectors, and 780 projections, on two Intel® Xeon® 2.67 GHz X550 processors.

VOXSI utilizes the total mass attenuation coefficients, grouping all main three photon interaction mechanisms, photoelectric absorption, Rayleigh scatter, and Compton scatter, of the attenuating material that was assigned to the intersecting voxel. The polychromatic X-ray attenuation in each voxel is calculated by using the total mass attenuation coefficient as a function of the photon energy and the radiological path length, according to Lambert-Beer’s attenuation law [31], [32]. The radiological path length through the voxel grid is calculated according to the algorithm described in [32]. The mass attenuation coefficients of a specific material can be calculated from the photon cross section database (XCOM version 3.1) [33]. Bowtie filters of different materials, described by their focal spot distance, thickness, and curvature, can be added to the CT simulation model.

To correct for cupping artifacts in polychromatic CT simulations, a basic beam hardening correction algorithm can be applied on the raw projection data. First, the algorithm assumes that the original simulation geometry consists of liquid water. Then, ideal projections (ignoring beam hardening) and spectrum hardened projections are calculated using the polychromatic X-ray spectrum, the detector energy response, and the radiological path length through the simulated water geometry. Finally, the ratio between the hardened and the ideal projection data was used to correct the projection data of the original simulation geometry for beam hardening artifacts. Because the raw projection data can be extracted from VOXSI, the user has the ability to apply their own beam hardening correction algorithm in external software.

An open source platform for image reconstruction, the ASTRA toolbox was implemented in VOXSI to reconstruct the simulated raw projection data [34]. The fan-beam projection data can be reconstructed with filtered back projection (FBP), a simultaneous iterative reconstruction technique (SIRT), a simultaneous algebraic reconstruction technique (SART), and a conjugate gradient least squares (CGLS) algorithm. These reconstruction techniques are based on a planar projector, while the projection data is simulated with either a curved or a planar projector. Therefore, the fan-beam projection data is converted to parallel-beam projection data by applying a rebinning procedure. Poisson noise can be added in the projection domain to simulate CT scans acquired at a lower imaging dose.

Single-energy CT and dual-energy CT image analysis toolkit

VOXSI has an implemented SECT and DECT analysis toolkit. The SECT analysis toolkit allows the user to calculate various image metrics (e.g. minimum, maximum, and mean Hounsfield Units ±1 standard deviation) in regions of interest defined by contours delineated by the implemented delineation tools. The metrics are exportable in Excel format (Windows operating systems only). The more extended DECT analysis toolkit (GUI) allows the user to calculate RED images according to Saito’s method, Zeff images according to both Landry’s method (2013) and Saito and Sagara’s method, PMIs following the principles in Yu et al., and material maps according to method described in Landry et al. (2011) [36], [37], [38], [39], [40]. Simulated low energy and high energy image data of the Gammex RMI 467 tissue characterization phantom (Gammex, Middleton, WI) inserts could be used to obtain the fit parameters.

Saito’s method calculates RED images based on a linear fit procedure, using energy-subtracted CT numbers , calculated by a weighted subtraction of and [35]. In Eqs. (1), (2), a, b, and are fit parameters.

Landry’s method calculates Zeff⁎ images based on a fit procedure [36]. Eq. (3) shows the fitting model between the ratio of reduced CT numbers , for or , which involves several fit parameters; A, B, C, D, E, F, m, and n.

Quantitative CT validation

To validate VOXSI’s simulation algorithm quantitatively, the RMI 467 phantom with its certified tissue mimicking inserts was scanned at three X-ray tube potentials (80, 120, and 140 kVp) with a SOMATOM Confidence® RT Pro (Siemens Healthineers, Forchheim, Germany) CT scanner. These SECT scans were all acquired with a constant CT dose index (32 cm) of 40 mGy and were reconstructed with the FBP algorithm (B30 kernel, Siemens terminology) in a 512 × 512 image matrix with a voxel dimension of 0.7 × 0.7 × 1.0 mm3. The specifications of the X-ray tube, the spectral filtration, and the detector energy response curve provided by the manufacturer were used to build the CT simulation model in VOXSI.

The acquired CT image data of the RMI 467 phantom was loaded in VOXSI to build a model of the phantom’s central slice using the region growing delineation tool†. This model is independent from the CT scanner. First, the phantom bulk and its inserts were delineated, and based on the delineation, VOXSI automatically calculates a unique material map phantom. Subsequently, the model of the RMI 467 phantom was simulated for the same three X-ray tube potentials and reconstructed with the FBP algorithm in the same matrix size and voxel dimensions as the acquired CT scans. The noise levels measured in the RMI 467 phantom bulk of the real CT scans were adopted in the VOXSI CT simulations of the RMI 467 phantom by adding realistic Poisson noise in the projection domain by giving an mAs value (Supplementary material 1). A visual comparison between a simulated and a measured RMI 467 phantom at 120 kVp is added to Supplementary material 2.

Circular regions of interest within the tissue mimicking inserts were defined on the acquired and the simulated CT scans to calculate the mean Hounsfield Unit (HU) values and their standard deviations for every tissue mimicking insert and every X-ray tube potential.

XCAT phantom dual-energy CT simulations

VOXSI’s RMI 467 phantom model was used to perform 80/140 kVp DECT simulations. Because the X-ray tube current (mAs) is non-equally divided over the dose of the two energy scans in a clinical DECT scanner, and must be known by the simulation algorithm to determine the noise level, an extra dual-spiral DECT scan (80/140 kVp) with a total CT dose index (32 cm) of 40 mGy was acquired. The RMI 467 phantom DECT simulations were then used to acquire the RED and Zeff calibration fits (Fig. 3, left panels).

Fig. 3

The relative electron density (RED) and effective atomic number (Zeff) calibration curves fitted on the simulated image data of the RMI 467 phantom (left panels) (±1 SD). In the Zeff calibration curve, only 8 of the 11 data points are clearly visible because 3 points closely overlap with another data point. The RED and Zeff images of the simulated XCAT phantom are depicted in the middle panels. The relative differences between the calculated images and the ground truth are shown in the right panels.

An XCAT phantom of the abdominal region was generated with a voxel dimension of 0.7 × 0.7 × 1.0 mm3. Eight different ICRU tissues, from VOXSI’s customizable material database, were assigned in XCAT’s voxel space [40]. 80/140 kVp CT imaging, with FBP reconstructions, of the XCAT phantom was then simulated with the known scanner characteristics.

Results

Quantitative CT validation: RMI 467 phantom

The overall difference between the simulated and the measured mean HUs were below 17 HU for all soft tissue mimicking inserts and below 48 HU for all osseous tissue mimicking inserts Figure 2. The maximum HU differences between the measured and the simulated CT numbers were obtained at the lowest X-ray tube potential of 80 kVp and in the densest osseous tissue mimicking inserts such as ‘CaCO3–50%’ (1.53 g/cm3) and ‘Cortical Bone’ (1.82 g/cm3). Increasing the X-ray tube potential to 120 kVp and 140 kVp reduced this maximum error in osseous tissue mimicking inserts to 10 HU and 6 HU respectively.

Fig. 2

Comparison of the measured and simulated CT numbers averaged in a circular region of interest (±1 standard deviation (SD)), in Hounsfield Units, for polychromatic 80, 120, and 140 kVp CT X-ray spectra of the RMI 467 soft tissue mimicking inserts (left panel) and the osseous tissue mimicking inserts (middle panel). The right panel plots the absolute differences between the measured and simulated CT numbers for all RMI 467 phantom inserts.

Dual-energy CT simulations: XCAT phantom

Large RED errors occurred around the body and near tissue transitions with a large density change. In the soft tissues, the RED error was smaller (±4%) but increases in regions subject to CT artifacts, i.e. beam hardening. The Zeff errors have a large variation in the soft tissue range (±12%), but a smaller variation in the osseous tissues (±6%). A consistent RED and Zeff overestimation of the lung tissue is visible in Fig. 3 (right panels).

Discussion

An SECT and DECT scenario simulator and image analysis toolkit, was presented for the simulation of voxelized volumes of interest in radiology and radiotherapy. The software, which will be released for general research use in the future, provides multiple degrees of freedom in the X-ray spectra, geometry, geometrical setup, detector energy response, reconstruction algorithms, and the image analysis. The simulation algorithm accuracy was quantitatively evaluated against a RMI 467 phantom with multiple tissue mimicking inserts that was scanned with a third generation CT scanner at 80, 120, and 140 kVp. The versatility of the software was demonstrated by an elaborated research example investigating the accuracy of the RED and Zeff determination algorithms in a 3D humanoid thoracic XCAT phantom with known tissue characteristics.

Based on the results presented in quantitative CT evaluation (Fig. 2), we expect that the relatively small differences (±15 HU) in CT numbers were not only originating from the accuracy of the CT simulation algorithm and the reconstruction algorithm, but also from the accuracy of the reported atomic compositions of the tissue mimicking inserts. The uncertainties of the atomic compositions of the RMI 467 phantom are not reported in their certificate. Additionally, the small differences in CT numbers could also be explained by the fact that the beam hardening correction algorithm implemented in the software is different from the undisclosed beam hardening correction algorithm in the commercial CT reconstruction software. Eventually, we expect the two different beam hardening correction implementations to be the main cause of discrepancies between the measured and the simulated CT numbers, mainly for the water insert where the Hounsfield Units should be around 0 HU, according to Hounsfield’s definition. However, for radiotherapy treatment planning with MV photon beams, 4–10% uncertainty in RED would result in no more than a 2% error in dose [41], [42], [43], [44].

In Fig. 3, the maximum RED errors occur around the body contour and near the regions with large density changes (±12%): e.g. trachea, lung, and bones, because of the boundary smoothing caused by the FBP image reconstruction algorithm. The Zeff calibration curve shows that a small change in the . will result in a large shift in Zeff value. This gives rise to noisy Zeff images, and therefore larger errors between the calculated Zeff image and the ground truth. The largest relative difference in Zeff is found in lung tissue because the Zeff of ICRU inflated lung tissue has a very low RED but a normal Zeff in the soft tissue range.

The CT simulation software can provide researchers with simulated raw projection data, which is more difficult to obtain from clinical CT scanners, to investigate e.g. material decomposition algorithm physical quantity extraction, SECT RED calculation, or projection based (metal) artifact correction algorithms, to name just a few applications [44], [45], [46], [47]. Currently, the ASTRA toolbox and the MATLAB implementation of FBP were implemented in the software. As an alternative, the open-source reconstruction toolkit (RTK) can be used to reconstruct the simulated projections [48]. RTK is mostly used to reconstruct cone-beam CT projections, but can also be used to reconstruct plane detector fan-beam CT projections, where no rebinning of fan projections to parallel projections is required. The ASTRA reconstructions in the section: “Quantitative CT validation” were repeated with RTK, no noteworthy differences in CT numbers were found.

The image analysis toolkit in the software provides tools to calculate image metrics for the SECT scans, and tools to perform a DECT analysis. A potential use of the CT scenario generator is to find optimal X-ray spectra, filtration materials or thicknesses for a DECT-based tissue assignment in patient-like geometries, given a user specified geometry (e.g. XCAT phantom) and detector energy response, which is useful for MC dose calculations in megavoltage photon, proton, and brachytherapy treatments. Different tissue assignment methods could be compared directly because the software still has the originally assigned materials in voxel space. In a real patient the absolute truth is missing and therefore only relative differences can be shown between various SECT or DECT based tissue assignment methods, e.g. comparisons of radiotherapy dose distributions [14], [49].

The presented software can also assess different RED/Zeff quantification algorithms and tissue assignment methodologies. Currently, the latter are based on individual pixel values: low and high energy CT HUs, RED, and Zeff, which have three major disadvantages. First, surrounding pixel values or body site information are not taken into account. Second, a material with comparable tissue characteristic could be wrongly assigned due to an increased noise level and third, no knowledge of the anatomy is considered. Using the flexibility of the XCAT phantom and other voxelized geometries, alternative methods making use of artificial intelligence e.g. deep learning algorithms, could be developed to improve tissue segmentation considering 3D voxel space and CT artifacts. The user is not limited to the analysis algorithms provided in the GUI. The simulation data can be exported as MATLAB or DICOM files which can be further analyzed in external software packages.

Table 1

The fit coefficients obtained for the RED and Zeff fit with the corresponding coefficients of determination. The average residuals for the RED fit and the Zeff fit were 0.01 and = 5.38E-4, respectively.

Fit	A	B	C	D	E	F	m	n	α	R2
Zeff	0.19	0.01	2.91E-5	0.19	0.01	2.91E-7	4.55	2.87	–	0.99
RED	1.04	0.98	–	–	–	–	–	–	0.90	0.99