A study of photon propagation in free-space based on hybrid radiosity-radiance theorem

Noncontact optical imaging has attracted increasing attention in recent years due to its significant advantages on detection sensitivity, spatial resolution, image quality and system simplicity compared with contact measurement. However, photon transport simulation in free-space is still an extremely challenging topic for the complexity of the optical system. For this purpose, this paper proposes an analytical model for photon propagation in free-space based on hybrid radiosity-radiance theorem (HRRT). It combines Lambert’s cosine law and the radiance theorem to handle the influence of the complicated lens and to simplify the photon transport process in the optical system. The performance of the proposed model is evaluated and validated with numerical simulations and physical experiments. Qualitative comparison results of flux distribution at the detector are presented. In particular, error analysis demonstrates the feasibility and potential of the proposed model for simulating photon propagation in free-space. 2009 Optical Society of America OCIS codes: (110.2970) Image detection systems; (170.3880) Medical and biological imaging; (170.5280) Photon migration; (170.7050) Turbid media. References and links 1. R. Weissleder, and U. Mahmood, “Molecular imaging,” Radiology 219, 316-333 (2001). 2. B. W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt. 6, 432440 (2001). 3. R. Weissleder, and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9, 123-128 (2003). 4. 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, 313-320 (2005). 5. 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, 1659-1670 (2000). 6. J. Tian, J. Bai, X. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27, 48-57 (2008). 7. R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. 28, 1701-1703 (2003). 8. R. B. Schulz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging 23, 492-500 (2004). 9. N. Deliolanis, T. Lasser, D. Hyde, A. Soubert, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett. 32, 382-384 (2007). 10. V. Ntziachristos, E. E. Graves, R. B. Schultz, and J. Ripoll, “Fluorescence molecular tomography: New detection schemes for acquiring high information content measurements,” IEEE International Symposium on Biomedical Imaging (ISBI 2004) 2, 1475-1478 (2004). 11. R. B. Schultz, J. Peter, and W. Semmler, “Comparison of noncontact and fiber-based fluorescence-mediated tomography,” Opt. Lett. 31, 769-771 (2006). (C) 2009 OSA 31 August 2009 / Vol. 17, No. 18 / OPTICS EXPRESS 16266 #111632 $15.00 USD Received 20 May 2009; revised 6 Aug 2009; accepted 23 Aug 2009; published 28 Aug 2009 12. H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V. Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic surface extraction,” Appl. Opt. 46, 3617-3627 (2007). 13. W. 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, 6756-6771 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-18-6756. 14. J. H. Lee, S. Kim, and Y. T. Kim, “Finite element method for diffusive light propagation in index-mismatched media,” Opt. Express 12, 1727-1740 (2004), http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-8-1727. 15. Y. Lv, J. Tian, H. Li, J. Luo, W. Cong, G. Wang, and D. Kumar, “Modeling the forward problem based on the adaptive FEMs framework in bioluminescence tomography,” Proc. SPIE 6318, 631801 (2006). 16. C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16, 20317-20333 (2008), http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-25-20317. 17. L. V. Wang, S. L. Jacques, L. Q. Zheng, “MCML-Monte Carlo modeling of photon transport in multi-layered tissues,” Comput. Met. Prog. Biomed. 47, 131-146 (1995). 18. D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159-169 (2002), http://www.opticsinfobase.org/abstract.cfm?URI=oe-10-3-159. 19. E. Margallo-Balbas, and P. J. French, “Shape based Monte Carlo code for light transport in complex heterogeneous tissues,” Opt. Express 15, 14086-14098 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-21-14086. 20. J. Ripoll, R. B. Schulz, and V. Ntziachristos, “Free-space propagation of diffuse light: theory and experiments,” Phys. Rev. Lett. 91, 103901-1-4 (2003). 21. J. Ripoll, and V. Ntziachristos, “Imaging scattering media from a distance: theory and applications of noncontact optical tomography,” Mod. Phys. Lett. B 18, 1-29 (2004). 22. R. B. Schulz, J. Peter, W. Semmler, C. D’Andrea, G. Valentini, and R. Cubeddu, “Quantifiability and image quality in noncontact fluorescence tomography,” Proc. SPIE 5859, 58590Z-1-8 (2005). 23. J. Yao, G. Hu, and J. Bai, “Modeling and validation of light propagation in free-space for noncontact nearinfrared fluorescent tomography,” J. Infrared Millim. W. 27, 330-332,369 (2008). 24. Q. Tian, Y. Liao, and L.Sun, Engineering optics (Tsinghua University Press, Beijing, 2004). 25. J. Zhang, X. Fang, Infrared Physics (Xidian University Press, Xi’an, 2004). 26. 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, 767-783 (1996). 27. J. Piley, 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, 462-467 (2000), http://www.opticsinfobase.org/abstract.cfm?URI=oe-7-13-462. 28. M. Aggarwal, and N. Ahuja, “A pupil-centric model of image formation,” Int. J. Comput. Vis. 48, 195-214 (2002). 29. M. Aggarwal, and N. Ahuja, “A new imaging model,” in Proceeding of IEEE Conference on Computer Vision (Institute of Electrical and Electronics Engineers, New York, 2001), pp. 82-89. 30. H. Li, J. Tian, F. Zhu, W. X. Cong, L. V. Wang, 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, 1029-1038 (2004). 31. D. Qin, H. Zhao, Y. Tanikawa, and F. Gao, “Experimental determination of optical properties in turbid medium by TCSPC technique,” Proc. SPIE 6434, 64342E (2007).


Introduction
With the development of the optical marker technique, optical molecular imaging technology, especially optical tomography, has been more and more widely used in many domains due to its reasonable spatial and temporal resolution and affordable cost [1,2,3,4].Optical tomography is a noninvasive imaging technology which aims to image the optical properties or abnormities of biological tissue from measurement of the transmitted light on the surface of the body [5].Among these various optical tomography modalities, diffuse optical tomography (DOT), fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT) are three main categories and have become research focuses over the past years for their excellent performance [6].Noncontact measurement could facilitate the use of large detector arrays at complete projection angles that are well-suited for tomography application.Furthermore, such an imaging strategy eliminates the use of restricted geometries, the need for individual fibers and matching fluids that are in contact with the imaging object [7,8,9], which are supposed to increase detection sensitivity, improve the spatial resolution and avoid the attenuation within matching fluids [10,11,12].In biomedical research, noncontact optical tomography has become a valuable tool for the noninvasive detection area.
In noncontact optical tomography, the external measured photons that escape from the turbid media surface are all detected with a high performance CCD camera.Thus, photon propagation both in turbid media and free-space must be modeled for such a noncontact optical tomography.Many numerical algorithms [13,14,15,16] have been proposed to simulate photon propagation in turbid media and performed high efficiency and good applicability, but they do not work when photons escape from the media surface and propagate in free-space.Monte Carlo method, which is a gold-standard method to dynamically trace various behaviors of photon transporting process [17,18,19], can handle the photon propagation problem both in media and free-space.However, it is time-consuming and needs to run a large number of photons to obtain a practicable result.An effective freespace photon transport model has been successfully proposed by J. Ripoll's group, who firstly presents and demonstrates the possibility of realizing qualitative noncontact optical tomography [7,20,21].However, the application of this model is limited in real environment for lacking of the consideration over perspective and depth-of-field effects of the objective in a complicated optical system [22,23].Based on the proposed model, R.B. Schulz [22] and J. Bai [23] put forward two new improved and simplified models respectively, using the perspective projection method through replacing the objective with a virtual pinhole.Nevertheless, the improved models neglect the depth-of-field effects and only work in the small aperture objective environment.Therefore, for the purpose of modeling this noncontact scheme, to develop a general and effective model for free-space photon propagation in a complicated optical system is becoming an urgent problem to be solved.
Based on the consideration of aforementioned problems, a hybrid radiosity-radiance theorem (HRRT) model is proposed and developed in this paper to simulate the photon transport process in free-space.Compared with available models, HRRT model takes a complicated optical system into account, including optical layout and objective effects analysis such as perspective effects, image aberrations and depth-of-field effects.In the HRRT model, Lambert's cosine law (LCL) that derived from the radiosity theory on the basis of Lambertian source theory combined with the radiance theorem, plays the leading role in calculating photon flux on CCD camera.Therein the radiance theorem embraces the analysis of perspective effects.On the other hand, the introduction of the virtual detector plane could contribute to eliminate the influence of image aberrations and depth-of-field effects of lens and establish a specific relationship between the media surface and CCD camera.The paper is organized as follows.Section 2 depicts the proposed HRRT model.In Section 3, the performance of the model is validated with numerical simulation and physical experiments involving difform phantom surface.In the end, relevant issues are discussed and conclusion is provided.

The Lambert's cosine law
Once the turbid media surface flux density has been obtained, the free-space propagation of photon flux radiated from media surface S to receiving plane A can be modeled.Assuming that each differential surface dS centered at a point r constitutes a new Lambertian source [20,21], it would act as a radiation source and irradiate to the surrounding space isotropically, that is, the radiance of the source is constant and independent of the solid angle but varies in different positions [24,25].Thus, the micro-unit power of a differential receiving area dA centered at d r can be obtained through the following equation according to the literature [20,24,25]: where ( ) d r,r ξ is a visibility factor that discards the surface points invisible from the receiving plane; cos s θ is defined as a cosine value of the angle between the differential surface normal and the directional vector from r to d r ; cos d θ accounts for the detector orientation with respect to the line of sight and can be equally defined as a cosine value of the angle between the differential detector normal and the directional vector from d r to r ; dS and dA are the area of the differential surface and receiving unit, respectively.
Eq. ( 1) is distinctively named LCL in this paper, which presents a formula derived from radiosity theory on the basis of the Lambertian source theory [26,27].Thus, as a principal part of the proposed HRRT model, LCL acts as a main theoretical basis to depict the energy transformation of diffuse photons when they are propagating in free-space.

The radiance theorem
The radiance theorem is an important component of the proposed HRRT model, which demonstrates that the radiance is conserved through a loss-less optical system.According to the Conservation Law of Basal Radiance and literature [25], the radiance of image can be calculated via that of object in an optical system.An ideal optical system, which neglects the loss caused by the reflection and absorption of lens system, is shown in Fig. 1.Under the assumptions that AB is an ideal lens system, dS is a differential area at the object plane perpendicular to the optical axis and dS ′ is that at the image plane, dΩ and d ′ Ω are solid angles of differential object and image area with respect to the ideal lens, r is the distance from the differential object area to the center of lens system and r′ is that of the differential image area, θ accounts for the angle of the line of sight with respect to the optical axis, the relationship between radiance of object L and that of image L′ can be obtained: 2) is a general description of the radiance theorem.And it is shown that the radiance of image is equal to that of object when an ideal optical system is taken into account.Since the irradiance radiated from a differential area dS in direction s within a solid angle dΩ is given by cos E L d θ = Ω .Thus, irradiance E at the lens radiated from the differential object area can be calculated using radiance L on the basis of the definition of solid angle where A is the effective area of the lens.Similarly, irradiance from the image area satisfies: where E′ is the irradiance at the lens radiated from the differential image area.Substituting Eqs.
(2) and (3) into Eq.( 4), the following relationship of irradiance between object and image differential area can be obtained: 2 where t is magnification ratio of the lens system given by / t u ν = ; u denotes the distance along the optical axis between a object plane and the lens system, ν denotes the perpendicular distance between the focused image and the lens system.Using the knowledge of geometric optics and magnification ratio, some useful equations can be derived and one of them has played an important part in the derivation of Eq. ( 5):   dS t dS ′ = (7) Based on the definition of radiant flux, we can obtain the radiant power of differential object area that irradiates to the lens system as: P EdS = (8) In the same way, the power received from the corresponding image area is: 5), ( 7) and ( 8) into Eq.( 9), we can obtain an expression as follows: P P ′ = (10) Equation ( 10) is a derived formula of the radiance theorem (DRT), which depicts the energy transformation when light passes through a lens system.The derivation of radiance theorem reveals the law of conservation of energy when photon propagates in free-space involving an ideal optical system.

Hybrid radiosity-radiance theorem model and numerical implementation
An imaging objective in an ideal optical system typically is not a single-lens, but a lens combination with finite thickness, and Gauss had postulated that this lens combination can be modeled as a Gaussian thick-lens model [28].Furthermore, the lens system in this model can be replaced by two equivalent refracting surfaces called principal planes.The detector position v is measured with respect to one of the principal planes and the object distance u is to the other [28,29], as shown in Fig. 2(a).In addition, they satisfy the lens law given by 1/ 1/ 1/ u f ν + = .On the basis of geometric and photometric mapping, Gaussian thick-lens model can be equivalently reduced to a thin lens one located at the first principal plane and the distance from the thin-lens to the detector changes to ν ′ [29], as shown in Fig. 2(b).Thus, the simplified detector position ν ′ satisfies the following equation: where h and f are the total thickness and the focus of the lens system, respectively.Based on the aforementioned simplification theory of an imaging objective, the highperformance CCD camera adopted in optical imaging experiments can be equivalently modeled with a simple optical system consisting of only two parts: a simplified thin-lens and a detector plane.The experimental configuration for the simplified noncontact measurement in real experiments is shown in Fig. 3. Here, a concept of virtual detector plane, whose position and size are only determined by the position information of the optical system and the size of the lens system, is introduced to be used as a transit station for combining LCL with DRT.This means that the photon flux is firstly registered at the virtual detector plane using LCL, which depicts the energy changes when photon transfers in free-space, and then mapped to the detector plane on the basis of DRT.Furthermore, the virtual detector plane and the detector plane are said to be conjugated in geometric optics, that is, they satisfy the lens law mentioned above.In other words, any point at the virtual detector plane must be exactly imaged at the detector plane.For the purpose of pursuing photon flux at the detector plane, we will rewrite Eq. ( 1) to obtain the value of its conjugated position at the virtual detector plane firstly: where the meanings of the parameters are the same as that in Eq. ( 1) except that the points at the detector plane are replaced by their conjugated points at the virtual detector plane.Using the geometric knowledge, the point d r at the detector plane can be determined by the following equation: where d is the distance between the point d r at the detector plane and its conjugated point vd r at the virtual detector plane, and can be determined by ; s is a unit direction vector along the line of sight (see Fig. 3).Eq. ( 12 where u′ is the object distance of the simplified optical system and can be calculated using the lens law when the imaging distance is determined; θ is the angle between the line of sight and the optical axis, as shown in Fig. 3.According to the DRT derived above, photon flux of all points at the detector plane is equal to that of the corresponding conjugated points at the virtual detector plane.Therefore, a generalized formula can be obtained using Eqs.( 10), ( 12) and ( 14) synchronously: where cos d θ has the same value as cos vd θ and is a visibility factor that discards the surface points invisible from the lens and depends mainly on the parameters of the lens system configured in the optical system.Important to recall that the visibility factor in Eq. ( 15) is intrinsically different from the one presented in Eq. (1) that does not take the lens factor into account and only applies to the case of no lens system.Furthermore, Eq. ( 15) can be conveniently rewritten to yield the following equation: Here, a function that accounts for photon propagation in free-space from the media surface to the detector plane consisting of a lens system has been introduced: Equations.( 16) and ( 17) are a generalized description of HRRT model for photon propagation in free-space.Compared with the published method in literature [20], the proposed HRRT model avoids the use of the NA (numerical aperture) dependent function which is very difficult to determine.Furthermore, the HRRT model improves the visibility factor which has taken the influence of the lens system into account and also decreases the influence of perspective effects via the introduction of the radiance theorem.Thus, the flux distribution at the detector plane in a complicated optical system can be obtained directly by the proposed HRRT model.
To sum up, the implementation framework for the aforementioned HRRT algorithm is shown as follows: Algorithm HRRT algorithm for modeling photon propagation in free-space Require: Obtain the information of optical system and surface flux distribution.
1: Compute the central position and the size of virtual detector plane using Eq. ( 14).

Experiments and results
To validate the HRRT model developed in this paper, several verification experiments were conducted in this section, including numerical simulation experiment based on Monte Carlo method and physical experiment using a lens-coupled CCD system.Important to recall is that all the media surface flux utilized in the HRRT model was obtained by Molecular Optical Simulation Environment, which is a Monte Carlo method based simulation platform [30] and can be accessed freely on the web site: http://www.mosetm.net/.For the sake of valid analysis, terms of mean error (ME) and root-mean-square error (RMSE) are introduced to scale their discrepancy in this section.In mathematics, the quanta of ME and RMSE can be defined as ∑ respectively, where i e is the ith error between the calculated result by HRRT and experimental one.And the RMSE can be better to characterize the randomicity of Monte Carlo method.

Numerical simulation experiments
Two simulation experiments involving difform phantom surface were designed and conducted to evaluate the performance of the proposed HRRT model.The phantoms adopted in the simulations were designed to be homogeneous and assigned the same optical properties: absorption coefficient .In these experiments, a commercial software named TracePro (Lambda Research Corporation, Littleton, MA), which is based on the Monte Carlo method to simulate the light transport processes under various conditions, is used to obtain the simulation experimental results.And its simulation results will be compared with the calculated results using the proposed HRRT model.Moreover, it is important to note that the detector plane overlaps with the virtual detector plane, that is, the term f approaches infinity and t equals to 1 in Eq. ( 17).

Flat surface simulation
The experimental configuration for the flat surface simulation is shown in Fig. 4(a).Firstly, a cylindrical phantom was adopted in the simulation, which was located at ( ) 0, 0, 0 and of 10 mm height, 31.9 mm diameter.Then, a cylindrical source of 0.01 mm height and 0.05 mm diameter, which approximates to a point source with a total power of 0.491 Watts , was placed at the center of the bottom plane of the phantom.In order to register the photons escaped from the phantom surface, a rectangle detector of 31.9 mm length and 2 mm width was placed at ( ) 0, 0, 5.2 , which was parallel to the top plane with a perpendicular distance of 0.2 mm .Photon flux distribution at the detector obtained by HRRT and TracePro TM are shown in Fig. 4(b) and 4(c), respectively.
Results for detector y position: 0.0 was examined and the flux normalized to its maximum value versus the detector x position is shown in Fig. 5(a).The plus lines represent the simulation experimental results and the asterisks show the calculations, and their ME and RMSE are about 5.61% and 7.69%.Better to evaluate the comparison between HRRT and TracePro TM , we used the mean filtering method to achieve smooth curves with the ME and RMSE being about 4.18% and 3.93%, as shown in Fig. 5(b).As it is expected, the calculated results based on the HRRT model are in good agreement with the experimental ones obtained by TracePro TM .

Curved surface simulation
To further test the performance of the HRRT model, a curved phantom surface was adopted in the second simulation experiment.In the simulation, we placed two sources in a cylindrical phantom and a detector parallel to the generatrix of the cylinder, as shown in Fig. 6     In these two simulation experiments, similar tendency and good agreement between the calculated results and the simulation experimental data were obtained, as seen from Fig. 4-7.Moreover, the experimental data obtained from TracePro TM were sparse for the statistical characteristic of the Monte Carlo method.However, the calculated results using HRRT model seem much smoother because the photons escaped from the phantom surface would act as a new Lambertian source.The comparison results indicated that the proposed HRRT model can be comparable to the Monte Carlo simulation.Furthermore, it can obtain much smoother photon flux distribution and avoid high time complexity.

Physical experiments
Besides the above numerical simulation experiments, two physical experiments were also carried out for demonstrating the accuracy of the proposed algorithm in this section.Two different shape physical phantoms, including a cubic phantom and a cylindrical phantom, were designed and adopted.As mentioned above, the ME and RMSE between the calculated results and experimental ones were used to scale the performance of the proposed HRRT algorithm.

Cubic phantom experiment
In the optical imaging experiment, a cubic phantom of 45 mm side length was designed and utilized.The phantom is made from nylon, and one small hole of 5.9 mm diameter and 23.5 mm depth from the top surface was drilled in the phantom to emplace the light source, as shown in Fig. 8. Compound solutions from a red luminescent light stick (Glowproducts, Canada) were injected into the phantom to serve as the internal light source.In luminescent light stick, the peroxide and ester compound solutions provide ATP and other co-factors for fluorescent dye to generate luminescence light whose central wavelength located about 650 nm [16].After injecting the light source, the free volume of the small hole was filled with the same material solid stick to avoid leakage of the light and then the cubic phantom was mounted on a 360° rotation stage.By rotating the phantom three 90°, luminescence light was registered by a scientific liquid-water-cooled back-illuminated CCD camera (Princeton Instruments/Acton 2048B, Roper Scientific, Trenton, NJ) which was coupled with an optical lens subsystem (Nikon Nikkor Micro 55 mm , f/2.8, Japan).For each side measurement, the experiments were performed in a totally dark environment and a 2048 2048 × pixels and 16 bits dynamic range image with 13. 5 13.5 m m µ µ × sized pixel was recorded with an exposure time of 120 s .In this section, four captured images were used to reveal the validity of the proposed HRRT model, as shown in Fig. 9(a)-9(d).As seen from the figures, a small inconsecutive region appeared at the bottom of the image area, which is induced by the obstruction of phantom holder.
The calculation setup for HRRT model is shown in Fig. 8(c), and Fig. 8(b) gives the schematic diagram of four different perspectives.The same length cube, consisting of 274626 nodes and 49152 surface elements, was adopted in the numerical calculation with a cylindrical light source located at its center.The optical parameters used in the calculation were measured at the wavelength about 660 nm by a time-correlated single photon counting system [31] and listed as follows: the absorption coefficient   Results of HRRT for four different perspectives are shown in Fig. 9(e)-9(h).Similar tendency and good agreement were observed in the cases examined.Particularly, Fig. 10(a)-10(d) depict the comparative curves of the results calculated by the proposed model and measured using the CCD camera.Although there were some little differences in 90° and 270° directions, which were intrinsically caused by the insufficient smoothness of experimental phantom surface, their ME and RMSE had still shown the conformability from the aspect of quantitative error analysis, as listed in Table 1.As it is expected, the HRRT worked well in simulating the photon propagation in free-space in complicated optical imaging system.

Cylindrical phantom experiment
In the cylindrical phantom experiment, a homogeneous phantom of 30 mm height and 15 mm radius was utilized to perform the physical experiment; the phantom is also made from nylon and has the same optical properties, and a luminescence light source of 2 mm height and 1 mm radius was embedded into the phantom with its center at ( ) 8, 0, 0 , as shown in Fig. 11  (a) and (b).In addition, the same optical system as described in the cubic phantom experiment was adopted to register luminescence photons.
We further illustrated the ability of our proposed HRRT model when the lens system and a curved phantom surface were taken into account.Similarly, results for detector z positions: 0.0    and the asterisks show the theoretical calculations.Moreover, another group of calculated results obtained by the published method in the literature [20] were also compared with the experimental results, as shown in Fig. 12(b).The ME and RMSE between the calculated results and experimental data are listed in Table 2. Good agreement between the theoretical calculations, which are obtained by both the HRRT model and the published method, and the experimental data measured using a CCD camera was observed.The comparison results enhanced that the HRRT model developed in this paper can be comparable to the published one.

Discussion and conclusion
Noncontact optical tomography technique has been rapidly developed and broadly used in various biomedical domains due to its unique advantages on resolution, specificity, directness, non-injury and acceptability.Such an imaging strategy could eliminate the need for individual fibers that are in contact with tissue, restricted geometries and matching fluids.Thus, it could significantly improve the experimental procedures and realize to visualize physiological and pathological processes in vivo at the cellular and molecular levels.Therefore, the study of photon propagation in free-space is necessary and urgent because it helps not only to lay foundation for the noncontact measurement theory but also to model the experimental configuration effectively in real environment.However, an effective model to study the freespace photon propagation has not been developed and still been a challenging and urgent issue to make a breakthrough; especially a complicate optical system is to be considered, such as complicated lens subsystem.In this paper, a novel hybrid radiosity-radiance theorem model has been presented to solve photon transport problem in free-space.In comparison with the TracePro TM , the published model and the physical phantom experiments, HRRT model has the following novel potentials.Firstly, some effective models and software are available, but these strategies are not perfect enough to solve the photon propagation problem in free-space.TracePro TM is based on the Monte Carlo method to simulate the photon propagation processes and can be considered as a standard to evaluate other analytic models.However, it is time-consuming and needs to simulate a large number of photons to obtain a relatively smooth result.The analytic model published in literature [20] has high efficiency but neglects the influence of perspective and depth-of-field effects.The proposed HRRT model in the paper can be comparable to TracePro TM and the published model in some degree and takes the influence of lens into account, such as the perspective and depth-of-field effects and image aberrations.Furthermore, the proposed model provides a probability to conquer the obstacle of quantitative analysis due to its consideration of the complicacy of optical system.Based on a precise result of qualitative analysis, quantitative results could be obtained through taking the transmittance factor into account.
Secondly, when the optical imaging experiment is conducted, not all the light escaped from phantom surface can reach the whole entrance pupil of the optical lens subsystem.If an incident ray bundle is partially blocked by the phantom surface profile, vignetting effect in optics appears.In the proposed HRRT model, the introduction of virtual detector plane and a discretizing method for it can provide a practical method to handle this vignetting phenomenon.Furthermore, it is also available to solve the depth-of-field problem brought by the optical lens subsystem for its specific conjugated relationship established with the detector plane.In addition, the new calculation method of the visibility factor can effectively reduce the image aberration effects caused by optical lens system.Thus, this method could model the photon propagation in free-space more factually and effectively.
Finally, the developed model in this paper has provided a way to realize the reconstruction of three-dimensional surface energy distribution from two-dimensional gray scale images obtained using a CCD camera.It is well known that the inverse problem of optical tomography is to reconstruct the abnormities or optical properties of biological tissue.However, all the reconstruction procedures usually begin with the precondition that the surface flux distribution of the tissue-like models is given, and then reconstruct the targets using the related reconstructed algorithms.Therefore, the reconstruction of three-dimensional surface energy distribution from two-dimensional gray scale images must be modeled so as to obtain more perfect reconstructed results.However, it is still an awkward problem to achieve this distinctive mapping relationship due to its energy mappings rather than three-dimensional shape reconstruction.Fortunately, the proposed model in this paper would provide theoretical grounds for achieving the knotty mapping through the inverse procedure of the HRRT model.And the corresponding work is undergoing.
In summary, we have presented a novel model to simulate photon propagation in freespace for complicated optical system.Furthermore, the correctness and effectiveness of the proposed algorithm have been demonstrated with the numerical simulation and physical experiments.The preliminary results show that the HRRT model is of great potential in the noncontact optical tomography technology.Our further study will concentrate on the quantitative analysis of the photon transport in free-space and the reconstruction of threedimensional surface energy distribution from two-dimensional gray scale images obtained using a CCD camera.

Fig. 2 .
Fig. 2. Schematic diagram for simplification of the Gaussian thick-lens model to a thin-lens model.(a) Gaussian thick-lens model; (b) Thin-lens model.
) illustrates that the point d r can be determined by the point vd r moving a distance d in the direction s .Using the formula of the lens law given by 1obtain the relationship between the point d r and its conjugated point vd r as follows:

3 :
Discretize the phantom surface and obtain the number of differential surface unit m .4: for i = 1 to n do the discretized form of Eqs. (
(a).Photon flux distribution at the detector obtained by HRRT and TracePro TM are given in Fig. 6(b) and 6(c), respectively.Furthermore, results for detector z position: 5.2 d z mm = was examined and Fig. 7(a) plots the flux normalized to its maximum value versus the detector y position.Their ME and RMSE are about 4.28% and 6.14%, which were lower than the results in the first simulation.The filtered curves were also shown in Fig. 7(b), with a much lower value of the ME and RMSE being about 1.45% and 1.70%.

Fig. 4 .
Fig. 4. Numerical simulation experiment for the flat surface.(a) Phantom configuration; (b) and (c) Photon flux density at the detector obtained by HRRT and TracePro TM .

Fig. 5 .
Fig. 5. Comparison between the calculated results using HRRT model and the experimental data obtained from TracePro TM at the position yd = 0.0mm at the detector.(a) Original curves; (b) Filtered curves based on (a).

Fig. 6 .
Fig. 6.Numerical simulation experiment for the curved surface.(a) Phantom configuration; (b) and (c) Photon flux density at the detector obtained by HRRT and TracePro TM .

Fig. 7 .
Fig. 7. Comparison between the calculated results using HRRT model and the experimental data obtained from TracePro TM at the position zd = 5.2mm at the detector.(a) Original curves; (b) Filtered curves based on (a).

Fig. 8 .
Fig. 8. Physical experiment for cubic phantom.(a) Physical phantom with one light source; (b) Schematic diagram of four different perspectives; (c) Schematic diagram of numerical calculation phantom.
4.0 mm are shown in Fig. 12. Fig. 12(a) describes the results of HRRT model and the cylindrical phantom experiment.The solid lines represent the experimental data

Fig. 11 .
Fig. 11.Physical experiment for cylinder phantom.(a) Physical phantom with one light source; (b) Schematic diagram of numerical calculation phantom.

Fig. 12 .
Fig. 12.Comparison between the calculated and experimental results at the detector positions: zd = 0.0, 2.0, 4.0mm.(a) Comparison between HRRT and experimental results; (b) Comparison between the published method and experimental results

Table 1 .
Error comparisons between the calculated and experimental results of four perspectives

Table 2 .
Error comparisons between calculated results and experimental data