of the delta-P1 approximation for recovery of optical absorption, scattering, and asymmetry coefficients in turbid media.

We introduce a robust method to recover optical absorption, reduced scattering, and single-scattering asymmetry coefﬁcients (cid:1)(cid:2) a , (cid:2)(cid:3) s , g 1 (cid:4) of inﬁnite turbid media over a range of (cid:1)(cid:2)(cid:3) s (cid:1) (cid:2) a (cid:4) spanning 3 orders of magnitude. This is accomplished through the spatially resolved measurement of irradiance at source– detector separations spanning 0.25–8 transport mean free paths (cid:1) l * (cid:4) . These measurements are rapidly processed by a multistaged nonlinear optimization algorithm in which the measured irradiances are compared with predictions given by the (cid:5) - P 1 variant of the diffusion approximation to the Boltzmann transport equation. The ability of the (cid:5) - P 1 model to accurately describe radiative transport within media of arbitrary albedo and on spatial scales comparable to l * is the key element enabling the separation of g 1 from (cid:2)(cid:3) s . © 2004 Optical Society of America OCIS codes: 100.3190, 170.3660, 170.5280, 170.7050, 290.1990, 290.4210.


Introduction
An understanding of photon migration within turbid media provides the basis for a wide variety of diffuse optical spectroscopy and tomography techniques. 1 Such techniques are enabled through the resolution of an appropriate inverse problem. One common approach to resolving such inverse problems is to use a radiative transport model to provide predictions for measured data based on an assumed set of optical properties. Such a model is then used in conjunction with a nonlinear optimization algorithm to determine the set of optical absorption and scattering properties that results in a radiative transport prediction that best matches the experimental data. In photon migration applications, the most commonly used radiative transport approximation is the standard diffusion model ͑SDA͒. 2,3 Although use of the SDA has proven quite fruitful in this respect, its limitations are significant. Specifically, the accuracy of radiative transport predictions provided by the SDA is seriously compromised when the radiant field displays even a moderate degree of angular asymmetry. Thus appropriate use of the SDA is restricted to highly scattering media ͓͑Ј s ͞ a ͒ տ 10͔, over length scales much larger than the transport mean free path l* ϵ ͑ a ϩ Ј s ͒ Ϫ1 , and locations remote from collimated sources or interfaces possessing significant refractive index mismatch. 4 Several approaches have been developed to extend photon migration methods for optical property determination in tissues over a broad range of reduced albedo ͓aЈ ϵ Ј s ͑͞ a ϩ Ј s ͔͒ [5][6][7] and in situations where source-detector ͑s-d͒ separations comparable with l* are utilized. 2,8 -13 Such conditions represent situations for which the SDA provides inaccurate radiative transport predictions. Whereas diffusion-based approaches have been successful in determining a and Ј s for samples ͑Ј s ͞ a ͒ տ 10 ͑Refs. 4 and 5͒, such approaches will always be inadequate for more highly absorbing media. Chemometric and neural network approaches have been applied to media with optical properties in the range 0.2 Յ ͑Ј s ͞ a ͒ Յ 1000. 6,7,14 However, these approaches are not based on a phenomenological forward model of radiative transport and have provided mixed results. As a result, the conditions under which such methods fail are not always clear. The approach that we investigate here is the use of an improved analytic radiative transport model, the ␦-P 1 approximation, to enable optical property determination over a broadened range of optical properties that accommodates both highly scattering and highly absorbing media.
In previous work 15 we have shown that the ␦-P 1 approximation provides accurate radiative transport predictions in both highly absorbing and scattering turbid media, as well as on length scales approaching l*͞4. Analysis of turbid media on spatial scales comparable with a transport mean free path is of interest because photons detected at these locations have typically undergone only a limited number of scattering events. As a result, the character of the light field depends on not only the absorption and reduced scattering coefficients but also the moments of the singlescattering phase function, i.e., g 1 , g 2 , and so on. 10,16,17 Thus implementation of radiative transport models such as the ␦-P 1 approximation may provide the opportunity for recovery of not only optical absorption ͑ a ͒ and reduced scattering ͑Ј s ͒ coefficients but also specific additional descriptors of the single-scattering phase function. In this regard, recent research by Jones and Yamada, 8 Bevilacqua et al. 9,10 and Charvet et al. 13 has demonstrated that spatially resolved diffuse reflectance measurements of highly scattering media collected at s-d separations comparable to l* can be used in conjunction with Monte Carlo simulations to determine not only a and Ј s but also important parameters related to the first two moments of the single-scattering phase function. In this paper we investigate the use of analytic predictions provided by the ␦-P 1 approximation to the Boltzmann transport equation together with a multistaged gradient-based nonlinear optimization algorithm to recover a , Ј s , and g 1 in turbid media of arbitrary albedo from spatially resolved photon migration measurements.

A. Radiative Transport Model
As noted earlier, we have previously demonstrated the ability of the ␦-P 1 approximation to provide excellent radiative transport predictions within an infinite medium down to length scales approaching l*͞4 over a broad range of reduced albedo aЈ. 15 The key element enabling such predictions is the decomposition of the light field into collimated and diffuse components. The separate consideration of collimated ͑using a ␦ function͒ and diffuse ͑using a P 1 approxi-mation͒ light components allows this analytic treatment to be sensitive to both first and second moments of the single-scattering phase function, g 1 and g 2 . 18 Our focus is in providing predictions for the measured radial irradiance resulting from a spherical source of finite radius embedded within an infinite turbid medium, as depicted in Fig. 1. Here we review elements of the derivation of the ␦-P 1 approximation relevant to the development of our multistaged algorithm.
The basis of the ␦-P 1 approximation lies within the Boltzmann transport equation ͑BTE͒ that describes the angular photon flux or radiance L͑r, ͒, which represents the rate of photon arrival at position r traveling in direction per unit area and unit solid angle. The time independent form of the BTE is where t ϭ a ϩ s is the total attenuation coefficient, s is the scattering coefficient, a is the absorption coefficient, S is the volumetric source, and p͑Ј 3 ͒ is the single-scattering phase function that describes the probability of scattering from direction Ј to . Derivation of the governing equations within the ␦-P 1 approximation proceeds in a manner identical to the derivation of the SDA, with two important modifications. These modifications involve the addition of a collimated contribution, in the form of a ␦ function, to both the phase function and the radiance approximation. This is the key element that enables the separation of g 1 from Ј s . Specifically, the phase function p͑Ј 3 ͒ employed in the ␦-P 1 approximation is where f is the fraction of light scattered directly forward and g* is a scattering asymmetry coefficient. Matching the first and second moments of Eq. ͑2͒ to that of the Henyey-Greenstein phase function fixes the parameters f and g* as Note that in Eq. ͑2͒ passage to the limit f 3 0 reduces the ␦-P 1 phase function to the Eddington phase function approximation used in the SDA. Substitution of the ␦-P 1 phase function ͓Eq. ͑2͔͒ into the BTE ͓Eq. ͑1͔͒ results in an equation identical to the BTE with transformed optical parameters 18 where * s ϭ s ͑1 Ϫ f ͒, * t ϭ * s ϩ a , and p*͑Ј 3 ͒ ϭ ͑1͞4͓͒1 ϩ 3g*͑ ⅐ Ј͔͒. This transformed BTE demonstrates that the natural optical parameters in the ␦-P 1 approximation are not a , Ј s , and g 1 but rather a , * s , and g*. These three parameters provide a representation of the optical properties of the medium complete to the second-order angular moment and can be derived from the first-and secondorder similarity relations. 9,10 Thus the optimization algorithm described below employs this natural set of variables. Other papers that consider similar issues in radiative transport have instead used an equivalent parameter set, namely, a , Ј s , and ␥ ϵ ͑1 Ϫ g 2 ͒͑͞1 Ϫ g 1 ͒. 9,12,16 The second modification employed in the derivation of the ␦-P 1 approximation is the addition of a collimated term to the radiance approximation where d ͑r͒ is the diffuse fluence rate, j͑r͒ is the radiant flux, P͑r, ͒ is the irradiance provided by the collimated source, and 0 is the direction of collimated light emitted by the source. Removal of the ␦-function term in Eq. ͑5͒ reduces the radiance approximation to that used in the SDA. For an infinite geometry illuminated with a spherical source of finite size, as depicted in Fig. 1, the specific expressions for the diffuse fluence rate d ͑r͒ and radiant flux j͑r͒ within the ␦-P 1 approximation are given by 15 where In Eqs. ͑6͒-͑10͒, P 0 is the power of the optical source, r 0 is the radius of the optical source, E 1 ͑ z͒ is the exponential integral defined as ͐ z ϱ ͓exp͑Ϫt͒͞t͔dt, and eff is the effective attenuation coefficient ϵ ͓3 a ͑ a ϩ Ј s ͔͒ 1͞2 .
The addition of a collimated contribution to both the phase function and radiance approximations provides better accommodation for those photons that are highly forward scattered and allows the resulting radiance prediction to be more accurate in the description of highly asymmetric light fields. These collimated terms relieve a substantial burden from the first-order term in the Legendre series to fully accommodate the angular asymmetry of the light field and scattering phase function. The resulting ␦-P 1 approximation describes a transport process with an asymmetry coefficient g* that is lower than the original, g 1 , thereby providing a more accurate representation of the radiance in highly forwardscattering media.
Using the expressions for d ͑r͒ and j͑r͒ shown in Eqs. ͑6͒-͑10͒, we can predict the radial irradiance I r measured at a s-d separation r using an optical fiber with acceptance angle with where the collimated fluence rate c ͑r͒ is given by The acceptance angle is given by ϭ sin Ϫ1 ͑N͞n͒, where N is the numerical aperture of the optical fiber and n the refractive index of the medium. The collimated fluence rate c ͑r͒ is an element particular to the ␦-P 1 approximation and is not present in SDAbased models.

B. Measurements
The results shown utilize published irradiance data measured within an infinite intralipid phantom ͑n ϭ 1.33͒ at ϭ 674 nm at s-d separations from 0.25l* to 8l* by use of an optical fiber with numerical aperture N ϭ 0.37 oriented radially towards a spherical source of radius r 0 ϭ 0.4 mm. 15 The optical properties and corresponding values for ͑Ј s ͞ a ͒ for each of the five phantoms measured are shown in Table 1. The number of irradiance measurements taken from the five phantoms ranged from 90 to 106 with s-d separations spanning 2 to 45.4 mm.

C. Algorithm Development
Using predictions for the measured radial irradiance specified by Eqs. ͑6͒-͑12͒, we developed a procedure to determine the optical properties of the medium in which the experimental measurements were performed. To develop an algorithm that can be used for both strongly absorbing and scattering media we found it essential to obtain first a rough estimate for l*. With this estimate for l*, intervals of s-d separation are identified that are information-rich with respect to the particular optical parameter of interest. Specifically, measurements at s-d separations comparable with l* are dominated by collimated and minimally scattered light and contain important information relating to g 1 . By contrast, measurements at s-d separations significantly larger than l* are dominated by multiply scattered light. Such measurements should be minimally sensitive to g 1 but should instead contain information related to a and Ј s . 10 This fact suggests that a multistaged optimization algorithm that processes data in specific intervals of s-d separation may prove advantageous for the recovery of a , Ј s , and g 1 . Our development of these ideas resulted in a three-stage algorithm, the flow of which is depicted in Fig. 2.
In each stage of the algorithm, measurements from a specific interval of s-d separations are supplied to a constrained Levenberg-Marquardt ͑LM͒ algorithm ͑MATLAB, MathWorks Inc., Natick, Massachusetts͒ that minimizes a sum of squares, 2 , with respect to the 'natural' parameters of the ␦-P 1 approximation, i.e., where M is the number of measurements, I m ͑r i ͒ is the irradiance measured at location r ϭ r i , I p is the irradiance as predicted by Eqs. ͑6͒-͑12͒, and i is the standard deviation of the irradiance measurement at location r ϭ r i . Estimates for a , Ј s , and g 1 are then obtained via algebraic manipulation from the optimized parameters, a , * s , and g*. To ensure full sampling of the parameter space and to avoid convergence to local minima, we executed thirty trials of the LM algorithm using initial guesses of the optical properties selected randomly from the intervals: a ʦ ͓10 Ϫ4 , 10 Ϫ1 ͔ mm Ϫ1 , s ʦ ͓0.03, 3͔ mm Ϫ1 , and g 1 ʦ ͓0.6, 0.99͔. Since we typically consider highly forward-scattering media, we use this information to our advantage and constrain the single-scattering asymmetry factor to the interval g 1 ʦ ͓0.6, 0.99͔ when optimizing over all parameters. When determining the optimal value of g 1 in the final stage of the algorithm, we expand the interval to g 1 ʦ ͓0.0, 0.99͔. Of the converged trials, the recovered parameter set with the lowest 2 value is selected.
Stage 1 of the algorithm utilizes the entire measurement set consisting of 96 -106 distinct measurements taken at s-d separations in the range 2.0 -45.4 mm to estimate l*. Thus while the parameters a , * s , and g* are all determined by use of the 2 figure of merit, only the value for l* derived from these estimates is used to identify the interval of s-d  separations for which the measured data is supplied to Stage 2 of the algorithm.
In Stage 2, we use the estimate for l* obtained in Stage 1 to supply radial irradiance measurements acquired at s-d separations r Ն 1.5l* to the LM algorithm to determine a and Ј s . Only data distal from the source are considered in this stage, as measurements in this region result from multiply scattered photons. Thus these measurements should be insensitive to the value of g 1 but instead contain information critical to the extraction of a and Ј s . The values of a , * s , and g* obtained from this stage are used to determine and fix the values of a and Ј s . The value of l* computed from the estimates of a and Ј s is again used to identify the interval of s-d separations for which the measured data is supplied to Stage 3 of the algorithm.
In Stage 3, we determine g 1 by supplying the measured data at s-d separations r Ն 0.5l*, constrain Ј s and a to the values recovered in Stage 2, and perform a single parameter optimization for g 1 . The inclusion of irradiance measurements proximal to the laser source provides measurements at s-d separations dominated by collimated and minimally scattered light essential for the extraction of g 1 .

A. Algorithm Performance
The estimated values of a , Ј s , and l* recovered in Stage 1 of the algorithm are provided in Table 2. These results represent the optical properties that would be recovered by a single-stage inversion algorithm that processes all of the measured data and uses radiative transport predictions provided by the ␦-P 1 approximation. Errors exceeding 40% are seen in the estimates for a and Ј s , indicating that an inversion algorithm that uses the entire measurement set to recover all three coefficients simultaneously gives poor results and is of limited use. However, while the individual optical property estimates are poor, the estimates for l* computed from the individual optical properties are slightly better and accurate to within 28% for media characterized by 0.33 Շ ͑Ј s ͞ a ͒ Շ 300.
Using the l* estimate provided by Stage 1, we supply only irradiance data at s-d separations r Ն 1.5l* to the LM algorithm in Stage 2 in which we determine a and Ј s . The values for a and Ј s , along with the corresponding value for l* recovered in Stage 2, are shown in the left portion of Table 3. The results demonstrate that a , Ј s , and l* are recovered with relative errors no worse than Ϯ16%, Ϯ23%, and Ϯ13%, respectively.
Fixing these values for a and Ј s and providing irradiance data located at s-d separations r Ն 0.5l* enables the recovery of g 1 in Stage 3. These results are shown in the final two columns of Table 3 and reveal relative errors in g 1 that are no worse than Ϯ39%. Repetition of Stages 2 and 3 of the algorithm, making use of the slightly improved l* estimates derived in Stage 2, rather than those initially found in Stage 1, produces minimal alteration in the optical property estimates. Figure 3 displays the measured data for the turbid media with ͑Ј s ͞ a ͒ ϭ 290 and 3.1 along with the ␦-P 1 model predictions, corresponding to the recovered set of optical properties. Moreover, supplying the algorithm with the same set of irradiance data mixed with an additional 2% random error results in a negligible change in the recovered optical properties. This is a further indi- cator of the robustness of this algorithm to potential sources of noise in the measured data. The recovery of a and Ј s with errors no worse 16% and 23%, respectively, over a range of ͑Ј s ͞ a ͒ spanning nearly 3 orders of magnitude is an impressive result. By comparison, diffusion-based approaches have only demonstrated recovery of these optical properties over a range of ͑Ј s ͞ a ͒ spanning no more than 2 orders of magnitude. 4,5 Specifically, Fantini et al. 4 developed an inversion algorithm that processed frequency-domain photon migration measurements and used a SDA-based radiative transport model to determine the optical properties of semi-infinite turbid phantoms in the range 10 Ͻ ͑Ј s ͞ a ͒ Ͻ 615. The relative errors of the resulting a and Ј s estimates were reported to be no worse than Ϯ4% and Ϯ15%, respectively. Kienle and Patterson 5 used a similar radiative transport model and processed spatially resolved reflectance measurements using a continuous optical source to recover a and Ј s of turbid phantoms with properties in the range 13 Ͻ ͑Ј s ͞ a ͒ Ͻ 300. This approach resulted in estimates for a and Ј s that were no worse than Ϯ28% and Ϯ14%, respectively. Thus the results of the approach shown here that use the ␦-P 1 rather than the standard diffusion approximation compare favorably with previous studies using continuous optical sources, are able to span a much larger interval of ͑Ј s ͞ a ͒. It is important to note that SDA-based approaches are not capable of spanning a larger range of optical properties, as the radiative transport predictions are erroneous for media of high absorption. 4,15,16 Moreover, because the SDA provides identical predictions for media of identical Ј s independent of g 1 , recovery of g 1 is not tractable by use of SDA-based approaches.
It is also instructive to compare our results with studies using approaches that do not use a radiative transport model. For example, Kienle et al. 14 developed a neural network algorithm to analyze spatially resolved reflectance measurements on large phantoms with optical properties spanning 6 Ͻ ͑Ј s ͞ a ͒ Ͻ 904. They reported the recovery of a and Ј s with relative errors that were no worse than Ϯ23% and Ϯ5%, respectively. Pfefer et al. 7 have investigated a variety of partial least-squares ͑PLS͒ and neural-network ͑NN͒ approaches to determine optical properties from spatially resolved reflectance measurements of highly attenuating phantoms over the range of 0.2 Ͻ ͑Ј s ͞ a ͒ Ͻ 25. For the most accurate of the PLS approaches they reported that relative errors in a and Ј s can be no worse than Ϯ41% and Ϯ36%, respectively, while the most successful of the NN approaches result in relative errors of approximately Ϯ41% and Ϯ29%, respectively. However, for a given set of measurements, the most accurate PLS͞NN approach is not known a priori, and recovery of optical coefficients from phantoms with properties close to the borders of the ͑Ј s ͞ a ͒ range often exhibited errors on the order of 100 -300%. Again the results of the ␦-P 1 approach that we have presented compare quite favorably with these other studies and have validity in highly absorbing media extending to ͑Ј s ͞ a ͒ Ϸ 0.33.
Beyond the recovery of a and Ј s , our algorithm provides very promising results for the recovery of g 1 . Close examination of Table 3 suggests that the accuracy with which g 1 is recovered is correlated with ͑Ј s ͞ a ͒. For moderately and strongly scattering media, i.e., ͑Ј s ͞ a ͒ Ն 3, g 1 is recovered to within Ϯ11%. This result compares quite favorably with the studies of Jones and Yamada 8 and Bevilacqua and Depeursinge, 10 who employed Monte Carlo-based approaches to recover g 1 in highly scattering media with an accuracy of Շ15%. However, the recovery of g 1 in more strongly absorbing media does not fare as well and produces much larger errors of ϩ39% and Ϫ21% for ͑Ј s ͞ a ͒ ϭ 1.0 and 0.33, respectively. Our investigations reveal that the separation of g 1 from Ј s requires a data set containing measurements of collimated and minimally scattered photons close to the source in addition to multiply scattered photons distal from the source. While data proximal to the source is critical for the separation of g 1 from Ј s , the data distal from the source is equally critical for an accurate determination of Ј s . The successful determination of Ј s relies on measurements taken sufficiently far from the source so that the light field is essentially diffuse. Our studies of the light transport in this situation indicate that increases in absorption ͓i.e., reduction in ͑Ј s ͞ a ͔͒ displace the production of a diffuse light field dominated by multiply scattered photons to larger radial positions relative to l*. 20,21 As a result, the accurate determination of optical properties for media with low ͑Ј s ͞ a ͒ may require measurements at s-d separations larger than 8l*, especially for the accurate determination of g 1 .

B. Effect of Measurement Number
The number of data points required for accurate determination of the optical properties was studied also for both practical and algorithmic reasons. Accu- Fig. 3. Normalized irradiance versus s-d separation for samples with ͑ s Ј͞ a ͒ ϭ 290 and 3.1. Symbols represent measured data ͑N ϭ 0.37͒, whereas curves represent ␦-P 1 model prediction corresponding to the recovered set of optical properties. The normalized irradiance I ϭ 4I r ͞P* 0 t 2 , P 0 being the power of the source.
rate recovery of optical properties with a minimal number of measurement points is experimentally advantageous and may also reduce the computation time. With respect to the optimization algorithm, we sought to explore how the terrain of the 2 hyperspace is altered when a reduced number of measurements is available. We examined the effect of systematic reduction of the number of measurements by supplying roughly every second, fourth, or eighth measurement of the original data to be sent to the algorithm. This approach provided for the removal of individual data points while preserving the range of s-d separations measured and permitted the use of the multistage algorithm without modification. Figure 4 provides the relative errors in the predicted optical properties as a function of the number of measurements used for each of the five phantoms. The results show that the errors are minimally affected even with an eightfold reduction in the number of measurements supplied to the algorithm. The sole exception is the case of ͑Ј s ͞ a ͒ ϭ 1, where a slight degradation in the recovery of a and Ј s is seen as the number of measurements is reduced. These results demonstrate that as few as 13-15 measurements can be supplied to the algorithm without a significant degradation in its performance.
Notably, the reduction in measurement number did not result in a comparable reduction in the execution time of the algorithm. When 90 -106 distinct measurements were provided to the optimization algorithm, execution required ϳ15 min on a 1.6-GHz Athlon MP processor, whereas data sets containing only 13-15 measurements required ϳ12 min. Thus an eightfold reduction in the number of data points resulted in only a 20% reduction in execution time. Only a minimal reduction in execution time is achieved, because although a reduced number of measurements led to a reduced time per iteration, the number of iterations necessary for convergence increased. Figure 5 illustrates this by displaying the trajectory to the recovered optical parameters in Stage 1 of the algorithm in the case of ͑Ј s ͞ a ͒ ϭ 14, with an initial guess of ͑ a ϭ 0.0004, * s ϭ 0.02, g* ϭ 0.43͒. When every data point is supplied to the algorithm, 22 iterations were required to converge to ͑ a ϭ 0.0075, * s ϭ 0.2704, g* ϭ 0.375͒, whereas when only every eighth data point was supplied to the algorithm, 159 iterations were required to converge to ͑ a ϭ 0.0074, * s ϭ 0.2777, g* ϭ 0.375͒. Thus the reduced number of measurements likely results in a less precise computation of the differential information necessary for the algorithm to navigate through the 2 topology and results in a more circuitous route to the global minimum.

Conclusions
We have demonstrated the efficacy of using a multistaged optimization algorithm that incorporates radiative transport predictions provided by the ␦-P 1 approximation to the Boltzmann transport equation to determine optical properties of infinite turbid media from spatially resolved photon migration measurements. With as few as 13 photon migration measurements taken at source detector separations in the range 0.25l* to 8l*, this algorithm can recover a , Ј s , and l* for media 0.33 Շ ͑Ј s ͞ a ͒ Շ 300, a range of optical properties spanning nearly 3 orders of magnitude. The exceptional range of applicability of this algorithm to successfully determine these optical properties distinguishes this approach from those that use radiative transport predictions that use the standard diffusion approximation or other computational techniques such as neural networks or chemometrics.
Besides the recovery of a , Ј s , and l* for highly scattering and absorbing media, our algorithm exhibits the capability of determining the single-scattering asymmetry parameter g 1 quite well for moderately to highly scattering media ͑Ј s ͞ a ͒ Ն 3. These results compare favorably with prior Monte Carlo-based efforts. Our studies of radiative transport in moderately to highly absorbing turbid media suggest that Fig. 4. Errors in the recovered optical coefficients for each phantom as a function of the number of measurements utilized for l* ͑E͒, a ͑{͒, s ͑‚͒, and g* ͑ƒ͒. Fig. 5. Route to the converged optical properties in Stage 1 for the case ͑Ј s ͞ a ͒ ϭ 14 with an initial guess of ͓ a , * s , g*͔ ϭ ͓0.0004, 0.02, 0.43͔ when every measurement ͑F͒ is used and when every eighth measurement ͑ᮀ͒. The symbols denote the optical properties found after every iteration. A total of 22 iterations were necessary when every measurement was used, whereas 159 iterations were necessary when only every eighth measurement was used.
the recovery of g 1 may also be possible if measurements made at larger s-d separations are available.
The principal elements responsible for the success of our implementation are the ␦-P 1 radiative transport model that accommodates the collimated component of the light field and the use of a multistaged optimization algorithm that appropriately identifies the data necessary for the recovery of the appropriate optical properties at each stage. The first element provides an accurate description of irradiance measurements on spatial scales comparable with a transport mean free path, while the second suggests the proper spatial segregation of the measurements necessary to enable the successful resolution of the inverse problem.
Clearly, much remains to be explored and improved with this approach. The fact that the ␦-P 1 approximation yields excellent radiative transport predictions over a broad range of optical properties prompted our initial focus to investigate its use in recovery of optical properties over an increased range of ͑Ј s ͞ a ͒. While this has clearly proven to be successful, the separation of the collimated and diffuse light components by the ␦-P 1 approximation also enables the determination of the single-scattering asymmetry coefficient g 1 , albeit successfully over a more narrow range of ͑Ј s ͞ a ͒. A focused examination of the accuracy and sensitivity with which g 1 can be determined when media of identical ͑Ј s ͞ a ͒ but different values of g 1 are examined is warranted and currently underway. Moreover, the recent development of radiative transport predictions within the ␦-P 1 approximation for semi-infinite media 20 promises to enable the development of a similar algorithmic approach to determine a , Ј s , and g 1 of semiinfinite media within the context of an analytic radiative transport model.