Research Papers: General

Influence of the phase function in generalized diffuse reflectance models: review of current formalisms and novel observations

[+] Author Affiliations
Katherine W. Calabro

Boston University, Department of Biomedical Engineering, 44 Cummington Street, Boston, Massachusetts 02215

Synopsys Inc., 377 Simarano Drive, Marlborough, Massachusetts 01752

Irving J. Bigio

Boston University, Department of Biomedical Engineering, 44 Cummington Street, Boston, Massachusetts 02215

Boston University, Department of Electrical and Computer Engineering, 8 St. Mary’s Street, Boston, Massachusetts 02215

J. Biomed. Opt. 19(7), 075005 (Jul 15, 2014). doi:10.1117/1.JBO.19.7.075005
History: Received January 31, 2014; Revised April 27, 2014; Accepted May 27, 2014
Text Size: A A A

Open Access Open Access

Abstract.  Diffuse reflectance spectroscopy, which has been demonstrated as a noninvasive diagnostic technique, relies on quantitative models for extracting optical property values from turbid media, such as biological tissues. We review and compare reflectance models that have been published, and we test similar models over a much wider range of measurement parameters than previously published, with specific focus on the effects of the scattering phase function and the source-detector distance. It has previously been shown that the dependence of a forward reflectance model on the scattering phase function can be described more accurately using a variable, γ, which is a more predictive variable for reflectance than the traditional anisotropy factor, g. We show that variations in the reflectance model due to the phase function are strongly dependent on the source-detector separation, and we identify a dimensionless scattering distance at which reflectance is insensitive to the phase function. Further, we evaluate how variations in the phase function and source-detector separation affect the accuracy of inverse property extraction. By simultaneously fitting two or more reflectance spectra, measured at different source-detector separations, we also demonstrate that an estimate of γ can be extracted, in addition to the reduced scattering and absorption coefficients.

Figures in this Article

Diffuse reflectance spectroscopy (DRS), which encompasses techniques with more specific designations, including elastic scattering spectroscopy and light scattering spectroscopy, has demonstrated high value in recent years relevant to the goal of diagnosing disease in vivo1,2 in a wide range of organs, such as brain,3 cervix,4 skin,5 colon,6 breast,7 esophagus,8 pancreas,9 and oral cavity.10 Most commonly, DRS entails the transmission of broadband light [anywhere from near-UV through visible to near-IR (NIR)] through a fiberoptic probe to a turbid medium; after propagating through the sample, a fraction of the light is collected at the surface as reflectance, at a specified distance from the source entry point. Figure 1 illustrates a commonly employed probe geometry, with one fiber used for illumination and a separate fiber for collection. The measured reflectance spectrum can then be analyzed with appropriate mathematical models to extract the optical properties of the tissue, such as the reduced scattering and absorption coefficients (μs and μa, respectively). These optical properties are directly determined by the cell morphology, the extracellular matrix, and the biochemistry and vascularization of the tissue, which undergo a predictable set of changes during the progression of disease.1 Thus, optical property differences are directly related to changes in the structural and biochemical tissue properties, and consequently, quantification of these optical properties can be used not only to inform about the tissue disease status, but also provide important contextual information about the physiological state of the tissue.

Graphic Jump LocationF1 :

Schematic of the two-fiber reflectance geometry.

The key, however, to reliable optical property quantification and disease classification using DRS is in the use of an accurate mathematical model to describe reflectance as a function of the optical properties of the tissue. Development of such reflectance models has been addressed using a variety of analytical and empirical approaches, which are sensitive to a range of parameters, including μs and μa, and the geometric specifications of the measurement probe. The description of reflectance as a function of optical property values was originally modeled using diffusion theory based on the diffusion approximation to the transport equation. 11 However, the use of diffusion theory is generally restricted to the conditions of the diffusion approximation: (1) scattering dominates absorption; and (2) the distance between the illumination and collection points is large compared to the mean-free-path for scattering. Except in the NIR spectral region, condition 1 is not generally met in tissue due to strong absorption by hemoglobin; and condition 2 (large source-detector separation) is often a limitation when accessing internal tissues or when highly localized tissue properties are sought. To avoid these restrictions, alternative models based on modifications to the diffusion theory, Monte Carlo simulations, or empirical formalisms from observation have been developed by various authors over the past 15 years.1224 As with diffusion theory, however, these modified models only account for the relationship between reflectance and the scattering and absorption coefficients and generally do not account for other factors that are important for describing light reflectance under nondiffuse conditions. The importance of the scattering phase function (angular scattering probability distribution) is one such factor that is chronically overlooked, having been addressed for the incorporation into forward models by only a handful of authors.3,19,25,26 Further, most of these models have been developed for a single illumination-collection geometry, focusing either on very small source-detector separations (500μm) or on larger separations (of the order of several millimeters, approaching the diffusion regime), with less attention paid to intermediate distances.

In this work, we present a generalized formalism that describes diffuse reflectance models for two-fiber collection geometries, such as that shown in Fig. 1, over a wide range of input parameters. This formalism provides an intuitive understanding of the trends concerning how reflectance is affected by tissue optical properties and by measurement geometries. In particular, we focus on the effects of the scattering phase function and the source-detector separation, ρ, both of which are known to have strong influence on reflectance. The description of this generalized formalism builds upon previous work, which has explored expressions for reflectance as it relates these factors. We review and highlight how these previous observations all contribute and can be comprehensively understood within the scope of a single framework.

In addition, we present several novel and significant observations that further enhance the understanding of the influence of the phase function and source-detector separation on reflectance models. These include the identification of a dimensionless scattering distance (within the nondiffuse regime) at which reflectance is insensitive to phase function, the use of a unique error-estimation approach to evaluate the accuracy of inverse property extraction (obtaining μs and μa estimates from in vivo reflectance measurements), and the demonstration of a simple method that allows for extraction of phase function information (in addition to μs and μa) when analyzing reflectance measurements at two or more different source-detector separations.

Taken as a whole, the trends and observations described in this work are intended to serve as a resource to better understand, from an intuitive perspective, how various input parameters can be expected to influence reflectance measurements, and the factors in probe geometry or experimental design that can enable the best performance for accurate retrieval of the desired tissue properties. Thankfully, with the recent advancement of graphical processing units (GPUs) for parallel computing, custom reflectance models can now be developed for a wide range of measurement geometries and tissue optical properties (including phase function) with relative ease.

We begin with a comprehensive review of diffuse reflectance modeling and of phase function formalisms and trends specific to the optical characterization of tissue in Sec. 2. This is followed by a description of the Monte Carlo code and methods used for the investigation of reflectance relationships along with a description of the novel error-estimation approach for investigating the accuracy of inverse-model property extraction in Sec. 3. Section 4 presents the generalized reflectance framework that describes the relationships among reflectance values, optical property values (with particular emphasis on phase function), and source-detector separation. Inverse-model error analysis results are also presented and discussed. Previous contributions by others are clearly identified, as well as our own novel advancements. We conclude with key insights and recommendations concerning experimental design to minimize errors due to phase function uncertainty.

Modeling Reflectance

In the context of DRS, a forward model is any function that quantitatively describes the relationships between measured reflectance and the multiple input parameters that will influence reflectance, such as the tissue optical properties and the geometry specifications of the optical measurement. The primary motivation for understanding these reflectance relationships is not for use in forward direction calculations, but rather for use in solving the inverse problem (calculating tissue optical properties from measured reflectance). The relationship between forward and inverse models is expressed generically in Eq. (1), where R represents measured reflectance, p(θ) is the scattering phase function, and G represents all probe geometry parameters. Display Formula

Forward:R=f[μs,μa,p(θ),G]Inverse:[μs,μa,p(θ)]=f1(R,G).(1)

Forward model

The forward model can be described in various ways either with analytical/empirical closed-form equations, or computationally using Monte Carlo simulations or look-up tables. Analytical forms include the diffusion theory equation and its many modified versions.11,13,16,26 Monte Carlo based reflectance models calculate reflectance using a single simulation, making use of similarity relationships to calculate reflectance for a large range of optical coefficients.24 Look-up tables, as their name implies, allow for calculation of reflectance by looking up a value based on a set of input optical parameters within a table of data collected from either phantoms or Monte Carlo simulated data.27 Empirical models use experimental data (again, either from phantoms or Monte Carlo simulated data) and fit it to a proposed closed-form equation that describes the data well.1723,28 Several of these empirical models,23,28 including one developed by our group,22 describe the relationship between reflectance and optical property coefficients (μs and μa) using the framework of Beer’s Law, where reflectance under conditions of negligible absorption, R0, is scaled based on the absorption coefficient, μa, and the average photon path length, L. Scattering-only reflectance, R0, is modeled to be a function of μs, and L is modeled to be a function of μs and μa. Display Formula

R=R0eμaL.(2)
The relationships describing R0 and L as functions of μs and μa contain constant coefficient values that are unique for a given optical probe geometry and tissue phase function, such that the model can be customized for a specific measurement probe.

We note that, regardless of the specific form of the reflectance model (empirical, analytical, look-up table, etc.), as long as the description of each forward model represents the same reflectance relationships (reflectance versus optical property values) under the same geometry conditions, their performance in the inverse direction will be equivalent. Therefore, in this paper, we explore the dependence of reflectance models on the phase function and probe geometry from a general perspective. Relationships between reflectance and optical property values will be illustrated in a graphical format—reflectance as a function of μs or μa—with the understanding that this information can be converted into one of the other model formats as needed. (Empirical equations and lookup table formats provide the most flexibility for this purpose.) These graphical reflectance relationships form the basis of the generalized reflectance model formalism that is presented in Sec. 4.

Inverse model

Because reflectance values are not uniquely associated with a single combination of scattering and absorption properties, using any model in the inverse direction is a fundamentally ill-posed problem and requires use of wavelength-dependent relationships for μs and μa to reduce the dimensionality of the problem. In models of tissue, the spectral shape of μs is often based on a power law,14,15,23,27 and the absorption coefficient is often determined by a linear combination of concentrations and extinction coefficients for oxyhemoglobin and deoxyhemoglobin [εHbO(λ) and εHb(λ), respectively] sometimes with the addition of a correction for a vessel-packing factor, Ccorr, to account for the fact that blood is confined within blood vessels rather than uniformly distributed.29,30 Thus, the tissue parameters can be represented as follows: Display Formula

μs=a(λλ0)bμa(λ)=f1Ccorr(r,λ)μa,blood(f2,λ),(3)
where λ0 is a normalization wavelength. Here, μa,blood, which is the absorption coefficient of whole blood, and Ccorr are calculated as30Display Formula
μa,blood(λ)=cHb,blood[f2εHbO(λ)+(1f2)εHb(λ)]Ccorr(r,λ)={1exp[2μa,blood(λ)r]2μa,blood(λ)r},(4)
where cHb,blood is the concentration of hemoglobin in whole blood. Other contributions to the tissue absorption coefficient (e.g., beta carotene for adipose tissue in the visible, or water for any tissue in the NIR) can be added to μa as appropriate for the specific tissue and wavelength range. The remaining five coefficients include a and b, representing the scattering coefficient at a normalization wavelength, λ0, and the scattering power-law exponent, related to average scatterer size, respectively, and f1, f2, and r, representing the blood volume fraction (%), the hemoglobin oxygen saturation (%), and the average blood vessel radius, respectively. These wavelength relationships for μs and μa [Eqs. (3) and (4)] can then be substituted into any desired forward model. Measured reflectance spectra are fit to the model using a nonlinear least-squares curve fitting method, such as the Levenberg-Marquardt method, with a, b, f1, f2, and r as the fitting coefficients. Thus, each measured spectrum can be characterized by these five physiological parameters. While these wavelength-dependent expressions for μs and μa are specific to tissue, by selecting alternative wavelength-dependent expressions, reflectance models can be used in the inverse direction to extract optical properties of any turbid medium.

Phase Functions

The phase function, p(θ), of a scattering sample specifies the probability distribution of scattering angles and is typically normalized to 1, integrated over all solid angles. To generalize the phase function, the anisotropy factor, g, is defined as the average cosine of the scattering angle, cos(θ). In tissue, g ranges, typically, between 0.75 and 0.99.31 The scattering coefficient, μs, and the anisotropy factor, g, are often combined into a single parameter, known as the reduced scattering coefficient: μs=μs(1g). This similarity relationship is based on the assumed condition that observed optical measurements are equivalent for any combination of g and μs that results in the same μs.32

Interest in the influence of the details of the phase function on reflectance measurements was originally reported by Mourant et al., Bevilacqua and Depeursinge, and Kienle et al., who demonstrated that reflectance is dependent on the specific form of phase function for reflectance measured close to the source (<1mm source-detector separation in tissue), thus negating the applicability of the reduced scattering coefficient similarity condition for short distances.19,33,34 Bevilacqua and Depeursinge further explored the relationship between phase function and reflectance by considering higher moments of the phase function and by developing a new similarity relationship applicable for small source-detector separations.3,19 More recently, Kanick et al. have also utilized this similarity function for single-fiber source-detector geometries.35 In the derivation of the new similarity relationship, the phase function is expanded into a series of Legendre polynomials, where gn are the n’th-order Legendre moments. The first-order moment, g1, is equivalent to the familiar anisotropy factor, g. To calculate its value from any arbitrary phase function, the following equation can be used:36Display Formula

g=g1=θ=0πp(θ)sinθcosθdθ,(5)
Similarly, using the second Legendre polynomial, (3cos2θ1)/2, the second moment, g2, can be calculated as Display Formula
g2=θ=0πp(θ)sinθ[12(3cos2θ1)]dθ.(6)
The calculations of g1 and g2 from Eqs. (5) and (6) assume that the phase function pi(θ) is already normalized, such that θ=0πpi(θ)sinθdθ=1. From these two moments, the new similarity relationship and parameter is defined Display Formula
γ=(1g2)(1g1).(7)
In the nondiffuse regime (i.e., when scattered light is detected after experiencing only a few scattering events), this additional parameter, γ, in addition to μs, helps to more accurately describe the scattering properties of a medium.19 Under these conditions, observed optical measurements will be equivalent for any combination of g, μs, and g2, which result in the same values for μs and γ. From a physical perspective, γ is representative of the relative contribution of near-backward (large-angle) scattering in a phase function. The consideration of γ in quantitative reflectance models has been limited to an empirical model by Bevilacqua and Depeursinge, which employs multiple source-detector separations for use with radially resolved analysis,3,19 and more recently, a model for single-fiber reflectance measurement geometries by Kanick et al.35 Extending the work of these two groups, the generalized reflectance formalism presented in Sec. 4 examines the role of the parameter γ in reflectance models for two-fiber geometries over a broad range of separations.

The use of γ is applicable to any phase function representation, including the commonly adopted Henyey-Greenstein (HG), modified HG (MHG), and more rigorous Mie theory formalisms, as described below.

Henyey-Greenstein and modified Henyey-Greenstein

The HG phase function is the most widely used function in the bio-optics community due to its simplified representation of tissue scattering.37 It is defined by the following equation: Display Formula

pHG(cosθ)=14π1gHG2(1+gHG22gHGcosθ)3/2.(8)
By selecting a given g value, all other Legendre moments are automatically fixed: gn,HG=g1n. Thus, the possible values for γ range from 1.0 to 2.0, as constrained by g1<1.0.

From the few experimental studies on true effective phase functions in tissue, it has been observed that the HG phase function underestimates high-angle backward scattering.34,38,39 To address this drawback, a modified version of the HG function was proposed that combines the standard HG function with an added Rayleigh scattering component, which contributes additional high-angle scattering. Known simply as the MHG phase function, it is defined as Display Formula

pMHG(cosθ)=βpHG(cosθ)+(1β)34πcos2θ,(9)
where β is the fractional contribution of the standard HG phase function, and conversely, (1β) is the fractional contribution of the Rayleigh scattering component.19 The first and second moments are easily calculated from the following: Display Formula
g1,MHG=βgHGg2,MHG=βgHG2+25(1β),(10)
where gHG is the g value used to construct the HG phase function contribution from Eq. (6). With the added flexibility of the MHG phase function, both g and γ can be controlled individually. Values of γ, however, are limited to the range of 1.0 to 2.0 due to mathematical constraints of the HG and MHG phase functions.

Mie theory phase function

A more rigorous estimation of the phase function in tissue can be calculated from Mie theory by approximating tissue as a distribution of discretely sized spherical particles. By combining the individual phase functions from each sphere size, weighted by relative concentrations of different sizes, an average phase function can be constructed to represent the bulk medium. One proposed approach to selecting the appropriate sizes and distribution of particles to best approximate tissue is to use a fractal distribution, based on observations that refractive index variations measured by phase-contrast microscopy scale according to a power law.4042 The number density, Ni, of each particle size di is then defined as Display Formula

Ni=(d0/di)α,(11)
where d0 is the representative size of a spherical volume containing the particles, and α is the fractal volume dimension that relates to the ratio of small versus large particles.40 A large α value corresponds to a greater ratio of small to large particles.

For each discrete particle size, Mie theory is used to calculate the scattering cross-section, σs(di), and the anisotropy factor, g(di), using values for the medium and particle indices of refraction that are representative of cytoplasm (nm) and cellular organelles (np), respectively. Schmitt and Kumar calculated average values of nm=1.352 and np=1.39 to 1.45.42 Using the assumption that waves scattered by individual particles add linearly, the anisotropy factor is defined as Display Formula

g=i=1mdiασs(di)g(di)i=1mdiασs(di),(12)
where m is the total number of discrete particle sizes in the distribution.43,44 The combined phase function, p(θ), can be calculated from a linear combination of the scattering amplitude, S11 (as defined by Bohren and Huffman in their formulation of Mie theory45), scaled by the number density of spheres. Display Formula
p(θ)=i=1mdiαS11(θ,di)i=1mdiα.(13)
Calculation of g1 and g2 can then be performed using Eqs. (5) and (6).

The use of Mie theory and the fractal distribution model of scattering sizes offers a unique opportunity to examine the relationship between g and γ. Throughout this paper, calculations are performed using refractive index values and particle diameter values from the literature:42,43np=1.42, nm=1.352, d=[10,100,150,400,600,800nm,1,3,10,20μm], and a wavelength of λ=600nm. The fractal dimension, α, is varied between 2 and 6, each value of α corresponding to a unique pair of g and γ. Different selections for np, nm, λ, and the discrete particle sizes will yield slightly different numerical results for the relationship between γ and g, but, in general, the trends remain the same.

Figure 2 demonstrates that the relationship between g and γ is nonlinear and that γ increases more rapidly for large g values, especially when g>0.8. This is especially important since most tissues have g values in this higher range. Further, the values of γ in this region are generally >2.0 (dependent on the values of np, nm, d, and λ), which is the upper limit of γ for the HG and MHG phase functions, indicating that a phase function constructed from Mie theory is more appropriate when modeling tissue. Values of γ in tissue have been found to range from 1.7 to 2.2, which highlights the importance of using a phase function capable of γ values >2.3 Bevilacqua and Depeursinge have previously analyzed the space of possible combinations of values for g1, g2, and γ for various phase functions, noting the larger range of possible combinations for Mie phase functions as compared to HG and MHG.19 The analysis here differs by limiting the Mie phase functions to those constructed using the fractal particle distribution model, providing a simplified understanding of the relationship between g and γ that is specifically representative of tissue.

Graphic Jump LocationF2 :

Relationship between the anisotropy values, g and γ (for the specific parameters np=1.42, nm=1.352, λ=600nm), as calculated from Mie theory using a fractal distribution of scattering particle sizes (d=[10,100,150,400,600,800nm,1,3,10,20μm]) with the fractal dimension, α, varied between 2 and 6. Differences in refractive index values and/or particle sizes will provide different numerical results, but the generalized shape of the curve will remain consistent.

For the purposes of comparison, Fig. 3 illustrates constructed functions for the HG, MHG, and Mie phase function descriptions. All three functions are for g=0.7. However, the HG function has γ=1.7, while the MHG and Mie functions have values of γ=1.4, which can be seen in the differences in high-angle scattering between the HG function and the MHG and Mie functions. Note also that, despite having the same γ value, the MHG and Mie phase functions are not identical. Differences between the two result from, and can be observed in variations in the higher-order Legendre moments (gn).

Graphic Jump LocationF3 :

Comparison of a Henyey-Greenstein (HG), a modified Henyey-Greenstein (MHG), and a Mie phase function, all with g=0.7, but with different γ values (HG: γ=1.7, MHG and Mie: γ=1.4).

It is important to note that these phase functions are presented using the standard representation that specifies the distribution of light intensity scattered as a function of angle [p(θ)], i.e., the relative power per unit solid angle that is observed at a given angle in a steady-state system. Alternatively, the phase function can be expressed as a probability distribution of scattering angles [p(θ)], i.e., the probability that a photon will scatter at a given angle. The difference reflects the inherent field versus photon representation issue common in optics; the two forms are related as p(θj)=p(θj)sinθj. In the literature, phase functions have been reported in both ways, without specifying which form is being used, which often leads to confusion.46,47 The angular representation, p(θ), needs to be used when sampling scattering direction changes in Monte Carlo simulations.

Monte Carlo Simulations

“Experimental” data for the work in this paper were obtained from Monte Carlo simulations. Computational simulation allows for a controlled and exact approach to the comparison of various models and input parameters, which would be more difficult experimentally using phantoms. A version of Monte Carlo for multilayered tissue48 adapted for use on a GPU was used for increased speed.49 This code was modified to launch and collect photons within the area and numerical aperture of user-defined fibers to model the two-fiber probe geometry (one illumination and one collection fiber), as illustrated in Fig. 1. Additionally, the code was modified to provide a choice in phase function, allowing the user to select phg, pmhg, pmie, or any other user-defined phase function at run time.

Monte Carlo simulations were run in sets, each having identical geometry and phase function input parameters. Table 1 reports the values of all input parameters used for the simulations. In each set, simulations were run for a set of combinations of reduced scattering coefficient values, μs, and absorption coefficient values, μa. For all simulations, constant values for the refractive index of the tissue, the refractive index of the fibers, the numerical aperture of the fibers, and the fiber diameter were used. All simulations were terminated upon the collection of 10,000 photons, such that the variances of the simulations were equivalent at 1% according to Poisson statistics. Because of the equivalence of variances, and their small 1% value, error bars are not included on plots of results.

Table Grahic Jump Location
Table 1Input parameters for Monte Carlo simulations. Parameters listed with an asterisk (*) are constant for all simulations.

To examine the effects of source-detector separation, sets of simulations were run for a range of values, ρ. And for exploration of the effects of phase function, simulations were run assuming an MHG phase function form using various combinations of g and γ. All combinations of the two variables were constructed, but due to the mathematical limitations of the form of the formula for MHG, phase function combinations with γ=1.9 and g=0.75 and 0.85 are not possible and, thus, could not be simulated and evaluated. Mie theory–based phase functions were also considered, using constant refractive index values, particle sizes d, and wavelength λ, while varying the fractal dimension, α. Resulting anisotropy values of g and γ were calculated using Eqs. (5), (6), (7), and (13). The values were selected such that the values of γ were spaced equally between 1.2 and 2.3; the range of γ was chosen to correspond to physiologically relevant anisotropy values between g=0.55 and 0.96 (see Fig. 2). The particle sizes and refractive index values are representative of tissue.42,43

While only two phase function models (MHG and Mie) were chosen here for analysis, note that regardless of the mathematical form of the phase function, the γ similarity relationship dictates that reflectance values will be equal for equal values of γ, and thus, the details of the higher-order moments of various models is not significant when considering reflectance model trends. Nonetheless, we note that only Mie theory can be invoked for values of γ>2.0.

Inverse Modeling: A Different Approach to Error Analysis

As previously demonstrated in the literature,19,26,35 reflectance models depend on phase function, and hence, the choice of a model to use in the inverse direction requires knowledge of the tissue phase function. However, in most cases, this information is not known a priori and, consequently, assumptions must be made. An incorrect assumption will affect the accuracy of extracted scattering and absorption properties from measured reflectance; we wish to quantify the degree to which an incorrect assumption affects these extracted values. Unfortunately, because there is currently no gold standard for measuring the optical and physiological parameters from point measurements in tissue, it is not possible to experimentally test the effects of model differences. Therefore, the error analysis presented in this work uses computationally constructed reflectance data, which allows for a more controlled and exact approach to the comparison of various models and input parameters.

This analysis first involves computationally constructing a set of reflectance spectra that are representative of the true model and tissue properties. Then the assumed reflectance model, which is different from the true model due to inaccurate assumptions about tissue phase function, is used in the inverse direction to extract the fitted optical property values. The differences between these extracted values and the actual property values, thus, provide a measure of error that represents the potential uncertainty in in vivo measurements. We note moreover, that this technique is novel in its approach to examining error in terms of extracted physiological values from reflectance spectra, as opposed to reporting error exclusively of the optical coefficient values.

It is difficult to define a single estimate of error because the errors in extracted quantities are directly related to μs and μa themselves. Hence, errors are calculated for multiple reflectance spectra, each with different scattering and absorption properties. A set of wavelength-dependent μa and μs spectra are first calculated using Eqs. (3) and (4), over a wavelength range from 400 to 800nm at a resolution of 1 nm (for a total of 400 data points). Values for the five physiological fitting variables (a, b, f1, f2, and r) were selected such that the resulting μa and μs spectra would span the ranges characteristic of tissue. For scattering, four spectra were selected to represent overall low scattering (a=8, b=1), medium scattering (a=16, b=1), high scattering (a=24, b=1), plus one with a large slope representing a high density of smaller scatterers (a=16, b=2). The reference wavelength λ0 was selected to be 600 nm. For absorption, spectra were selected for a physically representative range of absorption values (blood volume fraction, f1=0.5, 2, and 5%), all with a moderate hemoglobin oxygen saturation level (f2=75%). For all cases, the average blood vessel radius is held constant at r=10μm which is representative of capillaries. Using all combinations of the scattering and absorption conditions, a set of 12 unique scattering and absorption spectra were used for testing.

Computationally constructed reflectance spectra, representing measured data, are calculated using the test scattering (μs) and absorption (μa) spectra according to Eq. (1), where the forward reflectance model is specific to a single phase function [p(θ)] and measurement geometry (G). Reflectance values are calculated using the forward model via a lookup table composed of Monte Carlo simulated reflectance data specific to the unique pair of p(θ) and G. A bilinear interpolation algorithm (available in MATLAB®) was used to calculate reflectance values for combinations of μs and μa not explicitly listed in the table. Noise of 1% was added to the spectra to better represent experimentally acquired data.

As previously mentioned, when using reflectance models in the inverse direction, the least-squares fitting algorithm is insensitive to the specific form of the forward model as long as the modeled reflectance is negligibly different over the set of μa and μs values. Hence, the choice of model format to use in the inverse direction can be based primarily on preference. For flexibility, we have chosen to use an interpolation lookup table method (with the same lookup table format used to construct the reflectance spectra). This is easily implemented in MATLAB® using the lsqcurvefit toolbox (using the Levenberg Marquardt optimization algorithm), as it can be constructed with any user-defined computational function for the forward model regardless of functional form.

When calculating optical property extraction error associated with inaccuracies in phase function assumptions, reflectance spectra were calculated for all combinations of μs and μa spectra listed above, yielding a set of 12 reflectance spectra. The assumed reflectance model, representing an inaccurate choice of phase function, is then used in the inverse direction to fit each of the 12 true reflectance spectra and extract the fitting coefficients (a, b, f1, f2, and r) and, consequently, the μs and μa spectra. Errors in extracted optical and physiological values are reported as an average across all 12 combinations on a percentage basis for the physiological parameters (a, b, f1, f2, and r), and as the mean percentage error across the entire wavelength range for the optical coefficients (μa and μs).

Phase Function Dependence

Many reflectance models, specifically empirical models describing reflectance at small source-detector separations, have assumed an HG phase function with an anisotropy value g=0.9.17,18,20,22 As introduced in Sec. 2.2, the MHG and the Mie phase functions have been shown in the literature to be more representative of true scattering in tissue as a consequence of the higher probability for large-angle scattering.

To examine the impact of these phase functions that are more representative of tissue, specifically the influence of both g and γ, reflectance relationships were simulated for a range of g and γ values. The MHG phase function is user-friendly for this purpose, since it conveniently enables independent variation of g and γ. Phase functions constructed using values of g and γ described in Sec. 3.1 were examined using a source-detector separation of 250 μm. Figure 4 illustrates the resulting reflectance relationships when each of the phase functions is used for simulation.

Graphic Jump LocationF4 :

Reflectance relationship plots for a probe geometry with ρ=250μm, and different MHG phase functions. Colors and symbols correspond to different γ values and line styles correspond to different g values. Note that plots for the combinations g=0.95 and γ=1.7 and 1.9 are not present due to the mathematical limitations of the MGH phase function with these values.

In much of what follows, reflectance models are reported graphically in a set of two plots, as illustrated in Fig. 4. The left plot illustrates the relationship between the relative reflectance and the reduced scattering coefficient, μs, in the absence of absorption, and the right plot illustrates the relationship between reflectance and the absorption coefficient, μa. In the plot of reflectance versus μa, reflectance values are normalized to 1 at μa=0. This highlights the shape differences rather than the overall amplitude differences in reflectance, as caused by differences in μs. The two subplots correspond to the two components of the empirical model format inspired by Beer’s Law, presented in Eq. (2); the first represents Eq. (2) when μa=0 such that RREL=R0(μs), and the second represents Eq. (2) normalized such that RNORM=R/R0=eμaL. This graphical representation is useful for identifying trends, as it helps to distinguish the dependences of scattering and absorption separately.

Most importantly, Fig. 4 reveals that when γ is held constant, the influence of g is remarkably small. Conversely, when g is held constant, reflectance values vary significantly with γ. At μs=5cm1 and μa=0, reflectance for γ=1.3 is more than three times larger than for γ=1.9. Differences due to γ in the reflectance versus μa relationships are more moderate, but still show a significant dependence. This observation for a two-fiber geometry reflectance relationship follows from previous work by both Bevilacqua et al. and Kanick et al. (for radially resolved multifiber geometries and for single-fiber geometries) in that the critical phase function variable to consider when building reflectance models (for short source-detector distances) is γ, and not g.3,19,35 We further emphasize not only that γ is more influential than g, but also that the influence of g is limited; phase functions with different g values, but the same γ value, have reflectance relationships that are essentially similar.

Because the MHG phase function is restricted to γ values <2, and yet Mie phase functions constructed to represent tissue have possible γ values >2, we also examined how Mie phase functions, which are more representative of tissue, influence reflectance relationships for an even wider range of γ. These results are presented in Fig. 5. As compared to the reflectance relationships built with the MHG phase function, the differences due to γ in the reflectance versus μa relationships are qualitatively similar. The differences in reflectance for the reflectance versus μs relationships, however, have increased such that, at μs=5cm1, reflectance for γ=1.2 is nearly six times larger than for γ=2.3. Further, as γ increases, the differences in reflectance become more substantial. For example, when considering g values most representative of tissue (g>0.85), the differences in reflectance values of γ between 1.75 and 2.3 are comparatively greater than those for values between γ=1.2 and 1.75, even though, in the first case, g only ranges from 0.86 to 0.96, whereas in the second, g ranges from 0.55 to 0.86. This indicates that the impact of γ is even stronger for g values most commonly expected in tissue.

Graphic Jump LocationF5 :

Reflectance relationship plots for a probe geometry with ρ=250μm, for a range of Mie phase functions, each with unique g and γ values.

We offer an intuitive explanation for trends observed regarding the influence of γ on reflectance, which can be explained by the relative number of large-angle scattering events that are expected to occur. For the reflectance versus μs relationships, reflectance is smaller when γ is larger. At small source-detector-separations, such as that used in the above simulations (250 μm), for a photon to be collected, it can only scatter a modest number of times, and at least one of those events must be at a large angle such that the photon can turn around and reach the collector. Lower γ phase functions have larger probabilities for these backward-scattering events (see Fig. 3) and, therefore, allow for more photons to be redirected back toward the fibers and collected. Conversely, the low probability of backscattering characterized by a high-γ phase function will result in photons traveling far from the source and detector due to the predominance of forward scattering, and when the photon does experience a large-angle scattering event, it will likely be so far away from the detector that the probability for collection is low. This effect is amplified for small values of μs, since average path lengths are longer, and therefore, photons will propagate much farther from the collection fiber with a fewer number of scattering events. With fewer opportunities for large-angle scattering, fewer photons will reach the collector.

The influence of backscattering can also explain why reflectance is greater for high γ values in the reflectance versus μa relationships. As we have established, with a larger γ, the probability of high-angle scattering is lower, and the probability of forward scattering is larger. Many forward-scattering events will move the photon away from the detector very quickly. Therefore, the most likely scenario for a photon to be collected would then be to experience one high-angle scattering event in conjunction with a small number of forward-scattering events. When γ is lower, the photon experiences more scattering events at moderate scattering angles, which allows it to scatter more times while still staying near the collector, thus enhancing its probability for being collected. This larger number of scattering events corresponds to a longer average path length, and, as determined by Beer’s law, a longer traveled path length results in higher absorption and, thus, a lower reflectance.

Error analysis

Under practical conditions, when reflectance from a single collection fiber is being used for an inverse analysis to extract optical property values, assumptions about the value of γ and its corresponding reflectance model must be made since the exact value of γ is typically unknown. With most tissues of interest, even if the anisotropy value, g, is assumed to be >0.85, the range of γ is still quite large (γ=1.75 to 2.3), as is demonstrated in Fig. 2. Within this range, it is then logical to assume an average value of γ near 2.0. To estimate the effect that an inaccurate assumption will have on the accuracy of the extracted optical coefficient values, the novel error analysis method, described in Sec. 3.2, was employed. A set of 12 reflectance spectra were constructed from the μs and μa test spectra, using the reflectance model for γ=1.75. The reflectance model for γ=2.05 was then used in the inverse direction to extract the optical property values from the reflectance spectra. This represents a situation where the tissue being measured has a phase function with γ=1.75, but the inverse model incorrectly assumes a phase function with γ=2.05. This analysis was also performed assuming a large γ value, γ=2.3. Figures 6 and 7 illustrate the errors in optical coefficient values that result from this analysis for an example case with tissue physiological properties a=16, b=2, f1=5%, f2=75%, and r=10μm.

Graphic Jump LocationF6 :

(a) Fitting results of forward model to reflectance data built using a phase function with γ=1.75, but fit using a model for a phase function with γ=2.05. (Optical property spectra built assuming a=16, b=2, f1=5%, f2=75%, and r=10μm.) (b) The real and extracted optical coefficient spectra.

Graphic Jump LocationF7 :

(a) Forward model fitting results to reflectance data built using a phase function with γ=2.3, but fit using a model for a phase function with γ=2.05. (Optical coefficient spectra built assuming a=16, b=2, f1=5%, f2=75%, and r=10μm. (b) The real and extracted optical coefficient spectra.

When γ is underestimated, the μs and μa extracted values are higher than expected, and when γ is overestimated, the μs and μa extracted spectra have values lower than expected. The mean errors in optical and physiological values, averaged over both cases, are presented in Table 2. Note that these errors represent comparison of γ values with differences 0.3. In the extreme case, when comparing forward models with differences in γ>1.0, extraction errors in μa and μs were found to exceed 100%. But in realistic situations, assuming an average γ for the inverse fit is much more reasonable.

Table Grahic Jump Location
Table 2Average errors in the extracted physiological and optical property values by incorrectly assuming a model with γ=2.05, when the true tissue has γ=1.75 or γ=2.3.

We remind the reader that the errors reported here are the averages over the 12 combinations of μs and μa test spectra. However, the errors are also dependent on the specific μs and μa values being used. For scattering in particular, the low-scattering spectra (a=8cm1, b=1) will result in extraction errors in μs that are two to three times larger than for the high-scattering spectra (a=24cm1, b=1). Although partially due to the fact that percentage errors are greater for smaller values, this is primarily due to the greater differences in reflectance at smaller μs values from variations in γ (see Figs. 4 and 5). Because differences in reflectance as a function of μa, which are due to variations in γ, are much smaller than differences in reflectance as a function of μs, the errors in extracted μa, f1, f2, and r are also significantly smaller.

Note on goodness-of-fit

A new and remarkable observation from Figs. 6 and 7 is that, despite significant errors in the extracted optical properties, the reflectance spectra for the forward model fit are almost perfect! As compared to the fit of the reflectance spectra when the true and assumed values of γ match, the fits of the reflectance spectra when the true and assumed values of γ are different result in fittings of the reflectance model that are only slightly inferior. This is quantitatively verified by comparing the mean normalized residual of the fits: <RfitRtrue|/Rtrue> over all wavelengths. For matching γ, the average of the mean normalized residuals (over the 12 test spectra) is 0.0045, and for nonmatching γ, the average value is 0.0056, both of which are within the 1% noise that was added to the reflectance spectra.

The implication of this novel observation is that accuracy in the extracted optical coefficient values cannot be associated with the goodness-of-fit of the forward reflectance model to the data. This is of particular significance because the validity and suitability of a given reflectance model is often reported as being verified by the goodness-of-fit that it demonstrates when used to fit experimental reflectance data. In reality, the goodness-of-fit measure can only be used to verify a model if that model accounts for all sources of potential variability (phase function, source-detector separation, fiber diameter, fiber and probe material indices of refraction, etc.), and if the wavelength-dependent assumptions concerning the spectral shapes of μs and μa (when using wavelength-dependent fitting as is employed in this work) are well representative of the tissue being measured. In practice, it is prudent not to assume a significant level of certainty in these factors to make conclusive verifications based on goodness-of-fit alone.

While this observation regarding goodness-of-fit as a deceptive measure of model accuracy is most critically applicable to uncertainties in phase function, it is also applicable to conditions for which other model variables are inaccurate. For example, when the reflectance model being used to fit a measured spectrum is built using incorrect values for probe geometry details (e.g., numerical aperture, fiber diameter, etc.), similar inaccuracies in extracted physiological and optical properties will result despite excellent spectral fits.50 This suggests that building reflectance models that are unique to the specifications of individual measurement probes is important for maximal accuracy; the use of GPU-enabled Monte Carlo codes makes this much more feasible, especially when the data can be used directly in lookup tables for inverse model fitting.

Effects of Source-Detector Separation

To explore the influence of fiber separation, i.e., the distance between the source and detector fibers, reflectance relationships were constructed from Monte Carlo simulations for a number of source-detector separations, ρ=[250,500,1000,3000,2000,5000,10,000]μm. For each separation distance, several phase functions constructed from Mie theory were used. The separation distances were chosen to represent the nondiffuse regime (500μm), the diffuse regime (5mm), and the intermediate region (1, 2, and 3 mm). The intermediate region is of particular interest because reflectance models have not been extensively investigated for these source-detector separations. Thus, by examining all three regimes simultaneously, a more comprehensive understanding of how tissue optical properties influence reflectance can be attained.

As has been previously utilized by others,25,26 we employ a convenient method of investigating the combined effects of both source-detector separation and optical coefficient value in a single parameter for both scattering (μsρ) and absorption (μaρ). These values are referred to as dimensionless scattering and dimensionless absorption, respectively. Correspondingly, relative reflectance is presented in a dimensionless form as dimensionless reflectance (Rρ2). This makes it possible to investigate trends in reflectance relationships in terms of dimensionless variables without any loss of generality.

Figure 8(a) illustrates the dimensionless relationship between reflectance and scattering for all simulated values of ρ between 250 μm and 1 cm. Figure 8(b) illustrates the relationship between normalized reflectance (which is not explicitly presented as dimensionless since it represents the ratio of two dimensionless reflectance values such that the ρ2 values cancel) and dimensionless absorption; in this plot, data are provided for the case where dimensionless scattering (μsρ) is equal to 0.25. (Note that the relationship between normalized reflectance and dimensionless absorption is unique for each value of dimensionless scattering.) Given the values of μs and ρ used to generate the simulation data, only a limited number of combinations of μs and ρ yield a dimensionless scattering value of 0.25. Thus, the data in Fig. 8(b) are provided for only ρ=250μm, 500 μm, 1 mm, and 5 mm. The choice of presenting reflectance, μs, and μa in dimensionless units is for simplicity of normalization and to show consistency of trends over a large range of ρ values. It is important to note that when plotted for a fixed value of ρ using nondimensionless units, the shape of the curve remains the same; the dimensionless normalization simply scales the reflectance values.

Graphic Jump LocationF8 :

Dimensionless reflectance relationships for a range of source-detector separations. (a) Dimensionless reflectance versus dimensionless scattering. (b) Normalized reflectance versus dimensionless absorption for a dimensionless scattering value (μsρ) of 0.25. All data are for a phase function gamma value (γ) of 2.3.

The results presented in Fig. 8 also enable several important observations about reflectance relationships, illustrating the significant variations in reflectance over large dimensionless scattering and absorption ranges, especially for the reflectance versus scattering curve. In contrast to the nearly linear relation of relative reflectance versus scattering illustrated in Figs. 4 and 5, the plot in Fig. 8(a) demonstrates considerable nonlinearity. For low dimensionless scattering, the dimensionless reflectance increases as both scattering and source detector separation increase. But at a given distance, dimensionless reflectance reverses its relationship with μsρ and begins to decrease as μsρ increases.

The observation of this scattering-dependent reflectance relationship is a key component of the generalized reflectance formalism. Several authors have observed this behavior experimentally and it is also predicted by diffusion theory.25,51,52 Further, Kumar and Schmitt related the curve peak to an ideal distance at which reflectance is only slightly dependent on μs, which they reported to be between 2 and 5 mm, depending on the range of μs values of interest.51 For a typical tissue value of μs of 10cm1, this range of separation distances corresponds to dimensionless scattering values between 2 and 5, which correlates well with the peak of the dimensionless reflectance versus dimensionless scattering curve. It is at this peak value that moderate changes in μs will result in only minimal change in reflectance, suggesting that reflectance is insensitive to μs at source detector separations between 2 and 5 mm.

A novel and significant observation results from extending this analysis to consider variations in phase function. Figure 9 illustrates the dimensionless reflectance versus dimensionless scattering relationships for a range of phase function γ values (1.2, 1.75, and 2.3). The critical observation in this figure is the existence of an isosbestic point, where reflectance is not dependent on γ at a dimensionless scattering value of 0.7. (Note that we adopt the term “isosbestic” loosely, with the intent of describing a condition for which one measurement parameter is independent of another.) While reflectance is known to be insensitive to phase function in the diffusion regime (μsρ>10), it is an unexpected revelation that there is an additional singular dimensionless scattering distance, well below the diffusion regime, at which this insensitivity to phase function also exists. Because this isosbestic point is part of the generalized reflectance model, it is applicable to any turbid medium (i.e., independent of tissue model), and is not dependent on the fiber probe (including fiber separation, which is inherently part of dimensionless scattering).

Graphic Jump LocationF9 :

Dimensionless reflectance versus dimensionless scattering for a range of phase function γ values, plus diffusion theory. The three curves with varying γ values all intersect at a μsρ value of 0.7, although diffusion theory does not.

For low dimensionless scattering values, a high γ corresponds to lower reflectance, as was observed in Figs. 4 and 5 (where μsρ ranged from 0.125 to 0.75). As dimensionless scattering increases, however, a point is reached at which reflectance is insensitive to the value of γ. And at yet larger dimensionless scattering values, a larger value of γ results in a larger reflectance. As noted earlier, we hypothesize that this is because, in tissues with low γ, photons with more backward scattering and less forward scattering are more likely to remain close to the source. With high γ, more forward scattering and fewer backscattering events allow a photon to travel farther from the source before experiencing the high-angle scattering event that will direct it back toward the surface for collection. At large scattering and source-detector separations, this means that the proportion of photons collected at larger source-detector separations to photons collected at smaller separations is higher for large γ, leading to higher reflectance.

For comparison, the diffusion theory equation was used to calculate dimensionless reflectance as a function of dimensionless scattering, and is also included in Fig. 9 (dashed line). Recalling that diffusion theory is not dependent on phase function, above the isosbestic point the diffusion theory plot best approximates the γ=1.2 curve, and below the isosbestic point diffusion theory best approximates the γ=2.3 curve. This same relationship between diffusion theory and Monte ○Carlo simulations with different anisotropy g values was recently observed by Kim et al. for values of μs and ρ corresponding to dimensionless scattering from 0 to 1,52 which is consistent with our findings. At large dimensionless scattering values, the similarity between diffusion theory and the γ=1.2 plot is consistent with the fact that diffusion theory is descriptive of isotropic scattering; since the value γ=1.2 corresponds to g=0.55, the low γ curve (in this range) is closest to the isotropic case (g=0 and γ=1).

We now plot normalized reflectance versus dimensionless scattering, which provides different observations compared with the previously illustrated dimensionless reflectance versus dimensionless scattering relationships in Figs. 8 and 9. Figure 10(a) presents the relationships over a large range of dimensionless absorption values from 0 to 10, and for five values of dimensionless scattering between 0.1 and 8.0. We note that as dimensionless scattering increases, the plots relating normalized reflectance and dimensionless absorption do not monotonically progress. For values of dimensionless scattering below 2.5, normalized reflectance is larger for larger values of μsρ, as was seen in the relationships illustrated in Figs. 4 and 5. However, for dimensionless scattering values >2.5, normalized reflectance is smaller for larger values of μsρ. To better visualize this nonmonotonic relation, Fig. 10(b) plots normalized reflectance as a function of dimensionless scattering for a single value of dimensionless absorption, μaρ=4.0.

Graphic Jump LocationF10 :

(a) Normalized reflectance versus dimensionless absorption for five dimensionless scattering values between 0.1 and 8.0. (b) Normalized reflectance versus dimensionless scattering for a constant value of dimensionless absorption of 4.0.

This relationship was explored by Mourant et al., who reported on an optimized source-detector separation for minimizing the effects of scattering.53 In contrast to the analysis of Kumar and Schmitt,51 who identified a source-detector separation that reduced variance in reflectance, Mourant et al. examined the source-detector separation at which variation in photon path length is minimized (enabling application of Beer’s law for determining the absorption coefficient), which they experimentally determined to be 1.7mm. That observation is supported by the results presented in Fig. 10(b), which illustrates that for dimensionless scattering values between 2 and 3, the variation in normalized reflectance (which is directly related to path length) is minimal. When considering average values of μs in tissue, this range of dimensionless scattering values is consistent with the reported 1.7 mm source-detector separation distance for which normalized reflectance is independent of scattering (at a given μa).

Like the dimensionless reflectance versus dimensionless scattering relationships, the normalized reflectance versus dimensionless absorption relationships are also dependent on phase function. Figure 11 presents the percent difference between normalized reflectance values with a phase function value γ=1.2 and normalized reflectance values with a phase function value γ=2.3. The percent difference is illustrated for a range of dimensionless scattering values between 0.1 and 8.0. As expected, when absorption is zero, the error in normalized reflectance is zero, since, by definition, normalized reflectance is equal to 1.0 when μa=0. When dimensionless scattering is low (corresponding to low μs and/or ρ), the error in normalized reflectance due to phase function is maximal, reaching as much as 110%. As dimensionless scattering increases and approaches the diffuse regime, the error in reflectance due to phase function variation reduces significantly. However, even at a dimensionless scattering value of 8.0, the maximal reflectance error due to phase function is still over 20%, which is consistent with the known breakdown of the diffusion approximation for large values of absorption. The importance of this observation is that the phase function is an important variable for reflectance even at relatively large source-detector separations.

Graphic Jump LocationF11 :

Percent difference in normalized reflectance values between simulations when the phase function parameter γ=1.2 and γ=2.3. Percent differences are illustrated for a range of dimensionless scattering values between 0.1 and 8.0.

Error analysis at larger source-detector separations

In Sec. 4.1.1, the errors in physiological and optical fitting parameters that result from inaccurate assumptions about the phase function were examined for a probe geometry with ρ=250μm. Here, we analyze the errors that result from inaccurate phase function assumptions for larger distances: ρ=500μm, 1 mm, 2 mm, 5 mm, and 1 cm. The same analysis procedure that was described in Sec. 3.2 and implemented in Sec. 4.1.1 is used here. The errors for all tested reflectance spectra are presented in Table 3, representing the expected errors due to incorrect assumption about the average value of γ. Errors are smallest for ρ=1cm due to the minimal differences in the reflectance versus optical coefficient relationships as a function of phase function at this distance, and are largest (specifically for μs) when ρ=250μm since the dependence of phase function is greatest at small source-detector separations. However, errors do not decrease monotonically as source-detector separation increases. Rather, for μs, the error initially decreases when ρ=500μm before increasing to a local maximum when ρ=2mm and decreasing thereafter as ρ increases to 1 cm. The local error minimum at ρ=500μm is the result of the isosbestic point in the reflectance versus dimensionless scattering relationships, because variations in reflectance due to γ are at a minimum [Fig. 8(a)]. At ρ=500μm, the reflectance isosbestic point falls directly in the middle of the physiologically relevant range of μs for tissue (10cm1). With minimal variation in reflectance due to γ at this point, the inverse modeling errors that result from uncertainty in the tissue phase function are also at a minimum. This indicates that if the phase function of a tissue cannot be accurately estimated (as is typically the case), yet practical constraints do not permit large source-detector separations to take advantage of diffusion theory, it is advantageous to design a probe geometry with a source-detector separation close to 500 μm as a means to minimize errors due to phase function uncertainty.

Table Grahic Jump Location
Table 3Average errors resulting from inaccurate assumptions about the phase function of tissue when performing inverse calculations due to incorrect average value of γ used. Reported as averages over all test cases of γ wavelength dependence, and μs and μa test spectra.
Extracting Phase Function Information

While extracting μs and μa values from tissue has generally been of primary interest, the variation in reflectance as a function of γ presents an encouraging opportunity to extract additional information about the tissue using DRS. Given that γ relates to the relative backscattering contributions of a tissue phase function, which is dictated by the size distribution of scatterers, the ability to quantify γ can provide knowledge about the tissue microstructure. A handful of other authors have also examined this opportunity3,13,54,55 for other measurement geometries (e.g., single-fiber probes that both illuminate and collect light) and have observed that there is not enough information in a single reflectance spectrum to uniquely identify μs, μa, and phase function information simultaneously, and therefore, multiple measurements are necessary. In particular, Kanick et.al. have demonstrated that the value of γ can be extracted from multiple reflectance measurements using multiple single-fiber probes of different diameters.35 Extending this to the two-fiber collection geometry, collecting measurements at various source-detector separations is a natural choice, especially since the reflectance relationships vary significantly as a function of ρ.

The implementation of this approach involves measuring reflectance spectra at multiple distances and simultaneously fitting the multiple reflectance spectra to a reflectance model in two dimensions, with both ρ and λ being independent variables. In addition to the original inverse fitting coefficients described in Sec. 2.1.2, γ then becomes the sixth fitting coefficient. The reflectance model to which the data are fit is a set of three-dimensional lookup tables for which relative reflectance is defined as a function of μa, μs, and γ. Each three-dimensional lookup table is unique to a single source-detector separation geometry. Experimental reflectance spectra are fit to their respective probe geometry reflectance lookup tables, as would be done with only one reflectance measurement being analyzed, but a least-squares algorithm performs iterative calculations on all data simultaneously, optimizing the physiological fitting coefficients, including γ, for the entire set of data.

Using the μs and μa spectra described in Sec. 3.2, reflectance spectra were constructed for a sample value of γ and for a number of different source-detector separations. Various pairs of reflectance spectra, representing different pairs of source-detector separations, were then analyzed in the inverse direction using the above lookup table approach. In all cases, the extracted optical coefficient values, including the value of γ, were accurately retrieved with errors <0.25%. This indicates that with only two separate reflectance spectra and with the differences in source-detector separations as small as 250 μm (when the two reflectance spectra are measured using ρ=250 and 500 μm), not only can errors in the extracted optical properties due to uncertainty in phase function be significantly reduced, but an additional optical parameter can be gained (that of the phase function parameter gamma, γ), which can serve to better characterize and differentiate tissue types.

The work reported here constitutes a broad and critical examination of the combined effects of phase function and source-detector separation on the development of diffuse reflectance models for the purpose of extracting estimates about the optical and physiological properties of small volumes of tissue. We have reviewed previous advancements in the understanding of phase function and source-detector separation reflectance relationships, while also contributing a number of novel observations which are combined into a single and comprehensive formalism. Unlike previous developments for diffuse reflectance modeling (for which results were presented in forms specific to a given model construct), this formalism provides a generalized understanding of the ways that measured reflectance can be affected by a wide range of conditions. The key insights of this formalism are as follows:

  • Especially at small source-detector separations, reflectance is dependent not only on μs and μa, but also on the specific form of the phase function of the tissue. The phase function similarity variable γ, and not the more commonly used anisotropy factor g, is the key parameter to consider when uniquely characterizing the optical properties of a tissue. Reflectance values can vary by a factor of 3, depending on the value of γ. Consequently, reflectance relationships, which form the basis of reflectance models, are also dependent on the value of γ.
  • The effects of the optical property values (μs and μa) and the source-detector separation distance (ρ) on reflectance can be combined into dimensionless parameters, μsρ and μaρ. Reflectance can then be uniquely described by μsρ, μaρ, and γ.
  • Reflectance models are also strongly dependent on the source-detector separation. The relationship between reflectance and the reduced scattering coefficient, μs, is particularly interesting when described generally using dimensionless reflectance (Rρ2) and scattering (μsρ) parameters. At small dimensionless scattering values, the relationship between reflectance and μsρ is quasi-linear, with reflectance increasing with μsρ. As dimensionless scattering increases, the relationship becomes increasingly nonlinear. Reflectance increases with μsρ for smaller values before reaching a maximum, and then decreases with μsρ at larger values. The peak in reflectance versus dimensionless scattering exists at approximately μsρ=5. These relationships provide a general understanding of reflectance over the wide range of conditions that span the nondiffuse regime (small μsρ), the diffuse regime (large μsρ), and the rarely investigated intermediate regime.

The novel observations and contributions to this formalism presented in this work include the following:

  • The effects of inaccurate phase function assumptions can be explored using an error estimate approach that considers the errors not only in μs and μa, but also of physiologically related parameters during wavelength-dependent spectral fitting of reflectance measurements. Forward-model curve fits to reflectance spectra can be remarkably good, even when the models invoke inaccurate phase functions, indicating that goodness-of-fit is not an adequate measure of the suitability of a reflectance model for extracting optical properties.
  • When dimensionless scattering is small, reflectance is higher for lower values of γ, and when dimensionless scattering is high, reflectance is higher when γ is higher. At an intermediate value of dimensionless scattering, there is an isosbestic point at which reflectance is insensitive to both the g and γ parameters of the tissue phase function. This point is found at approximately μsρ=0.7 and can be used to guide design of reflectance probe geometries as a way to reduce errors resulting from uncertainty in the tissue phase function. To do so, the source-detector separation distance can be selected such that the isosbestic point falls in the middle of the range of expected measured μs values. Further, and perhaps unexpectedly, the range of the dimensionless scattering plot for μsρ=3 to 6 still exhibits significant phase function dependence despite being close to what is normally considered the diffuse regime.
  • By simultaneously fitting two or more reflectance spectra, measured using different source-detector separation distances, an estimate of the phase function value, γ, can be extracted in addition to the scattering and absorption coefficient results, which also reduces the error in extracted μs and μa. A straightforward means of achieving this is to use a three-dimensional lookup table that describes reflectance as a function of μs, μa, and γ.

The presentation of this generalized reflectance formalism is not intended to provide numerical model constructs that can be used directly in experimental diffuse reflectance studies. Rather, our intent has been to provide an intuitive understanding of the trends that form the foundations of useful reflectance models. According to the specific needs of an application, probe geometries should be designed to limit the errors associated with factors such as phase function; reflectance values should be collected as a function of μs, μa, and γ (ideally using Monte Carlo simulations such that input parameters can be most accurately controlled), and those data should be used to construct a probe-specific reflectance model that can be used in the inverse direction to extract accurate tissue optical and physiological properties.

This work was supported in part by the NIN/NCI (U54 CA104677), the Boston University Photonics Center, the NSF Graduate Research Fellowship, and the Center for Integration of Medicine and Innovative Technology (CIMIT).

Brown  J. Q. et al., “Advances in quantitative UV-visible spectroscopy for clinical and pre-clinical application in cancer,” Curr. Opin. Biotechnol.. 20, (1 ), 119 –131 (2009). 0958-1669 CrossRef
Bigio  I. J., Mourant  J. R., “Ultraviolet and visible spectroscopies for tissue diagnostics: fluorescence spectroscopy and elastic-scattering spectroscopy,” Phys. Med. Biol.. 42, (5 ), 803 –814 (1997). 0031-9155 CrossRef
Bevilacqua  F. et al., “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt.. 38, (22 ), 4939 –4950 (1999). 0003-6935 CrossRef
Weber  C. R. et al., “Model-based analysis of reflectance and fluorescence spectra for in vivo detection of cervical dysplasia and cancer,” J. Biomed. Opt.. 13, (6 ), 064016  (2008). 1083-3668 CrossRef
Salomatina  E. et al., “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt.. 11, (6 ), 064026  (2006). 1083-3668 CrossRef
Rodriguez-Diaz  E. et al., “Spectral classifier design with ensemble classifiers and misclassification-rejection: application to elastic-scattering spectroscopy for detection of colonic neoplasia.,” J. Biomed. Opt.. 16, (6 ), 067009  (2011). 1083-3668 CrossRef
Palmer  G. M. et al., “Monte Carlo-based inverse model for calculating tissue optical properties. Part II: application to breast cancer diagnosis,” Appl. Opt.. 45, (5 ), 1072 –1078 (2006). 0003-6935 CrossRef
Georgakoudi  I. et al., “Fluorescence, reflectance, and light-scattering spectroscopy for evaluating dysplasia in patients with Barrett’s esophagus,” Gastroenterology. 120, (7 ), 1620 –1629 (2001). 0016-5085 CrossRef
Wilson  R. H. et al., “Optical spectroscopy detects histological hallmarks of pancreatic cancer,” Opt. Express. 17, (20 ), 17502 –17516 (2009). 1094-4087 CrossRef
Schwarz  R. A. et al., “Autofluorescence and diffuse reflectance spectroscopy of oral epithelial tissue using a depth-sensitive fiber-optic probe,” Appl. Opt.. 47, (6 ), 825 –834 (2008). 0003-6935 CrossRef
Kienle  A., Patterson  M. S., “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A. 14, (1 ), 246 –254 (1997). 0740-3232 CrossRef
Zonios  G. et al., “Diffuse reflectance spectroscopy of human adenomatous colon polyps in vivo,” Appl. Opt.. 38, (31 ), 6628 –6637 (1999). 0003-6935 CrossRef
Hayakawa  C. K. et al., “Use of the delta-P1 approximation for recovery of optical absorption, scattering, and asymmetry coefficients in turbid media,” Appl. Opt.. 43, (24 ), 4677 –4684 (2004). 0003-6935 CrossRef
Sun  J. et al., “Influence of fiber optic probe geometry on the applicability of inverse models of tissue reflectance spectroscopy: computational models and experimental measurements,” Appl. Opt.. 45, (31 ), 8152 –8162 (2006). 0003-6935 CrossRef
Finlay  J. C., Foster  T. H., “Hemoglobin oxygen saturations in phantoms and in vivo from measurements of steady-state diffuse reflectance at a single, short source-detector separation,” Med. Phys.. 31, (7 ), 1949 –1959 (2004). 0094-2405 CrossRef
Hull  E. L., Foster  T. H., “Steady-state reflectance spectroscopy in the P3 approximation,” J. Opt. Soc. Ame. A. 18, (3 ), 584 –599 (2001). 1084-7529 CrossRef
Zonios  G., Dimou  A., “Modeling diffuse reflectance from semi-infinite turbid media: application to the study of skin optical properties,” Opt. Express. 14, (19 ), 8661 –8674 (2006). 1094-4087 CrossRef
Johns  M. et al., “Determination of reduced scattering coefficient of biological tissue from a needle-like probe,” Opt. Express. 13, (13 ), 4828 –4842 (2005). 1094-4087 CrossRef
Bevilacqua  F., Depeursinge  C., “Monte Carlo study of diffuse reflectance at source–detector separations close to one transport mean free path,” J. Opt. Soc. Am. A. 16, (12 ), 2935 –2945 (1999). 0740-3232 CrossRef
Amelink  A., Sterenborg  H. J., “Measurement of the local optical properties of turbid media by differential path-length spectroscopy,” Appl. Opt.. 43, (15 ), 3048 –3054 (2004). 0003-6935 CrossRef
Kanick  S. C. et al., “Measurement of the reduced scattering coefficient of turbid media using single fiber reflectance spectroscopy: fiber diameter and phase function dependence,” Biomed. Opt. Express. 2, (6 ), 1687 –1702 (2011). 2156-7085 CrossRef
Reif  R., A’Amar  O., Bigio  I. J., “Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media,” Appl. Opt.. 46, (29 ), 7317 –7328 (2007). 0003-6935 CrossRef
Bargo  P. R. et al., “In vivo determination of optical properties of normal and tumor tissue with white light reflectance and an empirical light transport model during endoscopy,” J. Biomed. Opt.. 10, (3 ), 034018  (2005). 1083-3668 CrossRef
Palmer  G. M., Ramanujam  N., “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: theory and validation on synthetic phantoms,” Appl. Opt.. 45, (5 ), 1062 –1071 (2006). 0003-6935 CrossRef
Kanick  S. C. et al., “Scattering phase function spectrum makes reflectance spectrum measured from Intralipid phantoms and tissue sensitive to the device detection geometry,” Opt. Express. 3, (5 ), 1086 –1100 (2012). 1094-4087 CrossRef
Vitkin  E. et al., “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun.. 2, (587 ), 1 –8 (2011). 2041-1723 CrossRef
Rajaram  N., Nguyen  T. H., Tunnell  J. W., “Lookup table–based inverse model for determining optical properties of turbid media,” J. Biomed. Opt.. 13, (5 ), 050501  (2008). 1083-3668 CrossRef
Kanick  S. C., Sterenborg  H. J. C. M., Amelink  A., “Empirical model description of photon path length for differential path length spectroscopy: combined effect of scattering and absorption,” J. Biomed. Opt.. 13, (6 ), 064042  (2008). 1083-3668 CrossRef
Talsma  A., Chance  B., Graaff  R., “Corrections for inhomogeneities in biological tissue caused by blood vessels,” J. Opt. Soc. Am. A. 18, (4 ), 932 –939 (2001). 0740-3232 CrossRef
Svaasand  L. O. et al., “Therapeutic response during pulsed laser treatment of port-wine stains: dependence on vessel diameter and depth in dermis,” Laser Med. Sci.. 10, (4 ), 235 –243 (1995). 1435-604X CrossRef
Cheong  W. F., Prahl  S. A., Welch  A. J., “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron.. 26, (12 ), 2166 –2185 (1990). 0018-9197 CrossRef
Graaff  R., “Similarity relations for anisotropic scattering in absorbing media,” Opt. Eng.. 32, (2 ), 244  (1993). 0091-3286 CrossRef
Kienle  A., Forster  F. K., Hibst  R., “Influence of the phase function on determination of the optical