In the first of our two papers dealing with light diffusion in -layered turbid media we studied the light propagation in the steady-state domain.1 In this second part, we concentrate on solutions of the time-resolved diffusion equation in both the time and frequency domains.
In general, measurements in the steady-state domain, for which the considered turbid medium is illuminated by a continuously emitting source and the remitted or transmitted light is detected steadily, are relatively inexpensive and simple to handle.2 Contrarily, time-resolved measurements are more complex and expensive, but they contain more information on the investigated tissues.3 Time-resolved measurements are performed in the time and frequency domains. In the time domain, a short laser pulse of the order of picoseconds is incident onto the turbid medium and the remitted or transmitted light is detected in a time-resolved manner.4, 5 In the frequency domain, the incident source is intensity modulated in the range of and the demodulation and phase shift of the remitted or transmitted light is detected.6
In the literature, experiments in all three measurement regimes were reported assuming that the investigated turbid media can be considered homogeneous or heterogeneous.7 In this paper, we consider layered geometries, a special case of heterogeneous turbid media. Many parts of the human body can be approximated as layered turbid media, e.g., the head, the arm, and the leg. For the mathematical description of light propagation in such large tissue volumes the diffusion equation is most often applied. The literature on the solution of the layered diffusion equation is reviewed in the accompanying paper.1
In this paper, we present the derived solutions of the time-resolved reflectance and transmittance from an -layered turbid medium having a finite thick ’th layer, and the time-resolved reflectance from an -layered turbid medium having an infinitely thick ’th layer. The solutions in the frequency domain are obtained by using the Fourier transform technique, as reported in the accompanying paper.1 Different methods are applied for calculation of the inverse Fourier transform to compare their performance and accuracy. The solutions in the time domain are obtained by applying the FFT (fast Fourier transform) to the solutions in the frequency domain. Additionally, we derive an analytical formula for the reflectance from a two-layered turbid medium having an infinitely thick second layer for the case of equal reduced scattering coefficients and refractive indices. Finally, the obtained solutions are compared to Monte Carlo simulations.
In this section, the derivation of the solutions of the -layered diffusion equation in the frequency and time domains is presented. We start with the diffusion equation in the time domain:is the fluence rate in real space, is the absorption coefficient, the diffusion coefficient, and is the reduced scattering coefficient of the turbid medium. By solving Eq. 1 for a source in time and space Green’s function of the problem is obtained with
For an arbitrary source in time and space , the fluence rate can be calculated with the knowledge of Green’s function by applying
We propose two different methods to calculate . In the first method, Eq. 2 is solved using the Laplace and Fourier transform techniques. We first apply the Laplace transform with respect to the coordinate and the 2-D Fourier transform relative to the and coordinates to obtain an ordinary differential equation depending on the coordinate in the inverse Laplace and Fourier space. Then, the boundary conditions are applied relative to the coordinate. The next step is to solve the resulting equation system. Finally, the inverse 2-D Fourier and Laplace transforms are employed to obtain the time-resolved fluence rate in real space. Note that the Laplace transform with respect to the coordinate can, in principle, be substituted by an additional Fourier transform. Thus in this case also, the inverse Laplace transform must be replaced by an inverse Fourier transform.
The challenge of this method is that it is difficult to find analytical expressions for the inverse Laplace and Fourier transforms.8 However, we succeed in deriving a fast analytical solution for a two-layered turbid medium in the time domain having an infinitely thick bottom layer for the same reduced scattering coefficient and refractive index of both layers. The derivation of this solution together with that for an infinite geometry is presented in the appendix.
In this section, we concentrate on the second method. First, solutions of the fluence rate in the frequency domain are derived numerically using the Fourier formalism that was explained in the accompanying paper.1 From these results the solutions in the time domain are derived using an FFT. For measurements in the frequency domain, the intensity of the incident source is modulated sinusoidally in time:and , where is the frequency of the oscillating intensity. In Eq. 4, we omit the steady-state portion and the conjugate-complex term for the description of the source. Inserting Eq. 4 into Eq. 3 yields the fluence rate
Subsequently, the 3-D integration in the spatial domain is performed. By assuming a source in space, , we obtainis Green’s function in the frequency domain. We can derive by solving the following differential equation, which is obtained from the Fourier transform relative to the coordinate of Eq. 2: is the Laplace operator, and is the velocity of light in the considered medium. Thus, Green’s function in the time domain can be derived by calculating Green’s function in the frequency domain, i.e., solving Eq. 7, and by applying the inverse Fourier transform:
For a vanishing intensity oscillation we get from Eq. 7 the steady-state diffusion equation, which was solved for -layered turbid media in the accompanying paper.1 Thus, to obtain the solution of the -layered diffusion equation in the frequency domain assuming a -source at and :-layered steady-state diffusion equation derived in the accompanying study1 can be used simply by replacing with , where is the velocity of light in layer . In addition, the four methods1 presented to perform the numerical 2-D inverse Fourier transform as well as the approximate solutions in the steady-state domain can also be applied in the frequency domain.
The reflectance from the -layered turbid medium in the frequency domain is computed from the fluence rate byis1 the Fresnel reflection coefficient and . For large distances from the source it is sufficient to use the Fick’s law for calculation of the reflectance9
In the frequency domain, the intensity-modulated source causes that the reflectance measured at the detector position is also modulated by the same frequency but the oscillation is phase shifted and the modulation is reduced. The quantities of interest in the frequency domain are the phase angle between the source and the detected signal and the modulation :
In summary, to obtain the solutions for the fluence rate in the frequency domain, the solutions in the steady-state domain presented in the accompanying paper1 are used by replacing with . The solutions in the time domain are calculated by computing the solution in the frequency domain for many frequencies (e.g. 512) and by using10 an FFT. Note that because is a causal function in time, it can be shown using the Hilbert transform that it is sufficient to calculate only the real part or the imaginary part of . Then, is obtained by an FFT of two times the real part or two times the imaginary part of , accelerating the calculations by a factor of 2.
The solutions of the -layered diffusion equation are compared to Monte Carlo simulations. The Monte Carlo method is within its stochastic nature an exact solution of radiative transfer theory, whereas diffusion theory is an approximation of radiative transfer theory. Thus, the validity of the diffusion equation for the considered geometrical and optical parameters can be tested by comparison with Monte Carlo simulations. In addition to the explanations given in the accompanying paper,1 for the calculation of the Monte Carlo simulations in the time domain we use a variance reduction scheme. The absorption coefficient is set to zero during the simulations for all layers and the length of the photon paths in each layer is scored. Thereby, the time-resolved reflectance can be calculated using the Lambert-Beer law.11 The results in the frequency domain are obtained by Fourier transforming (via an FFT) the Monte Carlo data calculated in the time domain.
We, first, show a comparison of solutions of the diffusion equation obtained with different methods for computing the inverse Fourier transform and calculated with the formulas presented in the appendix. Second, the solutions of the diffusion equation are compared to Monte Carlo simulations.
Comparison of Different Methods for Solving the Diffusion Equation
In the accompanying paper,1 we presented four different methods for calculation of the inverse Fourier transform to obtain the steady-state fluence rate in real space. The accuracy of the computed solutions can be tested with these methods. All methods can be applied also in the frequency domain and, subsequently, for the calculation in the time domain. Figure 1 compares the results of two methods for calculating the 2-D inverse Fourier transform, the numerical integration of the Hankel transform (solid curve, method C in Ref. 1) and the use of the correlation for solving the Hankel transform (circles, method D in Ref. 1). Figure 1 shows the time-resolved reflectance from a two-layered turbid medium having an infinitely thick second layer at a distance .
Both curves are indistinguishable in the figure. Thus, Fig. 2 shows the relative difference between both curves. The differences are smaller than for all time values over an intensity range of more than six orders of magnitude. A high accuracy of the obtained time-resolved fluence rate is important; e.g., for solving the inverse problem, the robust reconstruction of the optical properties of the different layers.
In the appendix, an analytical solution of the two-layered semi-infinite diffusion equation having the same reduced scattering coefficient and refractive index in both layers is derived, see Eq. 41. Figure 3 compares the reflectance calculated with this formula (circles) to that computed with the numerical solution presented in the theory section using method D in Ref. 1 for the inverse Fourier calculations (solid curve). The distance between the incident beam and detection is .
The disagreement between the two curves are small, so that we also plot the relative differences (see Fig. 4 ). The relative differences are again smaller than for the time values shown. Clearly, they can even be reduced by laying more stress on the numerical calculations, which on the other hand increases the computation time. The advantage of the analytical formula is its velocity. When the FFT is used for evaluation [see Eq. 43], the calculation time is in the range of a millisecond.
Comparison of the Solutions of the Diffusion Equation with Monte Carlo Simulations
First, we calculated the reflectance from a four-layered turbid medium having an infinitely thick fourth layer at a distance of from the incident light source in the time domain. Figure 5 shows the comparison of the results obtained by the Monte Carlo method (circles) and by the diffusion solution (solid curve).
The two solutions agree well except for reflectance values at short times where it is known that diffusion theory is not valid. To show the small disagreement more distinctly the relative difference is depicted in Fig. 6 . For times longer than about , the relative differences are mainly due to the statistics of the Monte Carlo simulations. The systematic relative differences between the solution of the diffusion equation and the Monte Carlo simulation are smaller than 1%.
In the frequency domain, we compare the results obtained with both models for the same four-layered turbid medium that was considered in the time domain. The geometrical and optical properties of the turbid medium are given in the caption of Fig. 5. Figures 7 and 8 show the phase and the modulation for the reflectance at a distance of calculated by the Monte Carlo simulations (circles) and the solution of the diffusion equation (solid curve). The phase and modulation are drawn versus the modulation frequency . These data were directly calculated in the frequency domain in the case of the solution of the diffusion equation, whereas in the case of the Monte Carlo method simulations were performed in the time domain and the FFT was applied to obtain the phase and modulation for different modulation frequencies.
In general, good agreement between the two methods is obtained. The differences are larger for the phase, mainly caused by errors due to the diffusion approximations at small times (compare Fig. 5). The advantage in the time domain is that the photons that are remitted at early times and that do not fulfill the diffusion approximation can be separated from the photons at longer times that are well described by the diffusion approximation. Contrarily, in the frequency domain this separation is not possible.
Solutions of the -layered diffusion equation were derived for turbid media having a finite and infinite thick ’th layer and mismatched refractive indices. The final results in the frequency domain were obtained by a numerical Fourier transform of the analytical solutions of the fluence rate in the Fourier space. As explained in detail in our accompanying paper1 we employed four different methods to be able to validate the accuracy of the calculations. The results in the time domain were calculated by an FFT using the solutions in the frequency domain at many frequencies.
We also derived an analytical solution for the special case of a two-layered turbid medium having the same reduced scattering and refractive index in both layers, of which the bottom one is infinitely thick. By applying the FFT fast solutions in the range of could be achieved with a state-of-the-art computer using one processor.
The accuracy of the numerical and analytical solutions was found to be smaller than for the evaluated times and geometrical and optical parameters. This number is a compromise between the accuracy and velocity of the calculations.
The derived solutions were compared to Monte Carlo simulations in the time and frequency domains and similar agreement as for homogeneous geometries was found. Solutions of the diffusion equation in the frequency domain showed larger differences compared to those in the time domain for longer times due to the diffusion approximation as discussed in the results section.
In the future we plan to apply the derived solutions for investigating the inverse problem, the reconstruction of the optical properties for some or all involved layers. Finally, note that executable files of a computer code for the solutions of the -layered diffusion equation in the time and steady-state domains are available on the Internet.12
We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).
In this section, we present the derivation of an analytical solution of the diffusion equation for the reflectance from a two-layered turbid medium having the same refractive index and the same reduced scattering coefficient in both layers. To this end, we apply the Laplace and Fourier formalism, which was explained in Sec. 2. As an introduction, we first present the same formalism for deriving the well-known solution in the infinite medium.
Appendix A: Infinite Turbid Medium
The starting point is the time-dependent diffusion equation for a delta source in space and time [see Eq. 2]. By applying the Laplace transform,and is employed: is the fluence rate in the Laplace-Fourier space, and . Next, the boundary conditions are introduced. First, is zero infinitely far away from the source: must be considered similar as in the accompanying paper:1
The ordinary differential Eq. 18 can be solved using
With the preceding four boundary conditions, the unknown constants are determined resulting in the following solutions:
Inserting in Eq. 24 yields
Subsequently, these equations are two-dimensionally inverse Fourier and Laplace transformed to obtain . The order of the three integrations can be chosen arbitrarily. We first apply the inverse Laplace transform and then two times the inverse Fourier transform. By using following formulas for the Laplace transform:, one gets
Finally, the 2-D inverse Fourier transform is employed using the formula (and its corresponding equation with respect to and ):13
Appendix B: Two-Layered Turbid Medium
Similar to the last subsection we started using the Laplace and Fourier transform to obtain an ordinary partial differential equation in the variable . Then the boundary conditions for the two-layered turbid medium were applied. The resulting equation system was solved, as presented in the accompanying paper.1 There, we got for the fluence rate from a two-layered turbid medium having an infinitely thick second layer in the transformed space, assuming that and ,, , and is the velocity of light in both layers.
We make use of the following auxiliary equationand the geometric series in the last step. Inserting Eq. 34 into Eq. 33 yields35 it follows that
Subsequently, we perform the inverse transforms. Employing the Laplace transform formulas14is the Heaviside step function, is the modified Bessel function of the ’th order, and the inverse Fourier transform with respect to the coordinate (and accordingly to the coordinate)
To accelerate the calculations we derive an alternative formula for . We start with Eq. 36, the first part of which is inverse Laplace and Fourier transformed as in the preceding, resulting in the first term of given in Eq. 41. The second part of Eq. 36 is transformed as follows:
Then, Eq. 42 is inverse Fourier transformed with respect to the and coordinates, whereas the inverse transform with respect to the time variable is not analytically evaluated. For the calculation of the inverse Fourier transform, we apply the corresponding shifting rule [see Eq. 29] with . In the inverse Laplace transform with respect to the time variable we used , because the imaginary axis is within the region of convergence. Thus, we obtain43 is identical to the Fourier transform and can be evaluated efficiently with the FFT method, which results in calculation times of the order of .