Fast photon-boundary intersection computation for Monte Carlo simulation of photon migration

Abstract. Monte Carlo (MC) method is generally used as a “gold standard” technique to simulate photon transport in biomedical optics. However, it is quite time-consuming since abundant photon propagations need to be simulated in order to achieve an accurate result. In the case of complicated geometry, the computation speed is bound up with the calculation of the intersection between the photon transmission path and media boundary. The ray-triangle-based method is often used to calculate the photon-boundary intersection in the shape-based MC simulation for light propagation, but it is still relatively time-consuming. We present a fast way to determine the photon-boundary intersection. Triangle meshes are used to describe the boundary structure. A line segment instead of a ray is used to check if there exists a photon-boundary intersection, as the next location of the photon in light transports is determined by the step size. Results suggest that by simply replacing the conventional ray-triangle-based method with the proposed line segment-triangle-based method, the MC simulation for light propagation in the mouse model can be speeded up by more than 35%.


Introduction
Optical molecular imaging is developing rapidly as a promising technology in the molecular imaging field. The establishment of a mathematical model of light propagation in biological media is essential in optical tomography. Diffuse equation (DE) is commonly used as the mathematical model, however, solving DE needs certain boundary conditions. 1,2 A Monte Carlo (MC) method is relatively flexible, as it does not require an assumption on boundary conditions. MC in the optical field is basically a statistical method that allows for constructing a stochastic model to simulate a large amount of light propagation trajectories. Since first introduced by Wilson and Adam into the field of laser-tissue interactions, 3 it has been widely used to simulate the light propagation in media. [4][5][6][7][8][9] Much work has been done to improve the efficiency of a MC simulation for modeling light propagation in turbid media. [10][11][12][13][14] Based on the information generated by a single MC baseline simulation, the scaling methods were developed for a wide range of optical properties. 15, 16 Liebert et al. applied a similar concept to simulate the time-resolved fluorescence imaging in layered turbid media using the MC method. 17 However, a single MC process is still the most time-consuming part owing to a large number of photon trajectory calculations.
In order to improve the efficiency of a single MC simulation for light propagation, Wang et al. developed an MC modeling of photon transmission in multi-layered tissues (MCML) 18 in 1995, while Tinet et al. 10 proposed a semianalytical Monte Carlo method for time-resolved light propagation. Each interaction point of a photon that transports in a medium contributes directly to a detector area, thus fewer photons are required to produce a given accuracy but with complex calculations. Zolek et al. 11 developed an optimized MC algorithm by approximation of the logarithmic and trigonometric functions. Boas et al. proposed a voxel-based model (tMCimg) to improve the capability of modeling turbid structures. 12 However, Binzoni et al. pointed out that the voxel-based model (tMCimg) was difficult to incorporate precise boundary information. 13 Margallo-Balbas et al., therefore, explored a shape-based method (triMC3D) 14 and used triangle meshes to define the complex boundary structures. This shape-based method gives a more accurate description of the interface than the voxel-based approaches. However, in the case with complicated geometry, the computation speed is bound up with the calculation of the photon-boundary intersection. In order to judge whether the photon hits its boundary, the ray-trianglebased method proposed by Möller et al. 19 is often used to calculate the photon-boundary intersection. 14,20 Even though this method can save the memory required, it needs to calculate the intersection between the ray and every triangle; therefore it is quite time-consuming. In fact, an intersection test between a line segment and every triangle can be taken before the calculation of the photon-boundary intersection, as the next location of the photon in light transports is determined by the step size. Thus, the photon-boundary intersection needs to be calculated only if an intersection exists within a triangle. When there is more than one intersection, it means that there are many boundaries along the line segment. Then only the intersection with the minimum distance to the photon's current position is calculated. This method is more effective than the method that considers all photon-boundary intersections, and makes it possible to reduce computation time.
In this paper, a fast way to determine the photon-boundary intersection in the shape-based MC simulation is presented. Triangle meshes are used to describe the boundary structure. A line segment, instead of the ray 14,20 is used to check if a photon-boundary intersection exists, as the intersection will be calculated only when it does exist. Therefore, only a part of triangles have to be taken into account when calculating the photon-boundary intersection, which greatly reduces the computation time. The conventional method and the proposed method for the intersection test in MC simulation for light propagation are detailed first. Then, the validation of the proposed algorithm is conducted in a semi-infinite medium. The fluence distribution within the medium and the remitted flux of photons obtained from MC with the proposed method are compared with those obtained from the analytic solution of DE and the conventional MC method, respectively. Finally, several sets of MC simulations under different conditions are executed to test the computational speed of the proposed method.
The paper is organized with the conventional ray-trianglebased method and the proposed line segment-triangle-based methods described in Sec. 2. In Sec. 3, three sets of simulation studies are conducted to validate the algorithm and test the computational speed. The results are shown and discussed. In the last section, the conclusion is presented.

The Ray-Triangle-Based Method
A ray RðtÞ with origin O and normalized direction D is defined as: (1) where t is the unknown distance the ray travels. And a triangle is defined by three vertices V 0 , V 1 , and V 2 . A point, Tðu; vÞ, in the plane where the triangle lies is given as: where ðu; vÞ are the coordinates within the triangle which satisfy u ≥ 0, v ≥ 0 and u þ v ≤ 1. The intersection between the ray, RðtÞ, and the triangle, Tðu; vÞ, can be calculated from Eqs. (1) and (2) as follows: Reordering the formula, we can get 19 where D ¼ s∕μ t ·ñ μ , s is the free path of the photon,ñ μ is the direction vector of the photon migration and μ t is the total interaction coefficient equal to the sum of the absorption coefficient and the scattering coefficient. It means that t, the distance from the ray origin to the intersection point, and ðu; vÞ, the intersection coordinates, can be obtained by solving Eq. (4). 19 But this approach is very time-consuming. The ray-triangle intersection algorithm tries to determine if the ray intersects with the triangle. If the distance t is in the range of (0,1), the ray-triangle intersection is within the distance of the next free path of the photon. In addition, when ðu; vÞ satisfies u ≥ 0, v ≥ 0 and u þ v ≤ 1, the ray-triangle intersection is within the triangle. Only when the solution meets these two conditions above at the same time, the next free path of the photon will intersect with the triangle defined by the three vertices V 0 , V 1 and V 2 . In the raytriangle-based method, solving Eq. (4) is implemented as a translation and change of base of the ray origin and a unit triangle. 19 The procedure of the conventional ray-trianglebased method, when determining whether the photon hits the boundary in the MC simulation, was demonstrated in a flowchart, as shown in Fig. 1. Fig. 1 Flowchart of the ray-triangle-based method.

Line-Segment Triangle-Based Method
However, the ray-triangle-based method has to take all triangles into account when calculating the intersection, which increases the computation time. In fact, when we refer to the path of the photon, it is a line segment rather than a ray, since every step of the photon migration is limited by the step size, as shown in Fig. 2. By taking advantage of this, the proposed method uses the line segment instead of the ray to compute the intersection. Therefore, the proposed method first checks whether a line segment-triangle intersection exists. Only when the intersection exists, can the calculation of the intersection continue.
The standard equation of the plane in which the triangle V 0 V 1 V 2 lies is given by where ðA; B; CÞ is the unit normal vector of the triangle and DD is the perpendicular distance from the origin to the plane. Before the start of photon transport simulation, all the standard equations of the planes in which the triangles lie are solved and stored. By incorporating Eqs. (1) and (5), the following equation is obtained when t is in the range of (0, 1), the line segment intersects with the plane in which the triangle lies.
If an intersection does exist, the next step is to determine whether the intersection is within the triangle. As illustrated in Fig. 3, supposing there is a triangular prism with the triangle that describes the media boundary being the base and three rectangular sides perpendicular to the triangle, the problem is transformed to determine whether the point of intersection is inside the triangular prism. If the point of intersection is within the triangular prism, the intersection is also within the triangle in the plane. As shown in Fig. 3, two line segments with the same origin O intersect with the plane in which a triangle (the blue part) lies respectively. The points of intersection are P1 and P2, and P2 is the point where photon hits the boundary.
The equations of three rectangular sides are defined respectively as: where ðS 0 ; S 1 ; S 2 Þ, ðS 4 ; S 5 ; S 6 Þ, ðS 8 ; S 9 ; S 10 Þ are the unit normal vector of the three rectangular sides respectively, and S 3 , S 7 , S 11 are the perpendicular distances from the origin to each plane. Coordinates of the median point, ðx 0 ; y 0 ; z 0 Þ, which are inside the triangle, can be obtained from the three vertices. The intersection point, ðx i ; y i ; z i Þ, can be solved by substituting t into Eq. (1). If conditions in Eq. (8) are all met, the intersection is inside the triangular prism.
The procedure of the proposed line segment-trianglebased method in calculating the photon-boundary intersection in MC simulation was illustrated in a flowchart, as shown in Fig. 4.
To sum up, the key concept of the line segment-triangle method is based on the idea that the step size can reduce the possible triangles intersected with the photon trajectory. Besides, the proposed algorithm can work out each intersection coordinates with fewer calculations. Thus the computational consumption of the proposed method can be greatly reduced.
The MC model with the proposed method is coded in C programming language and compiled with Microsoft Visual C++ 6.0. For comparison, both the MC methods were performed on the same PC workstation (2.00 GHz  Intel Xeon processor, 24 GB RAM) and used the same compiler. The analytical solution of the DE is calculated with MATLAB 7.6.0.324 (R2008a, The MathWorks). The Henyey-Greenstein (HG) phase function is employed for sampling scattering directions in the simulations. 20 Fresnel's equations are used to calculate the reflections and refractions at the boundaries where the refractive index mismatch occurs. These internal functions are just the same as in the conventional MC method.

Validation of the Algorithm in a Semi-infinite Medium
In this work, we simulated the photon propagation in biological medium using the well-established MC procedures. 9,10,12,14,21-23 Two sets of simulation studies were conducted to validate the proposed algorithm in a semi-infinite medium in this section. The first set of simulation studies was to compare the results obtained from the proposed method in MC simulation with those obtained from the analytic solution of DE. Owing to the diffusion approximation, 1,24-26 the simulations were executed in a relatively ideal model available for DE to yield an analytic solution accurate enough to be a validation standard. 12,21,22 The second set was to validate the proposed algorithm by comparing the results with those obtained from the conventional MC method. Compared with DE, the main advantage of the MC method is that it does not require an extra simplifying assumption. Thus this makes it capable of dealing with any situation. Considering the time consumed, the second set of simulations were conducted in a homogenous medium of the same geometry and the same source settings as that used in the first set of simulations. The optical properties were set according to that of biological tissues. There were 10 7 photons launched in both MC simulation sets.

Comparison to the analytical solution of DE
The structure of the medium used in the simulations was a cubic model whose surface was described by 1794 triangles generated by a commercially available software COMSOL Multiphysics 3.5 (COMSOL, Inc., Burlington, MA). The medium had a dimension of 60 × 60 × 40 mm 3 with its center positioned at ð0; 0; 0Þ mm. The optical parameters of the homogeneous medium were as follows: absorption coefficient μ a ¼ 0.005 mm −1 , scattering coefficient μ s ¼ 1.00 mm −1 , relative refractive index n ¼ 1.33, anisotropy factor g ¼ 0.1. The ambient medium is the air with refractive index n ¼ 1. A collimated point source was positioned at ð0; 0; 20Þ mm with the direction vector of ð0; 0; −1Þ.
The medium was approximately semi-infinite as the source was adequately far from the edges. The fluence distribution within the medium and the remitted flux of photons of the proposed MC method and the DE method are compared in Fig. 5. A decent agreement between both methods is observed. In Fig. 5(b), the mismatch near the source is due to the isotropic point source hypothesis within the DE. The slight discrepancy near the edge at 30 mm results from the influence of boundary conditions within the DE for the differences will narrow when enlarging the medium size. 12,21

Comparison to the results obtained from the conventional MC method
In this section, the comparison between the results from the MC simulation with the proposed method and with the conventional method was made. The comparison was conducted under the same settings as that used in Sec. 3.1.1, such as geometry size, pseudorandom generator, source settings and photon number. The optical parameters of the biological tissue were as follows: absorption coefficient μ a ¼ 0.005 mm −1 , scattering coefficient μ s ¼ 1.00 mm −1 , relative refractive index n ¼ 1.33, anisotropy factor g ¼ 0.9. The ambient medium was the air with refractive index n ¼ 1.
The fluence distribution within the medium and the remitted flux of photons of the proposed method and the conventional method are compared in Fig. 6. A good agreement between both methods is observed. Figure 6(a) shows the contour plot in the plane y ¼ 30 mm of the fluence distribution within the medium obtained from each method. The contour interval is 5 dB. The slight discrepancy is due to the number of photons used in the simulations since the MC method is a statistical model based on the stochastic nature of light propagation. Figure 6(b) compares the remitted flux of photons on the surface along the distance from the source. It should be noted that since both methods use surface triangles as detectors in the simulations, the remitted flux value for a target point is derived from the nearby triangles through interpolation. This helps to smooth the data to a certain degree.

Comparisons of Computation Time in Complex 3-D Geometry
In the following experiments, we compared the computation time of the proposed method and the conventional method. The three-dimensional (3-D) geometric structure of a mouse model obtained from an in-vivo experiment was used in this simulation. The torso part from the oxter to the abdomen of the mouse model with a height of 40 mm was chosen as the region to be investigated. The geometry was discretized in COMSOL Multiphysics using triangle meshes, as shown in Fig. 7. The optical parameters of the mouse model were listed in Table 1. The surrounding medium was still air with refractive index n ¼ 1. A collimated point source was positioned at ð10; 0; 10Þ mm with the direction vector of ð−1; 0; 0Þ. In addition to the geometry and optical characteristics of the medium, the total number of the triangle meshes describing the boundary of the mouse model also had a great impact on the computational speed. In our study, we compared the computation time in four cases with different numbers of triangle meshes when simulating  The results are shown in Table 2. We also compared the computational time of the conventional and proposed methods with four different numbers of simulated photons: PhotNum ¼ ð10 4 ; 10 5 ; 10 6 ; 5 × 10 6 Þ, when 6084 triangle meshes were used. The results are shown in Table 3. As shown in Tables 2 and 3, the computation time can be decreased by 35% to 45% by simply replacing the conventional method with the proposed method in MC simulations for light propagation in the mouse model. The proposed line segment-triangle-based method is equivalent to the conventional ray-triangle-based method in mathematics, which guarantees accuracy but the former is faster than the latter in the shape-based MC simulation for light propagation.

Conclusion
In this paper, we proposed a fast photon-boundary method in MC simulation of light propagation in complicated geometry. The proposed method uses the line segment instead of the ray in the conventional method to check if there exists a photon-boundary intersection. The MC algorithm using the proposed method was validated in a semi-infinite medium by comparing the results with those from the analytical solution of the DE and with those obtained from the MC simulation using the conventional method. Finally, the computation time of MC simulations for light propagation in complex 3-D geometry using the proposed method and the conventional method are compared, under the conditions of four different numbers of triangle meshes and four different sets of simulated photons. The results indicate that the MC simulation with the proposed method for intersection test is faster than that with the conventional method. In our future work, the tetrahedral-based model introduced in both Refs. 22 and 23 will be considered for its higher efficiency in the intersection test. In addition, hardware acceleration based on graphics processing units (GPU) 19,21,27 is gradually applied in MC simulation. The proposed method may accelerate the MC simulation of photon propagation more significantly when implemented in the GPU. In the future, we will also try to apply the presented method in GPU-based acceleration for MC simulation.