Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units

Abstract

We report a parallel Monte Carlo algorithm accelerated by graphics processing units (GPU) for modeling time-resolved photon migration in arbitrary 3D turbid media. By taking advantage of the massively parallel threads and low-memory latency, this algorithm allows many photons to be simulated simultaneously in a GPU. To further improve the computational efficiency, we explored two parallel random number generators (RNG), including a floating-point-only RNG based on a chaotic lattice. An efficient scheme for boundary reflection was implemented, along with the functions for time-resolved imaging. For a homogeneous semi-infinite medium, good agreement was observed between the simulation output and the analytical solution from the diffusion theory. The code was implemented with CUDA programming language, and benchmarked under various parameters, such as thread number, selection of RNG and memory access pattern. With a low-cost graphics card, this algorithm has demonstrated an acceleration ratio above 300 when using 1792 parallel threads over conventional CPU computation. The acceleration ratio drops to 75 when using atomic operations. These results render the GPU-based Monte Carlo simulation a practical solution for data analysis in a wide range of diffuse optical imaging applications, such as human brain or small-animal imaging.

Fig. 1

Block diagram of the parallel Monte Carlo simulation for photon migration. The curved dashed lines indicate read-only global memory access and curved solid line for read/write access..

Fig. 2

Determination of the reflection interface between the medium (shaded) and air (clear) voxels: (a) case with 2 intersections and (b) case with 3 intersections.

Fig. 3

Serial correlation of the logistic-lattice (N=5) based random number generator.

Fig. 4

Ratios between the missing and total accumulation events for regions >3 voxels away from the source at various threads and scattering coefficients.

Fig. 5

The comparisons between parallel Monte Carlo algorithm (MCX), tMCimg and the diffusion model for a semi-infinite medium: (a) the time courses at voxel (30,14,9), (b) the contour plots for t=0.1 to 2.1 ns with 0.5 ns step along plane y=30, (c) radial distribution of continuous-wave (CW) solution on the interface, (d) CW fluence contour plot (10 dB spacing) along plane y=30, (e) comparisons between atomic and non-atomic solutions near the source, and (f) domain diagram showing where the results were extracted. In (a), (c) and (e), medium refraction index is 1.37 for simulations with boundary reflections.

Fig. 6

The (a) continuous-wave and (b) time-resolved solutions (Media 1) of photon migration in an MRI head anatomy (transparent overlay). The color map depicts the logarithmic fluence values.

Fig. 7

Simulation speed (photon/ms) with various thread numbers and simulation parameters for (a) semi-infinite medium and (b) brain MRI atlas. “MT” - Mersenne-Twister RNG; “LL5” - Logistic-lattice of size 5; “fast” - linked with CUDA's fast math library; “atomic” – with atomic operations (otherwise, without). The acceleration ratio compared with CPU implementation is marked along the z-axis.