## 1.

## 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 semi-analytical 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-triangle-based 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-triangle-based 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.

## 2.

## Method

## 2.1.

### The Ray-Triangle-Based Method

A ray $R(t)$ with origin $O$ and normalized direction $D$ is defined as:

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\ge 0$, $v\ge 0$ and $u+v\le 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}

## (4)

$$[\begin{array}{ccc}-D,& {V}_{1}-{V}_{0},& {V}_{2}-{V}_{0}\end{array}]\left[\begin{array}{l}t\\ u\\ v\end{array}\right]=O-{V}_{0},$$^{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\ge 0$, $v\ge 0$ and $u+v\le 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 ray-triangle-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-triangle-based method, when determining whether the photon hits the boundary in the MC simulation, was demonstrated in a flowchart, as shown in Fig. 1.

## 2.2.

### 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:

## (7)

$${S}_{0}\xb7x+{S}_{1}\xb7y+{S}_{2}\xb7z+{S}_{3}=0\phantom{\rule{0ex}{0ex}}{S}_{4}\xb7x+{S}_{5}\xb7y+{S}_{6}\xb7z+{S}_{7}=0\phantom{\rule{0ex}{0ex}}{S}_{8}\xb7x+{S}_{9}\xb7y+{S}_{10}\xb7z+{S}_{11}=0,$$## (8)

$$({S}_{0}\xb7{x}_{i}+{S}_{1}\xb7{y}_{i}+{S}_{2}\xb7{z}_{i}+{S}_{3})\phantom{\rule{0ex}{0ex}}\xb7({S}_{0}\xb7{x}_{0}+{S}_{1}\xb7{y}_{0}+{S}_{2}\xb7{z}_{0}+{S}_{3})>0\phantom{\rule{0ex}{0ex}}({S}_{4}\xb7{x}_{i}+{S}_{5}\xb7{y}_{i}+{S}_{6}\xb7{z}_{i}+{S}_{7})\phantom{\rule{0ex}{0ex}}\xb7({S}_{4}\xb7{x}_{0}+{S}_{5}\xb7{y}_{0}+{S}_{6}\xb7{z}_{0}+{S}_{7})>0\phantom{\rule{0ex}{0ex}}({S}_{8}\xb7{x}_{i}+{S}_{9}\xb7{y}_{i}+{S}_{10}\xb7{z}_{i}+{S}_{11})\phantom{\rule{0ex}{0ex}}\xb7({S}_{8}\xb7{x}_{0}+{S}_{9}\xb7{y}_{0}+{S}_{10}\xb7{z}_{0}+{S}_{11})>0.$$The procedure of the proposed line segment-triangle-based 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.

## 3.

## Results and Discussion

## 3.1.

### 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}22.^{–}^{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}25.^{–}^{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.

## 3.1.1.

#### 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\times 60\times 40\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{3}$ with its center positioned at $(0,0,0)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. The optical parameters of the homogeneous medium were as follows: absorption coefficient ${\mu}_{a}=0.005\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, scattering coefficient ${\mu}_{s}=\phantom{\rule{0ex}{0ex}}1.00\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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}

## 3.1.2.

#### 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 ${\mu}_{a}=\phantom{\rule{0ex}{0ex}}0.005\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$, scattering coefficient ${\mu}_{s}=1.00\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## 3.2.

### 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)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 with ${10}^{4}$ photons: $\mathrm{TriNum}=(4952,5164,6084,13602)$. 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: $\mathrm{PhotNum}=({10}^{4},{10}^{5},{10}^{6},5\times {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.

## Table 1

The optical parameters of the mouse model.

μa | μs | n | G |
---|---|---|---|

0.03 mm−1 | 1.00 mm−1 | 1.33 | 0.9 |

## Table 2

Comparison of computation time of Monte Carlo simulations with different numbers of triangle meshes and the same number of photons (104).

Number of triangle meshes | Computation time of the conventional method (s) | Computation time of the proposed method (s) | Percentage improvement |
---|---|---|---|

4952 | 282 | 181 | 35.82% |

5164 | 302 | 183 | 39.40% |

6084 | 357 | 227 | 36.41% |

13602 | 856 | 532 | 37.85% |

## Table 3

Comparison of computation time of Monte Carlo simulations with different number of photons and the same number of triangle meshes (6084).

Number of photons | Computation time of the conventional method (s) | Computation time of the proposed method (s) | Percentage improvement |
---|---|---|---|

104 | 357 | 227 | 36.41% |

105 | 3913 | 2457 | 37.21% |

106 | 39279 | 24274 | 38.20% |

5×106 | 221150 | 121489 | 45.06% |

## 4.

## 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.

## Acknowledgments

This work is supported by the National Basic Research Program of China (973) under Grant No. 2011CB707701; the National Natural Science Foundation of China under Grant No. 81071191, 81271617; the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114; National Science and technology support program under Grant No. 2012BAI23B00.

## References

*in vivo*,” Med. Phys. 19(4), 879–888 (1992).MPHYA60094-2405 http://dx.doi.org/10.1118/1.596777 Google Scholar

## Biography

**Xiaofen Zhao** received her BS in biomedical engineering from Jilin University, Jilin, China, in 2010. She is currently a postgraduate student in the Biomedical Engineering Department at Tsinghua University, Beijing, China. Her research interest is in fluorescence molecular tomography for small animal imaging.

**Hongyan Liu** received her BS in biomedical engineering from Tsinghua University, Beijing, China, in 2011. Her research interest consists of fluorescence molecular tomography for small animal imaging.

**Bin Zhang** received his MS in mechanical and electronic engineering from University of Science and Technology of China in 2007. From 2007 to 2009, he was a research assistant in Shenzhen Institute of Advanced Integration Technology, Chinese Academy of Sciences. He is now a PhD candidate in the Department of Biomedical Engineering, Tsinghua University, Beijing, China. His research interest centers on fluorescence molecular tomography for small animals.

**Fei Liu** received her BS in biomedical engineering from Zhejiang University, Zhejiang, China, in 2008. She is currently a PhD candidate in the Biomedical Engineering Department, Tsinghua University, Beijing, China. Her research interest is in fluorescence molecular tomography for small animal imaging.

**Jianwen Luo** received his BS and PhD (with honors) in biomedical engineering from Tsinghua University, Beijing, China, in 2000 and 2005, respectively. He was a postdoctoral research scientist from 2005 to 2009, and an associate research scientist from 2009 to 2011 in the Department of Biomedical Engineering at Columbia University, New York. He became a professor at Tsinghua University in 2011, and was enrolled in the Thousand Young Talents Program in 2012. He serves as an Advisory Editorial Board member of the Journal of Ultrasound in Medicine, Associate Editor of Medical Physics, and Faculty Member of Faculty of 1000. His research interest is in biomedical imaging, including ultrasound and elasticity imaging.

**Jing Bai** received her MS and PhD from Drexel University, Philadelphia, Pennsylvania, in 1983 and 1985, respectively. From 1985 to 1987, she was a research associate and assistant professor at the Biomedical Engineering and Science Institute, Drexel University. In 1988, 1991, and 2000, she became an associate professor, professor, and Cheung Kong chair professor at Biomedical Engineering Department of Tsinghua University, Beijing, China. From 1997, she has been an associate editor for *IEEE Transactions on Information Technology in Biomedicine*. Her current research interests include mathematical modeling and simulation of the cardiovascular system, optimization of cardiac assist devices, medical ultrasound, telemedicine, home health care network and home monitoring devices, and infrared imaging. She has authored or coauthored 10 books and more than 300 journal papers.