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 = 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.(l) can be reduced to the following form

In deriving the above equation, we used the fact that

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₀ 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 τ=μ_tD as their ratio

where a=μ_s/μ_t is the single-scattering albedo and χ₀ = cosθ₀ is determined by the angle of acceptance (2θ₀) 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⁵ or larger. The calculated signals are compared to the measured ones to solve the inverse problem.

 

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.

Download Full Size | PDF

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

 

Fig. 2 (a) A schematic of the reflectometer system; PM: prism, PD: photodiode, S: sample. (b) The spatial filtering setup for measurement of the forward transmitted light signal I_f; L: lamp, M: monochromater, C: chopper, M1/M2: concave mirrors, LF: longpass filter, SA: sample assembly, SL: slit, PA: preamplifier, LIA: lock-in amplifier, PC: computer.

Download Full Size | PDF

 

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

Download Full Size | PDF

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₀ 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θ₀ 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 θ₀=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 τ=μ_tD > 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² 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⁵ 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.

 

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.

Download Full Size | PDF

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θ₀ 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%.

 

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.

Download Full Size | PDF

 

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.

Download Full Size | PDF

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.

 

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.

Download Full Size | PDF

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

 

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.

Download Full Size | PDF

 

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

Download Full Size | PDF

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.

 

Fig. 10 The 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.

Download Full Size | PDF

 

Fig. 11 The bar graphs of the relative changes in the optical parameters when the refractive index n_r or diffuse reflectance R_d or diffuse transmittance T_d is changed from the input values of n_r=1.3619, R_d=47.4%, T_d=30.4% at λ=550nm and n_r=1.3516, R_d=8.82%, T_d=16.3% at λ=1450nm to the signal MC code. The measured data were from the same sample used in Fig. 10 calculations. Some error bars of μs are too small to be seen.

Download Full Size | PDF

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.

Table 1. The optical parameters of intralipid solutions at λ=633nm

View Table

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.