A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm

The lack of a primary method for determination of optical parameters remains a significant barrier in optical study of turbid media. We present a complete system of experimental setups and Monte Carlo modeling tools for fast and accurate solution of the inverse problem from the measured signals of homogeneous turbid samples. The calibration of the instrument and validation of the Monte Carlo modeling have been carried out to ensure the accuracy of the inverse solution. We applied this method to determine the optical parameters of turbid media of 10% intralipid between 550 and 940nm and 20% intralipid between 550 and 1630nm. ©2006 Optical Society of America OCIS codes: 290.7050 Turbid media; 170.3660 Light propagation in tissues. _____________________________________________________________________________________________ References and Links 1. H. C. van de Hulst, Multiple light scattering: tables, formulas, and applications, vol. 1 & 2. (Academic Press, New York, 1980). 2. M. Born, E. Wolf,A. B. Bhatia, Principles of optics: electromagnetic theory of propagation, interference and diffraction of light, 7th ed. (Cambridge University Press, Cambridge, England, 1999). 3. H. Ding, J. Q. Lu, K. M. Jacobs,X. H. Hu, "Determination of refractive indices of porcine skin tissues and intralipid at eight wavelengths between 325 and 1557nm," J. Opt. Soc. Am. A 22, 1151-1157 (2005). 4. L. G. Henyey,J. L. Greenstein, "Diffuse radiation in the galaxy," Astrophys J 93, 70-83 (1941). 5. S. L. Jacques, C. A. Alter,S. A. Prahl, "Angular dependence of HeNe laser light scattering by human dermis," Lasers in Life Sciences 1, 309-333 (1987). 6. V. G. Peters, D. R. Wyman, M. S. Patterson,G. L. Frank, "Optical properties of normal and diseased human breast tissues in the visible and near infrared," Phys. Med. Biol. 35, 1317-34 (1990). 7. S. A. Prahl, M. J. C. van Gemert,A. J. Welch, "Determining the optical properties of turbid media by using the adding-doubling method," Appl. Opt. 32, 559-568 (1993). 8. I. V. Yaroslavsky, A. N. Yaroslavsky, T. Goldbach,H.-J. Schwarzmaier, "Inverse hybrid technique for determining the optical properties of turbid media from integrating-sphere measurements," Appl. Opt. 35, 67976809 (1996). 9. C. K. Hayakawa, B. Y. Hill, J. S. You, F. Bevilacqua, J. Spanier,V. Venugopalan, "Use of the delta-P1 approximation for recovery of optical absorption, scattering, and asymmetry coefficients in turbid media," Appl Opt 43, 4677-84 (2004). 10. A. H. Hielscher, J. R. Mourant,I. J. Bigio, "Influence of particle size and concentration on the diffuse backscattering of polarized light from tissue phantoms and biological cell suspensions," Appl. Opt. 36, 125-136 (1997). #70524 $15.00 USD Received 1 May 2006; revised 28 June 2006; accepted 3 July 2006 (C) 2006 OSA 7 August 2006 / Vol. 14, No. 16 / OPTICS EXPRESS 7420 11. Y. Du, X. H. Hu, M. Cariveau, X. Ma, G. W. Kalmus,J. Q. Lu, "Optical properties of porcine skin dermis between 900 nm and 1500 nm," Phys. Med. Biol. 46, 167-81 (2001). 12. X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang,X. H. Hu, "Determination of Complex Refractive Index of Polystyrene Microspheres from 370 to 1610nm," Phys. Med. Biol. 48, 4165-4172 (2003). 13. A. Roos, "Interpretation of integrating sphere signal output for nonideal transmitting samples," Appl. Opt. 30, 468-474 (1991). 14. J. W. Pickering, S. A. Prahl, N. Vanwieringen, J. F. Beek, H. J. C. M. Sterenborg,M. J. C. Vangemert, "DoubleIntegrating-Sphere System for Measuring the Optical-Properties of Tissue," Appl. Opt. 32, 399-410 (1993). 15. L. Wang, S. L. Jacques,L. Zheng, "MCML--Monte Carlo modeling of light transport in multi-layered tissues," Comput Methods Programs Biomed 47, 131-46 (1995). 16. S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star,M. J. van Gemert, "Optical properties of Intralipid: a phantom medium for light propagation studies," Lasers Surg. Med. 12, 510-9 (1992). 17. S. T. Flock, "The Optical Properties of Tissues and Light Dosimetry." Hamilton, Ontario, Canada: McMaster University, 1988. 18. C. J. M. Moes, M. J. C. van Gemert, W. M. Star, J. P. A. Marijnissen,S. A. Prahl, "Measurements and calculations of the energy fluence rate in a scattering and absorbing phantom at 633 nm," Appl. Opt. 28,, 22922296 (1989). 19. H. G. van Staveren, C. J. M. Moes, J. van Marle, S. A. Prahl,M. J. C. van Gemert, "Light scattering in Intralipid-10% in the wavelength range of 400-1100 nanometers,," Appl. Opt. 30, 4507-4514 (1991). 20. S. L. Jacques, "Optical properties of "Intralipid", an aqueous suspension of lipid droplet," (Oregon Medical Laser Center, 1998), http://omlc.ogi.edu/spectra/intralipid/index.html. 21. C. Chen, J. Q. Lu,X. H. Hu, "OPDISM Optical Parameters Determined by Integrating Sphere Measurement," (Biomedical Laser Laboratory, 2006), http://bmlaser.physics.ecu.edu. 22. M. C. Wang,E. Guth, "On the theory of multiple scattering, particularly of charged particles," Phys. Rev. 84, 1093-1111 (1951). 23. A. A. Kokhanovsky, "Small-angle approximations of the radiative transfer theory," J. Phys. D: Appl. Phys. 30, 2837-2840 (1997). 24. Z. Song, K. Dong, X. H. Hu,J. Q. Lu, "Monte Carlo simulation of converging laser beams propagating in biological materials," Appl. Opt. 38, 2944-2949 (1999). 25. Y. Du, "Optical Properties of Porcine Dermis in the Near Infrared Region between 900nm and 1500nm." Greenville, NC: East Carolina University, 2000, pp. 123. 26. V. R. Weidner,J. J. Hsia, "Reflection properties of pressed polytetrafluoroethylene powder," J. Opt. Soc. Am. 71, 856-861 (1981). 27. H. Ding, J. Q. Lu, W. A. Wooden, P. J. Kragel,X. H. Hu, "Refractive Indices of Human Skin Tissues at Eight Wavelengths and Estimated Dispersion Relations between 300 and 1600nm," Phys. Med. Biol. 51 (2006). 28. G. Hale,M. Querry, "Optical constants of water in the 200nm to 200 micrometer wavelength region," Appl. Opt. 12, 555-563 (1973). 29. K. Li, J. Q. Lu, R. S. Brock, B. Yang,X. H. Hu, "Quantitative Modeling of Skin Images Using Parallel Monte Carlo Methods," in Advanced Biomedical and Clinical Diagnostic Systems III, T. Vo-Dinh, W. S. Grundfest, D. A. Benaron, and G. E. Cohn, eds., Proc. SPIE 5693,82-87 (2005). 30. X. Ma, J. Q. Lu,X. H. Hu, "Effect of surface roughness on determination of bulk tissue optical parameters," Opt. Lett. 28, 2204-6 (2003). 31. X. Ma, J. Q. Lu, H. Ding,X. H. Hu, "Bulk optical parameters of porcine skin dermis tissues at eight wavelengths from 325 to 1557nm," Opt. Lett. 30, 412-414 (2005). ____________________________________________________________________________________________


Introduction
Development of a primary method as the "gold standard" for determination of macroscopic optical parameters of turbid media is highly desired in the optical study of turbid media such as biological tissues.The volumetric nature of multiple light scattering in turbid medium makes such a task a challenging inverse problem because of the difficulties in accurate measurement and subsequent modeling.Optical parameters may be defined precisely on the basis of a local theory with the radiative transfer equation (RTE) that has been widely used for modeling light transportation in turbid media [1].Established on the energy conservation principle, the RTE employs the absorption coefficient μ a , scattering coefficient μ s and phase function p(s, s') as the characterizing parameters, where s and s' are unit vectors along the light propagation directions.The RTE, however, has to be supplemented by proper boundary conditions to form a well-posted boundary value problem for modeling of a specific optical system.A reasonable one appears to be given by the Fresnel formulae for plane electromagnetic waves [2] in which light transportation at an interface depends on the mismatch of the refractive index between the two neighboring media.Even though no analytical model exists for complex turbid media such as the human skin tissues, experimental data on coherence reflectance curves demonstrate a good agreement with the Fresnel formulae [3].Consequently, it became customary to use μ a , μ s , p(s, s') and the real refractive index n r for optical characterization of a turbid medium.Furthermore, the Henyey-Greenstein (HG) function p(cosα) [1,4], where cosα=s⋅s', has been widely used to represent the phase function in turbid media, including soft human tissues [5].The HG function is a single-parameter function determined by the mean value of cosα, i.e., g=<cosα>, where g is also called the anisotropy factor.With p(cosα) as the phase function, optical characterization of a turbid sample is thus reduced to the determination of 4 scalar parameters: (μ a , μ s , g, n r ).In general, these parameters are functions of wavelength and should provide distinctive spectral attributes as the signature of a medium.
A primary method for determining the above parameters from measured signals should possess these features: (1) an experimental component that can be used to acquire measured signals in a wide spectral region and calibrated with readily available standard materials; (2) a fast modeling component that can be validated to obtain calculated signals accurately; (3) a robust inverse algorithm based on a well-posted inverse problem.Several methods have been developed to determine the optical parameters within the framework of RTE and Fresnel formulae as we discussed above or a diffusion approximation that can take into account of the anisotropic nature of scattering [6][7][8][9][10].None of them, however, has been rigorously shown to have the above features as a primary method.Among these, the method with an integrating sphere to measure the diffuse reflectance R d , diffuse transmittance T d and a spatial filtering setup to measure the collimated transmittance T c has been demonstrated to be capable of acquiring measured signals of relative high signal-to-noise ratios in a wide spectral region with the Monte Carlo (MC) modeling [6,8].This capacity is important in determination of the wavelength dependence of optical parameters with a wide-band light source.The method of combining integrating sphere based measurements with MC modeling has been shown numerically to yield unique solutions of optical parameters within reasonable ranges [11,12].
Despite of the strong potential, however, the accuracy and robustness of the integrating sphere based method have not been established because of the lack of fast numerical tools that can be validated for accurate modeling of the complex configuration introduced by the sample, the holder and the sphere [13,14].In this report, we present a system of experimental setups and Monte Carlo modeling tools that have been developed over the last few years in our laboratory to obtain the full set of parameters (μ a , μ s , g, n r ).The Monte Carlo method has been proved to be capable of yielding exact solutions of RTE boundary value problems in [12,15] and by results presented in this report.Therefore, the method described here can be applied to characterize turbid samples over large ranges of the above parameters.As an application, we determined the optical parameters of 10% and 20% intralipid samples in the spectral region from 550 to 940nm and from 550 to 1630nm, respectively.Intralipid has been used widely as a phantom for biological tissues while significant discrepancy existed in the published values of the parameter, especially in g [10,[16][17][18][19][20].Our results provide the complete spectra of optical parameters of this important phantom material from visible to the short-wave near-infrared to the existing body of knowledge that can be used to compare the accuracy of different methods.They also show that the presented system can be used as a primary method to acquire and invert measured signals into the optical parameters of homogeneous turbid samples in a wide spectral region.Furthermore, data acquisition in this system requires only relative measurements and thus needs no light source of calibrated intensity.The MC modeling tools described in this report are provided as a public domain source code package to interested researchers in this community [21].We will describe the method in the next section followed by the presentation and discussion of the instrument calibration, code validation and intralipid results.

Theoretical and numerical methods
Assuming an axial symmetry, we introduce here an approximate expression of the forward transmittance T f based on the RTE solutions published in [22,23] that can be used for validation of a MC code built for calculated signals.Consider a semi-infinite layer of turbid medium occupying 0 ≤ z ≤ D with z as the axis of symmetry and z=0 as the entrance surface.Inside the medium, the time-independent and source-free RTE is given by 2 0 0 ( , ) ( , ) ( , ') ( , ')sin ' ' ' where L(r, s) is the light radiance, r = xx + yy + zz with (x, y, z) as the unit vectors of respective Cartesian coordinate axis, s=sinθcosφx + sinθsinφy + cosθz is the unit vector along the propagation direction of light and μ t = μ a + μ s is the attenuation coefficient.For the case of a z-propagating plane wave incident on the medium with axial symmetry, one find L(r, s) = L(z, χ) with χ = cosθ.If the phase function p(s, s') is substituted by the HG function p(cosα), Eq.( 1) can be reduced to the following form ( In deriving the above equation, we used the fact that 0 (cos ) (cos ) where P l (m) ( ( ) m l P χ ) is the l-th order (associated) Legendre polynomial and {w l } are the Legendre expansion coefficients of the HG function.
For a normalized HG function of the form it is known that [1] 4πw l = (2l+1)g l .( 6 ) For a plane wave representing a collimated incident beam, the boundary condition may be expressed as the following by neglecting the effect of backscattered light at z=0 0 (0, ) , ( , ) 0, 2 where I 0 is the incident irradiance and δ(χ-1) is the Dirac delta function.It has been shown that the boundary-value problem of Eqs. ( 2) and ( 7) can be solved inside the homogeneous medium layer as [22,23] By expanding the exponential factor within the series and substituting w l with Eq. ( 6), the above series can be written in the following form to speed up the convergence It is obvious that the first term in the solution represents the direct or unattenuated component of the incident beam and thus leads to the Beer-Lambert law.The above solution, however, is exact only for the cases of μ s = 0 or g=1 because otherwise the boundary condition at z=0 given by Eq. ( 7) is inaccurate due to the presence of backscattered light.It is thus expected that this solution should be limited to turbid media with strongly forward phase functions.To use these results for validation of MC calculations, we integrate the radiance over a solid angle for incidence and transmitted irradiances and then obtain the forward transmittance T f through a sample of optical thickness τ=μ t D as their ratio where a=μ s /μ t is the single-scattering albedo and χ 0 = cosθ 0 is determined by the angle of acceptance (2θ 0 ) from the detector behind a pinhole or slit in measuring T f .The upper limits of the summation series, J and L, are to be determined by the convergences of the series.The first term on the right-hand side of Eq. ( 10) is the collimated transmittance T c .
A signal MC code has been built on the basis of our previous ones used for calculation of optical signals [11,24].Briefly, the MC method is a statistical simulation method of photon transportation which provides light distribution data equivalent to the solution of a boundaryvalue problem consisting of RTE and Fresnel formulae.In the signal MC code, each photon in a collimated incident beam is tracked along its trajectory in a sample assembly of a three-layer structure in which the turbid sample is sandwiched between two glass plates.Before tracking starts, the life time or total pathlength of the tracked photon in the turbid sample is determined by a random number distributed according to μ a .The free pathlength of the photon between two consecutive scattering events is determined similarly by a random number distributed according to μ s while the direction of scattered photon is determined by two random numbers in accordance of the HG phase function of g [24].The tracking of photons stops when they exit the sample assembly which are sorted into different fractions as defined in Fig. 1(a) and divided by the incident photon number to obtain the calculated signals of (R d , T d , T f ) cal .The variances or uncertainties in the calculated signals were found to be less than 1% if the number of tracked photons were set at 7.8x10 5 or larger.The calculated signals are compared to the measured ones to solve the inverse problem.
It should be noted that the normal-hemispherical configuration assumed in the signal MC code for the calculated R d and T d , as depicted in Fig. 1(a), is different from the configuration encountered in the measurements of R d and T d with an integrating sphere.In the former case one assumes a normally incident beam with the exiting photons not returning to the sample assembly, while in the latter case the photons enter and exit repeatedly into the sample assembly from the integrating sphere in a hemispherical-hemispherical configuration after the first incidence.Since the light reflectance and transmittance at an interface given by the Fresnel formulae depend on the incident angle, it is clear that the two types of definitions of R d and T d need to be carefully compared to examine their consistency in various samples with different parameters and geometry.Moreover, the equations used for obtaining R d and T d from the detected signals, given by Eq. ( 13) and ( 14) in the nest section, were derived without consideration of the glass plates in the sample assembly [25].Therefore, the effect of glass plates in the sample assembly should be also investigated.
To study these effects, which are very difficult to be determined analytically or experimentally, we developed a sphere MC code in which photons are tracked as they propagate between the sample assembly and the integrating sphere, including the baffle, before collection by the detector or escape from the sphere through the entrance port.The geometric shapes and sizes of the integrating sphere and detector used in the sphere MC code are the same as those of the experimental system.The photon reflection by the inner wall of the sphere was treated with a Lambertian model of hemispherical diffuse reflection [26] in which the reflection is guided by the Fresnel formulae with a randomly chosen surface normal.The wall reflectance is obtained by interpolating the values provided by the vendor between 400 and 1700nm.The sphere MC code has been validated by comparing the calculated R d with the measured values of diffuse reflectance standards, which has no glass plates; and it allows the quantitative comparison of different definitions of R d and T d and the investigation of the glass plate effect.

Experimental methods and inverse algorithms
Determination of (μ a , μ s , g, n r ) for a homogeneous turbid sample can be separated into two parts according to the difficulty of data inversion.If T f or the forward transmitted light irradiance I f is dominated by the direct component, as shown in Eq. ( 9) and (10), n r and μ t can be obtained from the coherent reflectance R c and forward transmittance T f by straightforward fitting to the Fresnel formulae and Beer-Lambert law, respectively.By comparison, a (=μ s /μ t ) and g have to be inverted from the scattered signals of R d and T d by a much more complicated iteration process based on MC modeling within the framework of the RTE and Fresnel formulae.The following description of the methods is thus divided into two subsections.In all measurements, the incident light was modulated by a mechanical chopper and light signals were detected by either a Si (400 to 940nm) or InGaAs (910 to 1700nm) photodiode connected to a preamplifier and a lock-in amplifier (SR830, Stanford Research).A second photodiode was used to produce normalization signal for elimination of the light source intensity fluctuation.The modulation frequency was set to either 370Hz for fast acquisition of strong R c signals (with incident angle between 45° and 80°) or 17Hz for detection of weak signals of T f , R d and T d with a high-gain and narrow bandwidth preamplifier.

Determination of the real refractive index n r and attenuation coefficient μ t
Specular or coherent reflectance R c of a turbid media such as the intralipid has been shown to be a gradually increasing function of the incident angle θ i .Furthermore, the measured coherent reflectance curve R c (θ i ) exhibit a very good agreement with the Fresnel formulae if a complex refractive index n=n r +in i is appropriately assigned to the turbid sample which provide a straightforward method to determine n r [3,27].We designed and constructed an automated reflectometer to measure R c (θ i ) which consists of a right-angle prism of BK7 glass and a photodiode with an aperture to detect the specularly reflected light beam, as shown in Fig. 2(a).The sample-prism assembly and the photodiode were rotated and translated separately to measure R c as the incident angle θ i was increased from 45° to 80° with the turbid sample held against the prism base.The design details of the reflectometer have been published elsewhere [3].Fig. 3 The real refractive index of 10% and 20% intralipid versus wavelength and the fitted Cornu dispersion relations (lines).The symbols and error bars represent the mean and standard deviation of the index determined from 6 measurements at each wavelength with 3 for s-and 3 for p-polarized incident beam [3].Inset: a measured coherent reflectance curve (black symbols) of 20% intralipid with a p-polarized incident beam at λ=633nm and fitted by the Fresnel formulae (red line).
The reflectometer was calibrated with deionized water by comparing the measured real refractive index from the critical angle θ c with the published values [28].Moreover, the water calibration provided the incident light intensity for determination of R c for θ i > θ c in which the system errors such as the air-prism reflection loss and prism absorption were eliminated.The dispersion of the real refractive index of water does not exhibit abnormal behavior between 400 and 1700nm [28].Therefore, we expect that the water based turbid samples such as intralipid should have similar dispersions.The coherent reflectance curves were measured with laser beams at 5 wavelengths between 442 and 1310nm for the 10% intralipid samples and at 8 wavelengths between 325 and 1550nm for the 20% intralipid samples.The n r were determined from these curves by a least-squares nonlinear fitting to the Fresnel formulae [3].A Cornu dispersion dispersion relation was applied for interpolation of the polarization averaged n r values by the follow equations where the wavelength λ is in the unit of nm.The coefficients in the above equation were found to be A=1.0055,B=7829.0 and C=−22269.6 for the 10% intralipid and A=1.2364, B=1261.1 and C=−9498.1 for the 20% intralipid.Eq. ( 11) was used for obtaining the real refractive index n r of the intralipid samples at desired wavelength for the MC calculations of R d and T d .The wavelength dependence of the polarization averaged n r is presented in Fig. 3 together with typical measured and fitted coherent reflectance curves.
The attenuation coefficient μ t was estimated from the measured plot of the forward transmitted light signal I f against the sample thickness D. Experimentally, an incident beam of irradiance I 0 at wavelength λ was obtained from a vertical line source of lamp filament coupled with a monochromater and collimated with a spherical mirror, as shown in Fig. 2(b).The nearly collimated beam has an elliptical profile with the vertical diameter of 30mm and horizontal diameter of 10mm at the entrance surface of the sample holder.The diameter of the beam incident on the sample was reduced to 6.35mm in diameter with a circular aperture.Another spherical mirror was employed to focus the forward transmitted light through the sample within a cone angle of 2θ 0 into an image of the vertical line source and filtered by a 0.5mm wide slit.A 30W tungsten-halogen lamp coupled with a monochromater (CM110, CVI Laser) with a 2nm resolution was used as a tunable light source to scan wavelength between 400 and 1700nm with a step size of 30nm.The viewing angle of the detector behind the slit in the horizontal plane was measured as θ 0 =1.0°, which has been used in all RTE and MC calculations of T f or I f presented in this report.We estimated μ t from the measured I f -D data by a Beer-Lambert law based fitting given by 0 ( ) where S is a numerical factor representing surface effect such as reflection loss.For samples of strong attenuation, with τ=μ t D > 4, the measured I f −D data deviates significantly from a straight line on a semi-log plot.In these cases, the estimated μ t values were determined by the two data points of the smallest D values because these I f values are much larger than others at large D values and their statistical weights dominate in the least-squares fitting.The values of μ t and n r were used as input parameters for the inverse determination of the albedo a and anisotropy factor g from the measured signals of R d and T d .

Determination of the scattering coefficient μ s and anisotropy factor g
The diffuse reflectance R d and transmittance T d of a turbid sample between two glass plates were measured with an integrating sphere (IS-060-SF, Labsphere Inc.) of 152mm diameter and with three unplugged ports.The two opposing ports of 38.1mm diameter were used to install the sample holder and a port plug, both carrying an aperture of 6.35mm diameter for light beam collimation.A holder of photodiode with a 3mm diameter sensor and an attached preamplifier was placed in the detector port of 12.5mm diameter for light detection.Both surfaces of the sample holder and photodiode were be placed flush with the inside surface of the sphere and were painted with diffuse white reflectance coating (WRC-680, Labsphere, Inc.) except the areas exposing the sample and photodiode sensor.The sample holder consisting of two BK7 optical windows of 25mm diameter and 3mm thickness as the glass plates with one glued to one end of a tapped aluminum tube and the other glued to a threaded ring which can be translated inside the tube by rotation.The sample thickness D was determined by the rotation angle of the ring to a precision of 0.008mm.In the measurements of T f , R d , and T d , the intralipid sample was filled into the holder through a small side hole by a syringe needle.Care must be taken to fill the sample holder with no air bubbles trapped between the glass plates, especially when D < 0.1mm.
The same schemes of light source and electronic detection were used in the spatial filtering, as shown in Fig. 2(b), and the integrating sphere measurements.The integrating sphere with the sample assembly was oriented in three configurations in a sequence, as shown in Fig. 1(b), to acquire three light readings from the photodiode (P T , P R and P C ) as the wavelength was tuned from 400 to 1630nm.While the P T and P R signals were due to the diffuse reflection and transmission from the turbid sample, P C was obtained as a signal related to the incident light power by rotating the sphere by 20° from the P R orientation so that the incident beam falls completely on the inner surface of the sphere.A baffle of same material as the inner wall of the integrating sphere is built in the sphere to prevent the detector from receiving light directly from the area illuminated by the incident beam.By considering the stead-state light distribution in an integrating sphere as a result of many rounds of light "bouncing" between the inner surfaces of sphere and the sample assembly/detector, one can derive the relations between the detector signals and R d and T d as the followings [14,25] 0 0 (1 ) cos 20 cos 20 where A=4πR 2 is the total surface area of the sphere, f is the area ratio of the three ports in the sphere to A and A s is the circular area of the sample assembly exposed to the integrating sphere.The factor of cos20° in Eq. ( 14) is used to account for the decreased power into the sphere at the oblique orientation for the P C measurement if the incident beam illuminates fully the entrance aperture of the sphere.
We developed an automated algorithm for determination of optical parameters from the measured signals of R d and T d .The calculated R d and T d from the signal MC code are compared with the measured signals to generate a total squared error function δ defined as ( ) ( ) To solve the inverse problem, the MC calculation is iterated with different a and g and the metric for guiding this process is provided by the function δ.The initial values of a and g are either estimated from the values of R d +T d and R d /T d , respectively, at the first wavelength or the inverse solution at the previous wavelength.Then an updated value of δ is obtained from another MC calculation with varied a and g along fixed directions with fixed step sizes.The directions and step sizes are then adjusted according to the relative change in δ in the previous two runs and this process is repeated until δ < δ m , where δ m =4x10 -4 is a pre-determined value consistent with the total experimental error in the R d and T d measurements.At the end of this iteration process, the optical parameters (μ a , μ s , g) are properly determined for samples with small or moderate μ t .For samples of large μ t where the semi-log plot of measured I c -D data deviate significantly from a straight line, MC calculated values of I f -D data are compared with the measured one to check the consistency.If the difference between the measured and calculated I c -D data is significant, the inverse solution from R d and T d is rejected for poor consistency.

Results
All MC calculations were carried out by tracking 7.8x10 5 photons on a single processing element (Intel xeon CPU of 3.06GHz and 1GB memory) of the parallel computing cluster in our lab.At this number of tracked photons, the variances in the output were found to be less than 1% and thus negligible.The MC codes for calculated signals and for sphere simulation were written in Fortran 90 and compiled by an Intel compiler on a Linux platform.

Validation of the signal MC code
With a single layer of turbid medium of an optical thickness τ assumed in the signal MC code, we compared the calculated T f , defined within the cone angle of 2θ 0 as shown in Fig. 1(a), with the RTE solution expressed in Eq. (10).To make the optical configuration in the MC calculations close to the one assumed for the RTE solution with an incident light beam of infinite diameter, we tested and adopted a diameter of the normally incident beam at 3.2mm which is much larger than the layer thickness D ranging between 0.03 and 0.51mm.The region for tracking the photons in the MC calculations was also set as borderless in the x-y plane.In obtaining T f (τ) from Eq. ( 10) as the RTE solution, we checked the convergence of the series and determined the upper limits for the summations to be J=20 and L=35.The results of MC and RTE calculated T f (τ) are compared in Fig. 4 with μ t =31(mm -1 ).
It can be seen from these results that the MC calculated T f (τ) agrees well with the RTE results for the cases of large g (> 0.75) and small τ (< 10), where the RTE solutions are sufficiently accurate under the boundary condition described by Eq. (7).For other cases, the backscattered photons contribute significantly to the irradiance at the boundary surface z=0 and RTE solution of Eq. ( 10) is expected to become inaccurate.MC calculation of the diffusely reflected photon density at z=0 confirmed the above statement (data not shown here).These results provided the validation of the signal MC code.

Calibration of the integrating sphere setup and validation of the sphere MC code
We calibrated the integrating sphere setup by R d measurements with two diffuse reflectance standards of nominal R d at 10% and 20% (SRS-40-010 and SRS-20-010, Labsphere, Inc.).The same set of measured data was also used to validate the sphere MC code by comparison with the calculated R d .In the sphere MC code, the vendor supplied reflectance data were used as the input parameters at different wavelengths for the sample and the calculated R d was obtained as the number ratio of photons received by the detector to that of the incident photons.The measured and calculated R d are plotted in Fig. 5  wavelength λ and compared with the values interpolated from the vendor supplied data.These results provided the validation of the sphere MC code and demonstrated that the errors in the measured R d and T d , in comparison with the vendor supplied data, were less than ±4%.
Fig. 5 The measured and calculated diffuse reflectance R d of two reflectance standards versus wavelength.The solid lines are interpolated results from the vendor supplied data.
Fig. 6 The comparison of R d and T d calculated with the signal and sphere MC code for samples between two glass holder plates with different edge reflectivity R. The optical parameters for the samples were set as μ s =5mm -1 , μ a =0.1mm -1 , n r =1.30 and (a) g=0.0;(b) g=0.9.The symbols represent the signal MC code results while the lines represent the sphere MC results: R=100%: and dashed lines; R=50%: ∇ and dotted lines; R=0%: o and solid lines.

Investigations of the sample holder effect and consistency of the R d and T d definitions
The validated signal and sphere MC codes were used to study these effects by comparing the calculated values of R d and T d .We note here that Eqs. ( 13) and ( 14) were derived without consideration of the sample holder [25] and therefore their consistency with the calculated signals from the signal MC code needs to be examined.For this purpose, we used the sphere MC code to simulate the detected signals from turbid samples of different thickness within a sample holder of two identical BK7 glass windows.The simulated values of P C , P R and P T were calculated by tracking and sorting each photon from the normally incident beam as they interact with the sample assembly and integrating sphere before finally escape from the system or collected by the detector.The simulated values of R d and T d were then determined from Eq. ( 14) and (15).In both MC codes we set the diameter of the turbid sample at 18mm; diameter of the glass window at 18mm, thickness at 3mm and refractive index at 1.60 with three different edge reflectivity of R=0%, 50% and 100%.An incident beam of diameter of 3.18mm was assumed and the calculated results are presented in Fig. 6.Based on these data, we concluded that the measured R d and T d defined in Eq. ( 13) and ( 14) are consistent with the definitions based on the photon number ratios used in the signal MC code.However, the edge reflectivity of the glass plates or the "leaked" photons from the edge was shown to affect the measured signals significantly.This requires a careful design of sample holder used in the experiment for accurate measurement of R d and T d and modeling of exact configuration in the signal MC code to ensure the stability of the inverse solution.

Determination of the optical parameters (μ a , μ s , g) of 10% and 20% intralipid solutions
As an application of the method presented above, we determined the optical parameters of 20% intralipid samples (Fresenius Kabi Clayton, L.P., Clayton, NC; lot #: 1022848) between 400 and 1630nm with a step size of 30nm.To compare with the published values of optical parameters of intralipid samples, often determined with 10% concentration [10,[16][17][18][19][20], we have also determined the parameters of 10% intralipid by diluting the 20% solution with deionized water using a volume ratio of 1:1.First, we measured the forward transmitted light signal I f of different thickness D with the spatial filtering setup shown in Fig. 2(b) and estimated from these data the attenuation coefficient μ t as a function of wavelength using the Beer-Lambert law as discussed before.Two examples of I f -D data with 20% intralipid samples are presented in Fig. 7.The I f measurements were carried out with 4 samples of D=0.079, 0.119, 0.159, 0.199mm to estimate μ t for the 10% intralipid.
The integrating sphere setup shown in Fig. 1(b) was then used to acquire 3 sets of measured signals of R d and T d from 3 samples of D=0.159, 0.191, 0.222mm between 400 and 1630nm for the 20% intralipid and 3 samples of D=0.238, 0.318, 0.397mm for the 10% intralipid.The edges of the two BK7 glass windows, used as the sample holder, were coated with silver to achieve a reflectivity close to the R value of 100% used in the signal MC code.The refractive index n r and attenuation coefficient μ t of the sample and the refractive index of the BK7 glass were used as the input parameters in the signal MC code at each wavelength.With the inverse procedure described around the Eq. ( 15), we determined the optical parameters of the intralipid samples by comparing the calculated signals with the measured R d and T d .Once the optical parameters at the first wavelength were determined, it took about 8 iterations on average for the inverse procedure to converge with a total CPU time of about 80 seconds at each of the successive wavelengths.After the optical parameters (μ a , μ s , g) were inversely determined from the measured R d and T d , the MC calculated I f -D data was compared to the measured data to check the consistency of the results, as shown in Fig. 7.The optical parameters of the 20% intralipid samples between 550 and 1630nm were found to pass this consistency test and are presented in Fig. 8.The results of the 10% intralipid samples between 550 and 940nm are plotted in Fig. 9 together with the curves of μ s and g by fitting the experimental data reported in [16].
Fig. 8 The wavelength dependence of the optical parameters of 20% intralipid samples.The symbols and error bars are the means and standard deviations obtained from three measurements with samples of thee different thicknesses.The dashed lines are visual guides and the blue solid line is a power law fitting of μ s =Cλ -2.407 , where C is a constant.
Fig. 9 The wavelength dependence of the optical parameters of 10% intralipid samples.The symbols and error bars are the means and standard deviations obtained from three samples.The lines are curves fitted from the μ s (dash-dot) and g (dash-dot-dot) data in [16].

Discussion
The inverse algorithm and an early version of the signal MC code with a Mie theory based phase function have been used to determine the refractive index of polystyrene spheres from the measured data with the same integrating sphere, which demonstrated a good agreement with the published data in the visible region [12].Our experience in the inverse determination of optical parameters from these previous [11,12] and current studies reveals sensitive relations between T f and μ t , R d +T d and a, R d /T d and g.These relations ensure the fact that the determination of optical parameters (μ a , μ s , g), or equivalently (μ t , a, g), from the measured signals of R d , T d and T f (or I f ) is well-posted inverse problem with a unique and wellbehaved solution, at least in the ranges of parameters most likely encountered in the investigations of biological tissues.To demonstrate this point, we calculated the total squared error δ from the measured signals as defined in Eq. ( 15) at two wavelengths using the signal MC code.The choice of these two cases are based on the value of the single-scattering albedo a which reaches 99.6% for high turbidity of the 20% intralipid at λ=550nm and decreases to 77.8% at 1450nm.A total of 10,000 MC simulations were performed in each case along either directions of the a and g axes with a fixed step size from the location of the optimized parameters and the contour graphs of δ are presented in Fig. 9.These results indicate clearly the existence of only one minimum of δ within reasonably large ranges of a and g.Inverse algorithms are relatively easy to develop in these situations as long as the accuracy of modeling can be ensured.At the same two wavelengths, we have also tested the sensitivity of inversely determined optical parameters on the measured data of refractive index n r , diffuse reflectance R d , and diffuse transmittance T d from one sample of 20% intralipid.We determined the new values of the parameters with the same inverse procedure by changing one of the measured data from the values used as the input data to the signal MC code.The results are shown in Fig. 10.Since the total uncertainty in the n r was determined to be ±0.002[3], the errors in the optical parameters caused by uncertainty of n r can be neglected.On the other hand, μ a exhibits high sensitivity to the measured data of R d and T d in the sample of high turbidity at λ=550nm.Similarly, g is very sensitive to the measured data in the sample of moderate turbidity at λ=1450nm.In contrast, μ s is very stable relative to the errors in the measured data because scattering dominate the interaction in all cases.Based on these results, we conclude that the uncertainties in the optical parameters determined by our method should be about ±25% for μ a , ±2% for μ s and ±10% for g.
It should be pointed out here that the optical parameters inversely determined from the measured R d and T d between 400 and 530nm were found to fail the consistency test on the I f -D data for the 20% intralipid samples.We attributed this to the inability of using the Beer-Lambert law to estimate μ t from the measured I f in which contribution by collimated transmitted light is very small at these wavelengths.Because of the large optical thickness (τ > 5) of the samples at even the smallest D (between 0.040 and 0.080mm with an error of ±0.005mm), μ t values at these wavelengths were underestimated from these data points which led to the underestimated μ s at the end of inverse determination.This caused the MC calculated I c -D data to lie significantly below the measured one due to the reduced contribution to I f by the underestimated forward transmission.Inverse determination in these cases is possible without the Beer-Lambert law fitting, however, by comparison of MC calculated R d , T d and I f -D data with the measured signals using (μ a , μ s , g) as three free parameters.Research on an efficient algorithm is currently underway.
Taking advantage of the algorithm simplicity and versatility of MC modeling method, we were able to show that the definitions of the calculated signals of R d and T d are equivalent to the measured signals.We have also analyzed the effect of sample holder on the measured signals to improve the holder design that was critical in achieving accurate measurements.Taken all together, the accurate measurement and modeling of the measured signals provides the foundation required for robust inverse determination of optical parameters.Therefore, the system of reflectometer, integrating sphere and spatial filtering setups and related modeling tools described in this report provides a primary method for determination of the optical parameters (μ a , μ s , g, n r ), in a large dynamic range, of highly turbid samples.While a widely used HG function was assumed in this work, extension of our method to other forms of the phase function is expected to be straightforward by defining g as the first moment of the phase function with possibly additional parameters.This method can be further improved by combining the measurement of coherent reflectance with the integrating sphere measurement into a single instrument using multiple channels of lock-in detection and a high-brightness tunable light source.The availability of the validated modeling tools online and methodology affords investigators the opportunity to characterize various phantoms and calibrate other methods including those designed for in vivo parameter measurements with imaging cameras [10,29].
Several previous papers reported all or some of the three parameter of μ a , μ s and g of intralipid samples at 633nm or selected wavelengths between 400 and 1100nm determined from the collimated transmittance T f and spatially resolved fluence rate data or the diffuse reflectance data of intralipid and intralipid with added absorbers (ink) samples modeled by diffusion equation or Mie theory [10,[16][17][18][19][20].We first note that the wavelength dependence of μ t reported in [16] agrees well with our μ s (≈μ t ) data of 20% intralipid as shown in Fig. 8; both exhibit a power dependence on wavelength described by μ s =Cλ -p with p=2.33 between 500 and 1000nm in [16] and p=2.41 between 550 and 1630nm in our study.In addition, the peak of μ a at 1450nm shown in Fig. 8 is consistent with that of the water [28], indicating that the absorption of intralipid solution is mainly determined by its water component for λ>1300nm.We have also include the wavelength dependence of μ s and g, fitted from the data in [16], in the Fig. 9 for comparison with our results of the 10% intralipid.While the μ s data show relatively good agreement, the g data have large differences from the values determined in this work, with our g data much smaller.This may be caused by the photons "leaking" out of side of the sample holder that was accounted for in our MC modeling.We further summarize the reported optical parameters of intralipid samples at 633nm in Table 1 for comparison with our results.It is obvious that large discrepancies exist between these results which may indicate that in additional to the differences caused by the measurement and modeling methods the material difference among intralipid samples could also play a role.One particular point that can be noted is that the absorption coefficient μ a determined in this work is more than one order of magnitude larger than most of the reported values.One possibility may lie in the small sample thickness used in our experiments which could place a lower limit of μ a detection at about 0.1mm.However, this hypothesis is not supported by the nearly 2-fold decrease in μ a between the 20% and 10% intralipid, as shown by the data in Figs. 8 and 9 and Table 1.Therefore, the lower limit of μ a detection by the integrating sphere based method is currently unknown and should be investigated in the future.In summary, we presented a primary method for determination of the optical parameters of (μ a , μ s , g, n r ) of a homogeneous turbid sample with smooth surfaces in a wide spectral region within the framework of radiative transfer theory and Fresnel formulae [1].This method consists of experimental setups that are relatively easy to construct and calibrate and fast modeling tools using the versatile MC codes that have been validated and are made available as the public domain source codes [21].We demonstrated the utility of this method by performing the inverse determination of optical parameters of a highly turbid medium of 10% and 20% intralipid between 550 and 940nm, 550 and 1630nm, respectively.Extension of this method to the cases of biological tissue samples with rough surfaces has been investigated [30,31] and is to be completed in the future by significant improvement in accuracy of surface roughness measurement, modeling tool validation and computing speed.

Fig. 1 (
Fig. 1 (a) Definitions of calculated signals.The dotted lines representing the "virtual" sphere in the signal MC code.(b) Three orientations of a three-port integrating sphere for acquisition of the detected signals of P R , P T and P C .SA: the three-layer sample assembly, B: baffle, PD: photodiode, PA preamplifier.

Fig. 4
Fig. 4 The calculated T f by RTE and the signal MC code versus the optical thickness τ with μ t =31(mm -1 ): (a) different g with a=96.8%;(b) different a with g=0.90.The solid red line represents the Beer-Lambert law and symbols (Δ) are the MC calculated collimated transmittance with unscattered photons.

Fig. 7
Fig. 7 The measured forward transmitted light signal I f versus the thickness D of 20% intralipid samples: (a) λ=550nm; (b) λ=1450nm.The solid red lines are the fitted results based on the Beer-Lambert law for estimating μ t and the dashed blue lines are calculated results from the sphere MC code with the inversely determined optical parameters.

Fig. 10 Fig. 11
Fig.10The contour graphs of the total squared error δ as a function of albedo a and anisotropy factor g obtained from the measured signals at two values of wavelength λ.The measured data were acquired from a 20% intralipid sample with thickness D=0.159mm.

Table 1
The optical parameters of intralipid solutions at λ=633nm f (D) curve + inside fluence rate/added absorber + diffusion model Method b: reflectance imaging + Mie model Method c: I f (D) curve + integrating sphere + MC model