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–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

However, a serious problem when using the FT [Eq. (1)] arises when $F(\omega )$ exhibits a slow algebraic decay for increasing angular modulation frequency $\omega$. As a result, one has to deal with a highly oscillatory integrand, which is known to be difficult to integrate numerically. For example, we consider which is the known infinite-space Green’s function of the frequency-domain DE ${\nabla }^{2}G(\mathbf{r},\omega )-{\kappa }^{2}G(\mathbf{r},\omega )=-\delta (\mathbf{r})/D$.⁶ 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 This result can also be confirmed numerically via evaluation of Eqs. (1) for (2) using the discrete FT (DFT). However, especially for small source–detector separations, an accurate and efficient reconstruction of $G(\mathbf{r},t)$ from its image $G(\mathbf{r},\omega )$ becomes a challenging task. For illustration purposes, Fig. 1 shows the time-resolved fluence for a relatively small source–detector separation of $r=0.1\mathrm{mm}$ compared to $1/{\mu }_s^{\prime }$. The smooth curve denotes the exact time-domain solution [Eq. (3)], whereas the oscillating curve corresponds to the numerical inversion of Eq. (2) using the DFT with $N={10}^4$ discrete frequencies.

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

where $B$ denotes the Bromwich path, which is typically given by a line parallel to the imaginary axis of the $s$-plane. In particular, the evaluation of Eq. (4) along the imaginary axis with $s=i\omega$ results in the Fourier integral [Eq. (1)]. Significant improvements in view of convergence, efficiency, and numerical stability can be achieved when evaluating Eq. (4) along a Hankel contour, which starts and ends at $-\infty$. The key point in this context is that the kernel function $\exp (st)$ becomes a damped wave due to $\mathrm{Re}(s)<0$. Very recently, different contours have been developed and analyzed concerning the numerical inversion of the LT, such as the parabola and hyperbola⁷ as well as a modified Talbot contour.⁸ For our considerations, we found it appropriate to employ the hyperbola contour in the following parameter form: where $s^{\prime }(\theta )=i\mu \cosh (\theta +i\varphi )$. Inserting Eq. (5) into Eq. (4) leads to the modified inverse Laplace integral where we have implicitly taken into account a real time-domain signal, which means that $f(t)=f^{*}(t)$, hence $F(s^{*})=F^{*}(s)$. Then, the application of the midpoint rule results in the following computable series form: where $s_k=s({\theta }_k)$ for ${\theta }_k=(k+1/2)h$ and $h$ is the uniform node spacing with $N$ being the number of nodes along the hyperbola in the upper half-plane $\mathrm{Re}(s)>0$. The remaining curve parameters $\mu$ and $\varphi$ as well as the node spacing $h$ depend on the problem under consideration. This means that if $F(s)$ is computationally less expensive or if the computational time is not the highest priority, one should use the optimized parameters $\mu =4.492075287N/t$, $\varphi =1.172104229$, and $h=1.081792140/N$. These values are in principle the same as those given in Ref. 7, whereas the number of digits has been increased. As an example, we repeat the above numerical experiment for the same optical and geometrical parameters. Instead of a ${10}^4$-point DFT, we evaluate Eq. (7) for $N=15$ (!). In Fig. 2, the solid line corresponds to the exact time-domain solution [Eq. (3)], whereas the filled dots denote the reconstruction of $G(\mathbf{r},t)$ from its image $G(\mathbf{r},s)=\exp (-\kappa r)/(4\pi Dr)$ with $\kappa =\sqrt{{\mu }_a/D+s/(Dc)}$. We have additionally computed the relative differences between the exact and reconstructed solutions, which are depicted in the inset. It can be furthermore reduced by using a slightly increased number of samples. We note that the LT approach works in the same manner even for very small source–detector separations, such as $r={10}^{-6}\mathrm{mm}$.

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⁷

where $\mathrm{\Lambda }=t_2/t_1$. The uniform node spacing becomes $h=A(\varphi )/N$, and the value for the angle $\varphi =1.09$ has been recovered by us as a result of numerical experiments. For the next comparison, we consider the time-resolved reflectance $R(\rho ,t)$ from a two-layered medium predicted by means of the DE. The first layer has a thickness of $L=2\mathrm{mm}$, whereas the second layer is assumed to be infinitely thick. The refractive indices of the first and second layers are denoted by $n_1$ and $n_2$, respectively. The surrounding nonscattering medium in $z<0$ has a refractive index of $n_0=1.0$. We consider the DE under the extrapolated boundary condition. In this context, we refer to the papers^4,9–13 dealing with the solution of the DE in layered media. The corresponding solution to be considered here has the form where $R(q,s)$ denotes the reflectance in the Hankel–Laplace space, which can be found in Ref. 4, and $J_0(x)$ is the Bessel function of the first kind. For verification purposes, we take into account the real-space Green’s function derived by Tualle et al.¹⁰ 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.¹⁰ 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\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¹⁴

where ${\kappa }_x=\sqrt{{\mu }_{ax}/D_x+s/(D_{x}c)}$ and ${\kappa }_m=\sqrt{{\mu }_{am}/D_m+s/(D_{m}c)}$. The subscripts $x$ and $m$ refer to the excitation and emission wavelengths of the light. The corresponding solution in the time domain can be recovered by means of the inverse LT and written in the following closed-form representation: where $\lambda =\sqrt{({\mu }_{am}-{\mu }_{ax})/(D_m-D_x)}$ and $\mathrm{erfc}(z)$ are the complementary error function. Note that, depending on the optical properties, $\lambda$ can be both real and imaginary. Of course, the final result is in all cases automatically a real number. Figure 4 displays the time-resolved fluorescence calculated for two different source–detector separations. The solid line corresponds to Eq. (11), whereas the symbols denote the reconstructed fluorescence by means of the discrete inverse LT [Eq. (7)], which has been considered for $N=15$. In this case, the solution to be inverted [Eq. (10)] is computationally less expensive, so we used a variable hyperbola contour. The agreement between the exact and the reconstructed fluorescence is excellent, where the relative differences are in the range of $\approx {10}^{-14}$.

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\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.¹⁵ For the same conditions as above but with a typical distance of $\rho =2\mathrm{mm}$ and a measurement dynamic range of six orders of magnitude,¹⁵ 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.¹⁶ In order to support the implementation of the presented theory, we provide a small MATLAB script on the Internet.¹⁷