JBO Letters

Application of the Laplace transform in time-domain optical spectroscopy and imaging

[+] Author Affiliations
André Liemert, Alwin Kienle

Institut für Lasertechnologien in der Medizin und Meßtechnik an der Universität Ulm, Helmholtzstr. 12, D-89081 Ulm, Germany

J. Biomed. Opt. 20(11), 110502 (Nov 18, 2015). doi:10.1117/1.JBO.20.11.110502
History: Received August 11, 2015; Accepted October 20, 2015
Text Size: A A A

Open Access Open Access

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,13 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 Display Formula

f(t)=12πF(ω)eiωtdω,ωR.(1)
However, a serious problem when using the FT [Eq. (1)] arises when F(ω) exhibits a slow algebraic decay for increasing angular modulation frequency ω. 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 Display Formula
G(r,ω)=exp(κr)4πDr,κ=μaD+iωDc,(2)
which is the known infinite-space Green’s function of the frequency-domain DE 2G(r,ω)κ2G(r,ω)=δ(r)/D.6 In this context, we have μa as the absorption coefficient and D=1/[3(μa+μs)] as the diffusion coefficient with μs being the reduced scattering coefficient. Furthermore, ω 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 by1,2Display Formula
G(r,t)=c(4πDct)3/2exp(μactr24Dct),t>0.(3)
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(r,t) from its image G(r,ω) becomes a challenging task. For illustration purposes, Fig. 1 shows the time-resolved fluence for a relatively small source–detector separation of r=0.1mm compared to 1/μs. 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=104 discrete frequencies.

Graphic Jump Location
Fig. 1
F1 :

Reconstruction of the time-resolved fluence via numerical inversion of the Fourier transform. The optical properties of the infinite medium are assumed to be μa=0.01mm1, μs=1mm1, and c=2.99×108/1.4ms1. The smooth curve denotes the exact time-domain solution [Eq. (3)], whereas the oscillating curve corresponds to the numerical inversion of Eq. (2) by means of the discrete FT.

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 Display Formula

f(t)=12πiBF(s)estds,sC,(4)
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ω 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 . The key point in this context is that the kernel function exp(st) becomes a damped wave due to Re(s)<0. Very recently, different contours have been developed and analyzed concerning the numerical inversion of the LT, such as the parabola and hyperbola7 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: Display Formula
s(θ)=μ+iμsinh(θ+iϕ),<θ<,(5)
where s(θ)=iμcosh(θ+iϕ). Inserting Eq. (5) into Eq. (4) leads to the modified inverse Laplace integral Display Formula
f(t)=1πIm(0F[s(θ)]es(θ)ts(θ)dθ),(6)
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: Display Formula
f(t)=hπIm[k=0N1F(sk)exp(skt)sk],(7)
where sk=s(θk) for θ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 Re(s)>0. The remaining curve parameters μ and ϕ 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 μ=4.492075287N/t, ϕ=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 104-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(r,t) from its image G(r,s)=exp(κr)/(4πDr) with κ=μ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=106mm.

Graphic Jump Location
Fig. 2
F2 :

Reconstruction of the time-resolved fluence by means of the inverse Laplace transform. The optical properties of the infinite medium are assumed to be μa=0.01mm1, μs=1mm1, and c=2.99×108/1.4ms1. The solid line corresponds to the exact time-domain solution [Eq. (3)], whereas the filled dots denote the reconstructed values by means of the series [Eq. (7)].

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(t1,t2). In this case, we use a fixed integration path in connection with the following time-independent parameters7Display Formula

μ=4πϕπ2A(ϕ)Nt2,A(ϕ)=arcosh[(π2ϕ)Λ+4ϕπ(4ϕπ)sinϕ],(8)
where Λ=t2/t1. The uniform node spacing becomes h=A(ϕ)/N, and the value for the angle ϕ=1.09 has been recovered by us as a result of numerical experiments. For the next comparison, we consider the time-resolved reflectance R(ρ,t) from a two-layered medium predicted by means of the DE. The first layer has a thickness of L=2mm, whereas the second layer is assumed to be infinitely thick. The refractive indices of the first and second layers are denoted by n1 and n2, respectively. The surrounding nonscattering medium in z<0 has a refractive index of n0=1.0. We consider the DE under the extrapolated boundary condition. In this context, we refer to the papers4,913 dealing with the solution of the DE in layered media. The corresponding solution to be considered here has the form Display Formula
R(ρ,t)=12π012πi[BR(q,s)estds]J0(qρ)qdq,(9)
where R(q,s) denotes the reflectance in the Hankel–Laplace space, which can be found in Ref. 4, and J0(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.10Figure 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×600 grid) can be evaluated in 0.2s 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.

Graphic Jump Location
Fig. 3
F3 :

Time-resolved reflectance from a two-layered medium with optical properties μa1=0.005mm1, μs1=1.2mm1, μa2=0.04mm1, μs2=0.8mm1, n1=1.4, and n2=1.6. The refractive index of the surrounding medium is set to n0=1.0, and the thickness of the first layer is given by L=2mm. The solid line is the reflectance obtained from Eq. (9), whereas the symbols correspond to the real-space Green’s function of Tualle et al.10

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 by14Display Formula

F(r,s)=μax4πDmDx1κx2κm2(eκmrreκxrr),(10)
where κx=μax/Dx+s/(Dxc) and κm=μam/Dm+s/(Dmc). 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: Display Formula
F(r,t)=c8πrμaxDmDxexp(μaxDmμamDxDmDxct)×[exp(λr)erfc(12rDmctλDmct)+exp(+λr)erfc(12rDmct+λDmct)exp(λr)erfc(12rDxctλDxct)exp(+λr)erfc(12rDxct+λDxct)],(11)
where λ=(μamμax)/(DmDx) and erfc(z) are the complementary error function. Note that, depending on the optical properties, λ 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 1014.

Graphic Jump Location
Fig. 4
F4 :

Time-resolved fluorescence caused by an isotropic point source embedded in an infinite medium. The optical properties are assumed to be μax=0.02mm1, μsx=1.5mm1, μam=0.01mm1, μsm=1.0mm1, and c=2.99×108/1.4ms1. The solid line corresponds to Eq. (11), whereas the symbols denote the reconstructed fluorescence by means of the series [Eq. (7)].

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 ρ=10mm. 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 ρ=2mm 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.

Patterson  M. S., , Chance  B., and Wilson  C., “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt.. 28, , 2331 –2336 (1989). 0003-6935 CrossRef
Martelli  F.  et al., Light Propagation Through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software. ,  SPIE Press ,  Bellingham, WA  (2010).
Arridge  S. R., , Cope  M., and Delpy  D.T., “The theoretical basis for the determination of optical pathlength in tissue: temporal and frequency analysis,” Phys. Med. Biol.. 37, , 1531 –1560 (1992). 0031-9155 CrossRef
Liemert  A., and Kienle  A., “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt.. 15, , 025002  (2010). 1083-3668 CrossRef
Liemert  A., and Kienle  A., “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep.. 3, , 2018  (2013). 2045-2322 CrossRef
Haskell  R. C.  et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A. 11, , 2727 –2741 (1994). 0740-3232 CrossRef
Weideman  J., and Trefethen  L., “Parabolic and hyperbolic contours for computing the Bromwich integral,” Math. Comput.. 76, , 1341 –1356 (2007). 0025-5718 CrossRef
Dingfelder  B., and Weideman  J. A. C., “An improved Talbot method for numerical Laplace transform inversion,” Numer. Algorithms. 68, , 167 –183 (2013). 1017-1398 CrossRef
Kienle  A.  et al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt.. 37, , 779 –791 (1998). 0003-6935 CrossRef
Tualle  J. M.  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). 0740-3232 CrossRef
Tualle  J. M.  et al., “Asymptotic behavior and inverse problem in layered scattering media,” J. Opt. Soc. Am. A. 21, , 24 –34 (2004). 0740-3232 CrossRef
Martelli  F.  et al., “Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfunction method,” Phys. Rev. E. 67, , 056623  (2003).CrossRef
Martelli  F.  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). 0031-9155 CrossRef
Farrell  T. J., , Patterson  M. S., “Diffusion modelling of fluorescence in tissue,” in Handbook of Biomedical Fluorescence. , , Mycek  M. A., and Pogue  B. P., Eds.,  CRC Press ,  Boca Raton, FL  (2003).
Pifferi  A.  et al., “Time-resolved diffuse reflectance using small source–detector separation and fast single-photon gating,” Phys. Rev. Lett.. 100, , 138101  (2008). 0031-9007 CrossRef
Patterson  M. S., and Pogue  B. W., “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt.. 33, , 1963 –1974 (1994). 0003-6935 CrossRef
© 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.

Citation

André Liemert and Alwin Kienle
"Application of the Laplace transform in time-domain optical spectroscopy and imaging", J. Biomed. Opt. 20(11), 110502 (Nov 18, 2015). ; http://dx.doi.org/10.1117/1.JBO.20.11.110502


Figures

Graphic Jump Location
Fig. 1
F1 :

Reconstruction of the time-resolved fluence via numerical inversion of the Fourier transform. The optical properties of the infinite medium are assumed to be μa=0.01mm1, μs=1mm1, and c=2.99×108/1.4ms1. The smooth curve denotes the exact time-domain solution [Eq. (3)], whereas the oscillating curve corresponds to the numerical inversion of Eq. (2) by means of the discrete FT.

Graphic Jump Location
Fig. 2
F2 :

Reconstruction of the time-resolved fluence by means of the inverse Laplace transform. The optical properties of the infinite medium are assumed to be μa=0.01mm1, μs=1mm1, and c=2.99×108/1.4ms1. The solid line corresponds to the exact time-domain solution [Eq. (3)], whereas the filled dots denote the reconstructed values by means of the series [Eq. (7)].

Graphic Jump Location
Fig. 3
F3 :

Time-resolved reflectance from a two-layered medium with optical properties μa1=0.005mm1, μs1=1.2mm1, μa2=0.04mm1, μs2=0.8mm1, n1=1.4, and n2=1.6. The refractive index of the surrounding medium is set to n0=1.0, and the thickness of the first layer is given by L=2mm. The solid line is the reflectance obtained from Eq. (9), whereas the symbols correspond to the real-space Green’s function of Tualle et al.10

Graphic Jump Location
Fig. 4
F4 :

Time-resolved fluorescence caused by an isotropic point source embedded in an infinite medium. The optical properties are assumed to be μax=0.02mm1, μsx=1.5mm1, μam=0.01mm1, μsm=1.0mm1, and c=2.99×108/1.4ms1. The solid line corresponds to Eq. (11), whereas the symbols denote the reconstructed fluorescence by means of the series [Eq. (7)].

Tables

References

Patterson  M. S., , Chance  B., and Wilson  C., “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt.. 28, , 2331 –2336 (1989). 0003-6935 CrossRef
Martelli  F.  et al., Light Propagation Through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software. ,  SPIE Press ,  Bellingham, WA  (2010).
Arridge  S. R., , Cope  M., and Delpy  D.T., “The theoretical basis for the determination of optical pathlength in tissue: temporal and frequency analysis,” Phys. Med. Biol.. 37, , 1531 –1560 (1992). 0031-9155 CrossRef
Liemert  A., and Kienle  A., “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt.. 15, , 025002  (2010). 1083-3668 CrossRef
Liemert  A., and Kienle  A., “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep.. 3, , 2018  (2013). 2045-2322 CrossRef
Haskell  R. C.  et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A. 11, , 2727 –2741 (1994). 0740-3232 CrossRef
Weideman  J., and Trefethen  L., “Parabolic and hyperbolic contours for computing the Bromwich integral,” Math. Comput.. 76, , 1341 –1356 (2007). 0025-5718 CrossRef
Dingfelder  B., and Weideman  J. A. C., “An improved Talbot method for numerical Laplace transform inversion,” Numer. Algorithms. 68, , 167 –183 (2013). 1017-1398 CrossRef
Kienle  A.  et al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt.. 37, , 779 –791 (1998). 0003-6935 CrossRef
Tualle  J. M.  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). 0740-3232 CrossRef
Tualle  J. M.  et al., “Asymptotic behavior and inverse problem in layered scattering media,” J. Opt. Soc. Am. A. 21, , 24 –34 (2004). 0740-3232 CrossRef
Martelli  F.  et al., “Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfunction method,” Phys. Rev. E. 67, , 056623  (2003).CrossRef
Martelli  F.  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). 0031-9155 CrossRef
Farrell  T. J., , Patterson  M. S., “Diffusion modelling of fluorescence in tissue,” in Handbook of Biomedical Fluorescence. , , Mycek  M. A., and Pogue  B. P., Eds.,  CRC Press ,  Boca Raton, FL  (2003).
Pifferi  A.  et al., “Time-resolved diffuse reflectance using small source–detector separation and fast single-photon gating,” Phys. Rev. Lett.. 100, , 138101  (2008). 0031-9007 CrossRef
Patterson  M. S., and Pogue  B. W., “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt.. 33, , 1963 –1974 (1994). 0003-6935 CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

PubMed Articles
Advertisement


 

  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.