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,12.–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 as6 In this context, we have as the absorption coefficient and as the diffusion coefficient with being the reduced scattering coefficient. Furthermore, is the angular modulation frequency and is the speed of light in the medium. The exact solution in the time domain is well known and given by1,2 Fig. 1 shows the time-resolved fluence for a relatively small source–detector separation of compared to . 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 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 LT7 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: 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 -point DFT, we evaluate Eq. (7) for (!). In Fig. 2, the solid line corresponds to the exact time-domain solution [Eq. (3)], whereas the filled dots denote the reconstruction of from its image with . 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 .
We now consider the case that 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 is needed for many time values in . In this case, we use a fixed integration path in connection with the following time-independent parameters74,910.11.12.–13 dealing with the solution of the DE in layered media. The corresponding solution to be considered here has the form 4, and 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.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 . 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 ( grid) can be evaluated in 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 by14Figure 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 . 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 .
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 . 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 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
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.