1.

Introduction

Elastic light scattering is ubiquitous in an almost unlimited number of both animate and inanimate media, and therefore strongly influences our visual perception of most objects. The common feature for all of these media is the presence of internal refractive index variations, which give rise to light diffraction that is macroscopically denoted as light scattering. A very important example of turbid media is human tissue. Thus, there is a potential biomedical benefit to the study of light scattering in tissue. By looking closely at light backscattering from tissue, it is possible to study its microstructure which is encoded in the subdiffusive scattering intensities. Our study aims at simplifying the decoding of these subdiffusive backscattering intensities.

The propagation of light in turbid, i.e., elastically scattering, media can be envisioned by a narrow beam of nonpolarized light vertically incident on the flat surface of a scattering medium with isotropic light propagation. In such a medium, the scattering coefficient ${\mu }_s$ is direction independent unlike for fibrous media. The envisioned beam of light will give rise to a spatially resolved reflectance (SRR) signal with radial symmetry and an intensity decline with increasing source detector separation (SDS) $\rho$. The slope of the decline as well as the overall backscattering intensities are dependent on the macroscopic scattering and absorption properties. See Fig. 1(a) for a large number of such SRR curves.

Fig. 1

Analytical forward calculations of (a) SRR and (b) SFDR for 1000 different scattering phase functions (see Table 1). The model calculations assume (a) a vertically incident Gaussian beam or (b) a spatially modulated projection, and a numerical detection aperture of one (see insets). The color for each of the 1000 curves relates to the $\sigma$ value of the underlying phase function. The parameter $\sigma$ will be introduced in Sec. 3. The vertical dashed lines indicate $\rho$ and $f$ values further considered in Figs. 2, 3, and 6.

On average, the further away photons are detected from the incident beam, the deeper they have penetrated into the turbid medium. In consequence, backscattering (i.e., reflectance) intensities directly adjacent to a narrow beam of incident light correspond on average to the most shallow penetration depth. These proximal reflectance intensities consist of both, deeply penetrating light with a large number of scattering interactions and shallow photons with superficial propagation paths. We denote the former contribution as diffusive background. Shallow photons which undergo a smaller number of scattering interactions will be referred to as subdiffusive light.

It is the ambition of this work to enhance both the understanding and the accurate quantification of this subdiffusive light. In particular, we are interested in how subdiffusive reflectance intensities relate to the nature of single scattering events and thus to the scattering phase function of scattering particles. The latter quantifies the probability distribution for light to be redirected into different scattering angles in each scattering interaction.

To most accurately predict subdiffusive light propagation in turbid media, Maxwell’s equations need to be solved requiring the precise microscopic and complex-valued refractive index distribution inside the media.¹ In most cases, like for biological tissue, this very detailed information is not available nor is it possible to computationally solve Maxwell’s equations for sections in the range of a cubic millimeter in a timely manner. Fortunately, the accuracy related to the solution of Maxwell’s equations is often not required for larger-sized turbid media, especially when near-field effects and dependent scattering of neighboring scattering particles are negligibly small.

In this study, we make use of the radiative transfer model¹ that states an advanced yet empirical and well-approved approximation to Maxwell’s theory. It treats light propagation in scattering media as a succession of independent and localized scattering and absorption events, rather than as the propagation of electromagnetic waves in a medium with continuous and complex refractive index variations. In order to solve the radiative transfer equation (RTE), a scattering phase function specifying single scattering events is required. For some scattering particles like spheres, ellipsoids or cylinders, this scattering phase function may be computed analytically using Mie theory or similar solutions to Maxwell’s theory. In most practical cases, scattering particles have a rather complex and mostly irregular shape with only roughly known refractive index distributions. This is especially true for cellular scattering in tissue, where each cell and cellular organelle varies in size and structure and would give rise to a different scattering phase function. This problem can be dealt with by averaging over an ensemble of single scatterers through numerical simulations or by goniometric phase function experiments on highly diluted scatterers.

Subdiffusive light is ultimately affected by the nature of single scattering events based on the small number of scattering interactions. Therefore, experimental quantification of subdiffusive reflectance offers the great opportunity to derive meaningful microstructural properties of a turbid medium. This is in contrast to the diffusive background and reflectance at larger SDSs which are both only weakly influenced by the scattering phase function.

The benefit of experimentally quantifying subdiffusive light relies on the availability of a disambiguate correspondence between the nature of single scattering and subdiffusive light intensities. Such a correspondence allows meaningful conclusions about a media’s microstructure based on its subdiffusive light intensities.

Bevilacqua and Depeursinge² have reported on the scattering phase function parameter $\gamma =(1-g_2)/(1-g_1)$ as a quantifier for subdiffusive reflectance. This parameter employs $g_1$ and $g_2$ constituting the first and second Legendre moment of the scattering phase function, respectively. Since its introduction, the parameter $\gamma$ was adopted in a great number of experimental and theoretical studies.^3–11 The number of reports does not only prove the usefulness of the parameter, but also the general interest and attributed benefit toward quantification of subdiffusive backscattering.

In two recent studies, we have focused on subdiffusive light in the spatial frequency regime.^10,11 In these studies, we demonstrated the obtainable microstructural contrast through quantification of $\gamma$¹¹ and revealed major ambiguities related to $\gamma$.¹⁰ We have reported that false or imprecise assumptions on a model phase function can lead to considerable uncertainties in the derived $\gamma$-values.¹⁰ This states a major limitation to the significance of derived $\gamma$ values as the actual scattering phase function of a turbid medium needs to be known for proper selection of a model phase function.

In Sec. 2, we will review the current approach for quantification of subdiffusive backscattering. In this context, we will demonstrate the inability of the scattering asymmetry parameter $g_1$ to independently predict subdiffusive backscattering and review the uncertainties related to the established parameter $\gamma$. In Sec. 3, we will build upon $\gamma$ making use of higher-order similarity relations and introduce a new subdiffusive quantifier. In doing so, we overcome some of the uncertainties related to $\gamma$ and demonstrate the accuracy of the new quantifier in Sec. 4. We finally conclude on the benefits and consequences of our work for the experimental study of subdiffusive backscattering in Sec. 5.

2.

Asymmetry $g_1$ and Phase Function Parameter $\gamma$

We start our analysis on subdiffusive backscattering by presenting a large number of analytical forward calculations. Figure 1(a) exemplifies SRR computed for a large number of scattering phase functions assuming a Gaussian beam of light with FWHM of $300\mu \mathrm{m}$ vertically incident on an ideally flat, semi-infinite and homogeneously scattering medium. The curves illustrate the backscattering intensity integrated over the entire detection half-sphere versus SDS $\rho$ from the point of incidence [see inset of Fig. 1(a)]. For all forward computations, the assumed reduced scattering intensity is ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and the absorption value was set to ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ with mean refractive indices $n=1.38$ and $n=1.00$ inside and outside of the turbid medium, respectively.

The large number of curves in Fig. 1(a) makes individual curves indistinguishable. Different coloring of the curves however clearly reveals the large phase function influence both at very small, but also at intermediate SDSs. We find backscattering intensities at the point of incidence ($\rho =0\mathrm{mm}$) to vary by more than a factor of ten, and even at $\rho =10\mathrm{mm}$ intensities are found to differ by up to 10% (not shown). At a SDS of $\rho \approx 0.8\mathrm{mm}$, we observe a local minimum of the phase function influence on SRR which has previously been reported by Calabro and Bigio.¹²

For Fig. 1(a), 1000 different phase functions were successively assumed and the details of the selected phase functions are stated in Table 1. The emphasis of the phase function selection lies at modeling of biological tissue while we deliberately added more diversity through Mie functions with high-refractive index mismatch and other more exotic phase function models. The same phase function selection presented in Table 1 will be used throughout this study. Please note that some of the outer refractive index values in the Mie function models are meant to roughly match that of water and do not equal the assumed mean refractive index value of the medium of $n=1.38$.

Table 1

Employed phase functions and corresponding parameter ranges for RTE-based analytical forward calculations presented throughout our study.

Similar to Fig. 1(a), Fig. 1(b) depicts the phase function influence on reflectance intensities in the spatial frequency domain for the same phase function range. The point source is in this case replaced by a spatially extended and laterally structured sinusoidal projection pattern [see inset of Fig. 1(b)].¹⁵ The plotted spatial frequency domain reflectance (SFDR) in Fig. 1(b) is a direct measure of the unitless amplitude retention in the backscattering signal. The advantage of SFDR over SRR is the ability to limit the light propagation depth through selection of a spatial frequency. High-spatial frequencies are more surface sensitive and thus the relative influence of scattering phase functions on SFDR can be observed to increase with spatial frequency [see Fig. 1(b)].

Our forward calculations are based on analytical solutions to the RTE derived by Liemert and Kienle.^16,17 These models assume semi-infinite and homogeneously scattering media with ideally flat surfaces neglecting roughness induced surface scattering or specular reflections. All analytical solutions are verified by Monte Carlo (MC) simulations and within noise limits provide identical yet much faster and noise-free computation of reflectance patterns. In consideration of computation time, the computation order $N$ for calculation of SRR and SFDR was set to $N=17$ and $N=25$, respectively.

The ambition in the search for subdiffusive quantifiers is a maximum reduction of the complexity of a scattering particle’s phase function $p(\theta )$ without loosing the information that is relevant to subdiffusive light. We denote $p$ as a function of the polar angle $\theta$ omitting the dependence on the azimuthal angle $\phi$. This states a frequent approximation that is especially valid for orientation randomized scatterers. For nonpolarized light scattering on homogeneous spheres, $p$ is independent of $\phi$.

Given the complexity of a real scattering phase function, it is impossible to exactly describe subdiffusive backscattering with only a single generalized parameter. However, a study by Turzhitsky et al.¹⁸ using principal component analysis could demonstrate accurate modeling of subdiffusive backscattering intensities for biologically relevant phase functions by very few principal components.

Several analytical studies have derived equations for computation of effective scattering functions corresponding to the angular light distribution after two or more scattering interactions.^19–22 These reports demonstrate not only the complexity in deriving such effective scattering functions, but give illustrative conclusions on how oscillatory features or strong forward peaks in $p(\theta )$ get canceled through multiple scattering. While these analytical models help to calculate and predict subdiffusive light intensities,²³ they do not hint at simple quantifiers for subdiffusive backscattering.

The notion of phase function dependent light propagation through independent scattering interactions is based on the radiative transfer model. This model is often applied through the use of MC simulations. Analytical solutions to the RTE express the scattering phase function in terms of Legendre moments $g_n$, with

Even though $g_1$ has strong influence on the effective light scattering intensities in the diffusion regime through modification of the parameter ${\mu }_s^{\prime }$, we find that it has almost no relevance as an independent subdiffusive parameter. Mourant et al.²⁴ have already demonstrated through MC analysis that large differences in subdiffusive backscattering exist for Mie and Henyey-Greenstein phase functions with equal asymmetries $g_1$. We seek to provide a more general conclusion on the relevance of $g_1$ as an independent parameter for subdiffusive backscattering. To this end we reuse the reflectance data presented in Fig. 1.

Figures 2(a) and 2(b) depict both SRR and SFDR for all 1000 phase functions against their corresponding asymmetry values $g_1$. SRR in Fig. 2(a) is given for the four SDSs $\rho =0,0.5,1.5$ and 3 mm and Fig. 2(b) provides reflectance intensities for four spatial frequency values with $f=0.1,0.2,0.4$ and $0.8{\mathrm{mm}}^{-1}$ (see dashed lines in Fig. 1). Every data point in Fig. 2 thus corresponds to reflectance at a certain radius (a) or frequency (b) for forward computations using one of the 1000 different scattering phase functions of Table 1.

Fig. 2

Analysis of subdiffusive reflectance intensities (a) for four different SDSs and (b) for four spatial frequencies versus the phase function asymmetry value $g_1$. The data equal that of the forward calculations presented in Fig. 1. Every data point corresponds to a different phase function listed in Table 1. The observable patterns in each point cloud relate to the different phase function models and their parameterization. Small correlation is found in either graph with respect to $g_1$. As a guide to the eye for noncolored print, data symbols in (b) for frequency $f=0.4{\mathrm{mm}}^{-1}$ are slightly enlarged.

Both SRR and SFDR data in Fig. 2 clearly demonstrate the inability of $g_1$ to quantify subdiffusive reflectance at any of the investigated SDSs (a) or spatial frequencies (b). This observation is based on the very weak correlation of each point cloud comprising 1000 data points with respect to the parameter $g_1$. Indeed, different phase functions, which deviate in their $g_1$ values by more than 0.8, are shown to produce similar reflectance intensities.

In spite of the demonstrated inability of $g_1$ to describe SRR and SFDR, it is worth mentioning that experimental quantification of $g_1$ has been shown to be possible in a backscattering geometry. Samatham and Jacques^25,26 have demonstrated sensitivity toward the asymmetry parameter using a reflectance mode confocal microscopy setup. The $g_1$-sensitivity of their system is achieved by the added depth control for detection owing to the confocality of their system.

The phase function parameter $\gamma$ was derived by use of similarity relations,² which is a concept that has been introduced by van de Hulst²⁷ and that was further developed by Wyman et al.²⁸ thereafter. The similarity relations are $N+1$ constraints for two sets of optical parameters (i.e., two sets of ${\mu }_a$, ${\mu }_s^{\prime }$ and $p(\theta )$) to produce the identical radiance in an unbounded medium within the context of a computational order $N$. In consequence, the exact solution (i.e., $N\to \infty$) requires an infinite number of similarity relations to be fulfilled, essentially forcing the first and second set of parameters to be identical. This provides the important insight that every phase function gives rise to its unique light propagation.

In first order approximation ($N=1$), the two similarity relations are given by²

Figure 3 plots our previously calculated SRR and SFDR data versus $\gamma$ and demonstrates its general significance in describing subdiffusive backscattering. For the three larger $\rho$ values in Fig. 3(a) strong correlation of SRR with $\gamma$ can be found. At $\rho =0$, this correlation is strongly reduced. In the spatial frequency domain [Fig. 3(b)], we find the $\gamma$-reflectance point clouds to be much less correlated for all four spatial frequencies. However, a comparison of Figs. 2 and 3 clearly reveals the superiority of $\gamma$ over $g_1$ in describing subdiffusive backscattering. Nevertheless, the strong dispersion of data points, especially in the spatial frequency domain [Fig. 3(b)], states a considerable drawback for the experimental quantification of subdiffusive reflectance.

Fig. 3

SRR for four different (a) SDSs and (b) SFDR for four spatial frequencies versus the phase function parameter $\gamma$. Both graphs rely on the data already presented in Figs. 1 and 2. Every data point corresponds to one of the phase functions listed in Table 1. $\gamma$ shows decent correlation with (a) SRR for three out of the four presented $\rho$ values while large uncertainties exists in the quantification of (b) SFDR for all four frequencies.

3.

Introducing the Subdiffusive Quantifier $\sigma$

In trying to more precisely describe subdiffusive backscattering, it appears as a natural step to make use of higher-order similarity relations. However, not only does a higher number of similarity parameters beyond $\gamma$ [i.e., $\delta =(1-g_3)/(1-g_1),\epsilon =(1-g_4)/(1-g_1),\dots$] increase the overall complexity,³⁰ but also do the added parameters carry little independent structural information. It is therefore our goal not to increase the number of subdiffusive parameters, but to combine all similarity parameters in a single yet less ambiguous quantifier to describe subdiffusive light propagation. We arrive at this goal by an empirical approach which builds upon the findings of Bevilacqua and Depeursinge and their established parameter $\gamma$.

As a first step, it is instructive to study typical photon paths in the subdiffusive scattering regime. Figures 4(a)–4(c) each contain a two-dimensional projection of 20 photon paths for different light propagation constraints. In all three graphs, the given photon paths correspond to MC simulations in a semi-infinite scattering medium ($n=1.38$, ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$). The simulation imitates tissue scattering by assuming a scattering phase function that corresponds to a fractal distribution of Mie scatterers with fractal dimension $\alpha =4.5$ (asymmetry $g_1=0.90$). Light is obliquely incident from air in the $x$–$z$-plane at an angle of 15 deg with photons entering the turbid media at $(x,y,z)=(0\mathrm{,}0,0)$, as indicated by the arrows in Figs. 4(a)–4(c). To mimic typical experimental detection geometries, we show only photon trajectories that exit the medium within an aperture angle of 10 deg at a SDS smaller than a set radius $r$.

Fig. 4

MC analysis of both subdiffusive and diffusive light propagation. The three columns in the upper part of the figure relate to MC simulations with different constraints for the detection radius $r$ and the allowed penetration depth $z_{\min }$. (a)–(c) 20 photon paths for each simulation exemplify the nature of each scattering regime. (d)–(f) Statistical evaluation of all scattering interactions of detected photons reveals (d)–(e) a strong bias toward high-angle scattering as compared to the expected angular scattering distribution. This bias is almost canceled in the diffusion regime [(c) and (f)]. (g) All simulated light propagation and detection constraints give rise to a more or less pronounced sensitivity toward subdiffusive and diffusive light.

To specifically select photons corresponding to different scattering regimes, the detection radius $r$ as well as the allowed photon penetration depth $z_{\min }$ differs among the three sets of photon paths. In Fig. 4(a), $r=0.1\mathrm{mm}$ with depth constraint $z_{\min }=-0.5\mathrm{mm}$. In Fig. 4(b), photons are detected within $r=0.3\mathrm{mm}$ and are allowed to travel to a depth of $z_{\min }=-1.5\mathrm{mm}$. Both dimensional constraints are lifted for the simulation corresponding to Fig. 4(c). Please note the different axes scaling of Figs. 4(a)–4(c).

The three different MC simulations where designed to correspond to light propagation with a minimum (a), an intermediate (b), and a maximum (c) share of diffusive light content. By diffusive light we refer to the unbiased and randomized photon paths whose scattering angles statistically reflect the probability distribution of the scattering phase function. Correspondingly, the distinct feature of subdiffusive light lies in the directional bias of its photon paths. This bias can be well studied in Figs. 4(a) and 4(b) which contain photon trajectories with mostly forward and large angle backward scattering. As the diffusive influence increases [Figs. 4(b) and 4(c)], there is an increasing number of paths featuring also intermediate scattering angles.

In Figs. 4(d)–4(f), we do a statistical analysis to confirm the visually perceived angular scattering bias. In this analysis we evaluate all scattering angles of detected photons for all three MC simulations. Figures 4(d) and 4(e) relate to the first two simulations and 4000 scattering interactions where evaluated each, corresponding to 970 and 395 detected photons, respectively. For Fig. 4(f), only 200 photon paths where evaluated corresponding to 46,250 scattering interactions in the unconstrained third simulation. In all three graphs, we display the ratio given by the number of encountered scattering angles $p_{\mathrm{MC}}$ for one of 20 angular bins divided by the expected occurrence of scattering events $p_{\mathrm{frac}}$ for each bin. $p_{\mathrm{frac}}$ corresponds to the angular scattering probability of the assumed fractal phase function scaled by the number of detected photons.

The statistics in Fig. 4(d) clearly reveal over-proportional scattering at high-scattering angles for the most constrained MC simulation. Indeed, high-angle backscattering is found to occur one hundred times more frequent than statistically expected from the scattering phase function. At the same time, intermediate scattering angles are under-represented due to the unlikeliness for photons to reach the detector within the boundary constraints after more randomized scattering interactions.

This directional bias is considerably decreased for the second simulation [Fig. 4(e)] due to a slightly stronger influence of more directionally randomized photon paths.

In the last simulation with neither lateral detection nor axial propagation constraint, the angular scattering bias is almost completely lost [Fig. 4(f)]. The incidence and detection geometry nevertheless sustain a small bias toward very high-angle scattering by allowing subdiffusive light to still be part of the detection signal.

Figure 4(g) sketches the relationship of the three MC simulations, indicating the difficulty to clearly distinguish between the diffusive and subdiffusive regime due to the smooth transition of the relative weight of biased and unbiased photon paths.

Further MC analysis beyond that shown in Fig. 4 reveals that the studied angular scattering bias has some dependency on the source–detection geometry and is very sensitive to the implemented depth constraints. These constraints might seem as a random limitation but actually reflect the sampling characteristics in spatial frequency domain imaging. Deeply penetrating light is prevented from contributing to high-spatial frequency reflectance through a demodulation approach.¹⁵

Similar to Figs. 4(d)–4(f), we also observe a qualitatively comparable angular scattering bias for different scattering phase functions. However, we find the intensity of the angular scattering bias to scale with the phase function asymmetry $g_1$.

As a conclusion from this MC analysis, it seems reasonable that accurate quantification of subdiffusive backscattering should somehow account for the demonstrated high weight of scattering at large angles. We seek to use higher order Legendre moments for quantification of subdiffusive light. In this regard one needs to be aware of the change in sign of the Legendre polynomials $P_i$ at angles close to $\theta =180\mathrm{deg}$ for subsequent orders (see inset of Fig. 5 for illustrations of $P_2$, $P_3$, and $P_4$ with positive function values represented by solid blue and negative values by dashed red curves). Correspondingly, we find through a large number of forward calculations using different scattering phase functions and different computational orders that the effect of subsequent similarity parameters (i.e., $\gamma ,\delta ,\epsilon ,\dots$) on backscattering is also inverted. Furthermore, increasingly high-Legendre moments correspond to the high-frequency oscillations of scattering phase functions and thus have a diminishing impact on the overall subdiffusive backscattering intensities.

Fig. 5

Angular weight curves for the expression defined in Eq. 5 for different values $c$. The solid blue curve ($c=0$) corresponds to the angular weight assigned to $\gamma$ and the solid black curve ($c=0.5$) to that of the newly defined subdiffusive quantifier $\sigma$. The inset depicts the Legendre polynomials $P_2$, $P_3$, and $P_4$ with their change in sign for angles close to 180 deg (solid blue curves for positive, and dashed red curves for negative function values).

Based on these observations we propose the following formula to quantify subdiffusive backscattering

In order to understand the physical relevance of this empirical equation and to demonstrate its suitability, Eq. (5) can be rewritten as

We are particularly interested in the quantity marked as “angular weight” in Eq. (7). This quantity provides some insight on the relative weight of the entire expression toward the different phase function angles. We seek to compare the angular weight of the presented formula with our previous statistical analysis on angular scattering bias. Figure 5 displays the formula’s weight for different constants $c$. The blue colored dashed curve for $c=0$ corresponds to the angular weight of $\gamma$ displaying equal emphasis on the forward and backward scattering direction. This does not quite correspond to the angular weights obtained through our statistical analysis in Fig. 4. However, as $c$ is increased, the backscattering direction progressively gains in angular weight.

For $c>0$, we compare the angular weight curves in Fig. 5 with the statistically found angular weight histograms in Figs. 4(d) and 4(e). Accordingly, the observed high-angular weight of large angle scattering is reflected by Eq. (7) which also attributes less weight to intermediate scattering angles as desired. For absolute comparison, the scaling factor $1/(g_1-1)$ and the offset $1/(1+c)$ of Eq. (7) have to be considered as well. While $1/(g_1-1)$ scales the overall weight intensities, $1/(1+c)$ reduces the weight at all angles and may thus reduce the displayed relative weight of forward scattering depending on the actual phase function. Consequently, by using Eq. (5) with $c>0$ one may establish a rather good correspondence with the actual angular weights found for subdiffusive light. Moreover, the formula’s built-in scaling property of the asymmetry value reflects the previous finding of our MC analysis.

We have tested the significance of Eq. (5) on a large number of analytical^16,17 forward calculations for different scattering phase functions and various source–detector geometries. Resulting from these calculations, we find on average best quantification of subdiffusive backscattering for $c=0.5$ and thus introduce the parameter $\sigma$ with

4.

Verification and Accuracy of $\sigma$

As a first step to demonstrate the significance of $\sigma$, we make use of the calculated SRR and SFDR data of Sec. 2 which was previously presented with respect to $g_1$ and $\gamma$ in Figs. 2 and 3, respectively. The same data is now displayed versus $\sigma$ in Fig. 6 showing both SRR (a) and SFDR (b). By comparing the dispersion of the 1000 data points for all four SDSs in Fig. 6(a) to those in Fig. 3(a), we observe that $\sigma$ describes SRR intensities more uniquely for three out of the four displayed SDSs. At $\rho =0.5\mathrm{mm}$, SRR is well described by both $\gamma$ and $\sigma$, yet data points are slightly more widespread in the graph versus $\sigma$. A noticeable gain in accuracy using $\sigma$ is especially found for $\rho =0\mathrm{mm}$, but also for $\rho =1.5\mathrm{mm}$ and $\rho =3\mathrm{mm}$.

Fig. 6

Analytical forward calculations of (a) SRR for four different SDSs $\rho$ and (b) SFDR for four spatial frequencies $f$. The shown reflectance data is equal to that presented in Figs. 1–3 but is plotted against the new subdiffusive parameter $\sigma$ this time. Every data point corresponds to a different scattering phase function. $\sigma$ is observed to have a stronger correlation with subdiffusive SRR and SFDR than $\gamma$ when compared to Fig. 3.

Improvements in quantifying subdiffusive backscattering are most evident in the spatial frequency domain [compare Figs. 3(b) and 6(b)]. For all four displayed spatial frequencies $\sigma$ states a much more unambiguous quantifier than $\gamma$.

In order to add some quantitation to the predictive ability of $g_1$, $\gamma$, and $\sigma$, a fifth degree polynomial was fitted to all SRR and SFDR point clouds in Figs. 2, 3, and 6 (not shown). We then calculated the root mean square deviation (RMSD) for each subdiffusive quantifier with respect to this polynomial. A relative overall error metric is obtained for the spatial domain and spatial frequency domain by, respectively, averaging the RMSDs over the four SDSs or four frequencies and division by the parameter ranges in $g_1$, $\gamma$ or $\sigma$. Table 2 summarizes these averaged relative errors $\delta g_1$, $\delta \gamma$, and $\delta \sigma$. Accordingly, relative RMSDs in $\gamma$ are found to be 1.5 and 3.5 times as large as that of $\sigma$ for the spatial domain and spatial frequency domain, respectively. This error metric should, however, be considered solely as a rough indicator which is largely dependent on the phase function selection and the assumed source and detection geometry.

Table 2

Relative error metrics for g1, γ, and σ as averaged over the data of Figs. 2, 3, and 6.

To also demonstrate the significance of $\sigma$ beyond both the four SDSs and the four frequency values analyzed in Fig. 6, all SRR and SFDR curves in Fig. 1 of the original data were colored according to their corresponding $\sigma$ value (with darker curves corresponding to lower $\sigma$ values).

The perpendicular illumination and high-aperture detection assumed for the presented data in Fig. 6 have some similarity with the constraints found in reflected-light microscopy yet they do not represent typical diffuse optical imaging constraints. Oblique projection of light is a typical feature of experimental approaches toward SRR or SFDR. Detection is often camera based featuring a limited and mostly small detection aperture. In order to verify and thereby establish $\sigma$ as a meaningful parameter relevant also to experimental geometries, we performed extensive forward calculations of subdiffusive reflectance for four different geometries of the source, the detection and the sample.

The following computations were performed for SFDR, owing to the much higher-computation speed of analytical solutions in the spatial frequency domain. This allows to perform computations at a higher-computation order $N$ for better verification of $\sigma$. High-computation orders are especially required for narrow detection apertures and high-spatial frequencies where oscillatory phase function features start to have a meaningful influence on reflectance intensities. In our forward calculations we again set the macroscopic optical properties to ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.01{\mathrm{mm}}^{-1}$. The selected absorption value has little significance due to the in most cases negligible influence of absorption on subdiffusive light propagation.¹⁰ Through scaling of spatial frequencies, the computational results hold true also for higher or lower reduced scattering values and thus give a generally meaningful analysis for the subdiffusive regime. This is similarly true for the previously presented SRR data, where proper scaling of both $\rho$ and SRR corresponds to scaled values of ${\mu }_s^{\prime }$ and ${\mu }_a$.

In the following computations we consider the four spatial frequencies $f=0.1,0.2,0.4$ and $0.8{\mathrm{mm}}^{-1}$ for all phase functions listed in Table 1, as before.

Figures 7(a), 7(d), 7(g), and 7(j) sketch the modeled geometries of source, detection and the sample. In the first forward modeling (first row of Fig. 7) we assume oblique incidence at an angle of 30 deg and a detection aperture of 0.34 [see Fig. 7(a)]. The second geometry has oblique incidence at 20 deg and a smaller aperture of 0.17 [Fig. 7(d)]. In the third row of Fig. 7 we assume a two-layer geometry where the scattering layer is only $100\mu \mathrm{m}$ thick resting on a highly absorbing bottom layer with ${\mu }_a=10{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.01{\mathrm{mm}}^{-1}$. This geometry restricts light propagation to within the depth of $100\mu \mathrm{m}$ leading to an almost pure ballistic (i.e., strongly subdiffusive) reflectance content. For this geometry, projection is at 20 deg and the detection aperture is set to 1 [Fig. 7(g)]. Finally, Fig. 7(j) shows the last forward modeling which assumes the same source–detection and sample geometry as in Fig. 7(g) with the bottom layer optical properties set to ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$.

Fig. 7

Results of RTE-based analytical forward calculations of SFDR to demonstrate the accuracy of the phase function parameter $\sigma$. The sketches in the left column depict the assumed geometries for source, detection and sample for the four different sets of forward calculations. For each model geometry, two plots are given in the central and right column showing the computed reflectance values with respect to $\gamma$ and $\sigma$, respectively. Every data point corresponds to the forward calculation using one of the phase functions stated in Table 1 for one out of four spatial frequencies ($f=0.1,0.2,0.4$ and $0.8{\mathrm{mm}}^{-1}$).

The second and third column of Fig. 7 display the computed reflectance values for all model geometries plotted versus $\gamma$ and versus $\sigma$, respectively. Every data point again corresponds to a single forward calculation using one of the scattering phase functions listed in Table 1. The four point clouds in each graph thus correspond to the four different spatial frequencies.

We always find the computed reflectance values in Fig. 7 to be more scattered when plotted versus $\gamma$ (second column) as compared to the corresponding $\sigma$ plots (third column). This is confirmed by the relative error $\delta \gamma$ and $\delta \sigma$ presented in and averaged over each graph. Thus, for all four model geometries $\sigma$ states a more unambiguous parameter and shows stronger correlation with subdiffusive reflectance intensities.

With decreasing detection aperture and increasing spatial frequency the uncertainties related to both parameters are shown to increase. This is caused by the increased sensitivity toward small angular phase function features such as Mie resonances. Indeed, the most deviating points in the data graphs of Fig. 7 correspond to scattering functions with a rather distinct or oscillatory peak in backward direction. Such phase function behavior is unlikely to be found for biological tissue but observable for single Mie scatterers or “modified powers of cosines”² phase functions.

The third model geometry corresponds to the reflectance from a very thin slice of turbid medium and reveals an almost linear relationship of reflectance intensities versus $\sigma$. The small layer thickness states an even stronger constraint to the light propagation depth than the selected spatial frequencies and thus reflectance intensities are almost identical for all frequency values. We therefore show SFDR for only a single spatial frequency ($f=0.1{\mathrm{mm}}^{-1}$) in Figs. 7(h) and 7(i) representative for all four frequency values.

The rather unnatural bottom layer optical properties in the layered forward models [Fig. 7(g) and 7(j)] are meant to demonstrate general applicability of $\sigma$ for layered sample geometries with unmatched optical properties.

All source geometries in Fig. 7 feature oblique projection of spatial frequencies which induces phase shifts in the reflectance pattern. In spite of not illustrating these phase shifts we included them in our computational analysis. Similar to the depicted reflectance amplitudes, we observe that $\sigma$ relates to a much more unambiguous description of the phase than $\gamma$ for all studied geometries and frequencies.

The correctness of the presented analytical data is ensured by numerous verifications using MC simulations with phase function look-up tables. In addition, we ensure numerical stability of all data points by repeated forward computations at different computational orders.

5.

Conclusion

We demonstrated that the use of the asymmetry parameter $g_1$ as an independent subdiffusive quantifier is inappropriate. The established phase function parameter $\gamma$ is much more useful but also carries large uncertainties. Through analysis of subdiffusive photon paths we arrived at the new parameter $\sigma$ which builds upon $\gamma$ and makes use of higher order similarity relations. This empirically derived parameter greatly enhances the accuracy in quantifying subdiffusively backscattered light.

Our findings have important consequences for the experimental study of subdiffusive light. The parameter $\sigma$ mostly lifts the necessity of finding proper model phase functions to best mimic investigated media. This is especially useful for the study of biological tissue with mostly unknown scattering phase functions.

From our study, we conclude, that it is not the phase function model and its precise shape of angular scattering that matters, but rather the correct $\sigma$ regime. Our verification studies indicate a high significance of $\sigma$ independent of the source and detection geometry and for both the spatial domain and the spatial frequency domain. This hints at general applicability to various experimental approaches such as differential path spectroscopy,^31,32 single or multiple fiber reflectance spectroscopy,^7,33 SRR,³⁴ SFDR,^11,15 and to other modalities without polarization constraints.^35,36

As expected, we find the accuracy of $\sigma$ to decrease with decreasing SDS (i.e., increasing spatial frequency) as small angular phase function features become more relevant. The same is true for decreasing detection apertures.

A second benefit of $\sigma$ comes with the enhanced replaceability of phase functions in forward models. Within the uncertainties displayed in Figs. 6 and 7 there exists a general ability to replace phase functions for forward modeling of subdiffusive (and diffusive) reflectance. Heavily oscillating or extremely strong forward peaked scattering phase functions with slowly decreasing Legendre moments can produce numerical instabilities in analytical forward models. In MC simulations, both low phase function asymmetry values, as well as phase function invertibility boost computation speed. Through $\sigma$-equivalent phase function replacement, it is possible to lower the asymmetry value of a model function without altering subdiffusive backscattering intensities much. The possible extend of such asymmetry lowering can be understood from Fig. 2.

Our study revealed that subdiffusive backscattering intensities are mostly determined by the high-angle scattering probability of individual scatterers scaled by the asymmetry value of their scattering phase function. The parameter $\sigma$ reflects this angular weight distribution. As a derived version of the established parameter $\gamma$, $\sigma$, is likely to relate to similar physical and microsctructural interpretation.⁶

Similar to the measurement of ${\mu }_a$ and ${\mu }_s^{\prime }$, experimental quantification of $\sigma$ requires forward modeling to allow for fitting of optical properties and thus solving of the inverse problem. In biomedical optics, the Reynolds–McCormick function may be a model of choice for performing $\sigma$ forward calculations. Along with its simplicity it performs well in modeling $\sigma$ values below $\sigma \approx 1.26$. Very small $\sigma$ values even below $\sigma =0.2$ are possible if the second input parameter $g^{\prime }$ of the model function (see Table 1) is allowed to become negative. In the same sense, Henyey–Greenstein may be used for $\sigma$ values up to $\sigma \approx 1.11$. The “modified powers of cosines” function allows for modeling of the entire $\sigma$ range presented, but achieves this in part through rather unnatural probability distributions. Consequently, the $\sigma$ and $\gamma$ uncertainties for this phase function type are among the largest along with those for discrete Mie scatterers.

We conclude that every experimental geometry as well as every scattering regime corresponds to a different relative influence of high-order Legendre moments and thus to a different ideal value for $c$. The subdiffusive quantifier $\sigma$ along with the empirically selected constant $c=0.5$ provides an accurate compromise for many imaging constraints. This allows for meaningful comparison across different studies on subdiffusive backscattering.