In this paper, we consider the Laplace transform (LT) for solving different time-dependent photon migration problems occurring in the biomedical optics field. It is shown that the LT exhibits important advantages in view of accuracy, efficiency, and numerical stability compared to the classical approach that uses the Fourier transform (FT) to obtain time-dependent quantities from data in the frequency domain. For typical applications in tissue spectroscopy or imaging, a speed-up of up to several orders of magnitude can be accomplished by applying the LT for both numerical or analytical solution approaches.

Modeling of light propagation in scattering media, such as biological tissue, in the mesoscopic and macroscopic scales is commonly performed using the radiative transport equation (RTE) or its approximation, the diffusion equation (DE). Analytical solutions of these equations in the time domain are restricted to relative simple geometries,^{1}2.^{–}^{3} whereas for a series of applications, efficient analytical solutions are available in the frequency domain.^{4}^{,}^{5} Usually, the FT is applied to obtain the time-domain solutions from the corresponding solutions in the frequency domain. Similarly, in the case of numerical solutions of these equations, calculations in the frequency domain are more efficient than in the time domain. Again, commonly, the FT is applied to obtain solutions in the time domains.

In this paper, we show, using exemplary analytical solutions of the DE, that the application of the LT is considerably more efficient and accurate compared to the use of the FT for obtaining time-domain solutions. We present comparisons for the fluence in an infinitely extended medium, for the reflectance from a two-layered medium, and for the fluorescence in an infinitely extended medium.

Analytical solutions of the time-dependent DE are provided either by the separation of variables method or the integral transforms, such as the FT, which is defined as

## (1)

$$f(t)=\frac{1}{2\pi}{\int}_{-\infty}^{\infty}F(\omega ){\mathrm{e}}^{i\omega t}\text{\hspace{0.17em}}\mathrm{d}\omega ,\phantom{\rule[-0.0ex]{2em}{0.0ex}}\omega \in \mathbb{R}.$$## (2)

$$G(r,\omega )=\frac{\mathrm{exp}(-\kappa r)}{4\pi Dr},\phantom{\rule[-0.0ex]{2em}{0.0ex}}\kappa =\sqrt{\frac{{\mu}_{a}}{D}+\frac{i\omega}{Dc}},$$^{6}In this context, we have ${\mu}_{a}$ as the absorption coefficient and $D=1/[3({\mu}_{a}+{\mu}_{s}^{\prime})]$ as the diffusion coefficient with ${\mu}_{s}^{\prime}$ being the reduced scattering coefficient. Furthermore, $\omega $ is the angular modulation frequency and $c$ is the speed of light in the medium. The exact solution in the time domain is well known and given by

^{1}

^{,}

^{2}

## (3)

$$G(r,t)=\frac{c}{(4\pi Dct{)}^{3/2}}\mathrm{exp}(-{\mu}_{a}ct-\frac{{r}^{2}}{4Dct}),\phantom{\rule[-0.0ex]{2em}{0.0ex}}t>0.$$Indeed, the above result can be improved by increasing the number of nodes at the expense of the computational time. This circumstance, however, leads to serious problems in many situations of high practical importance. An example of this is the reconstruction of the optical properties in experiments, which involves the solution of the inverse problem. Here, the solutions of the DE must be evaluated many times, so the calculation time becomes a critical parameter. On the other hand, when applying the time-dependent RTE, the solution of the forward problem can already be critical, because the computation of the frequency-domain data is significantly more expensive than in the case of the DE. Fortunately, besides the FT, one can alternatively make use of the LT

## (4)

$$f(t)=\frac{1}{2\pi i}{\int}_{B}F(s){\mathrm{e}}^{st}\text{\hspace{0.17em}}\mathrm{d}s,\phantom{\rule[-0.0ex]{2em}{0.0ex}}s\in \mathbb{C},$$^{7}as well as a modified Talbot contour.

^{8}For our considerations, we found it appropriate to employ the hyperbola contour in the following parameter form:

## (5)

$$s(\theta )=\mu +i\mu \text{\hspace{0.17em}}\mathrm{sinh}(\theta +i\varphi ),\phantom{\rule[-0.0ex]{2em}{0.0ex}}-\infty <\theta <\infty ,$$## (6)

$$f(t)=\frac{1}{\pi}\mathrm{Im}({\int}_{0}^{\infty}F[s(\theta )]{\mathrm{e}}^{s(\theta )t}{s}^{\prime}(\theta )\mathrm{d}\theta ),$$## (7)

$$f(t)=\frac{h}{\pi}\mathrm{Im}[\sum _{k=0}^{N-1}F({s}_{k})\mathrm{exp}({s}_{k}t){s}_{k}^{\prime}],$$We now consider the case that $F(s)$ is an arbitrary complicated function that is costly to evaluate. Thus, the number of its evaluation should be kept as small as possible. At the same time, we assume that $f(t)$ is needed for many time values in $t\in ({t}_{1},{t}_{2})$. In this case, we use a fixed integration path in connection with the following time-independent parameters^{7}

## (8)

$$\mu =\frac{4\pi \varphi -{\pi}^{2}}{A(\varphi )}\frac{N}{{t}_{2}},\phantom{\rule{0ex}{0ex}}A(\varphi )=\mathrm{arcosh}\left[\frac{(\pi -2\varphi )\mathrm{\Lambda}+4\varphi -\pi}{(4\varphi -\pi )\mathrm{sin}\text{\hspace{0.17em}}\varphi}\right],$$^{4}

^{,}

^{9}10.11.12.

^{–}

^{13}dealing with the solution of the DE in layered media. The corresponding solution to be considered here has the form

## (9)

$$R(\rho ,t)=\frac{1}{2\pi}{\int}_{0}^{\infty}\frac{1}{2\pi i}[{\int}_{B}R(q,s){\mathrm{e}}^{st}\text{\hspace{0.17em}}\mathrm{d}s]{J}_{0}(q\rho )q\text{\hspace{0.17em}}\mathrm{d}q,$$^{10}Figure 3 displays the time-resolved reflectance from a two-layered medium evaluated for two radial distances relative to the normally incident light beam. The solid line is the reflectance obtained from Eq. (9), whereas the symbols correspond with the real-space Green’s function of Tualle et al.

^{10}Concerning the inverse LT, we consider a fixed hyperbola contour and truncate the series [Eq. (7)] at $N=30$. As a consequence, the numerical evaluation of the time-resolved reflectance can be performed quite efficiently. For example, in the case of a three-layered model, the reflectance for 20 radial distances and 600 time values ($20\times 600$ grid) can be evaluated in $\approx 0.2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ using a small MATLAB script. The computation time can be furthermore reduced when using, e.g., the C programming language instead of MATLAB. Figure 3 confirms the good agreement between the different solutions for both radial distances. We note that the evaluation of the time-resolved reflectance at null source–detector separation would be a difficult task when using the conventional FT.

As a final example, we consider the time-resolved fluorescence caused by an isotropic point light source embedded in an infinite scattering medium. The fluorescence in Laplace space is given by^{14}

## (10)

$$F(r,s)=\frac{{\mu}_{ax}}{4\pi {D}_{m}{D}_{x}}\frac{1}{{\kappa}_{x}^{2}-{\kappa}_{m}^{2}}(\frac{{\mathrm{e}}^{-{\kappa}_{m}r}}{r}-\frac{{\mathrm{e}}^{-{\kappa}_{x}r}}{r}),$$## (11)

$$F(r,t)=\frac{c}{8\pi r}\frac{{\mu}_{ax}}{{D}_{m}-{D}_{x}}\mathrm{exp}(-\frac{{\mu}_{ax}{D}_{m}-{\mu}_{am}{D}_{x}}{{D}_{m}-{D}_{x}}ct)\times [\mathrm{exp}(-\lambda r)\mathrm{erfc}(\frac{1}{2}\frac{r}{\sqrt{{D}_{m}ct}}-\lambda \sqrt{{D}_{m}ct})+\mathrm{exp}(+\lambda r)\mathrm{erfc}(\frac{1}{2}\frac{r}{\sqrt{{D}_{m}ct}}+\lambda \sqrt{{D}_{m}ct})-\mathrm{exp}(-\lambda r)\mathrm{erfc}(\frac{1}{2}\frac{r}{\sqrt{{D}_{x}ct}}-\lambda \sqrt{{D}_{x}ct})-\mathrm{exp}(+\lambda r)\mathrm{erfc}(\frac{1}{2}\frac{r}{\sqrt{{D}_{x}ct}}+\lambda \sqrt{{D}_{x}ct}\left)\right],$$In this paper, we pointed out a recently developed algorithm for accurate, efficient, and reliable inversion of the LT, which overcomes some severe drawbacks of the classical approach in optical spectroscopy and imaging using frequency-domain data in connection with the FT. Numerical experiments have been carried out in order to demonstrate the resulting advantages in view of accuracy and numerical stability compared to the conventional FT. Typical speed-ups in computation time can amount to orders of magnitude. For example, in a typical time-domain experiment, the reflectance or transmittance signal is measured over three orders of magnitude. With the classical approach using the FT, one needs about 400 frequency values to achieve a relative accuracy for all time values of at least 0.1% for typical optical properties of biological tissue at a radial distance of $\rho =10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$. Contrarily, by applying the LT as described in this study, 15 frequency values are sufficient, which results in a speed-up of approximately one order of magnitude. In addition, the relative accuracy of the time-domain data is much better for the LT approach compared to the classical FT method. As a further example, we considered the time-domain reflectance experiments at small distances, which enable imaging with higher resolution.^{15} For the same conditions as above but with a typical distance of $\rho =2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ and a measurement dynamic range of six orders of magnitude,^{15} 5000 frequency values are needed for the FT to achieve an accuracy of at least one per mill, whereas with the LT approach, 25 frequency values are sufficient, which results in a speed-up of about two orders of magnitude. Further, the presented algorithm can be implemented for reconstructing time-domain quantities from frequency-domain data not only for analytical but also for numerical solutions of the DE such as with the finite element method, which is often used in optical imaging, accelerating these calculations by the same factors. In addition, the proposed approach can be applied for analytical and numerical solutions of the RTE, where the above-mentioned speed-ups are even more important due to the longer calculation times needed for solving the RTE compared to the DE. In addition, we derived the fundamental solution for the time-resolved fluence of fluorescence light in closed-form representation. In the literature, analytical solutions to this problem are available only for the case of equal reduced scattering coefficients of the incident light and the fluorescence light.^{16} In order to support the implementation of the presented theory, we provide a small MATLAB script on the Internet.^{17}

## Acknowledgments

The authors would like to thank Jean-Michel Tualle for providing numerical data for the time-resolved reflectance from the two-layered medium. One of the authors (A.L.) gratefully acknowledges the financial support by the Carl Zeiss Foundation, Germany.

## References

M. S. Patterson, B. Chance and C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. 28, 2331–2336 (1989).APOPAI0003-6935http://dx.doi.org/10.1364/AO.28.002331Google Scholar

F. Martelli et al., Light Propagation Through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software, SPIE Press, Bellingham, WA (2010).Google Scholar

S. R. Arridge, M. Cope and D.T. Delpy, “The theoretical basis for the determination of optical pathlength in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560 (1992).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/37/7/005Google Scholar

A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt. 15, 025002 (2010).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3368682Google Scholar

A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep. 3, 2018 (2013).SRCEC32045-2322http://dx.doi.org/10.1038/srep02018Google Scholar

R. C. Haskell et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.11.002727Google Scholar

J. Weideman and L. Trefethen, “Parabolic and hyperbolic contours for computing the Bromwich integral,” Math. Comput. 76, 1341–1356 (2007).MCMPAF0025-5718http://dx.doi.org/10.1090/S0025-5718-07-01945-XGoogle Scholar

B. Dingfelder and J. A. C. Weideman, “An improved Talbot method for numerical Laplace transform inversion,” Numer. Algorithms 68, 167–183 (2013).NUALEG1017-1398http://dx.doi.org/10.1007/s11075-014-9895-zGoogle Scholar

A. Kienle et al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37, 779–791 (1998).APOPAI0003-6935http://dx.doi.org/10.1364/AO.37.000779Google Scholar

J. M. Tualle et al., “Real-space Green’s function calculation for the solution of the diffusion equation in stratified turbid media,” J. Opt. Soc. Am. A 17, 2046–2055 (2000).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.17.002046Google Scholar

J. M. Tualle et al., “Asymptotic behavior and inverse problem in layered scattering media,” J. Opt. Soc. Am. A 21, 24–34 (2004).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.21.000024Google Scholar

F. Martelli et al., “Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfunction method,” Phys. Rev. E 67, 056623 (2003).http://dx.doi.org/10.1103/PhysRevE.67.056623Google Scholar

F. Martelli et al., “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, 2827–2843 (2007).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/52/10/013Google Scholar

T. J. Farrell, M. S. Patterson, “Diffusion modelling of fluorescence in tissue,” in Handbook of Biomedical Fluorescence, , M. A. Mycek and B. P. Pogue, Eds., CRC Press, Boca Raton, FL (2003).Google Scholar

A. Pifferi et al., “Time-resolved diffuse reflectance using small source–detector separation and fast single-photon gating,” Phys. Rev. Lett. 100, 138101 (2008).PRLTAO0031-9007http://dx.doi.org/10.1103/PhysRevLett.100.138101Google Scholar

M. S. Patterson and B. W. Pogue, “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt. 33, 1963–1974 (1994).APOPAI0003-6935http://dx.doi.org/10.1364/AO.33.001963Google Scholar