22 May 2017 Modeling subdiffusive light scattering by incorporating the tissue phase function and detector numerical aperture
Author Affiliations +
To detect small-scale changes in tissue with optical techniques, small sampling volumes and, therefore, short source–detector separations are required. In this case, reflectance measurements are not adequately described by the diffusion approximation. Previous studies related subdiffusive reflectance to γ or σ , which parameterize the phase function. Recently, it was demonstrated that σ predicts subdiffusive reflectance better than γ , and that σ becomes less predictive for lower numerical apertures (NAs). We derive and evaluate the parameter R p NA , which incorporates the NA of the detector and the integral of the phase function over the NA in the backward and forward directions. Monte Carlo simulations are performed for overlapping source/detector geometries for a range of phase functions, reduced scattering coefficients, NAs, and source/detector diameters. R p NA improves prediction of the measured reflectance compared to γ and σ . It is, therefore, expected that R p NA will improve derivation of optical properties from subdiffusive measurements.



Light that has traveled through tissue carries information about tissue properties, such as structure and biochemical composition. The distance that light travels determines the scale on which the tissue is investigated, as optical properties are averaged over the sampling volume. To detect localized, small-scale changes, techniques with small sampling volumes are, therefore, appropriate. Consequently, short source–detector separations are required to ensure most detected photons will have scattered only once or a few times. In this letter, we will investigate the reflectance measured in this so-called subdiffusive regime, specifically for overlapping source–detector geometries, relevant to, e.g., single fiber spectroscopy measurements.

In the diffusive regime, photon direction is randomized and reflectance is independent of the exact shape of the phase function [p(θ), the probability distribution of scattering angles]. In contrast, in the subdiffusive regime, the photon direction is not fully randomized; therefore, measurements are sensitive to the phase function.1,2 Advances have been made to derive analytic models for subdiffusive reflectance.34.5 However, for overlapping source–detector geometries, challenges remain regarding incorporation of the tissue phase function.3,4

To model light transport, solutions to the radiative transport equation (RTE) involve expanding the radiance into a series of i spherical harmonics and the phase function into i Legendre polynomials. The latter are weighted with their moments gi, where g1 is commonly referred to as the scattering anisotropy. For i=1, the diffusion approximation to the RTE is obtained, with the similarity relation μs(1g1)=μs*(1g1*), which expresses that tissues with a different scattering coefficient μs and a different scattering anisotropy g1, but equal μs=μs(1g1) yield the same reflectance. However, reflectance at short source–detector separations is inadequately described for i=1. Improvement is possible by increasing i, giving additional similarity relations. For i=2, Bevilacqua and Depeursinge6 suggested an alternative form of these similarity relations using the parameter γ


Various studies have investigated describing subdiffusive reflectance using γ in combination with μs, μa, and the detector diameter ddet.78.9.10 However, recent work showed that a range of γ values (for the same μs, μa, and ddet) can result in the same reflectance.11,12 Therefore, Bodenschatz et al.13 incorporated more similarity relations into the parameter σ, employing all phase function moments


σ=i=2(c)i21gi1g1  ,
where c depends on both the measurement geometry and the scattering regime (empirically, c=0.5 was found to give the best results). The authors concluded that σ predicts reflectance better than γ; however, σ becomes less predictive for lower numerical apertures (NAs).

We consider subdiffusive reflectance as the sum of a diffusive component and a semiballistic component: the diffusive component depends on μsddet, and, for the semiballistic component, only photons that have experienced a single backscattering event in combination with an arbitrary number of small-angle forward scattering events are considered. To model this contribution, not all details of the phase function are needed. Rather, we consider the phase function within the NA of the detector. The NA characterizes the range of angles over which a detector can accept incoming photons. We propose a theoretically derived parameter RpNA, which is related to the NA of the detector and the integral of the phase function over the corresponding angular interval in the backward (pNA,backward) and forward (pNA,forward) directions, as


We show the theoretical derivation of RpNA and the improved prediction of subdiffusive reflectance based on Monte Carlo simulations by RpNA compared to σ and γ. RpNA, therefore, has the potential to improve models used in subdiffusive measurements to extract optical properties.


Theoretical Derivation of RpNA

We consider semiballistic photons that are backscattered once in combination with an arbitrary number of forward directed scattering events and assume that these photons will only be detected if subsequent scattering events occur at angles smaller or equal to the acceptance angle θNA (Fig. 1), where θNA=arcsin(NA/nsample) and nsample is the refractive index of the sample. We introduce the probabilities pNA,backward and pNA,forward as




The contribution of semiballistic photons to subdiffusive reflectance scales with the probability of a single backscattering event in the pathlength interval dl; this probability equals μs·pNA,backward·dl. To determine the effect of forward semiballistic scattering, we use the result obtained by Wang and Wilson,14 which gives the probability of a photon traveling a pathlength l while experiencing N forward scattering events



Fig. 1

Semiballistic scattering: a combination of any number of forward scattering events N in combination with a single backscatter event. We assume photons will only be detected if all scattering angles θ are smaller than or equal to the acceptance angle θNA of the detector. The scattering angle θ is determined with respect to the previous scattering direction as indicated by the dashed lines.


Summing Eq. (6) over all possible numbers of forward scattering events, N from 0 to , we obtain


probNA,forward(l)=epNA,forward·μsl·eμsl=  eμsl(1pNA,forward)

Thus, the combined probability of one backscatter event and an arbitrary number of forward scattering events, each occurring within θNA, is given by integrating over all possible pathlengths, 0 to





For our Monte Carlo simulations, we modified the software of Prahl et al.15 (which was the core programming later used in Monte Carlo model of steady-state light transport in multilayered tissues (MCML)16 software), to allow the use of arbitrary phase functions using the method of Zijp and ten Bosch.17 The modified code was benchmarked against standard MCML using the Henyey–Greenstein phase function. We modeled geometries where the source and detector geometry were equal and overlapping, where photons were launched from a location based on a uniform distribution across the source with an angle from a uniform angular distribution within the NA. We performed simulations using 15 modified Henyey–Greenstein (mHG),6 177 double Henyey–Greenstein (dbHG),18 8 modified power of cosines (MPC),6 and 46 Reynolds McCormick19 phase functions (RMC, which is equivalent to the Gegenbauer kernel phase function) employing the parameters specified in Table 1 and applying the restrictions g10.5 and g2<0.9 to exclude biologically unreasonable phase functions.

Table 1

Parameters employed in the selection of phase functions.

Phase functionParameters
mHG0.01gHG0.95, 10 linear steps
0.01α0.99, 10 linear steps
dbHG0.5α0.9, 3 linear steps
0.91α0.99, 5 linear steps
0.05gf0.95, 10 linear steps
0.50gb0.05, 5 linear steps
MPC0.01α0.99, 10 linear steps
0.01N10, 10 logarithmic steps
RMC0.01α2.5, 10 linear steps
0.01gR0.95 to 0.2·α, 10 linear steps

For each set of phase functions, simulations were performed for three values of μsddet:0.1 (μs=10  cm1, ddet=100  μm), 1 (μs=100  cm1, ddet=100  μm), and 9 (μs=100  cm1, ddet=900  μm) in combination with an NA of 0.22 or 0.5. An absorption coefficient of 0.1  cm1 and refractive indices inside and outside the sample of 1.35 and 1, respectively, were used.



Figure 2 shows the simulated reflectance versus RpNA, σ, and γ for an NA of 0.22 and 0.5, and three different values of μsddet (0.1, 1, and 9). Reflectance correlates with all three parameters for μsddet of 0.1 and 1. This correlation is far less pronounced for μsddet=9. To compare RpNA, σ, and γ, we determined the spread in each of these three parameters for a chosen reflectance (±10%) relative to the total range of the parameter (Table 2). For example, for μsddet=0.1 and a reflectance of 0.001 (±10%), γ ranges from 1.52 to 2.11. The total range of γ is 0.58 to 2.66, so the variability is calculated as 0.59/2.08=0.28. For all combinations of NA (0.22 and 0.5) and μsddet values (0.1 and 1), we determined these variability values for three reflectance values and found the variability to be lowest for RpNA compared to σ and γ.

Fig. 2

Simulated reflectance versus RpNA, σ, and γ for (a)–(c) NA=0.22 and (d)–(f) NA=0.5. Symbols indicate μsddet values, and colors indicate phase function types. Note the log scales for both the reflectance and RpNA.


Table 2

Variability of RpNA, σ, and γ for μs′ddet=0.1 and μs′ddet=1, defined as the spread in RpNA, σ, and γ values for a chosen reflectance (±10%) relative to the total range of each parameter.




We have theoretically derived the parameter RpNA to model subdiffusive light scattering for overlapping source–detector geometries, and we have shown with Monte Carlo simulations that the reflectance depends on RpNA. In comparison to σ and γ, RpNA improves prediction of the reflectance. Our findings indicate that the reflectance does not depend on the details of the phase function—knowledge of all Legendre moments is not required—but only on the magnitude of the phase function within the acceptance angle of the detector (in the backward and forward directions).

Currently, the parameter γ is widely used to model subdiffusive reflectance. Since RpNA improves prediction of the measured reflectance, incorporating RpNA into subdiffusive models is expected to improve their reliability and thereby the estimation of other optical properties, such as μs. This improvement is significant because subdiffusive measurements have the potential to detect small-scale tissue changes.

For higher values of μsddet, the reflectance becomes more diffusive and, therefore, depends less on RpNA, σ, and γ. The contribution of diffuse photons is a remaining challenge for describing subdiffusive scattering. In our simulations, we kept the absorption coefficient constant. For subdiffusive reflectance, path lengths are short; therefore, the effect of μa will be minor. While the absorption affects the diffuse contribution to the reflectance, it has a minor effect on semiballistic photons. To implement RpNA for subdiffusive measurements, models have to be developed incorporating RpNA, the diffusive reflectance, and absorption for specific measurement geometries, such as single fiber reflectance spectroscopy.



In conclusion, the theoretically derived parameter RpNA predicts subdiffusive light scattering for overlapping source–detector geometries. Consequently, the reflectance does not depend on the details of the entire phase function, but on the phase function within the acceptance angles of the detector. Since RpNA improves prediction of the measured reflectance compared to σ and γ, the use of RpNA is expected to improve derivation of optical properties from subdiffusive measurements.


No conflicts of interest, financial or otherwise, are declared by the authors.


This research was supported by the Netherlands Organization for Scientific Research (Technology Foundation STW, iMIT-PROSPECT Grant No. 12707) and the Dutch Cancer Society/Alpe d’HuZes (Project No. 2014-7009). We would like to thank the Organization Committee of the Spinoza Chair from the Academic Medical Center of the University of Amsterdam for funding a “Spinoza” lectureship for Steven L. Jacques.


1. J. R. Mourant et al., “Influence of the scattering phase function on light transport measurements in turbid media performed with small source-detector separations,” Opt. Lett. 21, 546–548 (1996).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.21.000546 Google Scholar

2. A. Kienle, F. K. Forster and R. Hibst, “Influence of the phase function on determination of the optical properties of biological tissue by spatially resolved reflectance,” Opt. Lett. 26, 1571–1573 (2001).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.26.001571 Google Scholar

3. E. Vitkin et al., “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun. 2, 587–604 (2012). http://dx.doi.org/10.1038/ncomms1599 Google Scholar

4. M. Jia et al., “Virtual-source diffusion approximation for enhanced near-field modeling of photon-migration in low-albedo medium,” Opt. Express 23, 1337–1352 (2015).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.23.001337 Google Scholar

5. D. Piao and S. Patel, “Simple empirical master-slave dual-source configuration within the diffusion approximation enhances modeling of spatially resolved diffuse reflectance at short-path and with low scattering from a semi-infinite homogeneous medium,” Appl. Opt. 56, 1447 (2017).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.56.001447 Google Scholar

6. F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source-detector separations close to one transport mean free path,” J. Opt. Soc. Am. A 16, 2935 (1999).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.16.002935 Google Scholar

7. S. C. Kanick, H. J. C. M. Sterenborg and A. Amelink, “Empirical model of the photon path length for a single fiber reflectance spectroscopy device,” Opt. Express 17, 860–871 (2009).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.17.000860 Google Scholar

8. U. A. Gamm et al., “Quantification of the reduced scattering coefficient and phase-function-dependent parameter γ of turbid media using multidiameter single fiber reflectance spectroscopy: experimental validation,” Opt. Lett. 37, 1838–1840 (2012).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.37.001838 Google Scholar

9. S. C. Kanick et al., “Sub-diffusive scattering parameter maps recovered using wide-field high-frequency structured light imaging,” Biomed. Opt. Express 5, 3376 (2014).BOEICL2156-7085 http://dx.doi.org/10.1364/BOE.5.003376 Google Scholar

10. F. van Leeuwen-van Zaane et al., “In vivo quantification of the scattering properties of tissue using multi-diameter single fiber reflectance spectroscopy,” Biomed. Opt. Express 4, 696–708 (2013).BOEICL2156-7085 http://dx.doi.org/10.1364/BOE.4.000696 Google Scholar

11. N. Bodenschatz et al., “Model-based analysis on the influence of spatial frequency selection in spatial frequency domain imaging,” Appl. Opt. 54, 6725 (2015).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.54.006725 Google Scholar

12. P. Naglic et al., “Estimation of optical properties by spatially resolved reflectance spectroscopy in the subdiffusive regime,” J. Biomed. Opt. 21, 095003 (2016).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.21.9.095003 Google Scholar

13. N. Bodenschatz et al., “Quantifying phase function influence in subdiffusively backscattered light,” J. Biomed. Opt. 21, 035002 (2016).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.21.3.035002 Google Scholar

14. R. K. Wang and M. J. Wilson, “Vertex/propagator model for least-scattered photons traversing a turbid medium,” J. Opt. Soc. Am. A 18, 224–231 (2001). http://dx.doi.org/10.1364/JOSAA.18.000224 Google Scholar

15. S. A. Prahl et al., “A Monte Carlo model of light propagation in tissue,” in Dosimetry of Laser Radiation in Medicine and Biology, , G. Mueller and D. Sliney, Eds., Vol. IS, pp. 102–111, SPIE Press, Bellingham, Washington (1989). Google Scholar

16. L. V. Wang, S. L. Jacques and L. Zheng, “MCML–Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47, 131–146 (1995).CMPBEK0169-2607 http://dx.doi.org/10.1016/0169-2607(95)01640-F Google Scholar

17. J. R. Zijp and J. J. ten Bosch, “Use of tabulated cumulative density functions to generate pseudorandom numbers obeying specific distributions for Monte Carlo simulations,” Appl. Opt. 33, 533 (1994).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.33.000533 Google Scholar

18. P. Saccomandi et al., “Estimation of anisotropy coefficient of swine pancreas, liver and muscle at 1064 nm based on goniometric technique,” J. Biophotonics 8(5), 422–428 (2015). http://dx.doi.org/10.1002/jbio.v8.5 Google Scholar

19. L. O. Reynolds and N. J. McCormick, “Approximate two-parameter phase function for light scattering,” J. Opt. Soc. Am. 70, 1206–1212 (1980).JOSAAH0030-3941 http://dx.doi.org/10.1364/JOSA.70.001206 Google Scholar

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Anouk L. Post, Anouk L. Post, Steven L. Jacques, Steven L. Jacques, Henricus J. C. M. Sterenborg, Henricus J. C. M. Sterenborg, Dirk J. Faber, Dirk J. Faber, Ton G. van Leeuwen, Ton G. van Leeuwen, } "Modeling subdiffusive light scattering by incorporating the tissue phase function and detector numerical aperture," Journal of Biomedical Optics 22(5), 050501 (22 May 2017). https://doi.org/10.1117/1.JBO.22.5.050501 . Submission: Received: 14 February 2017; Accepted: 1 May 2017
Received: 14 February 2017; Accepted: 1 May 2017; Published: 22 May 2017

Back to Top