1.

Introduction

The optical characterization of biological tissue is an important goal for many diagnostic and scientific applications. The multilayered property of most tissue types presents the need for a depth-resolved measurement in order to obtain layer-specific information about optical properties. A complete characterization of the optical properties in each layer would involve knowledge of the reduced scattering coefficient ( [TeX:] $\mu _s^{\bf *} $ ${\mu }_s^{*}$ ), absorption coefficient (μ_a), and the phase function. Perhaps the most reliable existing method for measuring optical properties is the integrating sphere approach. The integrating sphere technique allows for the measurement of [TeX:] $\mu _s^{\bf *} $ ${\mu }_s^{*}$ , μ_a, and the anisotropy factor (g). However, integrating sphere measurements require the samples to be optically thin and the results are subject to optical alterations due to the sample preparation. Furthermore, integrating sphere measurements are less robust for media with large values of g, which is typically the case for many types of soft tissue. Noninvasive measurements of the reduced scattering coefficient and absorption coefficient in intact tissue are most commonly accomplished with diffuse reflectance measurements. Although this method can assess the bulk optical properties of the sample ( [TeX:] $\mu _s^{\bf *} $ ${\mu }_s^{*}$ and μ_a in large tissue volumes), the long penetration depth of light in biological samples (several millimeters) prevents these measurements (e.g., diffuse optical tomography) from having an adequate depth resolution to accurately characterize most epithelial layers. Additionally, diffuse reflectance measurements cannot yield any information about the scattering phase function.

In recent years, several techniques have been demonstrated that can measure local optical properties of tissue. For example the exponential decrease in intensity due to scattering away from an incident beam can be measured with a confocal microscope and used to estimate the local reduced scattering coefficient.¹ However, this measurement does not yield information about the scattering phase function and is limited to a small field of view. Some examples of other highly depth-selective light scattering methods that have been used for the characterization of epithelial tissue include techniques involving polarization gating,^{2, 3, 4} spatial gating,⁵ and coherence gating.^{6, 7, 8} The calculation of optical properties from any of these methods involves the solution of an inverse problem in light scattering at sub-diffusion length scales, where the backscattering is dependent on properties of the scattering phase function in addition to the reduced scattering coefficient. These methods use simplifying assumptions because the large number of variables in more precise models exceeds the number of measured parameters. For example, Mie theory-based models typically require assumptions that determine the refractive index of the spheres and the surrounding medium as well as the size distribution of the spheres; in other words, some a priori knowledge of the phase function is required. Despite the lack of a robust tissue scattering model, optical spectroscopy of epithelial tissue has been demonstrated to have the potential for noninvasive cancer diagnostics,^{9, 10, 11, 12, 13, 14} an attractive alternative to the conventional biopsy and pathology approach. The promising results obtained in these studies underscore the importance of knowing the optical properties of the tissue under investigation. Therefore, the challenge for these novel diagnostic methods is to be able to convert light scattering measurements into optical properties of the scattering medium.

Low-coherence enhanced backscattering (LEBS) is a method that utilizes the self-interference of backscattered light under partial coherence illumination.^{15, 16} This reflectance method is not hindered by the measurement geometry limitations that are present in the integrating sphere technique and does not require any manipulation of the sample. LEBS is capable of measuring the reflectance distribution at small scattering distances (i.e., separations between incident and collected ray pairs well below the transport mean free path length or even the mean free path length) resulting in a high sensitivity to the phase function.¹⁷ Furthermore, the capability for spectroscopic measurements allow LEBS to measure optical properties such as μ_a with depth selectivity.¹⁸ In this work, we will present a methodology for how LEBS can be used to measure [TeX:] $\mu _s^{\bf *} $ ${\mu }_s^{*}$ and two parameters of the scattering phase function (g and m). LEBS spectroscopy has previously been shown to have promise for the detection of several types of epithelial cancers.^{14, 19, 20} These studies have shown that LEBS parameters such as the peak width, enhancement factor, and spectral dependence have diagnostic potential for cancer detection. The aim of this work is to relate the measured LEBS parameters (width, enhancement factor, and spectral dependence) to the optical properties of the scattering media. This is accomplished by applying a general light scattering model that is well suited for characterizing continuous random media.

The general layout is as follows: we will first review the methodology for performing numerical simulation of LEBS in Sec. 2. We will then discuss a theoretical light scattering model based on a weak scattering approximation and a general equation for the refractive index correlation function in Sec. 3. The dependencies between the optical properties of the scattering media and the LEBS peak properties are described in Sec. 4. In Sec. 5, we will present several approaches for solving the inverse problem based on this model. We then follow with an evaluation of each of these approaches and an error analysis in Sec. 6. Finally, in Sec. 7, we will demonstrate the validity of the approach with an experimental measurement of a tissue phantom composed of a mixture of microsphere suspensions and compare the results of the measured optical properties from LEBS with theoretically expected results. The conclusions are then summarized in Sec. 8.

2.

Monte Carlo Method for Simulation of LEBS

Monte Carlo simulations can be utilized to obtain reflectance distributions corresponding to varying optical properties. In this work, a publicly available Monte Carlo code²¹ was modified to use the phase function described in Sec. 3 [Eq. 4b]. The code was validated for the Henyey–Greenstein case against another publicly available code that utilizes the Henyey–Greenstein phase function.²² The two-dimensional probability distribution of backscatterd light resulting from an infinitely narrow collimated beam, p(r, φ), can be obtained from the Monte Carlo simulation.²³ In the absence of polarization (i.e., no φ-dependence), this distribution can be integrated to obtain a one-dimensional backscattering distribution P(r), where P(r) = ∫p(r, φ)rdφ = 2πrp(r). As required in a Monte Carlo simulation of LEBS, the backscattering distributions were collected at small angles (0 deg to 10 deg) in order to minimize the alteration to P(r) due to finite angular collection.²³ The effects of the finite angular collection were further reduced by storing the radial position of the final scattering event as opposed to the exiting radial position at the surface of the medium. The results of a Monte Carlo simulation were generalized to any value of [TeX:] $l_s^*$ $l_s^{*}$ ( [TeX:] $l_s^* = {1 / {\mu _s^* }}$ $l_s^{*}=1/{\mu }_s^{*}$ ) by rescaling all of the lengths used in the simulation, a well-known scaling property of Monte Carlo.^{24, 25} It is also convenient to present the backscattering distribution relative to [TeX:] $l_s^* $ $l_s^{*}$ in order to obtain convergence for all optical properties in the diffusion regime of [TeX:] $r \gg l_s^*$ $r\gg l_s^{*}$ . The absorption coefficient was assumed to be negligible and all exiting photon packets were assigned equal weights. Therefore, it is assumed that at least part of the LEBS peak measurement is obtained for a nonabsorbing wavelength range. A discussion of how absorption measurements can be obtained from an LEBS peak is presented elsewhere.¹⁸

The LEBS angular distribution can be calculated as the two-dimensional Fourier transform of the radial reflectance intensity distribution, [TeX:] $p(\vec r)$ $p\left(\vec{r}\right)$ :

Applying Eq. 1 to obtain the LEBS peak requires a two-dimensional Fourier transform. However, in cases of circularly symmetric peaks, a Henkel transform can be used. We implemented the fast Fourier transform algorithm in MATLAB to compute the Henkel transform, allowing for faster and more memory-efficient computation. We also applied interpolation to translate one-dimensional P(r) distributions into a two-dimensional grid according to the average position of each radial bin in P(r).³⁶ Similarly, the average position of the central bin in p(x, y) was calculated, while the positions of all other bins in p(x, y) were approximated to be the central point of the bin. The method of interpolation was validated by performing simulations with varying grid resolutions, but otherwise identical parameters. It was found that the error due to the absence of interpolation can greatly affect peak properties due to the sensitivity of the LEBS signal to small radial distances.

Earlier findings indicated that Monte Carlo simulations accurately predicted LEBS peak properties from microsphere phantoms when the output Monte Carlo data is multiplied by a factor of 0.5.¹⁷ The Monte Carlo data presented here were likewise scaled by the same factor in order to obtain excellent agreement with experimental results.

3.

Theoretical Model Based on the Whittle–Matérn Correlation Function

The Monte Carlo method provides an exact solution to the radiative transfer equation given that the correct phase function and geometry is implemented and an adequate number of rays are traced such that the signal is not affected by numerical noise. For the case of spherical particle scattering, the Mie theory provides an exact prediction of the phase function. However, it is arguable that spheres are not an adequate approximation to scattering in tissue. The commonly used Henyey–Greenstein phase function, adapted from the light scattering of interstellar dust particles,³⁷ offers an attractively simple model in which the first moment of the phase function is determined by the anisotropy factor g. The disadvantage of this phase function is that its higher order moments are pre-determined by the first moment as g _n = g ⁿ, where g ₁ = g = 〈cos θ〉. In reality, higher order moments for most tissue phase functions are very different from those of the Henyey–Greenstein phase function³⁴ not giving the phase function enough flexibility to accurately describe measurements from many types of biological samples.^{38, 39} Therefore, a more flexible phase function would present a more realistic characterization of tissue. Fortunately, LEBS is sensitive to higher order moments of the phase function than g ₁. For instance, two phase functions that have the same value of g but different shapes (e.g., Henyey–Greenstein phase function with g = 0.9 versus Mie phase function with g = 0.9), result in different LEBS peak widths.¹⁷ This is because the shape of the phase function (i.e., higher order moments), not just the anisotropy factor, plays an important role in determining the backscattering at small length-scales. We will therefore describe a general model that results in a two-parameter phase function where the shape of the phase function and g can be altered independently.

Recent work has shown that statistical models of the refractive index fluctuations can be well suited for tissue characterization.^{40, 41, 42} These models involve the characterization of the refractive index correlation function, which is then used to predict the scattering properties of the medium. In particular, the correlation function that results from a fractal distribution of length scales has been shown to be of relevance for tissue scattering.^{41, 43} A fractal distribution of sizes results in a power law decay of the refractive index correlation function and a power law dependence of the backscattering intensity as a function of wavelength (λ). The power law dependence of scattering intensity as a function of wavelength has been observed by several authors and often used to characterize tissue scattering.^{35, 38} Sheppard presented a model for scattering that utilizes the Born approximation and a general form of the refractive index correlation function.⁴³ Rogers then extended this model to incorporate the Henyey–Greenstein phase function and the phase function corresponding to the mass fractal scattering regime. The model also includes phase functions resulting from other commonly assumed types of refractive index correlations including the stretched exponential, exponential, and Gaussian.

This general model involves a minimal number of assumptions and is well suited for characterizing continuous turbid media such as tissue. The refractive index correlation function is written as a dependence on three sample properties: the variance of the refractive index (Δn ²), the correlation length (l _c), and an additional parameter that determines the functional form of the correlation function (m). A weak scattering (Born) approximation then results in a two-parameter phase function, where the parameters l _c and m determine the shape of the phase function.^{43, 44} It is possible to achieve variations in the shape of the phase function while maintaining the same value of g, resulting in a more flexible model than the Henyey–Greenstein phase function. This flexibility is important for the characterization of superficial tissue as superficial scattering determines the reflectance at small length scales and the phase function plays a significant role in this regime.

In the Born approximation, all scattering properties (e.g., the differential scattering cross section and its derivative quantities such as the scattering coefficient and the phase function) of a medium with a continuously varying refractive index are completely determined through a Fourier transform relationship that involves the refractive index correlation function B _n(r). Because optical refractive index n is a linear function of the local mass density (n = n _water + αρ, where ρ is the portion of tissue solids by volume and refractive index increment α ∼ 0.17 to 0.2), B _n(r) is proportional to the density correlation function. There are two assumptions to the model: one needs a comprehensive description of the refractive index (or density) correlation function pertinent to tissue and the first order Born approximation has to be valid.

First, the refractive index correlation function is assumed to have a general form that is described by the Whittle–Matérn family of correlation functions:⁴⁵

The second assumption is the validity of the Born approximation. It has been shown that in case of a medium with continuous refractive index fluctuations, the approximation is valid if the refractive index fluctuations are small, relative to the correlation length:⁴⁶

Fig. 1

Whittle–Matérn phase functions for m = 3/2 and varying g (a) and g = 0.9 and varying m (b). A higher value of g results in and increase in the probability of forward scattering and a decrease in the probability of backscattering. A smaller value of m results in an increase in both the probability of forward and backward scattering.

4.

Dependence of LEBS Peak Properties on the Optical Properties of the Scattering Medium

The LEBS peak can be parameterized into three measurements: the peak width, enhancement factor, and spectral slope. In this section, we will define the width of the peak as the full width at half maximum at a given wavelength (a second useful definition will be introduced later). The enhancement factor is defined as the height of the peak relative to an isotropic semi-infinite reflectance standard (Ocean Optics Inc.). The reflectance standard serves to normalize the measured angular profile by the total reflected intensity per unit solid angle and unit wavelength. The solid angle and wavelength are determined by the angular and spectral resolution, respectively, of the detection arm. The LEBS spectral slope is defined as the rate of change in the enhancement factor with respect to wavelength. These parameters have been found to have diagnostic value for the detection of colon, pancreatic, and skin cancer.^{20, 47, 48, 49} In this section, we will discuss the dependence of these parameters on the optical properties [TeX:] $l_s^* $ $l_s^{*}$ , g, and m. Several inverse models for the calculation of these optical properties are presented in Sec. 5.

The LEBS signal is related to the Fourier transform of the backscattering probability distribution [Eq. 1] and some general understanding about LEBS can be attained by observing the dependence of the alterations of P(r) due to alterations in the phase function. Figures 2a, 2b show backscattering distributions, [TeX:] $P(r/l_s^*)$ $P(r/l_s^{*})$ , for the phase functions displayed in Figs. 1a, 1b, respectively. We can observe that several features of the phase function correlate with features at small length scales of P(r). First, the probability of forward scattering is correlated with the initial peak in P(r) that is present for values of r near 0. As g increases, the forward scattering peak of the phase function becomes higher, resulting in an increased peak in P(r) at small r [Fig. 2a]. Additionally, the peak at small r also increases when m decreases [Fig. 2b], also correlating with the forward scattering peak in the phase function [Fig. 1b]. Note that the values shown in the P(r) curves are dimensionless and are cases where absorption is neglected. The measurement of absorption from LEBS is a separate topic that is described elsewhere.¹⁸ This results in ∫P(r)dr = 1 for any value of [TeX:] $l_s^*$ $l_s^{*}$ , constraining the values of P(r) at [TeX:] $r < l_s^*$ $r<l_s^{*}$ to be approximately inversely proportional to [TeX:] $l_s^*$ $l_s^{*}$ . In other words, a medium with a longer [TeX:] $l_s^*$ $l_s^{*}$ will have a broader and shorter P(r) with identical area.

The relationship between the backscattering distribution P(r) and the phase function F(θ) can be understood by observing how the phase function is sampled for small exit radii ( [TeX:] $r \ll l_s^*$ $r\ll l_s^{*}$ ). In other words, we can observe the average distribution of scattering angles per backscattered ray that exits at a given radius. We therefore define H(θ) as the angular distribution of all scattering events that results in light exiting at distance r [as compared to the phase function F(θ) that is the angular distribution of a single scattering event]. H(θ) describes how the phase function is being sampled at the collection radius of interest. Figure 2c shows H(θ) at three exit radii for the case of g = 0.9 and m = 3/2. Interestingly, the probability of scattering into forward and backward directions is increased by approximately 2 orders of magnitude for [TeX:] $r = {{l_s^* }/ {100}}$ $r=l_s^{*}/100$ as compared to side scattering into angles around 90 deg. We can infer that the phase function is primarily sampled in the forward and backward directions for small radial distance and forward scattering media. Therefore, the height of the backscattering peak at small values of r in Figs. 2a, 2b can be related to the height of the phase function in the forward and backward directions (Fig. 1).

The enhancement factor is defined as the height of the LEBS peak at θ = 0. This value can be calculated from Eq. 1 as ∫P(r)c(r)dr. Therefore, for [TeX:] $l_s^* > L_{{\rm sc}}$ $l_s^{*}>L_{\mathrm{sc}}$ [i.e., P(r) extends much farther than c(r)], we would expect the enhancement factor to be approximately proportional to [TeX:] $L_{{\rm sc}} \mu _s^* $ $L_{\mathrm{sc}}{\mu }_s^{*}$ . The deviation to this approximate proportionality would be due to fluctuations in P(r) at r < L _sc, but we can neglect these fluctuations as a first-order approximation. We can therefore estimate the enhancement factor to be proportional to [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ for [TeX:] $L_{{\rm sc}} \mu _s^* < 1$ $L_{\mathrm{sc}}{\mu }_s^{*}<1$ , such that

Fig. 2

(a) and (b) Backscattering distributions obtained from Monte Carlo simulations at small radial distances corresponding to the phase functions shown in Fig. 1. (c) Scattering angle distribution for rays at three exit radii ( [TeX:] $r = {{l_s^* } / {100}}$ $r=l_s^{*}/100$ , [TeX:] $r = {{l_s^* } /{10}}$ $r=l_s^{*}/10$ , and [TeX:] $r = l_s^*$ $r=l_s^{*}$ ) for a scattering media with g = 0.9 and m = 3/2. Note that the distributions in (c) differ from the actual phase function because each backscattered ray has an equal contribution (as opposed to each scattering event having an equal contribution). The distribution of scattering angles is heavily weighted toward forward and backward scattering events for [TeX:] ${r / {l_s^* }} \ll 1$ $r/l_s^{*}\ll 1$ , indicating that a large majority of rays collected at small radii undergo primarily forward and backward scattering.

The LEBS peak width (W) is also related to the alterations of p(r) at small length scales. The width of the two-dimensional Fourier transform of p(r, φ)c(r, φ) is inversely proportional to the width of p(r, φ) at small r as well as the width of c(r, φ). Assuming a circular aperture as a secondary source and defining α as the angular extent of the source, we can relate the width of c(r), defined as L _sc, to α with the van Cittert–Zernike theorem: α = λ/(2πL _sc).²⁶ As a result, α/W is a dimensionless parameter that only depends on [TeX:] $L_{{\rm sc}} \mu _s^* $ $L_{\mathrm{sc}}{\mu }_s^{*}$ and the phase function. For [TeX:] ${{l_s^*} / {L_{{\rm sc}} }} \gg 1$ $l_s^{*}/L_{\mathrm{sc}}\gg 1$ , α/W decreases with increasing g [Fig. 4a] and increases with increasing m [Fig. 4b]. The width of the LEBS peak is primarily sensitive to g and m for [TeX:] ${{l_s^* } / {L_{{\rm sc}} }} > 10$ $l_s^{*}/L_{\mathrm{sc}}>10$ and has a weak dependence on [TeX:] $l_s^*$ $l_s^{*}$ . This is again the regime where the coherence function is isolating the first sharp peak in P(r) that is located at [TeX:] $r \ll l_s^*$ $r\ll l_s^{*}$ .

The spectral LEBS measurements can be utilized to obtain additional information about the optical properties of the medium. The scattering coefficient, μ_s, and the anisotropy coefficient, g, can be calculated from the Whittle–Matérn differential scattering cross section by following the procedure detailed by Rogers ⁴⁴ for the scalar wave case (neglecting the dipole factor). We obtain the following equation for g:

For kl _c ≫ 1, Eq. 6 has three regimes of behavior depending on the value of m:

Similarly, we can also express the kl _c dependence for the scattering coefficient:

The most relevant range of m for tissue scattering is 1 < m < 2 because it includes fractal scattering with mass fractal dimensions between 2 and 3, the Henyey–Greenstein phase function (m = 3/2), as well as stretched exponential correlation functions. From Eqs. 7, 8, we can obtain a dependence for the reduced scattering coefficient for the tissue scattering regime of kl _c ≫ 1 and 1 < m < 2,

The proportionality factor f(m) in Eq. 9 will be an m dependent function, with no dependence on k, l _c, or Δn ². By using Eq. 9 in Eq. 5 and employing the van Cittert–Zernike theorem to determine the spatial coherence length as L _sc = (kα)⁻¹, we can obtain a spectral dependence of the LEBS enhancement factor.

We will first assume that C _E can be approximated to be a constant (about 0.2 for [TeX:] $L_{{\rm sc}} \mu _s^* \approx 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\approx 1$ ). Then, we will extend the analysis to include all values [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ . Assuming that Δn ² does not significantly vary with the wavelength (which is a generally safe assumption considering the wavelength variation of refractive index for most biologically relevant materials), m can be measured from the LEBS spectral slope. If a relatively narrow range of wavelengths are used, the slope approximates the tangent to the E(λ) dependence and can be equated to the derivative of E(λ) with respect to λ.

We will refer to the quantity [TeX:] $S_n = \frac{{d[E(\lambda)]}}{{d\lambda }}\frac{{\left\langle \lambda \right\rangle }}{{\left\langle E \right\rangle }}$ $S_n=\frac{d\left[E\right(\lambda \left)\right]}{d\lambda }\frac{〈\lambda 〉}{〈E〉}$ as the normalized spectral slope.

In obtaining Eq. 12, we have assumed the proportionality factor C _E does not vary with optical properties. However, this assumption is not valid for [TeX:] $L_{{\rm sc}} \mu _s^* \ll 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\ll 1$ , as shown in Fig. 3. Therefore, it is important to understand how C _E depends on optical properties in this regime and how this variation will impact the spectral slope. This is accomplished by performing a series of Monte Carlo simulations with the Whittle–Matérn phase function and varying kl _c and m for [TeX:] $L_{{\rm sc}} \mu _s^* \ll 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\ll 1$ . The resulting values of C _E are plotted as a function of kl _c in Fig. 5a in log–log scale. Note that for a constant m, kl _c is directly related to [TeX:] $\tilde g$ $\tilde{g}$ through Eq. 3. In log–log scale, C _E has linear dependence on kl _c for kl _c > 1. This linearity indicates a power law behavior, where the power is determined by the slope of the curve and varies with m. A deviation from the power law dependence is only seen at kl _c < 1, where the value of C _E approaches the case of isotropic scattering (we will call this value C ₀). We will approximate the power law dependence on kl _c with a 2nd order polynomial function of m:

Fig. 3

Enhancement factor dependence on optical properties. The enhancement factor appears to be approximately proportional to [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ and weakly dependent on g (a) and m (b). However, when the proportionality coefficient [ [TeX:] $C_E = {E / {({L_{{\rm sc}} \mu _s^* })}}$ $C_E=E/\left(L_{\mathrm{sc}}{\mu }_s^{*}\right)$ ] is plotted, variations with g (c) and m (d) become apparent. These variations are attributed to the variation in the sharp peak of P(r) at small radial distances.

The dependence of the normalized spectral slope from Eq. 15 is shown in Fig. 5d for three cases of [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ . When [TeX:] $L_{{\rm sc}} \mu _s^* \sim 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\sim 1$ , the normalized spectral slope is 2m − 3 because S ≈ 0. When [TeX:] $L_{{\rm sc}} \mu _s^* \ll 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\ll 1$ , the spectral slope has a quadratic dependence on m, but is well approximated by 2m − 3 in the range of 0.7 < m < 1.5. This range of m is within the mass fractal regime of the Whittle–Matérn correlation function.

The changes in the LEBS peak enhancement, width, and spectral slope can be summarized as being related to optical properties. The enhancement factor is approximately proportional to [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ . The width is related to both g and m. The width increases as g increases and decreases with increasing m. Finally, the spectral slope is related to 2m − 3 in the mass fractal regime where 0.7 < m < 1.5. The spectral slope is also related to m outside of this range, but this relationship also depends on a, b, and c, from Eq. 15, which in turn depend on [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ .

5.

Inverse Models for Optical Property Measurements

There are several options for building an inverse model that can translate measured properties of the LEBS peak into the optical properties (g, [TeX:] $l_s^*$ $l_s^{*}$ , and m) or physical properties (l _c, Δn ², and m) of the scattering medium. In the following discussion, we will present these methods for solving the inverse problem to calculate optical properties as the equations presented in the preceding sections and elsewhere⁴⁴ can be used to convert between optical and physical properties.

One class of inverse models for the prediction of LEBS relies on a predictive model of P(r). The prediction of P(r) at [TeX:] $r \ll l_s^*$ $r\ll l_s^{*}$ is made possible due to several factors that simplify the inverse problem. First, absorption is assumed to be negligible. This is justifiable because the LEBS signal can be measured at a wavelength where tissue absorption is much less than the scattering (i.e., λ > 600nm). Furthermore, the spectral signal due to hemoglobin absorption can be removed by fitting with the hemoglobin absorption spectrum. Without the presence of absorption, the only length scale that can alter P(r) is l _s. For a nonabsorbing medium, it is also known that in the diffusion regime where [TeX:] $r \gg l_s^*$ $r\gg l_s^{*}$ , [TeX:] $l_s^*$ $l_s^{*}$ is the parameter that determines the scattering signal. By scaling the radius as [TeX:] ${r / {l_s^* }}$ $r/l_s^{*}$ , we remove the l _s dependence and also obtain the conversion of the backscattering signal for varying properties of the phase function at [TeX:] $r \gg l_s^*$ $r\gg l_s^{*}$ . Thus, the inverse model must account for the variations of [TeX:] $P({r / {l_s^* }})$ $P\left(r/l_s^{*}\right)$ at small [TeX:] ${r / {l_s^* }}$ $r/l_s^{*}$ due to changes in the phase function. One convenient approach to predict the changes in the phase function is to quantify the difference between isotropic [TeX:] $P({r / {l_s^* }})$ $P\left(r/l_s^{*}\right)$ , where g = 0 and nonisotropic [TeX:] $P({r / {l_s^* }})$ $P\left(r/l_s^{*}\right)$ . We have recently demonstrated how this can be done by taking the difference between g = 0.9 and g = 0 backscattering for two values of m.⁵¹ The backscattering probability can be predicted according to:

The relationship between P(r) and the three optical properties [TeX:] $l_s^*$ $l_s^{*}$ , g, and m is unique. Therefore, in principle, it should be possible to obtain optical properties by using one of the P(r) predictive models described above to calculate the LEBS peak [with Eq. 1] and thereby fit optical properties such that the observed angular dependence of the LEBS peak has minimal error with the LEBS peak predicted by the model. However, in cases where [TeX:] $L_{{\rm sc}} < l_s^*$ $L_{\mathrm{sc}}<l_s^{*}$ , the relationship between the three optical properties and the shape of the peak becomes nonunique, due to a finite measurement uncertainty. Figure 6 shows three LEBS peaks that are simulated for three different combinations of optical properties and an L _sc of 25 μm. The shapes of the peaks are indistinguishable, suggesting that there can only be two reliable parameters when fitting the angular dependence of the peak: the height, and width. The height of the peak can be fit by adjusting [TeX:] $l_s^*$ $l_s^{*}$ , and the width of the peak can be fit by adjusting either g or m. Fortunately, m can be independently obtained from the spectral dependence based on the analysis described in Sec. 4. The measurement of m can be obtained from the normalized spectral slope according to Eq. 12. As described above, this expression for m is accurate for the range 0.7 < m < 1.5. If the resulting value of m is outside of this range, an iterative process can be implemented by first estimating a value of m based on the normalized spectral slope and then adjusting it after [TeX:] $l_s^*$ $l_s^{*}$ is calculated based on the fit with the angular dependence of the peak. In this case, m would be calculated by using the expression in Eq. 15.

Fig. 6

LEBS peaks predicted from the PCA model for three sets of optical properties and an L _sc of 25 μm. The three LEBS peaks are indistinguishable, implying that one of the three parameters must be fixed in order to obtain a unique solution. For a fixed m, the height of the peak can be fit by adjusting [TeX:] $l_s^*$ $l_s^{*}$ and the width of the peak can be fit by adjusting g.

Another approach to obtain optical properties from an LEBS peak is through the construction of a lookup table that relates peak properties such as the LEBS width, enhancement, and spectral slope to the optical properties of the scattering medium. Lookup tables have the advantage of being computationally fast because the steps involved in minimizing the error and performing other intermediate calculations, mainly the two-dimensional Fourier transform, are performed ahead of time during the construction of the table. There are three optical properties ( [TeX:] $l_s^*$ $l_s^{*}$ , g, and m) and three measured parameters [enhancement (E), width (W), and normalized spectral slope (S _n)] making the process analogous to solving a system of three equations. An example procedure for implementing a lookup table is described in Fig. 7. Here, the problem is simplified because m can be independently measured from S _n. Two lookup tables are therefore needed, including one for Eand one for W, each being a function of [TeX:] $l_s^*$ $l_s^{*}$ , g, and m. After m is measured from the spectral slope, each of the tables would be interpolated into two-dimensional matrices that correspond to the obtained value of m, as shown in Fig. 7a for E and Fig. 7b for W. The enhancement factor table is then used to calculate [TeX:] $l_s^*$ $l_s^{*}$ as a function of g. The values of W corresponding to the obtained values of [TeX:] $l_s^*$ $l_s^{*}$ and g are then obtained from the width table. Finally, the measured value of the width is used to obtain g from the W(g) dependence. The measured value of g is then used to calculate [TeX:] $l_s^*$ $l_s^{*}$ from the dependence of [TeX:] $l_s^*$ $l_s^{*}$ as a function of g obtained earlier. A block diagram of the procedure is presented in Fig. 7c.

Fig. 7

(a) Cross section of E lookup table at m = 1.3. (b) Cross section of W lookup table at m = 1.3. (c) Flow chart of lookup table procedure for calculating optical properties from an LEBS peak measurement.

6.

Error Analysis of Inverse Models for Measuring Optical Properties with LEBS

We evaluated the three methods presented in Sec. 5 for calculating optical properties from LEBS by employing the Monte Carlo method for simulating the LEBS signal. A Monte Carlo simulation was used to simulate P(r), as described in Sec. 3. The LEBS signal was then calculated with a Fourier transform according to Eq. 1. The coherence length, L _sc, was varied between 25 and 173 μm, and the anisotropy factor g was varied between 0.2 and 0.95. [TeX:] $l_s^*$ $l_s^{*}$ and m were maintained at 1000 μm and 1.3, respectively. The chosen values of [TeX:] $\,l_s^*$ $l_s^{*}$ and m are within a range of properties typically observed from tissue.^{40, 43} The simulated LEBS peaks were interpolated to an angular resolution of approximately 0.0115 deg. This resolution corresponds to the experimental system used for validation. Details on the experimental setup were presented in Sec. 7. It was assumed that m could be accurately measured from the spectral dependence using the analysis presented in Sec. 4, and therefore the error in the measurement of m was not evaluated. Figure 8a shows the resulting error in the obtained value of [TeX:] $l_s^*$ $l_s^{*}$ using the difference model of P(r) for varying values of g and L _sc and Fig. 8b shows the error in the obtained value of g. This model results in a slight overestimation of g and underestimation of [TeX:] $l_s^*$ $l_s^{*}$ . The principle component analysis (PCA) method results in a more accurate prediction of P(r) and therefore allows for more accurate measurements of optical properties. The error in [TeX:] $l_s^*$ $l_s^{*}$ is reduced to <3% [Fig. 8c] and accuracy of the measurement of g is also improved [Fig. 8d]. The coefficients in both of these predictive models were optimized for g ⩾ 0.6 in order to obtain optimal performance for tissue-relevant values of g causing some sacrifice in the accuracy of the measurement for small values of g. On the other hand, the measurement of [TeX:] $l_s^*$ $l_s^{*}$ is relatively insensitive to the value of g and remains accurate even for small values of g.

The iterative process of calculating the Fourier transform in Eq. 1 for varying optical properties to minimize the error between the predicted and measured LEBS peak results in a lengthy computation time for both the PΔ and PCA methods. By comparison, a lookup table stores several relevant measurements from the LEBS peak and is much faster for optical property calculations after the table has been generated. We implemented two variations of a lookup table algorithm. The first scheme stored the enhancement factor (E) and LEBS peak width (W), defined as the full width at half maximum, each being a function of [TeX:] $l_s^*$ $l_s^{*}$ , g, and m. Figure 9a shows the percentage error in [TeX:] $l_s^*$ $l_s^{*}$ for varying values of g and L _sc obtained with this method. The error is minimal for small values of L _sc and increases up to a maximum of approximately 4% at the longest L _sc of 173 μm. Similarly, a high degree of accuracy is achieved in the measurement of g, shown in Fig. 9b, for small L _sc. This is because an increase in L _sc results in the peak width approaching the behavior that is observed in CBS. In CBS the peak width is proportional to [TeX:] ${\lambda / {l_s^* }}$ $\lambda /l_s^{*}$ where the behavior is well described by the diffusion approximation.³² In this regime, the peak width is insensitive to scattering characteristics of the phase function such as the anisotropy factor. Our second implementation of the lookup table overcomes this limitation by using a modified definition of the peak width. The second table uses a width measurement which is defined as the LEBS intensity at 0.5 deg relative to the enhancement factor (W _RI = I _LEBS(0.5 deg)/I _LEBS(0 deg). This definition of the peak width is more convenient to measure in an in vivo probe configuration⁵² where a small number of fibers are used to sample the angular backscattering signal. Furthermore, it is known that the periphery of the enhanced backscattering peak is sensitive to shorter path lengths and lower order scattering events.⁵⁰ Thus, the limitation of approaching diffusive scattering at longer coherence lengths can be overcome by adopting this new definition of the peak width. Figures 9c, 9d show the error in [TeX:] $l_s^*$ $l_s^{*}$ and g, respectively, obtained with the lookup table method that utilizes W _RI. There is virtually no difference in error for varying L _sc and g.

7.

Experimental Validation of Optical Property Measurements

The experimental setup used to obtain LEBS measurements has been extensively described in earlier work.^{17, 18, 51} In brief, broadband continuous-wave light from a 450 W Xenon lamp (Oriel Instruments) is collimated and spatially filtered with a high power circular aperture (Lenox Laser) in a 4-f lens system. The aperture and focal length of the collimating lens determine the L _sc of the illumination beam. It was ensured that L _sc was significantly smaller than the beam diameter such that the LEBS peak properties were solely determined by the coherence function c(r) and the sample properties. The exact backscattered light is reflected from a plate beam splitter, passes through a liquid crystal tunable filter for wavelength selection, and imaged via a Fourier lens with a CCD camera. Five measurements were obtained from the phantom with an L _sc of 173 μm and at a wavelength of 680 nm. Unpolarized illumination was obtained by removing the polarizer in the illumination path. An image of the experimentally obtained peak is shown in Fig. 10. Figure 10a shows the LEBS measurement relative to the average reflectance standard intensity. Although the LEBS peak is clearly visible, there is unevenness in the baseline due to artifacts in the CCD camera and the incident angle of the illumination beam onto the sample. These experimental artifacts were measured by removing the pinhole from the 4-f system, resulting in an L _sc of less than 10 μm, and using crossed polarizers between the illumination and collection. The flat field measurement of these isolated artifacts is shown in Fig. 10b. The corrected LEBS peak was then obtained by dividing by the flat field and subtracting the baseline [Fig. 10c].

Fig. 10

Illustration of procedure for removing artifacts in the baseline from the LEBS peak. The uncorrected measurement is shown in (a). In (b), a flat field measurement is shown that is obtained by using a very short L _sc and crossed illumination and collection polarizers. (c) Shows the corrected LEBS peak, obtained by dividing the raw measurement in (a) by the flat field in (b) and subtracting the baseline.

In order to validate the measurement of optical properties with LEBS, we constructed a phantom with a known scattering phase function. This was accomplished by combining polystyrene microsphere suspensions of varying sizes. Mie scattering theory was used to predict the phase function from each individual particle size. The phase function of the mixture of particles was predicted with a weighted average of the differential scattering cross sections of each of the constituent microspheres within the suspension. A mixture of 10% by weight solutions of 0.20, 0.82, and 1.3 μm diameter polystyrene microspheres (Thermo Scientific, Fremont, California) was combined in relative volume proportions of 0.36:0.27:1, respectively. The volume proportions were obtained by minimizing the square error between the predicted phase function and the Whittle–Matérn phase function for m = 1.3 (mass fractal dimension D _mf = 2.6) and g = 0.9. The resulting phantom phase function had an anisotropy factor of 0.899 ± 0.002 (mean ± Std. Dev.) after accounting for the variability in size according to the manufacturer specifications. The suspension was then diluted with deionized water to the desired value of [TeX:] $l_s^*$ $l_s^{*}$ (668 ± 33 μm). The expected range in [TeX:] $l_s^*$ $l_s^{*}$ is due to the tolerance on the bead concentration that is provided by the manufacturer. These scattering properties are in the range of typically reported values for soft tissue. In principle, a large variety of phantoms with controlled optical properties can be constructed by combining microsphere suspensions in this manner. The validity of the Born approximation for the constructed phantom was verified by applying Eq. 3 with the largest microsphere diameter serving as an overestimate of the correlation length and the variance of the refractive index fluctuations being obtained from the known volume fraction of polystyrene. The result is that the value of Δn ²(kl _c)² < 0.06 ≪ 1 for the phantom, which justifies the validity of the Born approximation.

Figure 11 shows the results from the microsphere phantom described above. The phantom phase function is compared to the desired phase function in Fig. 11a. Although some oscillations can be seen in the phantom phase function, the fit is excellent over several orders of magnitude, including the forward and backward scattering directions. Figure 11b shows the average measured LEBS peak compared to the Monte Carlo simulated peak. We used the lookup table presented in Sec. 6 to calculate the [TeX:] $l_s^*$ $l_s^{*}$ and g of the phantom from the measured value of E and W _RI. An [TeX:] $l_s^*$ $l_s^{*}$ of 687 ± 13 μm (mean ± 95% CI) and a g of 0.91 ± 0.04 (mean ± 95% CI) were obtained from the lookup table. The confidence intervals for both [TeX:] $l_s^*$ $l_s^{*}$ and g include the predicted values of these optical properties thus indicating an accurate measurement. The obtained measurements of g and [TeX:] $l_s^*$ $l_s^{*}$ are shown as box plots in Figs. 11c, 11d, respectively. The expected range according to the manufacturer specifications and the obtained 95% confidence intervals are also included as horizontal dashed and dotted lines.

Fig. 11

Tissue phantom for measurement of optical properties. (a) Shows a comparison of the target Whittle–Matérn phase function (solid red line) and the fit obtained by combining microspheres of three different sizes (dashed blue line) for g = 0.9 and m = 1.3. (b) Shows the resulting LEBS peak obtained from the experimental measurement of the microsphere phantom. (c) and (d) Box plots of the g and [TeX:] $l_s^*$ $l_s^{*}$ values, respectively, obtained from the lookup table with an assumed value of 1.3 for m. The horizontal red dashed lines represent the expected range of values according to the manufacturer specifications. The horizontal blue dotted lines represent the upper and lower 95% confidence interval on the mean. The horizontal solid lines in the box plot represent successive quartiles, with the solid red line corresponding to the median.

Although the phase function obtained from a mixture of microspheres [Fig. 11a] is in agreement with the Whittle–Matérn phase function for m = 1.3 and g = 0.9, it can be misleading to assign the parameters l _c and m from the Whittle–Matérn phase function to the microsphere phantom without taking the wavelength dependence into account. The values of l _c and m can only be assigned by ensuring an accurate fit to the wavelength-dependent variations in the scattering cross section. For example, this can be done by calculating [TeX:] $\mu _s^*$ ${\mu }_s^{*}$ as a function of λ and obtaining m from Eq. 12 as well as fitting the spectral dependence of the scattering phase function. In the presented results, the minimization was not done for a range of wavelengths, therefore, we have limited the phantom measurements to [TeX:] $l_s^*$ $l_s^{*}$ and g for an assumed value of m = 1.3.

8.

Discussion and Conclusions

We have presented a theoretical framework for applying LEBS to measure optical properties from tissue. The model involves applying the Born approximation to the Whittle–Matérn correlation function. The resulting two-parameter phase function encompasses the often used Henyey–Greenstein phase function (m = 3/2) as well as the mass fractal scattering regime (m < 3/2). The anisotropy coefficient g determines the average cosine of the phase function while the parameter m can independently control the shape of the phase function. The relationships that link LEBS peak properties (e.g., E, W, and S _n) to the scattering optical properties of the medium (e.g., [TeX:] $l_s^*$ $l_s^{*}$ , g, and m) are determined by applying the Whittle–Matérn phase function in a series of numerical LEBS experiments. E/L _sc is approximately proportional to [TeX:] $\mu _s^*$ ${\mu }_s^{*}$ (Fig. 3), W/α primarily depends on g and m, and only weakly depends on [TeX:] $\mu _s^*$ ${\mu }_s^{*}$ (Fig. 4), and S _n is a direct measure of m (Fig. 5).

Fig. 4

Dependence of LEBS peak width on optical properties [TeX:] $l_s^*$ $l_s^{*}$ , g, and m. The inverse width decreases as g increases (a) and increases as m increases (b). The width of the peak is scaled to the angular extent of the source, which is related to the spatial coherence length: α = λ/(2πL _sc).

Fig. 5

Effect of varying proportionality coefficient on the normalized spectral slope. In (a), the proportionality coefficient C _E is plotted with kl _c in log–log scale for three values of m. The lines represent linear fits for kl _c > 1. (b) Shows a least squares fit for the coefficients a, b, and c from Eq. 13 along with the parameter C ₀. The coefficient values decay to zero as [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ approaches 1 and the normalized spectral slope therefore approaches 2m − 3. (c) Evaluation of the accuracy of the 2nd order polynomial at predicting the power S from Eq. 13. The value of C _Ethat is predicted according the Eq. 13 closely mimics the coefficient obtained from Monte Carlo data for kl _c > 1. (d) Graph of the normalized spectral slope versus m from Eq. 15. The normalized spectral slope has a quadratic behavior for small [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ but is well approximated by the 2m − 3 relationship in the range of 0.7 < m < 1.5. For [TeX:] $L_{{\rm sc}} \mu _s^* \sim 1$ $L_{\mathrm{sc}}{\mu }_s^{*}\sim 1$ , the normalized spectral slope exactly follows the 2m − 3 relationship.

As shown in Sec. 4, m can be directly measured from the spectral slope in most tissue-relevant cases according to m = (S _n + 3)/2. This measurement is accurate for the approximate range of 0.7 < m < 1.5 [Fig. 5d], which is a subset of the mass-fractal regime (with the mass fractal dimension D _mf = 2m). In other words, the mass fractal dimension can be measured from the LEBS spectral slope according to D _mf = S _n + 3 for the range 1.4 < D _mf < 3. This result is obtained by quantifying the value of the proportionality coefficient C _E between the enhancement factor and [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ for varying values of kl _c and m. For kl _c < 1, C _E is constant with a value of approximately 0.2. For kl _c > 1, C _E has a power law dependence on kl _c, with the power depending on m and [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ [Fig. 5a]. For a given value of [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ , the power on the wavelength dependence is shown to exclusively depend on m [Fig. 5c]. This dependence is well approximated as quadratic. However, the range of 0.7 < m < 1.5 is also well modeled with the linear 2m − 3 dependence [Fig. 5d]. The results presented in Fig. 5 show that the spectral dependence relates directly to the mass fractal dimension even for very short coherence lengths in which the enhancement factor is no longer proportional to the reduced scattering coefficient. Furthermore, Fig. 5b shows that the constants a, b, and c converge to zero for larger [TeX:] $L_{{\rm sc}} \mu _s^*$ $L_{\mathrm{sc}}{\mu }_s^{*}$ , meaning that the exponent S from Eq. 13 becomes 0 and the enhancement is again proportional to λ^{2m − 3}. This indicates that the spectral shape of the LEBS peak is not affected even for longer coherence lengths where the proportionality coefficient on the enhancement factor begins to decrease [Fig. 3b].

Three types of inverse models for calculating optical properties of a scattering medium are described in Sec. 5. The first two models predict the shape of the backscattering distribution at all length scales. Both of these models predict the deviations from isotropic scattering by modeling the difference between nonisotropic and isotropic backscattering distributions. The first method (PΔ model) utilizes only two distributions that are calculated as the difference between isotropic backscattering and the backscattering from a chosen optical property. Two optical properties are chosen for optimal accuracy in the tissue-relevant range and can predict the backscattering distribution for any g and m, according to Eq. 16. The second method captures the difference between the backscattering distribution from a given optical property and an isotropic medium by implementing a principle component analysis. This results in a model that uses slightly more stored data but can obtain better accuracy. The third type of prediction model utilizes a lookup table for the three measured parameters (W, E, and S _n). It is important to note that a good model for backscattering at small length scales is not sufficient to be able to obtain optical properties from an LEBS peak, despite the fact that Eq. 1 implies this. This is illustrated in Fig. 6, which shows that different combinations of [TeX:] $l_s^*$ $l_s^{*}$ , g, and m can result in LEBS peaks with nearly identical shapes. The reason for this is because p(r) is multiplied by c(r) in Eq. 1. Therefore, short coherence lengths limit the range of p(r) that is accurately obtained thereby making the inverse problem more difficult. There are only two effective independent parameters in a given LEBS peak line shape: the width and enhancement. This problem can be resolved by obtaining a third independent parameter, S _n, from the spectral dependence of the LEBS peak.

An error analysis for each of the described methods is presented in Sec. 6. The errors were calculated by obtaining LEBS peaks from Monte Carlo data, downsampling to the resolution of an existing LEBS instrument, and then obtaining the optical properties with one of the described methods by treating [TeX:] $l_s^*$ $l_s^{*}$ and g as unknowns. The LEBS peaks obtained from Monte Carlo do not contain spectral information and therefore the value of m was assumed to be known. In the experimental case, any error in calculating m from the spectral slope will propagate and also contribute to the error in the measurement of g because the LEBS peak width is sensitive to both m and g. Some error would also propagate to the measurement of [TeX:] $l_s^*$ $l_s^{*}$ , but this would be a minor influence because the LEBS enhancement weakly depends on m and g. Out of the three models that were described, the best performance is obtained with the lookup table for E, S _n, and W _RI [Figs. 9c, 9d]. The lookup table for E, S _n, and W had excellent performance for short coherence lengths, but the accuracy declined for longer spatial coherence lengths [Figs. 9a, 9b]. This is due to the fact that there is a transition from LEBS to enhanced backscattering (EBS) as the coherence length increases. In EBS, the peak width is inversely proportional to [TeX:] $l_s^*$ $l_s^{*}$ and there is no dependence of the width on g or other properties of the phase function. This is why the EBS measurement is well-approximated by the diffusion approximation.³² Interestingly, the results presented in this work suggest that it is still possible to measure properties of the phase function from EBS measurements obtained with coherent illumination by observing the periphery of the peak. Specifically, W _RI at longer L _sc contains the same information as W at shorter values of L _sc, as seen from the results in Fig. 9. Two lookup tables are tested that use two different definitions of the width. The other two presented approaches that utilize models of P(r) for predicting the LEBS peak had errors that did not strongly depend on L _sc. This is because these models result in a prediction of the entire LEBS peak line shape and include the periphery of the peak. The error in these models is due to inaccuracies in fitting the coefficients from Eq. 16 or Eq. 17, as well as limitations due to the simplicity of the models. For example, the PCA-based method is more complicated in that it utilizes an additional component, but results in improved accuracy for measuring g and [TeX:] $l_s^*$ $l_s^{*}$ (Fig. 8).

Fig. 8

Error analysis of optical property measurements that utilize predictive models of P(r). (a) and (b) show the error in [TeX:] $l_s^*$ $l_s^{*}$ and g, respectively, for difference model described by Eq. 16. (c) and (d) show the error in [TeX:] $l_s^*$ $l_s^{*}$ and g, respectively for the PCA method described by Eq. 17. The black dashed line indicates an ideal measurement. The PCA method results in less error than the difference method for tissue-relevant optical properties.

Fig. 9

Error analysis of lookup table methods for obtaining optical properties with LEBS. (a) and (b) show the error in [TeX:] $l_s^*$ $l_s^{*}$ and g, respectively, for the lookup table that utilizes the enhancement factor, spectral slope, and full width at half maximum of the LEBS peak. The error is minimal for [TeX:] $L_{{\rm sc}} \ll l_s^*$ $L_{\mathrm{sc}}\ll l_s^{*}$ and increases with increasing L _sc. When the modified definition of the LEBS width is used [W _RI = I _LEBS(0.5 deg)/I _LEBS(0 deg)], the error in [TeX:] $l_s^*$ $l_s^{*}$ and g [(c) and (d), respectively] improves and becomes independent of L _sc. The method that uses W _RI is able to obtain precise values of gand [TeX:] $l_s^*$ $l_s^{*}$ for the entire tested range.

As an experimental example, we develop a scattering phantom with a phase function that can closely mimic the Whittle–Matérn phase function. The phantom is composed of a mixture of microspheres of three sizes with a relative concentration that minimizes the error between the resulting phase function and the target Whittle–Matérn phase function [Fig. 11a]. The LEBS peak is measured from the constructed phantom and compared with the predicted Monte Carlo simulation [Fig. 11b]. The lookup table for E, S _n, and W _RI is used to measure the optical properties. The phantom optical properties were estimated to be within the range [TeX:] $l_s^*$ $l_s^{*}$ = 668 ± 33 μm and g = 0.899 ± 0.002, according to the tolerances provided by the manufacturer. The measured optical properties were [TeX:] $l_s^*$ $l_s^{*}$ = 687 ± 13 μm (mean ± 95% CI) and g = 0.91 ± 0.04 (mean ± 95% CI). The confidence intervals of the measurements overlap with the ranges estimated from the manufacturer specifications, suggesting an accurate measurement of these properties. A measurement of m was not obtained from the phantom because m and l _c were not utilized as fitting parameters in the minimization, as this would require a fit of the spectral dependence of the Mie scattering cross section and perhaps a larger number of sphere sizes in the mixture.

It is important to point out that the presented work does not treat any polarization effects. Unpolarized illumination is described along with a scalar-wave model. The resulting optical property measurement is an average over all polarization orientations. It has previously been demonstrated that EBS has the potential to measure properties of birefringent samples.⁵³ The polarization-dependent characterization of tissue can therefore potentially be assessed with LEBS given that the model and experiment is modified to incorporate polarization and orientation dependent features.

In summary, LEBS has several important advantages for measuring optical properties over other existing techniques. LEBS is a depth-selective method that allows for measurement of backscattering at short penetration depths and has the potential to be implemented in vivo. LEBS can be used to measure not only [TeX:] $\mu _s^*$ ${\mu }_s^{*}$ but also the shape of the phase function as characterized by g and m. Additionally, in weakly scattering media such as soft tissue where the Born approximation can be applied, the parameter m provides a measurement of the type of the refractive index correlation function. The values of the optical properties can be measured with a choice of the three presented inverse models (Secs. 5, 6). In Sec. 7, an experimental phantom with known properties of the phase function was described and [TeX:] $l_s^*$ $l_s^{*}$ and g measured from LEBS were shown to be in agreement with predicted values.