Light transport in turbid media with non-scattering, low-scattering and high absorption heterogeneities based on hybrid simplified spherical harmonics with radiosity model.

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 (SPN) 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 SPN 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 SPN 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.


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][2][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][5][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][8][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][16][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][16][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][20][21][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][20][21][22][23][24][25][26], the hybrid RTE-diffusion method [27][28][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][20][21][22][23][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 690nm) [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-, lowand 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.

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]: where [ For the void problem, the radiosity theory was used to model light propagation in the nonscattering 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 + ′ at r′ , if r and r′ 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 where a μ is the absorption coefficient of a non-scattering region, B is the interface between the non-scattering and scattering regions, r and r′ are points on interface B , θ and θ ′ are the angles between the surface normal at r and r′ , and the unit vector where j β varies with the order of 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: ψ is the piecewise linear nodal basis function. Then, the corresponding variational forms of the HSRM are derived as the following expression: The M , Ξ and F could be expressed as:

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 threedimensional 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: where f i is the value of the intensity at the ith 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.

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.  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 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  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 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 nonscattering regions and was more suitable for depicting light propagation in whole-body small animal imaging.

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, twodimensional 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.

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

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

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