Quantifying the properties of two-layer turbid media with frequency-domain diffuse reflectance

Noncontact, frequency-domain measurements of diffusely reﬂected light are used to quantify optical properties of two-layer tissuelike turbid media. The irradiating source is a sinusoidal intensity- modulated plane wave, with modulation frequencies ranging from 10 to 1500 MHz. Frequency-dependent phase and amplitude of diffusely reﬂected photon density waves are simultaneously ﬁtted to a diffusion-based two-layer model to quantify absorption ~m a ! and reduced scattering ~m s 9! parameters of each layer as well as the upper-layer thickness ~ l ! . Study results indicate that the optical properties of two-layer media can be determined with a percent accuracy of the order of 6 9% and 6 5% for m a and m s 9 , respectively. The accuracy of upper-layer thickness ~ l ! estimation is as good as 6 6% when optical properties of upper and lower layers are known. Optical property and layer thickness prediction accuracy degrade signiﬁcantly when more than three free parameters are extracted from data ﬁts. Problems with convergence are encountered when all ﬁve free parameters ~m a and m s 9 of upper and lower layers and thickness l ! must be deduced.


Introduction
Visible and near-infrared diffuse reflectance spectroscopy is often used for noninvasive, in vivo characterization of tissue optical properties. [1][2][3][4] Tissue optical properties, namely, absorption ͑ a ͒ and reduced scattering ͑ s Ј͒, are sensitive to the concentration of light-absorbing molecules and light-scattering structures. The primary tissue contributors to nearinfrared light absorption are generally assumed to be hemoglobin ͑both oxy and deoxy forms͒, water, fat, cytochromes, and melanin. 5,6 Important light scatterers are tissue structural elements that have a dimension comparable to the optical wavelength ͑e.g., ϳ0.6 -1 m͒, such as cellular components and fibrous materials in the extracellular matrix ͑e.g., collagen and elastin͒. Light-absorbing and light-scattering structures, in turn, are markers of tissue function and structure, such as delivery and utilization of ox-ygen, distribution of blood and water in tissue, as well as density of cells and collagen. 5,[7][8][9] The sensitivity of optical properties to tissue structural and biochemical elements may provide a means to assess physiologic processes such as malignant transformation, inflammation, or response to thermal injury. However, these events may alter bulk tissue optical properties minimally or induce substantial changes in highly localized regions. Consequently, accurate quantification of the magnitude and spatial extent of optical property variations is essential for sensitively observing and characterizing disease-related functional and structural changes.
Presently, most methods that are used to quantify tissue optical properties are based on the diffusion approximation to the Boltzmann transport equation. 1,10,11 Diffusion theory of light transport in tissue treats photons as particles, undergoing random elastic collisions between scattering sites. Net transport of photons is driven primarily by a photon density gradient. To obtain an explicit solution to the diffusion equation, the medium is often assumed to be homogeneous and semi-infinite in extent. [12][13][14] Measurements are usually performed in a diffuse reflectance configuration, with the source and detector oriented normal to the surface and separated by a distance .
Solutions for semi-infinite geometry in a reflectance configuration may be adequate to characterize tissues that are deep and relatively homogeneous, even though most tissues have well-demarcated layers with different optical properties. A case in point is the skin covering the deep muscles of the limbs or the glandular tissues of the breast. Effects of overlying tissue layers on measurements may be minimal, provided that the top layer is thin, i.e., l Ͻ 3 mm, the source detector separation is Ͼ Ͼl, and its properties are similar to the bottom layer. 15 As the top layer becomes thicker, use of a semi-infinite homogeneous model to regress data collected from layered tissues introduces uncertainties in the calculated optical properties. Typically, values obtained with semi-infinite model functions do not agree with optical values from either layer, 16,17 but tend to represent fractional contributions of both layers. This has an overall effect of diluting the intrinsic optical contrast between layers. Because optical contrast is an important feature that often signifies pathologic transformations and occurs in selective drug localization, 18 -20 precise measurement of optical contrast in tissue layers would enhance the photon migration diagnostic value.
Model accounting for layered heterogeneity may accurately preserve optical contrast originating from pathologic, injurious, and pharmacological processes that predominantly affect superficial tissues. In principle, this would allow selective characterization of a , s Ј, and the thickness of each layer from a single measurement set. Several investigators have examined layered diffusion-based models to account for tissue heterogeneity. [15][16][17]21,22 Schmitt et al., 21 Dayan et al., 22 and Kienle et al. 17 have presented analytical models for photon diffusion in two-layer turbid media. These approaches employ point or pencil sources that require one to solve the diffusion equation in three-dimensional space. Their results suggest that spatially resolved steady-state reflectance data are not sufficient to accurately calculate optical properties and thickness of layered media. Additional information, such as the phase delay from the frequency-domain or time decay profile from time-domain measurements, may be needed for complete characterization of layered structures.
We recently developed a diffusion-based model that assumes a planar source irradiating two-and threelayer turbid media. 6,23 By assuming a planar source, we can convert the three-dimensional partial differential equations into relatively simple, ordinary differential equations because all equations depend only on the z-spatial ͑depth͒ coordinate. A plane source can be assumed when the beam diameter is much greater than the light penetration depth of the media. This requirement to approximate a plane source imposes a limit on lateral resolution. However, the harmonically varying plane-wave model has several key advantages: ͑1͒ compared with spatially resolved reflectance, measurements are less susceptible to local variations on the tissue surface and thus yield more precise global values, ͑2͒ measurements can be performed in a noncontact reflectance mode, ͑3͒ signals at high modulation frequencies are highly sensitive to superficial tissues, ͑4͒ nonlinear fitting of experimental data can be performed rapidly because the solutions are analytic, and ͑5͒ solutions to the diffusion approximation with a plane-wave source can be readily extended to media with three or more layers.
In this study we use a plane-wave source and the corresponding analytical model to quantify optical properties of two-layer media. Frequency-domain phase and amplitude of diffuse photon density waves in two-layer phantoms are measured for modulation frequencies ranging from 10 to 1500 MHz. We employ tissue phantoms with well-defined optical properties and upper-layer thickness to evaluate the feasibility of using model functions to characterize two-layer media. Our goals of this study are to ͑1͒ assess the feasibility and accuracy of using frequency-domain planar photon density waves ͑PPDW's͒ to quantify optical properties and upperlayer thickness of two-layer turbid media and ͑2͒ determine the maximum number of parameters ͑ a1 , s1 Ј, a2 , s2 Ј, and l ͒ that can be extracted from frequency-domain amplitude and phase data with acceptable accuracy.

Theory
Diffusion equation solutions for light propagation in two-layer media have been presented by Svaasand et al. 6,23 The special case of PPDW propagation is shown in Fig. 1. Turbid structures consist of a planar layer of thickness l on top of a semi-infinite support. The upper-layer surface is bounded by air. The coordinate ͑ϩz͒ denotes direction normal to the air-media surface, pointing into the media. Coordinates ͑ x, y͒ are parallel to the surface. Light absorption ͑ a ͒ and reduced scattering ͓ s Ј ϭ ͑1 Ϫ g͒ s ͔ of the top and bottom layers are denoted as ͑ a1 , s1 Ј͒ and ͑ a2 , s2 Ј͒, respectively. The g coefficient is the Fig. 1. Schematic depicts the parameters of two-layer turbid media. The top layer is characterized by the absorption coefficient ͑ a1 ͒, reduced scattering coefficient ͑ s1 Ј͒, and thickness ͑l ͒. a2 and s2 Ј denote, respectively, the absorption and reduced scattering coefficients of the bottom layer. The refractive indices of both layers are n, and air is n 0 ϭ 1. The light source is an intensitymodulated plane wave. The source strength decays exponentially at rates that depend on layer optical properties. average cosine of the scattering angle. The layers are assumed to have the same refractive index n, giving a light speed of c ϭ ͓3 ϫ 10 11 ͑͞n mm͒s͔.
The time-dependent diffusion equation for light transport in turbid media is given by 10 -12 where is the light fluence rate ͑W mm Ϫ2 ͒, tr ϭ a ϩ s Ј ͑mm Ϫ1 ͒ is the transport coefficient, and q ͑W mm Ϫ3 ͒ is the source. For a plane source shown in Fig.  1, the time-dependent diffusion equation ͓Eq. ͑1͔͒ is effectively reduced to an equation with only z-spatial dependence, i.e., the fluence rate is independent of the ͑x, y͒ coordinates and is solely a function of z. Furthermore, for the sinusoidal intensity-modulated source at angular frequency , the fluence rate varies at the same frequency, having the form ϳ exp͑it͒. With a harmonically varying planar source, Eq. ͑1͒ can be rewritten as where ␦ c and c are, respectively, the effective complex penetration depth and the complex diffusion coefficient. These complex coefficients, ␦ c and c , are defined by The complex diffusion coefficient relates the diffuse flux j to the fluence rate , given by Fick's law: Inside the media, the strength of the coherent light source distribution is assumed to decay exponentially and is given by where P 0 is the initial source power transmitted through the media and the subscripts 1 and 2 refer to properties of the top and bottom layers, respectively. Boundary conditions at the layer-layer and mediaair interfaces are used to determine the unique solu-tion to Eq. ͑1͒. At the layer-layer interface, the fluence rate and diffuse flux j are continuous: At the media-air interface, the partial current boundary condition relates the diffuse flux to the fluence rate immediately below the surface and is given by 12 where R eff is the effective reflection coefficient, determined when the Fresnel reflection coefficient for unpolarized light is integrated over the hemisphere. Haskell et al. 12 have calculated R eff values of 0.431 and 0.493 for media with refractive indices of 1.33 and 1.40, respectively. For the source specified in Eqs. ͑6͒, the general solutions to Eq. ͑1͒ for each region are given by where the constants 1 and 2 are defined as To determine the three unknown constants A 1 , A 2 , and A 3 , the general solutions given by Eqs. ͑9͒ are substituted into the three boundary conditions ͓Eqs. ͑7͒ and ͑8͔͒ to yield a linear system of three equations and three unknowns ͑ A 1 , A 2 , and A 3 ͒. The linear system of equations can be expressed in matrix form as follows: where the ␣ ij and ␤ i terms are given by where ␦ c,i , c,i , i , and are defined in Eqs. ͑3͒, ͑4͒, ͑8͒, and ͑10͒, respectively, and the subscripts 1 and 2 specify the top or bottom layer, respectively. Solving Eq. ͑11͒ for the vector A, we obtain the unknown constants A 1 , A 2 , and A 3 . The matrix approach has the advantage that it can be readily extended to media of three or more layers. For example, media with three layers will have five unknown constants and five boundary conditions, giving a linear system of five unknowns and five equations. The diffuse reflection coefficient, defined as the ratio of the diffusely reflected flux to the incident flux, is a measureable quantity given by The absolute value ͉␥͉ and phase ␥ of the diffuse reflection coefficient yield, respectively, the amplitude attenuation and phase shift of the diffusely reflected flux relative to the incident flux. Amplitude attenuation and phase shift of diffusely reflected light measured at a range of modulation frequencies are simultaneously fit to the amplitude and phase functions of Eq. ͑14͒ to extract upper-layer thickness and optical properties of two-layer media.

A. Frequency-Domain Instrument
A broadband, frequency-domain photon migration ͑FDPM͒ instrument is used to measure the phase and amplitude of photon density waves for 201 frequencies over the range of 10 -1500 MHz. The FDPM instrument is depicted schematically in Fig. 2͑a͒ Fig. 2. Schematic outline of components and construction of the frequency-domain instrument that was used to measure amplitude and phase of photon density waves for frequencies ranging from 10 to 1500 MHz. ͑b͒ A two-compartment container was used to hold the top and bottom layer liquid phantoms; layers were separated by a thin, transparent polyvinyl film. Light emitted from the 100-m source fiber was expanded to form a 5-cm-diameter illumination spot. GPIB, general-purpose interface bus; APD, avalanche photodiode.
laser diode to produce sinusoidal intensitymodulated light. Light emitted from the laser diode is coupled into a graded-index, multimode optical fiber with a 22-deg acceptance angle ͑Oz Optics, Ltd., Ontario, Canada͒.
To generate plane waves, sinusoidal intensitymodulated light emitted from the source fiber is expanded and collimated to form a uniform 5-cm spot. The beam is used to irradiate the turbid two-layer media. A detecting fiber ͑FT-1.0-EMT, Thorlab, Inc., Newton, N.J.͒ is slanted at 25 deg relative to the surface normal ͓as shown in Fig. 2͑b͔͒ to selectively collect diffusely reflected light and reject specularly reflected light. Because a detecting fiber of 1 mm in diameter was placed in the path of the incident light, it cast a thin shadow on the 5-cm irradiation spot. However, the inhomogeneity from the shadow was small and had a negligible effect on the overall uniformity of the planar source.
Diffusely reflected light collected by the detecting fiber is focused onto the avalanche photodiode. The rf output of the avalanche photodiode is directed to the network analyzer receiver. The receiver heterodynes the rf signal with an internal reference to determine the phase and amplitude of the test measurement. FDPM instrument phase and amplitude uncertainties are less than Ϯ0.30 deg and Ϯ2.5%, respectively.

B. Two-Layer Phantoms
Two-layer phantoms were constructed from a nearinfrared absorbing dye ͑water-soluble nigrosin, Aldrich Chemical Company Inc., Milwaukee, Wis.͒ and light-scattering Intralipid ͑Intralipid 20%, Pharmacia, Inc., Clayton, N.C.͒. Optical properties of the top and bottom layers can be independently and precisely varied when the dye and Intralipid concentrations are altered. A specially designed container, Fig. 2͑b͒, was used to hold the two-layer liquid phantoms. The container had two compartments, with the top and bottom components separated by 50-mthick transparent polyvinyl film ͑Reynolds Film, Reynolds Metals Co., Richmond, Va.͒. Liquid-dye and Intralipid mixtures were placed in the bottom and top compartments to form two-layer media. Optical properties of the liquid phantoms ͑i.e., dye-Intralipid mixtures͒ were calculated from the known concentrations of the dye and Intralipid 24 and were determined independently from frequency-domain data performed in infinite geometry. Techniques for quantifying optical properties of turbid, infinite media are thoroughly described elsewhere. 11,25,26 The thickness of the top layer could be varied continuously and was determined with an accuracy of Ϯ10 m by use of a micrometer. Concentrations of dye and Intralipid were precisely varied to set the optical properties of each layer.
Five sets of data were collected. In each set, one parameter ͑ a1 , s1 Ј, a2 , s2 Ј, or l ͒ was varied whereas the other four parameters remained constant. To change the thickness of the top layer, the volume of the liquid phantom in the upper compartment was varied. To change the absorption of the top layer ͑ a1 ͒, a known volume ͑e.g., 1 Ϯ 0.01 ml͒ of the top layer was removed. We replaced the removed volume with an equal amount of liquid phantom that had the same optical scattering but significantly higher absorption ͑25ϫ͒. Using this procedure, we increased the absorption of the top layer while the scattering remained constant. The magnitude of this absorption increase can be calculated from the volumes and concentrations of ͑1͒ the top layer and ͑2͒ the phantoms used in the exchange. Scattering of the top layer was similarly incremented. In this case, the liquid phantom used in the exchange had the same absorption but significantly higher scattering ͑10ϫ͒ as compared with the optical properties of the top layer. Optical properties of the bottom layer were varied by similar procedures. Table 1 summarizes the optical properties and thickness of the two-layer media for the five data sets. Parameters that were varied and their ranges are listed for each data set.
To investigate the potential influence of the thin polyvinyl film, we recorded phase and amplitude data from two-layer phantoms that had identical optical properties in the top and bottom layers. The toplayer thickness was varied from 0 to 10 mm. FDPM phase and amplitude data obtained from the layered phantoms were compared with corresponding data from semi-infinite media with the same optical properties.
Phase and amplitude of PPDW's were measured on two-layer media for frequencies spanning 10 -1500 MHz. FDPM measurements were performed on semi-infinite ͑i.e., one-layer͒ phantoms of known optical properties at the beginning of each experiment so as to calibrate the FDPM instrument response. The calibration phantom had optical properties of a ϭ 0.010 mm Ϫ1 and s Ј ϭ 0.620 mm Ϫ1 .

C. Data Fitting and Analysis
In preparation for data fitting, raw phase and amplitude data were corrected for the FDPM instrument response, determined from calibration measurements. Instrument phase ͑ ir ͒ and amplitude ͑ A ir ͒ response were calculated as follows: where cal and A cal are, respectively, the phase and amplitude of calibration measurements, and ␥ and ͉␥͉ are, respectively, the theoretical phase and amplitude response calculated with Eq. ͑14͒ for the known optical properties of the calibration phantom. Calibration-corrected phase and amplitude data were simultaneously fit to the phase and amplitude of the diffuse reflectance functions ͓ ␥ and ͉␥͉ of Eq. ͑14͔͒ to extract properties of two-layer media. To perform simultaneous fitting of calibration-corrected phase and amplitude data, a Nelder-Mead ͑simplex͒ algorithm was used to minimize the least-squares merit function 27-29 : where i and A i are calibration-corrected phase and amplitude data, i and ͉␥͉,i are phase and amplitude uncertainties ͑standard deviations͒, ␥ and ͉␥͉ are phase and amplitude model functions ͓Eq. ͑14͔͒, and the subscript i denotes data at modulation frequency i . One or more parameters ͑ a1 , s1 Ј, a2 , s2 Ј, or l ͒ were iterated until a 2 minimum was found. The remaining variables were fixed to real values. The fit algorithm was terminated when one of the following conditions occurred: ͑1͒ the changes in 2 value and fit variable reached their tolerance of 1 and 0.5, respectively; or ͑2͒ the number of iterations exceeded 500 times the number of fitting parameters. Fits were performed with different initial guess values to avoid local minima and assess convergence. The simplex algorithm was computationally intensive and relatively slow, but was selected because it is simple to implement and is robust in finding minima.
Four fit strategies were performed on the data sets to extract various combinations of parameters. In the first case, only one parameter was fit ͑i.e., free floating͒ whereas the other four were fixed. The fit parameter was selected to be the same variable that was experimentally altered. For example, only thickness was fitted for data set 1, where thickness was experimentally varied. The other four properties were set to real values. In the second case, two parameters were fit whereas the remaining three properties were set to real values. One of the two fit variables was selected to be the same variable that was experimentally altered. All possible combinations of two-variable fits were performed for each data set, e.g., ͑l, a1 ͒, ͑l, s1 Ј͒, ͑l, s2 Ј͒, and so on. In the third case, three parameters were fit whereas the remaining two properties were set to true values. Again, one of the three fit variables was selected to be the same variable that was experimentally altered, and all possible combinations of three-variable fits were performed for each data set. In the fourth case, all five parameters ͑ a1 , s1 Ј, a2 , s2 Ј, and l ͒ were sought when each was allowed to float free.
To evaluate the feasibility of using the PPDW model to accurately quantify properties of two-layer media, we computed the absolute and percent difference between fit and real values. The percent difference was defined as ͑real-fit͒͑͞real͒ ϫ 100%. Absolute and percent difference data for each parameter were pooled together from each fit case, e.g., difference a data from the two-variable fit combinations were pooled from data sets 2 and 4 where a was experimentally varied. In the interest of clarity, a1 and a2 were considered jointly in the accuracy calculation, as were s1 Ј and s2 Ј. The standard deviations of pooled difference data for the various parameters, i.e., l, a , s Ј, were calculated for one-, two-, three-, and five-variable fits. These standard deviations were used to assess accuracy when we quantified optical properties.

Results
FDPM data acquired from two-layer phantoms with identical top-and bottom-layer optical properties were compared with corresponding semi-infinite homogeneous samples. Differences between layered and homogeneous samples were within our normal phase and amplitude measurement uncertainties, demonstrating that the thin transparent film did not significantly perturb two-layer FDPM data.
The graphs in Fig. 3 compare the measured phase ͑a͒ and amplitude ͑c͒ of PPDW's to the simulated phase ͑b͒ and amplitude ͑d͒. Measured phase and amplitude were performed on two-layer media with properties described by data set 1, where thickness was varied from 0.00 to 10.40 mm, absorption of both layers was the same, and scattering of the bottom layer was twice that of the top layer ͑ a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , a2 ϭ 0.0056 mm Ϫ1 and, s2 Ј ϭ 0.630 mm Ϫ1 ͒. We calculated the simulated data by substituting the real property values of twolayer media into the diffuse reflection coefficient of Eq. ͑14͒. Experimental data agreed well with simulated data, which predicted an increase in phase and amplitude as the less-turbid top layer becomes thicker. Effects of thickness on phase and amplitude were minimal when thickness exceeded 0.5␦ dc of the top layer ͑ϳ7 mm under this experimental con-dition͒, indicating that a significant fraction of detected density waves originates from the superficial region. ␦ dc is the dc penetration depth, i.e., the complex penetration depth ␦ c at frequency ϭ 0. ␦ c is decreased with modulation frequency, signifying that high-frequency PPDW's selectively probe superficial regions and are minimally perturbed by deep structures ͑ϳ␦ c ͒. An example of simultaneous phase and amplitude fits to two variables is shown in Fig. 4; the fit variables were l and s1 Ј. Fit and experimental data are in excellent agreement, with typical 2 values between 1 and 2 per degree of freedom.

Case 1: Fitting One Variable
Phase and amplitude data from the sets presented in Table 1 were simultaneously fit to model functions to calculate the experimentally varied parameter, e.g., data from set 2 were used to quantify a1 because it was varied in the experiment. Figure 5 plots the fit-derived versus real values for the five sets: ͑a͒ thickness l, ͑b͒ a1 , ͑c͒ s1 Ј, ͑d͒ a2 , and ͑e͒ s2 Ј. Fit values for different initial guesses are denoted by the symbols ࡗ and OE, and the solid lines are for expected values. Fit-derived and expected values were in excellent agreement when only one parameter was free floating. Furthermore, fit variables converged to the same value regardless of the initial guess values. Absolute and percent differences between fit and real values were calculated for each data set presented in Table 1. Data within the ranges of l ͑0.5-7.0 mm͒, a ͑0.005-0.036 mm Ϫ1 ͒, and s Ј ͑0.3-2.0 mm Ϫ1 ͒ were pooled together. Note that these ranges were subsets of the experimental ranges and were selected Fig. 3. Plots show response of PPDW phase and amplitude as the top-layer thickness ͑l ͒ is increased from 0.00 to 10.40 mm. Optical properties of the layers were a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , a2 ϭ 0.0056 mm Ϫ1 , and s2 Ј ϭ 0.630 mm Ϫ1 . ͑a͒ Measured phase is compared with ͑b͒ simulated phase, and ͑c͒ measured amplitude is compared with ͑d͒ simulated amplitude. Simulated phase and amplitude were calculated from known properties of the two-layer media. based on the following observations. At low values of l and a , the percent errors relative to the real values were high, even though the absolute errors were small. The percent errors for thickness were also large when l was greater than 7 mm. This is especially true for the multiple variable fit, but less so for the single variable fit. In the interest of comparison, however, the same ranges are selected for each of the three fitting cases. Standard deviations of the absolute and percent difference data were determined and summarized in column 1 of Table 2 for the specified ranges.

Case 2: Fitting Two Variables
Phase and amplitude from the five data sets were fit to quantify two properties of layered media. All possible combinations of two-variable fits were performed for each experimental set, with one of the two fit variables set to the experimentally varied parameter. A representative result of two-variable fitting is shown in Fig. 6, where phase and amplitude data from set 5 were analyzed to determine a2 and s2 Ј. In ͑a͒, fit-derived s2 Ј values for two different initial guesses ͑ࡗ, OE͒ are plotted against real s2 Ј values, which were experimentally varied. Fit-derived ͑ࡗ, OE͒ and expected ͑solid line͒ a2 values are compared with real s2 Ј values in ͑b͒. Recovered a2 and s2 Ј values are in good agreement with true properties of the bottom layer. Significantly, the calculated s2 Ј values accurately track experimental variations in the scattering parameter, whereas the fit-derived a2 values remain constant. Results demonstrate that properties of two-layer media can be independently and sensitively recovered. Fit values reached convergence for different initial guess values, analogous to the convergence result seen for one-variable fits. Absolute and percent differences between fit and real values were calculated for two-variable fits. Difference data from the all-experimental sets that were within the specified ranges were pooled together, and their standard deviations were calculated for each parameter ͑l, a , s Ј͒. In column 2 of Table 2 we summarize the uncertainties in quantifying two variables from phase and amplitude data for the specified ranges.

Case 3: Fitting Three Variables
FDPM data from the five experimental sets were fit to quantify three parameters of two-layer media. All possible combinations of three-variable fits were performed for each set, e.g., ͑l, a1 , s1 Ј͒, ͑l, a2 , s2 Ј͒, and so on. One of the three fit variables corresponded to the parameter that was experimentally varied. Figure 7 is a representative example of a three-variable ͑l, a1 , s1 Ј͒ fit for data from experimental set 1, where thickness was varied. Fit-derived values for l, a1 , and s1 Ј are plotted against real thickness in ͑a͒, ͑b͒, and ͑c͒, respectively. Fit values deviated significantly from expected values and often failed to converge for different initial guesses of the fit parameters. Convergence failures were due to premature termination of the fitting algorithm when the maximum allowable iterations were reached. Significantly, fit values converged for different initial guesses when the algorithm was allowed to have unlimited iterations, suggesting that 2 space was relatively flat or had local minima.
Absolute and percent differences between fit and real values were calculated for all three-variable fits in each set. Difference data from all sets that were within the specified ranges were pooled together, and standard deviations were determined for each parameter ͑l, a , s Ј͒. In column 3 of Table 2 we summarize the uncertainties in quantifying the properties of two-layer media for three-variable fits. Uncertainty values are applicable for the specified ranges. Fig. 4. Example showing agreement between experimental data and results of the fit: ͑a͒ experimental phase ͑s͒ ͑solid curves͒; ͑b͒ experimental amplitude ͑s͒. Experimental phase and amplitude are from data set 1, and the properties of the two-layer medium were l ϭ 1.55 mm, a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , a2 ϭ 0.0056 mm Ϫ1 , and s2 Ј ϭ 0.630 mm Ϫ1 . Fit data are from use of the simplex ͑Nelder-Mead͒ algorithm to simultaneously fit phase and amplitude data to extract l and s1 Ј.

Case 4: Fitting Five Variables
We attempted to extract five properties of two-layer media from phase and amplitude data. Figure 8 is a representative result for five-variable fitting of data set 4, where a2 was experimentally varied. In the interest of clarity, only one initial guess is shown for ͑l, a1 , s1 Ј, and s2 Ј͒ in ͑b͒ and ͑c͒. With fivevariable fits, the differences between fit and real values were unacceptably large, typically greater than 200%. Moreover, fit values failed to converge for different initial guesses even when the maximum number of iterations was unlimited, suggesting that the inversion problem for five-variable fitting of phase and amplitude data does not have a unique solution.

Discussion
In principle, frequency-domain PPDW phase and amplitude data can be used to quantify all five properties ͑ a1 , s1 Ј, a2 , s2 Ј, and l ͒ of two-layer media because Fig. 5. PPDW phase and amplitude data were acquired from two-layer media with properties specified in Table 1. Phase and amplitude measurements were fit to model functions ͓Eq. ͑14͔͒ to estimate one variable. Only the experimentally varied property was calculated, whereas other parameters were set to real values. The true value for each data set is specified on the x axis of the subplots. Results of one-variable fits are shown in subplots ͑a͒ through ͑e͒ for the five data sets of Table 1 ͓͑a͒ through ͑e͒ correspond to sets 1 through 5, respectively͔. Fit values for different initial guesses are specified on the subplots as initial guess 1 ͑ࡗ͒ and initial guess 2 ͑OE͒, and solid lines indicate expected values. the function for the diffuse reflection coefficient ͓Eq. ͑14͔͒ is unique. Experimentally, we were able to quantify at most three variables with acceptable accuracy. Although the standard for acceptable accuracy is application specific, an uncertainty of Ϯ10% may be reasonable for applications requiring quantitative assessments, whereas an uncertainty Ϯ20% is more appropriate for qualitative assessments.
For one-and two-variable fits, optical properties and thickness were measured with percent accuracies of the order of Ϯ a ϭ Ϯ10% ͑Ϯ0.002 mm Ϫ1 ͒, Ϯ s ϭ Ϯ5% ͑Ϯ0.060 mm Ϫ1 ͒, and Ϯl ϭ Ϯ12% ͑Ϯ0.5 mm͒ within the ranges specified in Table 2. In media where the top-layer thickness was greater than 50% of the dc penetration depth ͑␦ dc,1 ͒, fit-derived values for thickness l and bottom-layer optical properties a2 and s2 Ј were less accurate ͑versus media with thinner upper layers͒. An example of this effect is shown in Fig. 7. Fit values for optical properties of the bottom layer and thickness, in particular, deteriorated significantly as the top-layer thickness increased.
We can explain the deterioration in the accuracy of estimating a2 , s2 Ј, and l with thicker top layers in terms of fractional signal contributions. As the top layer extends into the sample, diffusely reflected light originating from the upper regions constitutes a larger fraction of the total detected signal. Consequently, when l Ն 0.5␦ dc of the upper layer, it is easier and more robust to quantify a1 and s1 Ј from FDPM data than it is to estimate l, a2 , and s2 Ј. However, for many layered biological structures such as skin and subcutaneous tissues, l is less than 0.5␦ dc , making estimates of both upper and lower properties tractable in real tissues.
Percent differences between the recovered and the actual values systematically increase with more freefloating parameters, as can be observed from Table 2. When the number of fit variables was increased to two, the percent uncertainties for l and a were two times higher than the uncertainties for one-variable fits. For three free-floating variables, uncertainties for l and a were, respectively, three and four times higher. Interestingly, the uncertainty for s Ј increased minimally with additional free-floating variables. This is likely due to the fact that highfrequency photon density waves are heavily dominated by scattering. A small change in s Ј dramatically perturbs the PPDW phase and ampli- Standard deviations were computed from experimental data that were within the following ranges: l ͑0.5-7.0 mm͒, a ͑0.005-0.036 mm Ϫ1 ͒, and s Ј ͑0.3-2.0 mm Ϫ1 ͒. Percent standard deviations were rounded to the nearest integer. Fig. 6. Representative results of two-variable fits for data set 5 are plotted in ͑a͒ and ͑b͒. Reduced scattering of the bottom layer ͑ s2 Ј͒ was experimentally increased, whereas the other properties were held constant ͑l ϭ 2.35 mm, a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , and a2 ϭ 0.0054 mm Ϫ1 ͒. FDPM amplitude and phase data were simultaneously fit to model functions to calculate optical properties ͑ a2 and s2 Ј͒ of the bottom layer, whereas the remaining parameters ͑ a1 , s1 Ј, and l ͒ were assigned to real values. ͑a͒ Fit s2 Ј values for two different initial guesses are plotted ͑ࡗ, OE͒ and expected values are shown as a solid line. ͑b͒ Fit a2 values ͑ࡗ, OE͒ and expected a2 values ͑solid line͒. a2 was kept constant in the experiment, and thus expected values appear as a horizontal line in ͑b͒.
tude as compared with equivalent variations in a or l; thus prediction of s Ј is more robust.
Practically, it was not feasible to fit four or more parameters because of the uncertainties in the phase and amplitude data. The inverse problem ͑i.e., data fitting͒ does not yield a unique solution when four or more variables are left freely floating. Given the finite measurement uncertainties and fit termination criteria, the number of potential solutions satisfying the termination criteria grows dramatically with more free variables. 30 Specifically, numerous combinations of the fit parameters ͑ a1 , s1 Ј, a2 , s2 Ј, and l ͒ yield similar phase and amplitude within the bounds of measurement errors. This is especially true for variables that minimally perturb the phase and amplitude, i.e., when the partial derivative of the phase and amplitude with respect to the variable is insignificant. Under these conditions, calculated values are likely to deviate from the true solution and will fail to converge regardless of the initial guess. Nevertheless, simultaneous estimation of five variables from FDPM data may be feasible with improvements in the accuracy of phase and amplitude measurements.
Alternatively, frequency-domain PPDW's could be used in conjunction with standard diffuse reflectance approaches that use point-source detector pairs to estimate all five parameters. For example, FDPM methods employing large source detector separations ͑Ͼ Ͼl ͒ have been shown to provide good estimates of lower-layer optical properties. 15 Thus complementary point-source measurements could yield a priori knowledge that helps constrain model fits and produce better estimates of l, a1 , and s1 Ј.
We performed simulations to determine the necessary instrument accuracy to fit five variables while maintaining the same level of accuracy seen for twovariable fits. Equation ͑14͒ was used to generate phase and amplitude data of two-layer media with the same optical properties listed in Table 1. Normally distributed noise of different magnitude was added to the data. Noise-containing phase and amplitude data were fit to model functions and the five variables were calculated. Results of these simulations show that the instrument phase and amplitude uncertainties should be approximately 0.005 deg and 0.01%, respectively, to fit five variables with comparable accuracy that is seen for two-variable fits. Even though these uncertainty values are presently unattainable for our broad-bandwidth device, dedicated instruments specified for discrete frequencies may eventually be able to achieve this level of performance. Alternatively, instruments that are capable of measuring phase and amplitude at higher Fig. 7. Representative results of three-variable fit for data set 1 are shown in ͑a͒-͑c͒. Thickness l was experimentally increased whereas the other four properties were kept constant ͑ a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , a2 ϭ 0.0056 mm Ϫ1 and s2 Ј ϭ 0.630 mm Ϫ1 ͒. FDPM amplitude and phase data were simultaneously fit to model functions to calculate three parameters ͑ a1 , s1 Ј, and l ͒ of the top layer. a2 and s2 Ј were set to real values. ͑a͒ Fit values of l for two different initial guesses ͑ࡗ, OE͒ and expected l values ͑solid line͒. ͑b͒ Fit a1 values with different initial guesses ͑ࡗ, OE͒ and expected a1 values ͑solid line͒. ͑c͒ Fit s1 Ј values for different initial guesses ͑ࡗ, OE͒ and expected s1 Ј values ͑solid line͒. frequencies, i.e., Ͼ2 GHz, may also improve the accuracy of recovering all five free parameters in twolayer systems. This is due to the increased scattering sensitivity and better spatial localization of high-frequency PPDW's. Thus, under certain conditions, the shallow penetrance of higher frequencies could be used to improve estimates of upper-layer properties, which, in turn, would enhance overall fitting performance.

Conclusions
We have demonstrated that a diffusion-based PPDW model can be used to describe light propagation in two-layer media. We have confirmed experimentally that three properties of two-layer media can be estimated from frequency-domain phase and amplitude measurements of PPDW's. For values within the specified ranges, the percent accuracy when l, a , and s Ј are estimated are, respectively, of the order of Ϯ19%, Ϯ21%, and Ϯ6% when three variables are fit. Accuracy of l and a estimates are significantly better for one-variable fits ͑Ϯl ϭ Ϯ6%, Ϯ a ϭ Ϯ5%͒. Furthermore, simulation studies suggest that all five properties can be estimated provided that the instrument phase and amplitude accuracies are significantly improved.
For media with optical properties comparable to tissues, 31 PPDW has a depth sensitivity in the 1-8-mm region ͑i.e., ϳ0.1-0.5␦ dc ͒. Structures at these depths are generally beyond the depth resolution limit of coherent methods ͑e.g., optical coherence tomography, confocal microscopy͒ and are too superficial for methods that employ large source detector separations. 14,15 Consequently, the PPDW approach is well suited to characterize macroscopic optical properties and layer thickness in regimes that are intermediate between those accessible with traditional high-resolution coherent and low-resolution diffuse techniques. Clinical applications that are amenable to PPDW measurements include measurements of optical properties, coagulation depth, and extent of necrosis following burn injury and the assessment of optical properties of epithelial structures during dysplastic and malignant transformation.  Representative results of five-variable fit for data set 4 are shown in ͑a͒-͑c͒. a2 was experimentally increased whereas the other four properties were kept constant ͑l ϭ 2.10 mm, a1 ϭ 0.0054 mm Ϫ1 , s1 Ј ϭ 0.310 mm Ϫ1 , and s2 Ј ϭ 0.630 mm Ϫ1 ͒. Amplitude and phase data were simultaneously fit to model functions to estimate all five parameters ͑ a1 , s1 Ј, a2 , s2 Ј, and l ͒. ͑a͒ Fit a2 values for different initial guesses ͑ࡗ, OE͒ and expected a2 values ͑solid line͒. ͑b͒ Thickness ͑l ͒ and s2 Ј. On the left y axis of ͑b͒, fit ͑ࡗ͒ and expected ͑solid line͒ thickness values are plotted against real a2 values. On the right y axis of ͑b͒, fit ͑s͒ and expected ͑dashed line͒ s2 Ј values are plotted versus real a2 values. ͑c͒ Properties of the top layer ͑ a1 and s1 Ј͒: fit ͑ࡗ͒ and expected ͑solid line͒ a1 on the left y axis; fit ͑s͒ and expected ͑dashed line͒ s1 Ј on the right y axis.