Original Research ARTICLE
Parallelized Monte Carlo Photon Transport Simulations for Arbitrary Multi-Angle Wide-Field Illumination in Optoacoustic Imaging
- 1School of Precision Instruments and Optoelectronics Engineering, Tianjin University, Tianjin, China
- 2Tianjin Key Laboratory of Biomedical Detecting Techniques and Instruments, Tianjin University, Tianjin, China
Optoacoustic imaging (OAI) or photoacoustic imaging can resolve the distribution of tissue chromophores and optical contrast agents deep inside tissue from the optoacoustic detection with multi-spectral illumination. The development of a fast and accurate modeling method for the photon transport in OAI is necessary to quantitatively evaluate the tissue optical parameters. This paper presents a parallelized Monte Carlo modeling method especially for OAI (MCOAI) to simulate photon transport in bio-tissues with arbitrary multi-angle wide-field illumination. The performance of the MCOAI method is verified by comparison with the graphics processing unit (GPU)-accelerated MCX method in the typical cases with the pencil beam and ring source illumination. The simulation results demonstrate the GPU-based MCOAI method has equivalent accuracy and significantly improved computation efficiency, compared to the MCX method. The simulations with cylindrical and hemispherical source illumination further illustrate that the MCOAI method can effectively implement three-dimensional photon transport simulation for typical illumination geometries of OAI systems. A cross-section of Digimouse is selected as a realistic heterogeneous phantom illuminated with six different light sources usually employed in OAI, in order to prove the necessity of establishing photon transport modeling in OAI for the quantitative visualization in deep tissues. The MCOAI method can provide a powerful tool to efficiently establish photon transport modeling with arbitrary illumination modes for OAI applications.
Optoacoustic imaging (OAI) or photoacoustic imaging is a hybrid imaging modality that can provide multi-scale optical visualization with high depth-to-resolution ratio in biological tissues [1–6]. The distribution of the initial acoustic pressure or absorbed light energy density P0, is generally reconstructed from optoacoustic signals detected by ultrasonic transducers in traditional OAI methods. The initial acoustic pressure P0 is the product of optical absorption coefficient and light fluence in the irradiated region [7, 8]. Since light fluence in the tissue is attenuated along the penetration depth, the result of traditional OAI, i.e., the distribution of P0, cannot directly reflect the optical properties of deep tissues, such as optical absorption coefficient [9, 10]. To accurately reveal the optical functional information of deep tissues, a precise photon transport modeling approach is necessary to be developed for offering assistance to quantitative OAI methods.
As the gold standard for predicting distribution of light fluence in the biological tissue, the Monte Carlo (MC) simulation with strong applicability is preferred for the precise modeling with multi-angle wide-field illuminations usually employed in OAI systems. MC simulation methods have been widely developed for optical imaging modalities [11–16]. The MCML and MCXYZ methods by Wang et al. have been used to study light propagation in multi-layered tissues [11, 12]. Several light source modes have been realized including the pencil beam and the finite beam with a fixed incident direction, as well as the isotropic beam. An accelerated MCML has been developed to reduce the computation time based on graphics processing unit (GPU) . The MMC and GPU-based MCX methods by Fang et al. provide more light source modes, containing the slit beam, the planar beam, and the cone beam, with the option of sinusoidal modulated amplitude [14, 15]. For the multi-angle wide-field illumination modes usually employed in OAI systems, the conventional MC methods can not simultaneously and randomly generate the initial angles and positions of emitted photons using the GPU framework. Thus, the number of initial angles and positions is limited due to the low computation efficiency of the central processing unit (CPU) framework. The limited number of initial positions produces great challenges for efficient light transport modeling of arbitrary multi-angle wide-field illumination modes in OAI systems, such as ring, cylindrical and hemispherical sources for typical OAI systems.
In this work, we present a MC approach based on GPU to quickly and accurately establish light transport modeling for OAI systems with arbitrary multi-angle wide-field illuminations, termed as MCOAI. The incident angles and positions of emitted photons can be simultaneously and randomly generated from their pre-set probability distributions in MCOAI, allowing its effective performance for simulating arbitrary light sources with GPU acceleration. By comparing with the MCX method, the accuracy of MCOAI method is verified in the simulation with the pencil beam illumination. Comparison results using the ring source mode in optoacoustic tomography indicate that the MCOAI method performs higher efficiency and greater precision than the MCX method. The proposed MCOAI method can also effectively achieve three-dimensional (3D) light transport simulation for cylindrical and hemispherical illumination modes usually employed in OAI systems. Finally, several illumination modes are set up for a cross-section of Digimouse, illustrating the necessity of MCOAI for the quantitative OAI.
Method and Model
The MCOAI method can achieve simultaneous generation of both initial angles and positions of emitted photons by random numbers from pre-set probability distributions, or random selection from customized tables in the GPU framework. Owing to the sufficient number of initial angles and positions generated in the GPU framework, the MCOAI method can effectively conduct the precise photon transport modeling of arbitrary 3D media illuminated by arbitrary multi-angle wide-field light sources. In Cartesian coordinates, all objectives are uniformly divided into voxels. Optical properties of the media, such as absorption coefficient μa, scattering coefficient μs, scattering anisotropy factor g and refractive index n, are assigned to each voxel. The incident beam is simulated by emitting photon packets, whose trajectories are calculated in the media. Probability models and pseudorandom number generators are utilized to determine the initial position and direction, the path length between two adjacent scattering events, the scattering angle, transmission or reflection, as well as termination of the photon packet. As photons propagate from one voxel to the next, they deposit some energy (“weight”) into the voxel, depending on absorption coefficient μa of the voxel. Photons are tracked until terminated in Russian roulette or escaping from the media boundary. The deposited energy is digitally accumulated in a 3D matrix, and is normalized to output power at the end of the simulation. Figure 1 illustrates the flowchart of MCOAI for the multi-angle wide-field illumination mode.
Figure 1. Diagram of the MCOAI simulation for photon transport. The multi-angle wide-field illumination mode is established through randomly generating initialized direction and position for each emitted photon.
In MCOAI, the random number is randomly generated in [0, 1]. For the ring source with the radius of r in the plane of z = z0, a random number R0 is used to obtain the longitude angle ϕ = 2πR0 uniformly distributed in [0, 2π]. After an initialized direction ϕ is randomly generated, the initialized position rs and the direction vector can be expressed as
For the cylindrical source with the radius of r and the height of H, a random number R1 is employed to produce the longitude angle ϕ = 2πR1, and another random number R2 is used to compute the elevation position h = HR2 randomly distributed in [0, H]. The initialized position rs and the direction vector can be expressed as
For the hemispherical source with the radius of r, two random numbers R3 and R4 are used to generate the longitude angle ϕ = 2πR3 and latitude angle θ = πR4/2, respectively. The initialized position rs and direction vector can be expressed as
For other types of source mode, the initial direction and position of emitted photon can be customized aforehand, and then the random extraction is adopted to randomly obtain these parameters based on the probability distributions. For light source models with uneven light intensity distribution, packet weights at different initialized positions can be prepared beforehand. The packet weight of emitted photon at a certain initialized position rs is W(rs), depending on the light intensity at the position rs.
A total of 108 photons are launched in the simulation of the proposed MCOAI method. One hundred and twenty-eight threads are set for parallel calculation. Four simulated samples are employed to validate the performance of the MCOAI method. The optical properties of the samples are set as the representative values in the optical window . The MCOAI and MCX methods are run on the same computer with an NVIDIA GTX 750 GPU with 64 GB memory.
1) Sample 1: A homogeneous cube: A cubic domain of 60 mm × 60 mm × 60 mm is simulated shown in Figure 2. The volume is from −3 to +3 cm in x and y axis, and from 0 to 6 cm in z axis, with the number of grids in each axis of 60. The pencil beam is located at rs = (30,30,0) with an incident direction vector = (0,0,1); the medium is homogeneous with the absorption coefficient μa = 0.01 mm−1, the scattering coefficient μs = 10 mm−1, the refractive index n = 1.37 and the anisotropy g = 0.9.
2) Sample 2: A homogeneous cylinder: A cylinder sample with the center located at (6, 6, 8) mm has a radius of 5 mm and a height of 6 mm. The voxel size is 50 μm × 50 μm × 50 μm to build the accurate curved boundary. The optical properties for the cylinder are μa = 0.1 mm−1, μs = 10 mm−1, n = 1.37 and g = 0.9. A ring source with the radius of 5 mm is located at the plane of z = 8 mm, shown in Figure 3A. Figure 4A illustrates a 3D cylindrical illumination distributed around the cylinder boundary.
3) Sample 3: A homogeneous hemisphere: A hemisphere sample centered at (6, 6, 6) mm has a radius of 5 mm. The voxel size and optical properties of the hemisphere are the same as those of the Sample 2. Figure 4C shows a hemispherical illumination distributed around the hemisphere surface.
4) Sample 4: Cross-section of Digimouse: A cross-section of the widely-used Digimouse atlas has been chosen as a more realistic tissue sample. g = 0.9 and n = 1.37 are set to be homogeneous in the sample. μa and μs of the sample are shown in Table 1. Several types of light sources are simulated, including the finite-width light-sheet illumination from one side (left, right, top, or bottom) and four sides, as well as ring light-sheet illumination.
Figure 2. Comparisons between the MCOAI and MCX methods using the pencil beam to illuminate a cube (Sample 1): (A) the fluence map along the plane of y = 30 mm from the MCOAI result, (B) the fluence map along the plane of y = 30 mm from the MCX result, and (C) the fluence contours with 10 dB spacing along the plane of y = 30 mm. ϕ is the light fluence.
Figure 3. Comparisons between the MCOAI and MCX methods using the ring source to illuminate the cylinder sample (Sample 2): (A) the ring source located at the plane of z = 8 mm, (B) the fluence map at the plane of z = 8 mm from the MCOAI result, (C) the fluence map at the plane of z = 8 mm from the MCX result of 360 sources, (D) the fluence map at the plane of z = 8 mm from the MCX result of 720 sources, (E–G) the enlarged details of the fluence map in (B–D). ϕ is the light fluence.
Figure 4. The illumination modes achieved by the MCOAI method: (A) the cylindrical source mode to illuminate the cylinder sample (Sample 2), (B) the fluence map at the plane of z = 9 mm from the MCOAI result, (C) the fluence map at the plane of y = 6 mm from the MCOAI result, and (D) the corresponding 3D view of the fluence map; (E) the hemispherical source mode to illuminate the hemisphere sample (Sample 3), (F) the fluence map at the plane of z = 6 mm from the MCOAI result, (G) the fluence map at the plane of y = 6 mm from the MCOAI result, and (H) the corresponding 3D view of the fluence map.
Results and Discussion
To evaluate the accuracy of the MCOAI method, a comparison with the MCX method is implemented using the pencil beam to illuminate Sample 1. The fluence maps at the plane of y = 30 mm from the MCOAI and MCX results are presented in Figures 2A,B, respectively. The fluence contours with 10 dB spacing are also displayed in Figure 2C. The MCOAI result matches well with the MCX result in both fluence maps and contours, presenting the accurate performance of the MCOAI method. To quantitatively assess the accuracy of MCOAI result in Figure 2, the root mean square error (RMSE) has been calculated to show the difference between the MCOAI result and the MCX result. The RMSE is defined as
where N is the total number of the image pixels, rMCX(n) and rMCOAI(n) are the values of fluence maps at the nth pixel in the MCX result and MCOAI result, respectively. The RMSE value between Figures 2A,B is 0.0008, indicating high consistency of the MCOAI result and the MCX result.
To demonstrate the fast and accurate modeling of the MCOAI method for the multi-angle wide-field illumination mode, a cylindrical sample (Sample 2) in Figure 3A has been chosen, with a ring source illumination usually employed in OAI. The calculation time of the MCOAI and MCX methods is illustrated in Table 2. The run time of MCOAI is 1.7 min. In the MCX method, amounts of pencil beam sources are set at the assigned positions around the cylinder boundary. The MCX simulations with 360 and 720 pencil beam sources (108 photons for each beam) to form the ring illumination take the calculation time of 40 and 80 min, respectively. A total of 108 photons are launched for the proposed MCOAI method. For each emitted photon, the initial angle and position are randomly generated. Thus, a total of 108 initial angles and positions are randomly generated. With the sufficient number of photons (or the number of randomly generated initial angles and positions) and GPU acceleration, the computation time of the proposed MCOAI method is much faster than the MCX method.
The fluence maps at the plane of z = 3 mm from the MCOAI and MCX results are shown in Figures 3B–D. The zoom-in images at the curved boundary are, respectively, illustrated in Figures 3E–G. Although the distributions of fluence maps by the three MC simulations are generally similar in Figures 3B–D, the detailed image of fluence at the boundary in Figure 3E is obviously smoother and more authentic than the images in Figures 3F,G. Since 108 photon packets are emitted at random initialized positions along the ring shape in the MCOAI method, the finite number of initialized positions set in the MCX method is much less than the number in the MCOAI method, resulting in the uneven fluence distributions at the boundary in the MCX results. The comparison results of the focused light illumination in Figure S1 are consistent with the results of the ring source in Figure 3. The standard deviation (SD) has been calculated to analyze the regional uniformity of the MCOAI and MCX results indicated by the black squares in Figures 3E–G. The SD is defined as
where r(n) is the value of fluence map at the nth pixel and rmean is the mean value in the black square. The SD values of the MCOAI result, the MCX result of 360 sources and the MCX result of 720 sources are 0.0142, 0.0828, and 0.0520, respectively. The lowest SD value of the MCOAI result indicates the highest uniformity. The higher fidelity and efficiency of the MCOAI method can be attributed to introducing the random generation approach for the incident angle and position of each emitted photon into the parallel computing framework. Therefore, MCOAI method is more appropriate to accurately and quickly establish multi-angle wide-field illumination modes for optoacoustic tomography.
Furthermore, two typical 3D source modes (Figures 4A,E) usually existed in OAI are illustrated by MCOAI method. For the cylindrical illumination mode, the fluence maps at the planes of z = 9 mm, y = 6 mm and 3D view are presented in Figures 4B–D. For the hemispherical illumination mode, the fluence maps at the planes of z = 6 mm, y = 6 mm and 3D view are displayed in Figures 4F–H. The boundary and interior distributions of fluence maps in two homogeneous media are smooth and authentic owing to the sufficient incident positions generated randomly. The results demonstrate the light transport modeling in 3D simulation with multi-angle wide-field illumination modes can be efficiently and accurately established by the proposed MCOAI method.
This MCOAI approach is further used to simulate various light illumination modes on a cross-section of Digimouse (Sample 4) shown in Figure 5. The optoacoustic tomography technology has been widely employed for the whole-body small-animal imaging by integrating lots of cross-section images. The cross-section of Digimouse containing several main organs listed in Table 1 is selected.
Figure 5. Simulation of the Sample 4 under different light illumination modes by the MCOAI method: (A) finite-width light-sheet illumination from left (a), right (b), top (c), and bottom (d) side, as well as four sides (abcd); (B) ring light-sheet illumination; (C) pseudo-color image of the absorption coefficient μa; (D) pseudo-color image of the scattering coefficient μs; (E) images of initial acoustic pressure P0 under the finite-width light-sheet illumination from each side as well as four sides, and under the ring light-sheet illumination; (F) images of fluence ϕ under the finite-width light-sheet illumination from each side as well as four sides, and under the ring light-sheet illumination; (G) the corresponding fluence contours (with 10 dB spacing). Scale bar, 2 mm.
The MCOAI method was further utilized to illustrate the necessity for quantitatively evaluating functional information of biological tissue from OAI results. Figures 5C,D present the grayscale images of μa and μs of Sample 4, indicating the heterogeneity in different organs. Figure 5E shows the images of initial acoustic pressure P of Sample 4 using six different illumination modes. As displayed in the first four columns of Figure 5E using one-sided illumination, the intensity distribution of targets and background in the P images plummet when moving away the light source. In the last two columns of Figure 5E, due to the illumination modes are both full-angle wide-field, the P images are closer to the μa images while still presenting obvious differences. The results show the P image cannot directly reflect accurate quantitative distributions of μa in the sample, affected by the light attenuation along tissue depth. The fluence maps under these six different illumination modes are shown in Figure 5F, and the contours with 10 dB spacing are also displayed in Figure 5G. The light fluence ϕ in deep tissues is attenuated along the penetration depth, thus the fluence distributions in the region near the light source are higher than those in the region far from the light source. Under the finite-width light-sheet illumination modes only from one side, the heterogeneous distribution of ϕ is extremely evident. Under the four-sided light-sheet and ring light-sheet illumination modes, the Sample 4 receives full-angle wide-field irradiation, the fluence maps ϕ are attenuated from boundary to center. The results demonstrate the MCOAI method should be utilized to provide the fluence distribution of ϕ, in order to quantitatively recover the distribution of optical properties.
To quantitatively evaluate the optical properties of the biological tissues, it becomes essential to develop an efficient and precise photon transport modeling method for multi-angle wide-field illumination modes usually employed in OAI systems. In this work, we present the MCOAI method to randomly generate the incident angles and positions of emitted photons based on parallel computing framework for efficiently simulating arbitrary wide-field light sources. The MCOAI method can offer the same accuracy in the simulation using pencil beam as the MCX method. However, the MCOAI method for the multi-angle wide-field light sources, such as the ring source, is capable of providing more accurate and authentic distributions of light fluence in a few tenths of the time, compared to the MCX method. The GPU framework has been utilized in both the MCX method and the MCOAI method for acceleration. For the pencil beam and planar beam in the MCX method, the initial angles and positions of emitted photons are generated in the GPU framework. However, for the multi-angle wide-field illumination modes usually employed in OAI systems, the MCX method can not simultaneously and randomly generate the initial angles and positions of emitted photons in the GPU framework. The proposed MCOAI method can achieve simultaneous generation of both initial angles and positions of emitted photons by random numbers in the GPU framework, greatly reducing the running time. If the optical properties of samples at the selected wavelength in visible or near-infrared spectrum can be determined, the light transport modeling at any selected wavelength can be achieved by the proposed MCOAI method. The proposed method has been demonstrated with typical illumination modes and uniform light intensity, which can be explored on user-defined illumination modes and uneven light intensity for further study. The MCOAI method is able to precisely conduct the light transport modeling for all typical OAI systems, having great potential of contributing to the further investigation of quantitative OAI methods. In this work, the MCOAI method focuses on providing an accurate light transport simulation for the forward modeling in quantitative OAI, not involving the inverse problems. The MCOAI method is designed to be an open source for user-friendly employment.
Data Availability Statement
Publicly available datasets were analyzed in this study. The data of this article will be made available by the authors to any qualified researcher.
JL and FG guided the work. TL did the simulation. TC, SM, SL, and XX gave the guide for analysis. JL and TL wrote the manuscript. All authors contributed to the article and approved the submitted version.
This work was supported by the grants from the National Natural Science Foundation of China (Grant Nos. 81771880, 81401453, 81671728, 81871393, and 81971656) and Tianjin Municipal Government of China (19JCQNJC12800).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2020.00283/full#supplementary-material
2. Liu Y, Nie L, Chen X. Photoacoustic molecular imaging: from multiscale biomedical applications towards early-stage theranostics. Trends Biotechnol. (2016) 34:420–33. doi: 10.1016/j.tibtech.2016.02.001
5. Lu T, Wang Y, Zhang S, Li J. Photoacoustic mesoscopy:pushing toward the depth limit in the high-resolution optical imaging for biomedical applications and clinical potentials. Instrumentation. (2016) 3:29–42.
9. Wang Y, He J, Li J, Lu T, Li Y, Ma W, et al. Toward whole-body quantitative photoacoustic tomography of small-animals with multi-angle light-sheet illuminations. Biomed Opt Express. (2017) 8:3778–95. doi: 10.1364/BOE.8.003778
Keywords: optoacoustic imaging, quantitative visualization, multi-angle wide-field illumination, Monte Carlo, graphics processing unit
Citation: Lu T, Li J, Chen T, Miao S, Li S, Xu X and Gao F (2020) Parallelized Monte Carlo Photon Transport Simulations for Arbitrary Multi-Angle Wide-Field Illumination in Optoacoustic Imaging. Front. Phys. 8:283. doi: 10.3389/fphy.2020.00283
Received: 15 April 2020; Accepted: 23 June 2020;
Published: 31 July 2020.
Edited by:Chao Tian, University of Science and Technology of China, China
Reviewed by:Zhen Yuan, University of Macau, China
Xiaojing Gong, Shenzhen Institutes of Advanced Technology (CAS), China
Copyright © 2020 Lu, Li, Chen, Miao, Li, Xu and Gao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jiao Li, firstname.lastname@example.org