## Abstract

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

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

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

where L(**r**, **s**) is the light radiance, **r** = x**x** + y**y** + z**z** 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.(l) can be reduced to the following form

In deriving the above equation, we used the fact that

$$\phantom{\rule{28.5em}{0ex}}=\sum _{l=0}^{\infty}{w}_{l}\left\{{P}_{l}\left(\mathrm{cos}\phantom{\rule{.2em}{0ex}}\theta \right){P}_{l}\left(\mathrm{cos}\phantom{\rule{.2em}{0ex}}\theta \text{'}\right)+2\sum _{m=1}^{l}\frac{\left(l-m\right)!}{\left(l+m\right)!}{P}_{l}^{m}\left(\mathrm{cos}\phantom{\rule{.2em}{0ex}}\theta \right){P}_{l}^{m}\left(\mathrm{cos}\phantom{\rule{.2em}{0ex}}\theta \text{'}\right)\mathrm{cos}\left[m\left(\phi -\phi \text{'}\right)\right]\right\},$$

which yields

where *P*_{l}
(*m*) (${P}_{l}^{m}$
(*χ*)) 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]

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

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 boundary-value 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.8×10^{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.

## 3. 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 If 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 nr 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].

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

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 3 mm 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]

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}=4×10^{-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.

## 4. Results

All MC calculations were carried out by tracking 7.8×10^{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 against the
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%.

#### 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 PC, 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–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].

### 5. 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 a well-posted inverse problem with a unique and well-behaved 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 measurable 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–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.

Method a: I_{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

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.

## Acknowledgement

We thank Xiaoyan Ma, R. Scott Brock and Yalin Ti for their helps.

## References and Links

**01. **H. C. van de Hulst, *Multiple light scattering: tables, formulas, and applications*, vol. 1 & 2. (Academic Press, New York, 1980).

**02. **M. Born, E. Wolf, and A. B. Bhatia, *Principles of optics: electromagnetic theory of propagation, interference and diffraction of light*, 7th ed. (Cambridge University Press, Cambridge, England, 1999).

**03. **H. Ding, J. Q. Lu, K. M. Jacobs, and 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). [CrossRef]

**04. **L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys J **93**, 70–83 (1941). [CrossRef]

**05. **S. L. Jacques, C. A. Alter, and S. A. Prahl, “Angular dependence of HeNe laser light scattering by human dermis,” Lasers in Life Sciences **1**, 309–333 (1987).

**06. **V. G. Peters, D. R. Wyman, M. S. Patterson, and 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). [CrossRef] [PubMed]

**07. **S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding-doubling method,” Appl. Opt. **32**, 559–568 (1993). [CrossRef] [PubMed]

**08. **I. V. Yaroslavsky, A. N. Yaroslavsky, T. Goldbach, and H.-J. Schwarzmaier, “Inverse hybrid technique for determining the optical properties of turbid media from integrating-sphere measurements,” Appl. Opt. **35**, 6797–6809 (1996). [CrossRef] [PubMed]

**09. **C. K. Hayakawa, B. Y. Hill, J. S. You, F. Bevilacqua, J. Spanier, and 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). [CrossRef] [PubMed]

**10. **A. H. Hielscher, J. R. Mourant, and 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). [CrossRef] [PubMed]

**11. **Y. Du, X. H. Hu, M. Cariveau, X. Ma, G. W. Kalmus, and J. Q. Lu, “Optical properties of porcine skin dermis between 900 nm and 1500 nm,” Phys. Med. Biol. **46**, 167–81 (2001). [CrossRef] [PubMed]

**12. **X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang, and X. H. Hu, “Determination of Complex Refractive Index of Polystyrene Microspheres from 370 to 1610nm,” Phys. Med. Biol. **48**, 4165–4172 (2003). [CrossRef]

**13. **A. Roos, “Interpretation of integrating sphere signal output for nonideal transmitting samples,” Appl. Opt. **30**, 468–474 (1991). [CrossRef] [PubMed]

**14. **J. W. Pickering, S. A. Prahl, N. Vanwieringen, J. F. Beek, H. J. C. M. Sterenborg, and M. J. C. Vangemert, “Double-Integrating-Sphere System for Measuring the Optical-Properties of Tissue,” Appl. Opt. **32**, 399–410 (1993). [CrossRef] [PubMed]

**15. **L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput Methods Programs Biomed **47**, 131–46 (1995). [CrossRef] [PubMed]

**16. **S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. van Gemert, “Optical properties of Intralipid: a phantom medium for light propagation studies,” Lasers Surg. Med. **12**, 510–9 (1992). [CrossRef] [PubMed]

**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, and S. A. Prahl, “Measurements and calculations of the energy fluence rate in a scattering and absorbing phantom at 633 nm,” Appl. Opt. **28**,, 2292–2296 (1989). [CrossRef] [PubMed]

**19. **H. G. van Staveren, C. J. M. Moes, J. van Marle, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nanometers,,” Appl. Opt. **30**, 4507–4514 (1991). [CrossRef] [PubMed]

**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, and X. H. Hu, “OPDISM - Optical Parameters Determined by Integrating Sphere Measurement,” (Biomedical Laser Laboratory, 2006), http://bmlaser.physics.ecu.edu.

**22. **M. C. Wang and E. Guth, “On the theory of multiple scattering, particularly of charged particles,” Phys. Rev. **84**, 1093–1111 (1951). [CrossRef]

**23. **A. A. Kokhanovsky, “Small-angle approximations of the radiative transfer theory,” J. Phys. D: Appl. Phys. **30**, 2837–2840 (1997). [CrossRef]

**24. **Z. Song, K. Dong, X. H. Hu, and J. Q. Lu, “Monte Carlo simulation of converging laser beams propagating in biological materials,” Appl. Opt. **38**, 2944–2949 (1999). [CrossRef]

**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 and J. J. Hsia, “Reflection properties of pressed polytetrafluoroethylene powder,” J. Opt. Soc. Am. **71**, 856–861 (1981). [CrossRef]

**27. **H. Ding, J. Q. Lu, W. A. Wooden, P. J. Kragel, and 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). [CrossRef] [PubMed]

**28. **G. Hale and M. Querry, “Optical constants of water in the 200nm to 200 micrometer wavelength region,” Appl. Opt. **12**, 555–563 (1973). [CrossRef] [PubMed]

**29. **K. Li, J. Q. Lu, R. S. Brock, B. Yang, and 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). [CrossRef]

**30. **X. Ma, J. Q. Lu, and X. H. Hu, “Effect of surface roughness on determination of bulk tissue optical parameters,” Opt. Lett. **28**, 2204–6 (2003). [CrossRef] [PubMed]

**31. **X. Ma, J. Q. Lu, H. Ding, and X. H. Hu, “Bulk optical parameters of porcine skin dermis tissues at eight wavelengths from 325 to 1557nm,” Opt. Lett. **30**, 412–414 (2005). [CrossRef] [PubMed]