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\pi i\u222bBF(s)est\u2009ds,s\u2208C,$(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\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 $\u2212\u221e$. 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 hyperbola^{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: Display Formula$s(\theta )=\mu +i\mu \u2009sinh(\theta +i\varphi ),\u2212\u221e<\theta <\u221e,$(5)
where $s\u2032(\theta )=i\mu \u2009cosh(\theta +i\varphi )$. Inserting Eq. (5) into Eq. (4) leads to the modified inverse Laplace integral Display Formula$f(t)=1\pi Im(\u222b0\u221eF[s(\theta )]es(\theta )ts\u2032(\theta )d\theta ),$(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\pi Im[\u2211k=0N\u22121F(sk)exp(skt)sk\u2032],$(7)
where $sk=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 $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 $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(\u2212\kappa r)/(4\pi Dr)$ with $\kappa =\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\u22126\u2009\u2009mm$.