Andreu Badal

Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit

PURPOSE It is a known fact that Monte Carlo simulations of radiation transport are computationally intensive and may require long computing times. The authors introduce a new paradigm for the acceleration of Monte Carlo simulations: The use of a graphics processing unit (GPU) as the main computing device instead of a central processing unit (CPU). METHODS A GPU-based Monte Carlo code that simulates photon transport in a voxelized geometry with the accurate physics models from PENELOPE has been developed using the CUDATM programming model (NVIDIA Corporation, Santa Clara, CA). RESULTS An outline of the new code and a sample x-ray imaging simulation with an anthropomorphic phantom are presented. A remarkable 27-fold speed up factor was obtained using a GPU compared to a single core CPU. CONCLUSIONS The reported results show that GPUs are currently a good alternative to CPUs for the simulation of radiation transport. Since the performance of GPUs is currently increasing at a faster pace than that of CPUs, the advantages of GPU-based software are likely to be more pronounced in the future.

A PENELOPE-based system for the automated Monte Carlo simulation of clinacs and voxelized geometries - application to far-from-axis fields

PURPOSE Two new codes, PENEASY and PENEASYLINAC, which automate the Monte Carlo simulation of Varian Clinacs of the 600, 1800, 2100, and 2300 series, together with their electron applicators and multileaf collimators, are introduced. The challenging case of a relatively small and far-from-axis field has been studied with these tools. METHODS PENEASY is a modular, general-purpose main program for the PENELOPE Monte Carlo system that includes various source models, tallies and variance-reduction techniques (VRT). The code includes a new geometry model that allows the superposition of voxels and objects limited by quadric surfaces. A variant of the VRT known as particle splitting, called fan splitting, is also introduced. PENEASYLINAC, in turn, automatically generates detailed geometry and configuration files to simulate linacs with PENEASY. These tools are applied to the generation of phase-space files, and of the corresponding absorbed dose distributions in water, for two 6 MV photon beams from a Varian Clinac 2100 C∕D: a 40 × 40 cm(2) centered field; and a 3 × 5 cm(2) field centered at (4.5, -11.5) cm from the beam central axis. This latter configuration implies the largest possible over-traveling values of two of the jaws. Simulation results for the depth dose and lateral profiles at various depths are compared, by using the gamma index, with experimental values obtained with a PTW 31002 ionization chamber. The contribution of several VRTs to the computing speed of the more demanding off-axis case is analyzed. RESULTS For the 40 × 40 cm(2) field, the percentages γ(1) and γ(1.2) of voxels with gamma indices (using 0.2 cm and 2% criteria) larger than unity and larger than 1.2 are 0.2% and 0%, respectively. For the 3 × 5 cm(2) field, γ(1) = 0%. These figures indicate an excellent agreement between simulation and experiment. The dose distribution for the off-axis case with voxels of 2.5 × 2.5 × 2.5 mm(3) and an average standard statistical uncertainty of 2% (1σ) is computed in 3.1 h on a single core of a 2.8 GHz Intel Core 2 Duo processor. This result is obtained with the optimal combination of the tested VRTs. In particular, fan splitting for the off-axis case accelerates execution by a factor of 240 with respect to standard particle splitting. CONCLUSIONS PENEASY and PENEASYLINAC can simulate the considered Varian Clinacs both in an accurate and efficient manner. Fan splitting is crucial to achieve simulation results for the off-axis field in an affordable amount of CPU time. Work to include Elekta linacs and to develop a graphical interface that will facilitate user input is underway.

Monte Carlo study of the effects of system geometry and antiscatter grids on cone‐beam CT scatter distributions

PURPOSE The proliferation of cone-beam CT (CBCT) has created interest in performance optimization, with x-ray scatter identified among the main limitations to image quality. CBCT often contends with elevated scatter, but the wide variety of imaging geometry in different CBCT configurations suggests that not all configurations are affected to the same extent. Graphics processing unit (GPU) accelerated Monte Carlo (MC) simulations are employed over a range of imaging geometries to elucidate the factors governing scatter characteristics, efficacy of antiscatter grids, guide system design, and augment development of scatter correction. METHODS A MC x-ray simulator implemented on GPU was accelerated by inclusion of variance reduction techniques (interaction splitting, forced scattering, and forced detection) and extended to include x-ray spectra and analytical models of antiscatter grids and flat-panel detectors. The simulator was applied to small animal (SA), musculoskeletal (MSK) extremity, otolaryngology (Head), breast, interventional C-arm, and on-board (kilovoltage) linear accelerator (Linac) imaging, with an axis-to-detector distance (ADD) of 5, 12, 22, 32, 60, and 50 cm, respectively. Each configuration was modeled with and without an antiscatter grid and with (i) an elliptical cylinder varying 70-280 mm in major axis; and (ii) digital murine and anthropomorphic models. The effects of scatter were evaluated in terms of the angular distribution of scatter incident upon the detector, scatter-to-primary ratio (SPR), artifact magnitude, contrast, contrast-to-noise ratio (CNR), and visual assessment. RESULTS Variance reduction yielded improvements in MC simulation efficiency ranging from ∼17-fold (for SA CBCT) to ∼35-fold (for Head and C-arm), with the most significant acceleration due to interaction splitting (∼6 to ∼10-fold increase in efficiency). The benefit of a more extended geometry was evident by virtue of a larger air gap-e.g., for a 16 cm diameter object, the SPR reduced from 1.5 for ADD = 12 cm (MSK geometry) to 1.1 for ADD = 22 cm (Head) and to 0.5 for ADD = 60 cm (C-arm). Grid efficiency was higher for configurations with shorter air gap due to a broader angular distribution of scattered photons-e.g., scatter rejection factor ∼0.8 for MSK geometry versus ∼0.65 for C-arm. Grids reduced cupping for all configurations but had limited improvement on scatter-induced streaks and resulted in a loss of CNR for the SA, Breast, and C-arm. Relative contribution of forward-directed scatter increased with a grid (e.g., Rayleigh scatter fraction increasing from ∼0.15 without a grid to ∼0.25 with a grid for the MSK configuration), resulting in scatter distributions with greater spatial variation (the form of which depended on grid orientation). CONCLUSIONS A fast MC simulator combining GPU acceleration with variance reduction provided a systematic examination of a range of CBCT configurations in relation to scatter, highlighting the magnitude and spatial uniformity of individual scatter components, illustrating tradeoffs in CNR and artifacts and identifying the system geometries for which grids are more beneficial (e.g., MSK) from those in which an extended geometry is the better defense (e.g., C-arm head imaging). Compact geometries with an antiscatter grid challenge assumptions of slowly varying scatter distributions due to increased contribution of Rayleigh scatter.

penMesh —Monte Carlo Radiation Transport Simulation in a Triangle Mesh Geometry

We have developed a general-purpose Monte Carlo simulation code, called penMesh, that combines the accuracy of the radiation transport physics subroutines from PENELOPE and the flexibility of a geometry based on triangle meshes. While the geometric models implemented in most general-purpose codes-such as PENELOPEs quadric geometry-impose some limitations in the shape of the objects that can be simulated, triangle meshes can be used to describe any free-form (arbitrary) object. Triangle meshes are extensively used in computer-aided design and computer graphics. We took advantage of the sophisticated tools already developed in these fields, such as an octree structure and an efficient ray-triangle intersection algorithm, to significantly accelerate the triangle mesh ray-tracing. A detailed description of the new simulation code and its ray-tracing algorithm is provided in this paper. Furthermore, we show how it can be readily used in medical imaging applications thanks to the detailed anatomical phantoms already available. In particular, we present a whole body radiography simulation using a triangulated version of the anthropomorphic NCAT phantom. An example simulation of scatter fraction measurements using a standardized abdomen and lumbar spine phantom, and a benchmark of the triangle mesh and quadric geometries in the ray-tracing of a mathematical breast model, are also presented to show some of the capabilities of penMesh.

Monte Carlo reference data sets for imaging research: Executive summary of the report of AAPM Research Committee Task Group 195

The use of Monte Carlo simulations in diagnostic medical imaging research is widespread due to its flexibility and ability to estimate quantities that are challenging to measure empirically. However, any new Monte Carlo simulation code needs to be validated before it can be used reliably. The type and degree of validation required depends on the goals of the research project, but, typically, such validation involves either comparison of simulation results to physical measurements or to previously published results obtained with established Monte Carlo codes. The former is complicated due to nuances of experimental conditions and uncertainty, while the latter is challenging due to typical graphical presentation and lack of simulation details in previous publications. In addition, entering the field of Monte Carlo simulations in general involves a steep learning curve. It is not a simple task to learn how to program and interpret a Monte Carlo simulation, even when using one of the publicly available code packages. This Task Group report provides a common reference for benchmarking Monte Carlo simulations across a range of Monte Carlo codes and simulation scenarios. In the report, all simulation conditions are provided for six different Monte Carlo simulation cases that involve common x-ray based imaging research areas. The results obtained for the six cases using four publicly available Monte Carlo software packages are included in tabular form. In addition to a full description of all simulation conditions and results, a discussion and comparison of results among the Monte Carlo packages and the lessons learned during the compilation of these results are included. This abridged version of the report includes only an introductory description of the six cases and a brief example of the results of one of the cases. This work provides an investigator the necessary information to benchmark his/her Monte Carlo simulation software against the reference cases included here before performing his/her own novel research. In addition, an investigator entering the field of Monte Carlo simulations can use these descriptions and results as a self-teaching tool to ensure that he/she is able to perform a specific simulation correctly. Finally, educators can assign these cases as learning projects as part of course objectives or training programs.

Spatiotemporal Monte Carlo transport methods in x-ray semiconductor detectors: Application to pulse-height spectroscopy in a-Se

PURPOSE The authors describe a detailed Monte Carlo (MC) method for the coupled transport of ionizing particles and charge carriers in amorphous selenium (a-Se) semiconductor x-ray detectors, and model the effect of statistical variations on the detected signal. METHODS A detailed transport code was developed for modeling the signal formation process in semiconductor x-ray detectors. The charge transport routines include three-dimensional spatial and temporal models of electron-hole pair transport taking into account recombination and trapping. Many electron-hole pairs are created simultaneously in bursts from energy deposition events. Carrier transport processes include drift due to external field and Coulombic interactions, and diffusion due to Brownian motion. RESULTS Pulse-height spectra (PHS) have been simulated with different transport conditions for a range of monoenergetic incident x-ray energies and mammography radiation beam qualities. Two methods for calculating Swank factors from simulated PHS are shown, one using the entire PHS distribution, and the other using the photopeak. The latter ignores contributions from Compton scattering and K-fluorescence. Comparisons differ by approximately 2% between experimental measurements and simulations. CONCLUSIONS The a-Se x-ray detector PHS responses simulated in this work include three-dimensional spatial and temporal transport of electron-hole pairs. These PHS were used to calculate the Swank factor and compare it with experimental measurements. The Swank factor was shown to be a function of x-ray energy and applied electric field. Trapping and recombination models are all shown to affect the Swank factor.

hybrid

The computational modeling of medical imaging systems often requires obtaining a large number of simulated images with low statistical uncertainty which translates into prohibitive computing times. We describe a novel hybrid approach for Monte Carlo simulations that maximizes utilization of CPUs and GPUs in modern workstations. We apply the method to the modeling of indirect x-ray detectors using a new and improved version of the code MANTIS, an open source software tool used for the Monte Carlo simulations of indirect x-ray imagers. We first describe a GPU implementation of the physics and geometry models in fastDETECT2 (the optical transport model) and a serial CPU version of the same code. We discuss its new features like on-the-fly column geometry and columnar crosstalk in relation to the MANTIS code, and point out areas where our model provides more flexibility for the modeling of realistic columnar structures in large area detectors. Second, we modify PENELOPE (the open source software package that handles the x-ray and electron transport in MANTIS) to allow direct output of location and energy deposited during x-ray and electron interactions occurring within the scintillator. This information is then handled by optical transport routines in fastDETECT2. A load balancer dynamically allocates optical transport showers to the GPU and CPU computing cores. Our hybridMANTIS approach achieves a significant speed-up factor of 627 when compared to MANTIS and of 35 when compared to the same code running only in a CPU instead of a GPU. Using hybridMANTIS, we successfully hide hours of optical transport time by running it in parallel with the x-ray and electron transport, thus shifting the computational bottleneck from optical tox-ray transport. The new code requires much less memory than MANTIS and, asa result, allows us to efficiently simulate large area detectors.

\scriptsize{\mathrm{MANTIS}}

A code for Monte Carlo simulations of radiation transport using a Graphics Processing Unit (GPU) is introduced. The code has been developed using the CUDATM programming model, an extension to the C language that allows the execution of general purpose computations on the new generation of GPUs from NVIDIA. The accurate Compton and Rayleigh interaction models and interaction mean free paths from the PENELOPE package, and a generic voxelized geometry model, have been implemented in the new code. The secondary particles generated by Compton, photoelectric and pair-production events are not transported. An ideal x-ray detector and a cone beam source can be defined to reproduce an imaging system and facilitate the simulations of medical imaging applications. A 24-fold speed up factor with the GPU compared to the CPU is reported for a radiographic projection of a detailed anthropomorphic female phantom. A description of the simulation algorithm and the technical implementation in the GPU are provided. This work shows that GPUs are already a good alternative to CPUs for Monte Carlo simulation of x-ray transport.

: a CPU?GPU Monte Carlo method for modeling indirect x-ray detectors with columnar scintillators

Purpose: To present a Monte Carlo(MC) x‐ray transport simulation code that uses a massively parallel graphics processing units (GPU) to generate medical images of voxelized patient models and estimate organ doses. Methods: We updated MC‐GPU, a CUDA‐based open‐source code for the simulation of radiographic projections and CT, to include the capability to tally 3D dose distributions. Since electrons are not transported, deposited doses are assumed to be equivalent to photon KERMA. The validity of this assumption was investigated by comparing it with the general‐purpose MC code PENELOPE in the simulation of a posterior‐anterior chest examination with 1010, 60‐keV x rays and a realistic female phantom with 1‐mm voxels. The stability of the simulation in a 480‐core GPU was also benchmarked using a sequential execution in a single CPU core. A utility to post‐process the voxel doses and report average and peak organ doses and effective dose for the procedure (and their uncertainty) has been developed. Results: The accuracy of the dose values generated by MC‐GPU was successfully validated. Average relative differences in the doses from PENELOPE and MC‐GPU executed in CPU and GPU were below 2%, in the expected range for statistical fluctuations. The simulation speeds were 46554 x rays/s with PENELOPE (electron transport disabled), 158216 x rays/s with MC‐GPU in a CPU, and 5124867 x rays/s in a GPU (simulation time in the GPU: 32.5 minutes, 32‐fold speedup compared to the CPU). The average and peak doses in breast and skin were estimated as 1.1 and 3.4 eV/g, and 1.7 and 15.6 eV/g per history respectively. Conclusion: The presented code can estimate organ doses 110 times faster than a general‐purpose MC code with comparable accuracy. A multi‐GPU execution of MC‐GPU has the potential to provide accurate organ doses in a CT scan in near real‐time.

Monte Carlo simulation of X-ray imaging using a graphics processing unit

We present a Monte Carlo (MC) simulation method for studying the signal formation process in amorphous Selenium (a-Se) imaging detectors for design validation and optimization of direct imaging systems. The assumptions and limitations of the proposed and previous models are examined. The PENELOPE subroutines for MC simulation of radiation transport are used to model incident x-ray photon and secondary electron interactions in the photoconductor. Our simulation model takes into account applied electric field, atomic properties of the photoconductor material, carrier trapping by impurities, and bimolecular recombination between drifting carriers. The particle interaction cross-sections for photons and electrons are generated for Se over the energy range of medical imaging applications. Since inelastic collisions of secondary electrons lead to the creation of electron-hole pairs in the photoconductor, the electron inelastic collision stopping power is compared for PENELOPEs Generalized Oscillator Strength model with the established EEDL and NIST ESTAR databases. Sample simulated particle tracks for photons and electrons in Se are presented, along with the energy deposition map. The PENEASY general-purpose main program is extended with custom transport subroutines to take into account generation and transport of electron-hole pairs in an electromagnetic field. The charge transport routines consider trapping and recombination, and the energy required to create a detectable electron-hole pair can be estimated from simulations. This modular simulation model is designed to model complete image formation.