## Abstract

Modeling light propagation in the whole body is essential and necessary for optical imaging. However, non-scattering, low-scattering and high absorption regions commonly exist in biological tissues, which lead to inaccuracy of the existing light transport models. In this paper, a novel hybrid light transport model that couples the simplified spherical harmonics approximation (SP_{N}) with the radiosity theory (HSRM) was presented, to accurately describe light transport in turbid media with non-scattering, low-scattering and high absorption heterogeneities. In the model, the radiosity theory was used to characterize the light transport in non-scattering regions and the SP_{N} was employed to handle the scattering problems, including subsets of low-scattering and high absorption. A Neumann source constructed by the light transport in the non-scattering region and formed at the interface between the non-scattering and scattering regions was superposed into the original light source, to couple the SP_{N} with the radiosity theory. The accuracy and effectiveness of the HSRM was first verified with both regular and digital mouse model based simulations and a physical phantom based experiment. The feasibility and applicability of the HSRM was then investigated by a broad range of optical properties. Lastly, the influence of depth of the light source on the model was also discussed. Primary results showed that the proposed model provided high performance for light transport in turbid media with non-scattering, low-scattering and high absorption heterogeneities.

© 2013 Optical Society of America

## 1. Introduction

Three-dimensional (3D) optical imaging in preclinical research has been rapidly developing and extensively applied in biomedical fields with the advantages of high sensitivity, simple operation and low cost [1–3]. The foundation of 3D optical imaging is developing accurate and feasible reconstruction algorithms to reflect the information of the light source inside biological tissues at the molecular level [4]. However, most of the 3D optical imaging reconstruction algorithms rely on accurate and feasible forward models of light propagation in biological tissues [4–6].

The radiative transfer equation (RTE) has been addressed to accurately model light propagation in a turbid medium [5, 6]. However, the RTE is hard to solve in complicated 3D irregular medium and needs extensive computational time [4, 6]. Therefore, several approximations of the RTE were applied to model the light transport in a turbid medium, such as the diffusion equation (DE) [4, 7–9], the discrete ordinates equation (S_{N}) [10, 11], the spherical harmonics equation (P_{N}) [12, 13], the phase approximation (PA) [14] and the simplified spherical harmonics approximation equation (SP_{N}) [15–17]. The high order approximations, such as the S_{N}, P_{N} and PA, can accurately model light propagation in a scattering medium, but the computational burden limits their application in complex heterogeneities. The DE, however, provides a high computational efficiency and is only valid in the high scattering region [18]. To balance computational accuracy and efficiency, SP_{N} has been recently presented for 3D optical imaging, which provides high accuracy for light transport in low scattering or high absorption regions [15–17]. Nevertheless, all of the highly effective approximations become invalid in non-scattering regions.

Previous studies have demonstrated that an inaccurate forward model and reconstructed images would be obtained if the non-scattering region was not well processed [19–22]. To address the problem, several solutions have been proposed to deal with light transport in the non-scattering region, such as the hybrid radiostiy-diffusion method (HRDM) [19–26], the hybrid RTE-diffusion method [27–29], and the hybrid Monte Carlo-diffusion method [30]. Because of the high computational efficiency, the HRDM has caught the researchers' eye and obtained successful applications in 3D optical imaging [19–24]. Recently, the HRDM has been successfully applied to the reconstruction of bioluminescence tomography in the application of gastric cancer detection [26]. However, because of the characteristics of diffusion equation, the HRDM limits its applications to specific problems in which the high scattering regions were coupled with non-scattering ones. In a biological body, low-scattering or high absorption regions always exist simultaneously with non- and high-scattering tissues. For example, the non-scattering stomach is partly surrounded by the low-scattering liver (in the spectrum of 690*nm*) [31, 32]. Thus, HRDM is not a perfect solution to whole-body small animal imaging.

To provide an ideal solution for whole-body small animal imaging, in which high-, low- and non-scattering or high absorption tissues simultaneously exist in the solving domain, a novel hybrid light transport model that couples SP_{N} with the radiosity theory was presented in this paper, called HSRM for short, which is the abbreviation of the hybrid SP_{N}-radiosity method. In the presented model, the radiosity theory was used to characterize the light transport in non-scattering regions and SP_{N} was selected to handle the scattering problems. The hybrid model was coupled by a Neumann source which was constructed by the light transport in the non-scattering region and formed at an interface between the non-scattering and scattering regions. The Neumann source was superposed onto the original light source to construct a new synthetic light source. The accuracy and effectiveness of the HSRM were first verified with both regular and digital mouse model based simulations and a physical phantom based experiment. The feasibility and applicability of the HSRM was then investigated by the simulations performed over a broad range of optical properties. Lastly, the influence of depth of the light source on the model was also discussed.

## 2. Methods

In the HSRM method, SP_{N} was used to depict the light transport in the scattering region. Based on RTE, the three dimensional SP_{N} was obtained after a series of inferential reasonings in the planar geometry with the spherical harmonics approximation [15]. The expression of the SP_{N} can be described as follows [15, 33]:

For the void problem, the radiosity theory was used to model light propagation in the non-scattering region. The hybrid model was coupled by a Neumann source which was constructed by the light transport in the non-scattering region and formed at an interface between the non-scattering and scattering regions. Therefore, on the interface between the non-scattering and scattering region, the inward directed light distribution at $r$ was created by the light current ${J}^{+}({r}^{\prime})$ at ${r}^{\prime}$, if $r$ and ${r}^{\prime}$ are visible mutual. As a result, the contribution of the light distribution from a non-scattering region to the scattering region could be recognized as the equivalent of a Neumann source. In this paper, we defined the equivalent source as [34]:

*N*and could be deduced according to Eq. (2) and Eq. (3). The specific expression for the third order is presented as follows:

At the interface between the non-scattering and scattering regions, an equivalent Neumann source was formed and superposed into the original light source. Thus, a new synthetic light source $Q$was constructed for the hybrid light transport model, with the following expression:

where $S$is the original light source in the scattering region, $q$ is the equivalent Neumann source on the interface between the non-scattering and scattering region, and $Q$is the new synthetic light source for the hybrid light transport model. Thus, substituting the original light source $S$ with the new synthetic light source $Q$ in Eq. (1), the HSRM can be modeled by the following expressions:Because of the excellent performance of the FEM, it has been widely applied in optical imaging [7, 16, 17, 20–22, 24–29]. The HSRM-based light transport model was also solved by the FEM. So the composite moments ${\phi}_{i}$ and the original light source $S$ can be approximated by the piecewise linear bases:

where ${\varphi}_{i,k}$ and ${s}_{i,k}$ are the values at the*k*th node, and ${\psi}_{k}(r)$ is the piecewise linear nodal basis function. Then, the corresponding variational forms of the HSRM are derived as the following expression:

After discrete processing, the finite element formulation of the HSRM can be obtained as the following matrix equation:

The$M$,$\Xi $and$F$could be expressed as:

## 3. Experiments and results

In this section, to validate the accuracy and the feasibility of the proposed HSRM model, three groups of experiments were conducted. Firstly, the accuracy of the HSRM model was three-dimensionally validated with three numerical simulations, including using both the regular geometries and a digital mouse model. To further verify the accuracy of the HSRM model, a cubic phantom based experiment was also conducted. Secondly, using the three-dimensional cylindrical model, we systematically investigated the applicability of the HSRM model with simulations performed over a broad range of optical properties. Lastly, a series of numerical simulations with different depth of the light source were carried out to observe the influence of the depth of the light source on the surface light flux distribution.

Previous studies have demonstrated that the third order of SP_{N} can obtain an ideal balance between accuracy and computation time [16, 17]. Thus, the third order of the HSRM was evaluated in the experimental section. In the simulations, all of the results were compared with those of MOSE, a Molecular Optical Simulation Environment software which is based on the Monte Carlo method used to simulate light transport in turbid media [35, 36], to validate the accuracy of the HSRM. In the physical phantom experiment, the result was compared with the measured light flux map. Furthermore, the HRDM was selected as a reference to illustrate the superiority of the HSRM. To quantitatively evaluate the difference between the calculation and the standard, the average relative error (*ARE*) was introduced as:

*f*is the value of the intensity at the

_{i}*i*th sample point;

*N*is the total number of sample points; the superscript

*std*denotes that the intensity was obtained from MOSE or measurements and

*cal*represents that the intensity was calculated by the HSRM or HRDM.

*ARE*was used to reflect the discrepancy between the calculation of HSRM or HRDM and the standard of MOSE or measurements. It means that the closer

*ARE*gets to zero, the better the performance of the calculated method.

#### 3.1 Accuracy verification experiments

In this subsection, three groups of verification experiments were conducted to validate the accuracy of the HSRM, including regular shape geometry and digital mouse model based simulations and a physical cubic phantom based imaging experiment.

### 3.1.1 3D regular geometry based simulations

Two regular shape geometries were designed to verify the accuracy of the HSRM, including the spherical and cylindrical geometries, as shown in Fig. 1(a)-1(b). The related geometrical parameters and optical properties are listed in Table 1, where Model NO.1 and NO.2 correspond to the sphere and cylinder based simulations respectively. To ensure the reliability of the simulated result of MOSE, a 10^{8} photon simulation was performed. Comparisons among the results of the HSRM, HRDM and MOSE are shown in Fig. 2, where Fig. 2(a)-2(c) show the compared results for Model NO.1 and Fig. 2(d)-2(f) are those for Model NO.2 at the heights of *z* = 0, 3.5, 7 *mm*. In Fig. 2, the solid red curves represent the results of MOSE; while the blue asterisk curves and the green circle ones show those of the HSRM and HRDM respectively. From Fig. 2, we found that there was a better agreement between the results of the HSRM and MOSE than those between the HRDM and MOSE. For Model NO.1, *ARE* between the HSRM and MOSE was 1.803%, which was much smaller than that between the HRDM and MOSE which was 8.068%. For Model NO.2, *ARE* between the HSRM and MOSE was 1.743%, while the one between the HRDM and MOSE was 4.346%.

Except for the fact that a better agreement between the HSRM and MOSE was obtained, with *ARE* less than 2%, some other interesting conclusions could also be addressed. Firstly, the HSRM performed much better than the HRDM when the high absorption or low scattering regions existed in the solving subject. Secondly, the HSRM exhibited good performance in the non-scattering regions, such as the detection points at about 0-15 and 45-60 in Fig. 2(a)-2(c) and points at about 0-10 and 45-50 in Fig. 2(d)-2(e) where the results were affected greatly by the non-scattering regions. Last but not least, although the HSRM had good agreement with MOSE, both of the results of the HSRM and MOSE seemed a little fluctuant, as shown in Fig. 2. This may have been induced by the stochastic nature of the Monte Carlo method and insufficient mesh discretization of the HSRM. The comparison results demonstrated the accuracy of the HSRM model and illustrated its superiority over the HRDM.

### 3.1.2 3D digital mouse based simulation

For the application of HSRM in *in vivo* studies, a physical model should be first constructed as follows. Firstly, the tissue regions should be properly segmented from the anatomical structure of small animals which can be scanned by micro-CT or MRI. Secondly, the tissue regions were divided into high scattering, low scattering, and non-scattering regions according to their optical properties. Accurately, the optical properties should be ideally measured *in vivo*. Considering the difficulty in obtaining *in vivo* measurements of the optical properties, the optical properties commonly adopted the calculated values according to an empirical formula presented in the literature [32]. To mimic light propagation *in vivo* at the highest extent, a digital mouse model was employed to verify the accuracy of HSRM in an irregular model in this subsection. The torso of a digital mouse atlas extracted from the CT and cryosection data was used to construct the digital mouse model, and the main organs were selected to mimic the heterogeneities, as shown in Fig. 1(c) [37]. The detailed optical properties of each organ are listed in Table 2 at the wavelength around 690 *nm* [32]. The stomach was considered as the non-scattering region in this simulation. To simulate the light source which expressed the information of the lesion, a solid sphere source with the radius of 1.6 *mm* was located beside the stomach with the coordinates of (30.4, 15.2, 26.4) *mm*. To obtain relatively smooth results, the torso of the digital mouse was discretized into 94738 tetrahedral elements and 16998 nodes.

Distribution of the HSRM, HRDM and MOSE were obtained and shown in Fig. 3(a)-3(c), which had been normalized to their maximum intensities respectively. From Fig. 3(a)-3(c), we found that both the results of HSRM and HRDM approached those of MOSE well. However, the HSRM performed much better, while the result of HRDM appears to be more divergent. Especially, at the boundary and the innermost of the density spot in the distribution map of the photon density, the result of HSRM was more approximate to that of MOSE than the HRDM, which might be caused by the diffuse nature of the diffusion equation used in the HDRM. To further observe the differences of the results, two profiles were chosen at the heights of *z* = 26.5 and 34.5 *mm*, as shown in Fig. 3(d) and 3(e), and their relative errors are listed in Table 3. In Table 3, *ARE* is the average relative error between the calculated results of the HSRM or HRDM and the simulated one of MOSE on the whole detection points; *LARE* was defined to depict the average relative error between the calculated results and the one of MOSE on the partial detection points where the photon density was highly affected by the low-scattering regions, such as the detection points 15 to 25 in Fig. 3(d) and 15 to 20 in Fig. 3(e).

From Fig. 3(d)-3(e) and *ARE* listed in Table 3, it is obvious that both results of HSRM and HRDM were close to those of MOSE. However, the profiles at the detection positions of 15 to 25 in Fig. 3(d) and 15 to 20 in Fig. 3(e) show that the HSRM performed better than the HRDM as shown by *LARE* in Table 3. Because the influence of the low-scattering region mainly focuses on the detection points of 15 to 25 in Fig. 3(d) and 15 to 20 in Fig. 3(e), the HSRM and HRDM both performed well on all of the detection points. However, the HSRM performed better on the selected partial influenced detection points. On the other hand, there are some fluctuations in all of the results of the HSRM, HRDM and MOSE, which might be caused by the fact that the mesh was derived from the irregularity of the body surface of the digital mouse. The digital mouse based simulation demonstrated that the HSRM performed better than the HRDM when the low-scattering region existed simultaneously with the non-scattering regions and was more suitable for depicting light propagation in whole-body small animal imaging.

### 3.1.3 Physical cube phantom based experiment

To further verify the performance of the HSRM, a physical cubic phantom made of nylon was designed for the experiment. The dimensions of the cubic phantom were 30 *mm* in width and 30 *mm* in height. Two holes of different sizes were drilled in the phantom as shown in Fig. 1(d). The first hole with a radius of 1 *mm* and a depth of 16 *mm* was injected with the fluorescent solution of 2 *mm* in height to simulate the light source. The center of the light source was at the coordinate of (12, 12, 15) *mm*. The second hole was employed to simulate the non-scattering region, with a radius of 3 *mm* and a depth of 18 *mm*. The cavity region of 13 *mm* in height and 3 *mm* in radius was constructed by blocking up the top of the second hole with a nylon rod of 5 *mm* in height. The peak wavelength of the light source was about 660 *nm*, so the measured optical properties of the nylon phantom were listed as: the absorption coefficient was 0.023 *mm*^{−1} and the scattering coefficient was 20 *mm*^{−1} [7]. After the phantom was injected with the light source and mounted on the animal holder, two-dimensional luminescent images were acquired by our prototype imaging system of ZKKS-Direct 3D (jointly developed by Guangzhou Zhongke Kaisheng Medical Technology Co., Ltd. and Xidian University).

Because the two holes were drilled in diagonal lines, symmetry exists in the distribution of the photon density on the surfaces of the cube. Similarly, all of the results were normalized to each of their own maximum intensity and are shown in Fig. 4(a)-4(c). Therein, Fig. 4(a) and (b) show the two adjacent images on the phantom surface calculated results of the HSRM and HRDM respectively. Figure 4(c) shows the results mapped from the prototype imaging system. To further characterize the differences, the comparison results for the detection positions at *z* = 15 and 11 *mm* were extracted and shown in Fig. 4(d) and 4(e). From Fig. 4, we found that the similar tendency and good agreement were obtained between the measurements and the calculated results of both the HSRM and HRDM. The average relative errors between the measurements and the calculated results were about 1.824% and 1.818% respectively. Because of the high scattering of the nylon phantom, the performance of the HSRM and HRDM were both approved. Thus, the comparison results indicated that the proposed HSRM performed as well as the HRDM for simulating light transport in the medium with high scattering and non-scattering regions simultaneously, and showed good agreement with the measurement of the physical experiment.

#### 3.2 Investigation of the optical proprieties

In this subsection, a group of experiments was conducted to investigate the applicability of the HSRM by a series of optical properties. In the investigation, the solving domain was divided into the high scattering, low scattering and non-scattering regions based on the ratio of the reduced scattering coefficient to the absorption coefficient, according to our comprehension of the description from the literature [19, 23–25, 27]. If ${{\mu}^{\prime}}_{s}/{\mu}_{a}\ge 10$, the region was defined as the high scattering one, ${{\mu}^{\prime}}_{s}/{\mu}_{a}<10$ as the low scattering one and ${{\mu}^{\prime}}_{s}\approx 0$ as the non-scattering one. At the same time, the tissue was also divided to the high and low absorption region according to the threshold ${\mu}_{a}=0.1$. If the absorption was larger than the threshold, the tissue was subject to the high absorption region; otherwise, it was the low absorption region. According to the definition, the scattering tissues were divided into the high-scattering and high-absorption (HSHA), the high-scattering and low-absorption (HSLA), the low-scattering and high-absorption (LSHA) and the low-scattering and low-absorption (LSLA) groups.

To investigate the applicability of the HSRM model, four groups of experiments were conducted. The model was a simple cylinder, as shown in Fig. 5(a) and the related geometrical parameters are listed Table 4. In all of the four groups of experiments, the optical properties of the void region were the same, where the absorption ${\mu}_{a}$ was 0.25*mm*^{−1} and the scattering ${\mu}_{s}$ was nearly 0 *mm*^{−1}, and the optical parameters of the scattering region T1 was changed in the different groups of experiments. In the first group of experiments, three experiments were performed with the HSHA types of optical properties under the same ratio of ${{\mu}^{\prime}}_{s}/{\mu}_{a}\text{=}12$, in which the optical parameters $({{\mu}^{\prime}}_{s},{\mu}_{a})$ of the scattering region T1 were respectively set as the following: (1.2, 0.1), (2.4, 0.2), and (3.6, 0.3). In the second group, the three were performed with the HSLA types of optical properties under the ratio of ${{\mu}^{\prime}}_{s}/{\mu}_{a}\text{=30}$, in which the optical properties were respectively set as: (0.6, 0.02), (1.5, 0.05), and (2.4, 0.08). For the third group with the LSHA types of optical properties and the fourth one with the LSLA types of optical properties, the ratio ${{\mu}^{\prime}}_{s}/{\mu}_{a}$ of the two groups was set to be 3 and 6 respectively, and the corresponding optical parameters of the two groups were respectively as: (0.36, 0.12), (0.72, 0.24), and (1.08, 0.36) and (0.18, 0.03), (0.3, 0.05), and (0.42, 0.07). Figure 6 shows the comparison profiles of the four groups of experiments. Figure 6(a)-6(c) show the comparison profiles of the first three experiments with the HSHA optical property at the height of *z* = 0 *mm*; Fig. 6(d)-6(f) are the second three experiments with the HSLA optical property at the height of *z* = 0 *mm*; Fig. 6(g)-6(i) and Fig. 6(j)-6(l) are respectively the third three with the LSHA optical property and the fourth three with the LSLA optical property at the height of *z* = 0 *mm*. *ARE* of all of the experiments are listed in Table 5.

From Fig. 6 and Table 5, some conclusions can be addressed. Firstly, the HSRM performed better than the HRDM for all of the experiments with the series of different optical properties. It meant that the HSRM could improve the accuracy and obtain a more extensive application than the HRDM. Secondly, the HSRM had better improvement compared with the HRDM under the condition of both the LSHA and LSLA types of optical properties. It meant that the HSRM method was more suitable to depict light transport in the low scattering region. Thirdly, the HSRM and the HRDM both performed well with the HSLA optical property. However, to balance accuracy and effectiveness, the HRDM was more suitable to be employed in the case of the HSLA optical property. Lastly, although the novel HSRM model could improve the accuracy compared with the HRDM, it may be broken down in some cases with the HSHA type of optical properties, such as the case shown in Fig. 6(c). For Fig. 6(c), the transmission light may be affected greatly by high absorption and high scattering when the light source was deep inside the object. Simultaneously, some unsatisfactory performances of the HSRM also existed, such as in Fig. 6(b) and 6(j), which may have been caused by the lower order of SP_{N} used in the evaluation. The performance may be improved with an increase in the order of SP_{N}. Overall, the HSRM would achieve a much larger scope of application than the HRDM in whole-body small animal imaging.

#### 3.3 Investigation of the depth of the source in light transport

To investigate the influence of the depth of the light source in light transport, three experiments were performed. The models used in this subsection were the same as those used in Subsection 3.2, as shown in Fig. 5(b), where the localization of the light source is described. The optical parameters of the scattering region were set to be ${\mu}_{a}\text{=}0.3m{m}^{-1}$ and ${\mu}_{s}\text{=}36m{m}^{-1}$. By changing the localization of the light source, three experiments were obtained. To simulate the different depths of the light source, the center of the light source for the three experiments was set to be (3, 0, 0), (4.5, 0, 0) and (6, 0, 0) *mm*, in which the distance from the source center to the phantom surface was 7, 5.5, 4 *mm* respectively, as shown in Fig. 5(b). Comparison of the results of the HSRM, HRDM and MOSE for all three experiments are shown in Fig. 7 and the quantitative comparisons are listed in Table 6.

From Fig. 7 and Table 6, we observed that the depth of the source had a great influence on light transport under the high-scattering and high-absorption condition. Some phenomenon could also be found that the improvement of the HSRM appeared with a decrease in the depth of the light source. This may have been induced by the high scattering regions, in which the mean free path was very small. However, the HSRM also had a better performance than the HRDM. It meant that the HSRM had a more extensive application and was more accurate in describing the light transport compared with the HRDM.

## 4. Discussion and conclusion

Constructing an accurate and feasible light transport model is an essential task for the forward and inverse problems of three-dimensional (3D) optical imaging. Due to the complexity and diversity of the biological tissues, an accurate and efficient light transport model is urgently needed. For the problem of light transport in the media where the non-scattering region existed simultaneously with the low-scattering or high-absorption regions, a novel hybrid simplified spherical harmonics approximation (SP_{N}) and radiosity theory model (HSRM) was proposed in this paper.

The non-scattering, low-scattering and high absorption regions commonly exist in most species that can be simultaneously imaged optically. Rigorous simulations and experiments have demonstrated that the proposed HSRM model obtained good accuracy and exhibited great superiority over the HRDM in dealing with light transport in the media where the non-scattering region existed simultaneously with the low-scattering or high-absorption regions. Concretely, the HSRM would obtain good agreement with the results of MOSE when the ratio of the reduced scattering coefficient to the absorption coefficient was larger than 2 and the absorption coefficient was no larger than a certain value (1 in our studies). However, there are still some weaknesses for the proposed HSRM model. First, the HSRM is not all inclusive. It cannot achieve an approved result for all types of optical properties. Sometimes the HSRM would break down, especially for the case that the absorption coefficient is too big, which has been demonstrated in our simulations. Secondly, although the HSRM obtained high accuracy and much more extensive applications compared with the HRDM, the computational cost also increased. To further improve the accuracy and extend the application, we need to increase the order of the SP_{N} used in the HSRM, which will further increase the burden of the computation. Lastly, the result of the HSRM was also affected by the depth of the light source. The lower order HSRM model may not work with the increase in depth.

Overall, the proposed HSRM has a good potential in the applications of whole-body animal imaging and exhibits a more extensive application in modeling light propagation. Our future work will focus on the application of the HSRM model in small animal studies and the acceleration of the HSRM model. The corresponding results will be reported later.

## Acknowledgments

This work is supported in part by the Program of the National Basic Research and Development Program of China (973) under Grant No. 2011CB707702, the National Natural Science Foundation of China under Grant Nos. 81090272, 81227901, 81101083, 81201137, the Open Research Project under Grant 20120101 from SKLMCCS, and the Fundamental Research Funds for the Central Universities.

## References and links

**1. **V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. **23**(3), 313–320 (2005). [CrossRef] [PubMed]

**2. **J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. **27**(5), 48–57 (2008). [CrossRef] [PubMed]

**3. **J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. **7**(7), 591–607 (2008). [CrossRef] [PubMed]

**4. **A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. **50**(4), R1–R43 (2005). [CrossRef] [PubMed]

**5. **S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. **25**(12), 123010 (2009). [CrossRef]

**6. **A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. **202**(1), 323–345 (2005). [CrossRef]

**7. **W. X. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express **13**(18), 6756–6771 (2005). [CrossRef] [PubMed]

**8. **Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express **14**(18), 8211–8223 (2006). [CrossRef] [PubMed]

**9. **V. A. Markel and J. C. Schotland, “Inverse scattering for the diffusion equation with general boundary conditions,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **64**(3), 035601 (2001). [CrossRef] [PubMed]

**10. **K. Peng, X. Gao, X. Qu, N. Ren, X. Chen, X. He, X. Wang, J. Liang, and J. Tian, “Graphics processing unit parallel accelerated solution of the discrete ordinates for photon transport in biological tissues,” Appl. Opt. **50**(21), 3808–3823 (2011). [CrossRef] [PubMed]

**11. **Z. Yuan, X.-H. Hu, and H. Jiang, “A higher order diffusion model for three-dimensional photon migration and image reconstruction in optical tomography,” Phys. Med. Biol. **54**(1), 67–88 (2009). [CrossRef] [PubMed]

**12. **S. Wright, M. Schweiger, and S. R. Arridge, “Reconstruction in optical tomography using the PN approximations,” Meas. Sci. Technol. **18**(1), 247–260 (2007). [CrossRef]

**13. **P. Surya Mohan, T. Tarvainen, M. Schweiger, A. Pulkkinen, and S. R. Arridge, “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys. **230**(19), 7364–7383 (2011). [CrossRef]

**14. **W. Cong, A. Cong, H. Shen, Y. Liu, and G. Wang, “Flux vector formulation for photon propagation in the biological tissue,” Opt. Lett. **32**(19), 2837–2839 (2007). [CrossRef] [PubMed]

**15. **A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. **220**(1), 441–470 (2006). [CrossRef]

**16. **K. Liu, Y. Lu, J. Tian, C. Qin, X. Yang, S. Zhu, X. Yang, Q. Gao, and D. Han, “Evaluation of the simplified spherical harmonics approximation in bioluminescence tomography through heterogeneous mouse models,” Opt. Express **18**(20), 20988–21002 (2010). [CrossRef] [PubMed]

**17. **Y. Lu, A. Douraghy, H. B. Machado, D. Stout, J. Tian, H. Herschman, and A. F. Chatziioannou, “Spectrally resolved bioluminescence tomography with the third-order simplified spherical harmonics approximation,” Phys. Med. Biol. **54**(21), 6477–6493 (2009). [CrossRef] [PubMed]

**18. **A. D. Klose, “The forward and inverse problem in tissue optics based on the radiative transfer equation: A brief review,” J. Quant. Spectrosc. Radiat. Transf. **111**(11), 1852–1853 (2010). [CrossRef] [PubMed]

**19. **H. Dehghani, D. T. Delpy, and S. R. Arridge, “Photon migration in non-scattering tissue and the effects on image reconstruction,” Phys. Med. Biol. **44**(12), 2897–2906 (1999). [CrossRef] [PubMed]

**20. **H. Dehghani, S. R. Arridge, M. Schweiger, and D. T. Delpy, “Optical tomography in the presence of void regions,” J. Opt. Soc. Am. A **17**(9), 1659–1670 (2000). [CrossRef] [PubMed]

**21. **J. Riley, H. Dehghani, M. Schweiger, S. R. Arridge, J. Ripoll, and M. Nieto-Vesperinas, “3D optical tomography in the presence of void regions,” Opt. Express **7**(13), 462–467 (2000). [CrossRef] [PubMed]

**22. **D. Yang, X. Chen, S. Ren, X. Qu, J. Tian, and J. Liang, “Influence investigation of a void region on modeling light propagation in a heterogeneous medium,” Appl. Opt. **52**(3), 400–408 (2013). [CrossRef] [PubMed]

**23. **M. Firbank, S. R. Arridge, M. Schweiger, and D. T. Delpy, “An investigation of light transport through scattering bodies with non-scattering regions,” Phys. Med. Biol. **41**(4), 767–783 (1996). [CrossRef] [PubMed]

**24. **S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light in scattering media: A direct method for domains with nonscattering regions,” Med. Phys. **27**(1), 252–264 (2000). [CrossRef] [PubMed]

**25. **J. H. Lee, S. Kim, and Y. T. Kim, “Modeling of diffuse-diffuse photon coupling via a nonscattering region: a comparative study,” Appl. Opt. **43**(18), 3640–3655 (2004). [CrossRef] [PubMed]

**26. **X. Chen, D. Yang, X. Qu, H. Hu, J. Liang, X. Gao, and J. Tian, “Comparisons of hybrid radiosity-diffusion model and diffusion equation for bioluminescence tomography in cavity cancer detection,” J. Biomed. Opt. **17**(6), 066015 (2012). [CrossRef] [PubMed]

**27. **T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol. **50**(20), 4913–4930 (2005). [CrossRef] [PubMed]

**28. **T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. P. Kaipio, “Hybrid radiative-transfer-diffusion model for optical tomography,” Appl. Opt. **44**(6), 876–886 (2005). [CrossRef] [PubMed]

**29. **D. Gorpas and S. Andersson-Engels, “Evaluation of a radiative transfer equation and diffusion approximation hybrid forward solver for fluorescence molecular imaging,” J. Biomed. Opt. **17**(12), 126010 (2012). [CrossRef] [PubMed]

**30. **T. Hayashi, Y. Kashio, and E. Okada, “Hybrid Monte Carlo-diffusion method for light propagation in tissue with a low-scattering region,” Appl. Opt. **42**(16), 2888–2896 (2003). [CrossRef] [PubMed]

**31. **V. Y. Soloviev, G. Zacharakis, G. Spiliopoulos, R. Favicchio, T. Correia, S. R. Arridge, and J. Ripoll, “Tomographic imaging with polarized light,” J. Opt. Soc. Am. A **29**(6), 980–988 (2012). [CrossRef] [PubMed]

**32. **G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. **50**(17), 4225–4241 (2005). [CrossRef] [PubMed]

**33. **Y. Lu, H. B. Machado, A. Douraghy, D. Stout, H. Herschman, and A. F. Chatziioannou, “Experimental bioluminescence tomography with fully parallel radiative-transfer-based reconstruction framework,” Opt. Express **17**(19), 16681–16695 (2009). [CrossRef] [PubMed]

**34. **J. Ripoll, R. B. Schulz, and V. Ntziachristos, “Free-space propagation of diffuse light: theory and experiments,” Phys. Rev. Lett. **91**(10), 103901 (2003). [CrossRef] [PubMed]

**35. **H. Li, J. Tian, F. P. Zhu, W. X. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method,” Acad. Radiol. **11**(9), 1029–1038 (2004). [CrossRef] [PubMed]

**36. **N. Ren, J. Liang, X. Qu, J. Li, B. Lu, and J. Tian, “GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues,” Opt. Express **18**(7), 6811–6823 (2010). [CrossRef] [PubMed]

**37. **B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. **52**(3), 577–587 (2007). [CrossRef] [PubMed]