## Abstract

The inverse problem in three-layered scattering media is investigated using simulations. Instead of using the common diffusion approximation, the light propagation is modeled using the spherical harmonics ($P_N$) approximation in the time domain. The optical properties are determined by fitting the $P_3$ approximation to solutions obtained by the $P_7$ approximation, representing an almost exact solution to the radiative transfer equation. Poisson noise is added to the data to simulate time-correlated single photon counting measurements. It is shown that, with simulated two-distance measurements, it is possible to derive the optical properties accurately, especially the absorption coefficient of the third layer and the reduced scattering coefficient of the upper layer. In the case of unknown layer thicknesses, solutions with different parameters but very similar reflectance curves are found, which can lead to larger errors in the identified optical properties. For known layer thicknesses, the retrieval of all properties improves. Altogether, the use of the spherical harmonics approximation enables identifying optical properties from time domain reflectance data much more accurately than using the diffusion equation.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Light is used for diagnostic purposes in a wide range of biomedical applications. The light transport in human tissue is e.g. employed to reveal cancer and monitor cerebral hemodynamics [1,2]. Clinically relevant parameters are deduced from optical properties of the tissue which are determined using a theoretical model for the light propagation. A common model used for biological tissue is the radiative transfer equation (RTE) [3,4]. Analytical solutions to the RTE can be obtained applying different approximations. The most widely used one is the diffusion approximation, which can be derived from the RTE [5,6]. The diffusion equation has been solved for infinite, half space [7], layered and other geometries [8], but has limitations for short times, small source-detector separations and large absorption.

Multi-layered solutions play an important role in medical applications and have been of ongoing interest. Dayan *et al.* [9] analyzed the time and spatial dependence of the reflectance employing a two-layered model with approximated solutions to the diffusion equation. Using steady state and frequency domain reflectance data, the determination of the optical properties was examined by Kienle *et al.* [10] for a two-layered medium with an exact solution to the diffusion equation. They were able to determine the optical properties, but also point out the limited accuracy of the diffusion equation. Alexandrakis *et al.* [11] tested the accuracy of the same two-layered diffusion model and also indicated the need for a more accurate solution to the RTE. A theoretical model for many-layered diffusive media in spatial frequency domain was given by Ripoll *et al.* [12] and used for fitting optical properties. Liao and Tseng [13] proposed an iterative algorithm for recovering the optical properties of two- and three-layered media using a frequency domain diffusion model.

The Fourier counterpart to the frequency domain are time domain measurements. They provide more information than merely measuring at a single frequency and thus are used for investigating layered structures, as well. Kienle *et al.* [14] studied the time domain approach in more detail and also performed *in vivo* measurements [15]. Different procedures were proposed to extract the optical properties from time domain measurements [16,17]. With single-distance, time resolved measurements, Steinbrink *et al.* [18] were able to determine changes in the absorption using a layered Monte-Carlo approach. Martelli *et al.* [19] used multi-distance, time domain measurements to extract optical properties from phantoms. They also fitted a two-layered diffusion model to the reflectance data. The theory was extended to three layers by Martelli *et al.* [20]. Time domain measurements with a small source-detector separation were performed by Pifferi *et al.* [21] to detect brain activity. Alerstam *et al.* [22] combined the injection and detection of light into a single fiber for time-of-flight spectroscopy. Many of the above works used the diffusion approximation to describe the light propagation, since it can be calculated fast and easily. Others solved the RTE numerically with Monte-Carlo simulations [23] which is computationally more expensive and incorporates statistical noise.

The spherical harmonics ($P_N$) method offers a better approximation than the diffusion theory to solve the RTE. Especially, for increasing $N$, the $P_N$ solution approaches the exact RTE solution. With the method of rotated reference frames [24,25], it has recently been derived for a slab with matched boundary conditions [26], a semi-infinite medium using mismatched boundary conditions and an improved treatment of the incident light beam [27,28] and layered geometries [29]. The work of Hull and Foster [30] considers the $P_3$ approximation with a simplified boundary condition. Gardner *et al.* [31] solve the $P_N$ equations with a Fourier decomposition method. Chin *et al.* [32] employed the $P_3$ approximation to calculate the steady state, angularly resolved radiance for the determination of the optical properties in the case of a spherically symmetric, infinite turbid medium. The incomplete $P_5$ approximation was used by Liu *et al.* [33] within a very similar setup.

In this work, the determination of the optical properties of three-layered media is investigated based on simulated time domain measurements. This is accomplished by fitting relevant optical parameters to simulated data using the $P_3$ approximation. It is shown that the determination of the optical properties of a three-layered medium is in principle possible and superior to a model using diffusion theory. Under certain conditions some challenges remain.

## 2. Method

In Fig. 1 a schematic view of the simulated time domain measurement is shown. A pulsed laser beam is incident on a laterally-infinite, three-layered medium. An infinitely short laser pulse is simulated with a Gaussian spatial beam profile

where the beam radius was chosen as $\rho _{\mathrm {w}}= {0.75}\,{\textrm{mm}}$. The first two layers have finite thicknesses $l_1$ and $l_2$, respectively. The third layer extends infinitely. All layers are assumed to have the same refractive index $n = 1.4$ and anisotropy factor $g = 0.8$. Above the top layer, air with $n = 1.0$ is presumed. Different refractive indices for the individual layers are not feasible with the current version of the employed model. Changes in the refractive index between the layers complicate the boundary conditions and impose the problem of additional reflections amid the boundaries. The layers may differ in their thicknesses, in the absorption coefficients $\mu _{\mathrm {a}i}$ and the reduced scattering coefficients $\mu ^\prime _{\mathrm {s}i}$ ($i \in \{1,2,3\}$). For a given distance $\rho$ from the source and a given time $t$, the reflectance at the surface $R(t,\rho )$ can be calculated using the analytical solutions derived by Liemert*et al.*[29]. This approach uses the $P_N$ method [24] to obtain the solution to RTE, where $N$ denotes the order of the approximation. The $P_N$ approximation approaches the exact solution for high orders $N$. Already for orders $N \geq 3$ it offers significant improvements over the diffusion approximation [29], especially at short times. For the given configuration, choosing $N=7$ results in a practically exact solution to the RTE.

Using multiple points in time $t_1$,…, $t_M$ a time dependent reflectance measurement can be simulated. These time domain curves depend on the optical properties of the layers. The goal is to recover the unknown optical properties from a time dependent measurement. Therefore, a curve fitting algorithm using the Ceres Solver library [34] was implemented. It can fit the reflectance curves in the time domain adapting the optical parameters.

To investigate the retrieval of the optical properties from time-resolved reflectance curves, 1000 combinations of different parameters were created by randomly choosing values within the ranges of typical biological tissue [35] shown in Table 1. For all data sets, a time-resolved reflectance curve at times ranging from ${5\times 10^{-11}}\,{\textrm{s}}$ to ${1.2\times 10^{-8}}\,{\textrm {s}}$, in steps of ${1\times 10^{-11}}\,{\textrm{s}}$, was calculated for two distances $\rho _1= {8}\,{\textrm {mm}}$ and $\rho _2={16}\,{\textrm{mm}}$ using the $P_7$ approximation. In time-correlated single photon counting measurements the number of photons per time bin is counted. This leads to a Poisson distribution in the measured counts per bin [36]. To simulate such measurements, the reflectance curves were scaled to have a certain maximum count rate $R_{\mathrm {max}}$ (typically $R_{\mathrm {max}}= 1 \times 10^{5}$ counts). To account for the statistical nature of the measurement process, Poisson noise was added to the curves. For each point in time, a random number was generated following the Poisson distribution with the expectation value equal to the scaled reflectance at that time, resulting in a noisy reflectance curve. In real experiments, effects of the measuring system like the impulse response function (IRF) need to be considered. Those are neglected in this study. It was not the goal of this work to study the influence of the IRF on the determination of the optical properties. In this sense, the ideal case without the influence of the measuring system was investigated. Convolution with the IRF may considerably deteriorate the results.

After the data were generated, the unknown parameters were fitted to the curves using the $P_3$ approximation. Although the data generation and the fitting both rely on the $P_N$ approximation, the systematic error due to the incomplete approximation is strongly reduced with the higher order $N=7$. For the given optical properties and geometry several time resolved reflectance curves were compared to Monte-Carlo-simulations which confirmed the negligible error on the $P_7$ approximation. For other configurations, like angularly resolved radiance comparison, higher orders may still be needed. Only the part of the reflectance curve starting at $0.8 R_{\mathrm {max}}$ before the maximum value up to $10^{-3} R_{\mathrm {max}}$ after the maximum was taken into account. The fit parameters are the optical properties and optionally the layer thicknesses and scaling factors $k$. As shown in Table 1, the starting values of the fit parameters were fixed within the allowed interval. The least square fitting algorithm was also restricted to vary these parameters only within the range given in Table 1 to reduce the total cost $\chi ^2$:

where $R_{s,j}$ is the noisy, simulated data and $R_{t,j}$ is the predicted value by the $P_3$-solution for a given data point. $M$ is the number of data points in the data set and $P$ is the number of fit parameters. The weights $w_j$ of the least squares are estimated to be the square root of the data points of the scaled, noisy curves, due to the Poisson statistic. They are explicitly taken from the noisy data to resemble measurements, where the weights also have to be estimated from the noisy signal. The properly chosen weights lead to a total cost $\chi ^2$ that has a value around unity, if the correct fit parameters are found.After fitting all 1000 different parameters sets, a post selection is performed. First, only problems where the algorithm converged within 100 iterations are selected. Second, the cost $\chi ^2$ of the curve fitting problem is considered. As stated above, the cost should be around 1. But some problems still have a cost which is $> 1.5$. This is a hint that the algorithm was not able to find the global minimum and those problems are excluded as well. If there are systematic deviations in the model, even for the global minimum, the cost $\chi ^2$ can be significantly larger than 1.5 (e.g. for the diffusion theory). Anyway, we opted for the $\chi ^2\leq 1.5$ criterion to identify the convergence to global minima, with the downside of possibly missing some useful information. For a model without systematic deviations this criterion can substantially improve the results by excluding fits that converged to a local minimum with very different reflectance curves. We will refer to these problems, where the fit algorithm converged and the scaled cost was $\leq 1.5$, as “good problems” or “good fits”. The evaluation of the errors in the fit parameters will consider only those “good problems”.

To evaluate the curve fitting problems, the absolute values $e_{\mathrm r}$ of the relative errors of the fit parameters of each problem $i$ were calculated in the following way:

where $p_{\mathrm r}$ is the true value and $p_{\mathrm f}$ is the found value of a parameter. Afterwards, the mean, the standard deviation $\sigma$ and the median of the $e_{\mathrm r}$ were taken over all “good problems”.## 3. Results and discussion

#### 3.1 Fits including layer thicknesses

Fitting both layer thicknesses, the absorption coefficients and reduced scattering coefficients, the parameters can be obtained reasonably well. From the 1000 fitting problems 828 fulfilled the criteria for “good fits” with $R_{\mathrm {max}}= 1\times 10^{5}$, the errors in the fit parameters are shown in Table 2. The reduced scattering coefficient of the first layer and the absorption coefficient of the third layer can be obtained with mean errors less than ${2}\,{\% }$, the other properties are more difficult to retrieve. Especially the layer thickness and the absorption coefficient of the second layer show large errors. It is also noticeable that the standard deviation can be quite large showing that some problems with very wrong parameter estimates exist. In practice, often only relative reflectance data are measured. For the absolute fit, the amplitude of the measured reflectance signal must be known. This can be difficult in practice since a reference must be measured to incorporate the energy of the laser pulse. Often, the energy is not constant over time which further complicates the measurement of an absolute reflectance signal. With relative measurements, only the shape of the signal is obtained (up to a scaling factor). Therefore, two scaling parameters were included in the fitting process (one for each distance).

Using these relative fits, the retrieval of the parameters deteriorated. The reduced scattering coefficient of the first layer and the absorption coefficient of the third layer can still be obtained with fair accuracy.

#### 3.2 Fits for fixed layer thicknesses

When the layer thicknesses are known *a priori*, the other parameters are obtained quite well in most cases for absolute fits. From the 1000 fitting problems 995 fulfilled the criteria for “good fits” with $R_{\mathrm {max}}= 1\times 10^{5}$ and the results are shown in Table 3. The parameters with the best results are still $\mu _{\mathrm {a} {3}}$ and $\mu _{\mathrm {s} {1}}^\prime$. If relative fits are performed, only the reduced scattering coefficient of the first and second layer and the absorption coefficient of third layer show very small relative errors for most fitting problems. In this case 899 problems directly converged as “good fits”.

All the simulations were also preformed with a lower level of noise, which improved the results. The corresponding data can be found in Appendix A.1, Table 5 and 6.

Since the diffusion approximation is widely used for describing light propagation in tissue, the same data were also fitted using a three-layered diffusion model [37]. Due to systematic discrepancy of the diffusion model to the $P_7$ approximation, only few of the fitting problems showed a $\chi ^2\leq 1.5$ for absolute fits. For relative fits, the majority of the fit problems can satisfy the “good fits” criteria, after all. Nevertheless, the fit parameters are obtained with a significantly lower accuracy. The data are shown in Appendix A.2, Table 7. The given source-detector separations limit the validity of the diffusion approximation within the given range of optical properties. For larger distance or longer times, the diffusion theory approximates the RTE more accurately. Even for those cases, the $P_3$ approximation is still superior, albeit not as substantial. As the shorter distances and times contain more information about the upper layers, switching to larger distances is not always advisable, since it can hinder the retrieval of the optical properties of the top layers.

The mean computation time for one inversion varied with the number of fit parameters. On a server with a total of 16 threads, fitting only the optical properties took around 4 minutes per problem using the $P_N$ approximation with $N=3$. With an explicit $P_3$ solution [38], the computation time can be further reduced by a factor of nearly 5. The diffusion approximation is still over an order of magnitude faster to compute.

The results indicate that determination of the optical properties with the diffusion equation is mainly limited by its inaccurate approximation of the RTE. With the $P_3$ solution, although still an approximation, the close determination of the optical properties is substantially hindered by the remaining level of noise in the data. This demonstrates that the $P_3$ approximation has a substantial advantage over the diffusion theory for the considered setup, since the leftover errors in the fit parameters are not significantly caused by a systematic discrepancy.

#### 3.3 Ambiguity problem

The above errors in the fitted parameters may show a large spread, which is indicated by the high values for the standard deviation as well as the difference between the mean and the median values. In Fig. 2 the distribution of the $e_{\mathrm r}$ in $\mu _{\mathrm {a} {3}}$ is illustrated in a histogram for relative fits with unknown layer thicknesses. Although the majority of the found absorption parameters has an absolute relative error clearly below 5 %, the largest error found is 42 %. Further investigation of the specific fitting problems that have large errors in the found parameters shows that a second set of parameters may exist that exhibit nearly the same reflectance curves (as shown in Fig. 3 and Table 4). This is mainly the case, if more fit parameters are involved (e.g. layer thicknesses and scaling parameters). Similar problems were already encountered by Pham *et al.* [39].

Another source of error is a low sensitivity of the reflectance data in some fit parameters. If e.g. the top layers are thick and highly absorbing and scattering, the properties of the lower layer have little influence on the obtained reflectance. Even for the case with the best results (known layer thicknesses and absolute reflectance curves) the errors in $\mu _{\mathrm {a} {3}}$ can reach up to 18.8 %.

## 4. Summary and conclusion

The determination of the optical properties of a three-layered turbid medium using simulated time domain measurements was investigated. Using the three-layered $P_3$-approximation, time domain reflectance data were fitted at two source-detector separations simultaneously. At large, it was possible to determine the optical properties, but not for all configurations with high accuracy. Specifically for unknown layer thicknesses and relative fits, the absorption coefficients of the first two layers and the reduced scattering coefficient of the third layer were difficult to retrieve with low uncertainty. If the layer thicknesses were known, the accuracy in the determination of the optical properties improved. In addition, it was shown that absolute reflectance measurements can further enhance the results. The absorption coefficient of the deepest layer is effectively identified with good accuracy. If the first two layers are merely surface layers, where one is not interested in the absorption coefficients, this can readily be exploited in many applications. One main caveat that was found during the investigation is the ambiguity of the time-resolved reflectance curves within the given level of noise. For some sets of parameters, a second set of parameters might be found that results in nearly the same reflectance curves at the considered distances. Recently, García *et al.* [40] presented a state-estimation instead of a curve fitting algorithm to approach this ambiguity problem by providing prior information. They used the diffusion approximation, but underline the generality of the inverse method. The $P_3$ approximation is an adequate substitution in this case. It does display much lower systematic approximation errors and is also suited in the regime were the diffusion approximation fails (short times and small distances from the source). For real life applications other challenges remain. The impulse response function of the measuring device has to be taken into account and a different geometry compared to the assumed layered model may lead to further discrepancies. Nevertheless, the $P_3$ approximation is a promising tool for time-resolved reflectance measurements and offers a large improvement over the diffusion theory.

## A. Appendix

## A.1 Results with lower noise

In addition to the above used $R_{\mathrm {max}}= 1\times 10^{5}$, the same reflectance curves were additionally scaled to $R_{\mathrm {max}} = 1\times 10^{6}$ before adding the Poisson noise. The reduced level of noise led to an improved determination of the optical properties. In Table 5 the fits including the layer thicknesses are shown, corresponding to Table 2. Table 6 contains the results for fixed layer thicknesses in comparison to Table 3. The low number of “good fits” in the case of relative fits with unknown layer thicknesses is the result of the fit algorithm converging to local minima of the cost function $\chi ^2$. Repeating the fit with different start values increases the number “good fits” significantly.

## A.2 Results with diffusion theory

The time resolved reflectance data were furthermore fitted using a three-layered diffusion model [37]. Only after several times changing the start values randomly, the majority of the relative fits had a $\chi ^2\leq 1.5$. As seen in Table 7, the optical properties are obtained with a significantly lower accuracy compared to the results obtained with the $P_3$ solution. The noise was the same as in the main part of the work ($R_{\mathrm {max}}= 1\times 10^{5}$), relative fits were performed.

## Funding

Deutsche Forschungsgemeinschaft (DFG); Carl-Zeiss-Stiftung.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. **73**(7), 076701 (2010). [CrossRef]

**2. **T. Myllylä, M. Harju, V. Korhonen, A. Bykov, V. Kiviniemi, and I. Meglinski, “Assessment of the dynamics of human glymphatic system by near-infrared spectroscopy,” J. Biophotonics **11**(8), e201700123 (2018). [CrossRef]

**3. **S. Chandrasekhar, * Radiative Transfer* (Dover Publications, 1960).

**4. **W. F. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. **26**(12), 2166–2185 (1990). [CrossRef]

**5. **A. Ishimaru, “Diffusion of light in turbid material,” Appl. Opt. **28**(12), 2210–2215 (1989). [CrossRef]

**6. **M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “Diffusion approximation revisited,” J. Opt. Soc. Am. A **26**(5), 1291–1300 (2009). [CrossRef]

**7. **R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A **11**(10), 2727–2741 (1994). [CrossRef]

**8. **S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. **37**(7), 1531–1560 (1992). [CrossRef]

**9. **I. Dayan, S. Havlin, and G. H. Weiss, “Photon Migration in a Two-layer Turbid Medium a Diffusion Analysis,” J. Mod. Opt. **39**(7), 1567–1582 (1992). [CrossRef]

**10. **A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. **37**(4), 779–791 (1998). [CrossRef]

**11. **G. Alexandrakis, T. J. Farrell, and M. S. Patterson, “Accuracy of the diffusion approximation in determining the optical properties of a two-layer turbid medium,” Appl. Opt. **37**(31), 7401–7409 (1998). [CrossRef]

**12. **J. Ripoll, V. Ntziachristos, J. P. Culver, D. N. Pattanayak, A. G. Yodh, and M. Nieto-Vesperinas, “Recovery of optical parameters in multiple-layered diffusive media: theory and experiments,” J. Opt. Soc. Am. A **18**(4), 821–830 (2001). [CrossRef]

**13. **Y.-K. Liao and S.-H. Tseng, “Reliable recovery of the optical properties of multi-layer turbid media by iteratively using a layered diffusion model at multiple source-detector separations,” Biomed. Opt. Express **5**(3), 975–989 (2014). [CrossRef]

**14. **A. Kienle, T. Glanzmann, G. Wagnières, and H. van den Bergh, “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt. **37**(28), 6852–6862 (1998). [CrossRef]

**15. **A. Kienle and T. Glanzmann, “In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,” Phys. Med. Biol. **44**(11), 2689–2702 (1999). [CrossRef]

**16. **F. Martelli, A. Sassaroli, Y. Yamada, and G. Zaccanti, “Method for measuring the diffusion coefficient of homogeneous and layered media,” Opt. Lett. **25**(20), 1508–1510 (2000). [CrossRef]

**17. **F. Martelli, S. D. Bianco, and G. Zaccanti, “Procedure for retrieving the optical properties of a two-layered medium from time-resolved reflectance measure ments,” Opt. Lett. **28**(14), 1236–1238 (2003). [CrossRef]

**18. **J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. **46**(3), 879–896 (2001). [CrossRef]

**19. **F. Martelli, S. D. Bianco, G. Zaccanti, A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, and R. Cubeddu, “Phantom validation and in vivo application of an inversion procedure for retrieving the optical properties of diffusive layered media from time-resolved reflectance measurements,” Opt. Lett. **29**(17), 2037–2039 (2004). [CrossRef]

**20. **F. Martelli, A. Sassaroli, S. D. Bianco, and G. Zaccanti, “Solution of the time-dependent diffusion equation for a three-layer medium: application to study photon migration through a simplified adult head model,” Phys. Med. Biol. **52**(10), 2827–2843 (2007). [CrossRef]

**21. **A. Pifferi, A. Torricelli, L. Spinelli, D. Contini, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, A. Dalla Mora, F. Zappa, and S. Cova, “Time-Resolved Diffuse Reflectance Using Small Source-Detector Separation and Fast Single-Photon Gating,” Phys. Rev. Lett. **100**(13), 138101 (2008). [CrossRef]

**22. **E. Alerstam, T. Svensson, S. Andersson-Engels, L. Spinelli, D. Contini, A. D. Mora, A. Tosi, F. Zappa, and A. Pifferi, “Single-fiber diffuse optical time-of-flight spectroscopy,” Opt. Lett. **37**(14), 2877–2879 (2012). [CrossRef]

**23. **A. H. Hielscher, H. Liu, B. Chance, F. K. Tittel, and S. L. Jacques, “Time-resolved photon emission from layered turbid media,” Appl. Opt. **35**(4), 719–728 (1996). [CrossRef]

**24. **V. A. Markel, “Modified spherical harmonics method for solving the radiative transport equation,” Waves Random Media **14**(1), L13–L19 (2004). [CrossRef]

**25. **G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiative transport equation in rotated reference frames,” J. Phys. A: Math. Gen. **39**(1), 115–137 (2006). [CrossRef]

**26. **M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “The Green’s function for the radiative transport equation in the slab geometry,” J. Phys. A: Math. Theor. **43**(6), 065402 (2010). [CrossRef]

**27. **A. Liemert and A. Kienle, “Light transport in three-dimensional semi-infinite scattering media,” J. Opt. Soc. Am. A **29**(7), 1475–1481 (2012). [CrossRef]

**28. **A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep. **3**(1), 2018 (2013). [CrossRef]

**29. **A. Liemert, D. Reitzle, and A. Kienle, “Analytical solutions of the radiative transport equation for turbid and fluorescent layered media,” Sci. Rep. **7**(1), 3819 (2017). [CrossRef]

**30. **E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the ${\textrm {P}}_{{3}}$ approximation,” J. Opt. Soc. Am. A **18**(3), 584–599 (2001). [CrossRef]

**31. **A. R. Gardner, A. D. Kim, and V. Venugopalan, “Radiative transport produced by oblique illumination of turbid media with collimated beams,” Phys. Rev. E **87**(6), 063308 (2013). [CrossRef]

**32. **L. C. L. Chin, A. E. Worthington, W. M. Whelan, and I. A. Vitkin, “Determination of the optical properties of turbid media using relative interstitial radiance measurements: Monte Carlo study, experimental validation, and sensitivity analysis,” J. Biomed. Opt. **12**(6), 064027 (2007). [CrossRef]

**33. **L. Liu, W. Wan, J. Li, H. Zhao, and F. Gao, “Simultaneous recovery of a full set of optical properties in turbid media using incomplete P5 approximation to CW radiance,” Opt. Lett. **43**(17), 4188–4191 (2018). [CrossRef]

**34. **S. AgarwalK. Mierle, *et al.*, “Ceres Solver,” http://ceres-solver.org.

**35. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**(11), R37–R61 (2013). [CrossRef]

**36. **D. V. O’Connor and D. Phillips, * Time-correlated Single Photon Counting* (Academic Press, 1984).

**37. **A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt. **15**(2), 025002 (2010). [CrossRef]

**38. **A. Liemert and A. Kienle, “Explicit solutions of the radiative transport equation in the P3 approximation,” Med. Phys. **41**(11), 111916 (2014). [CrossRef]

**39. **T. H. Pham, T. Spott, L. O. Svaasand, and B. J. Tromberg, “Quantifying the properties of two-layer turbid media with frequency-domain diffuse reflectance,” Appl. Opt. **39**(25), 4733–4745 (2000). [CrossRef]

**40. **H. García, G. Baez, and J. Pomarico, “Simultaneous retrieval of optical and geometrical parameters of multilayered turbid media via state-estimation algorithms,” Biomed. Opt. Express **9**(8), 3953–3973 (2018). [CrossRef]