Many different parts of the human body exhibit a layered structure, e.g., the head or the extremities. To calculate the light propagation in these layered tissues, which is important for the application of light in therapeutical and diagnostic medicine, fast and accurate forward solutions are required.
Usually, it is assumed that the radiative transfer theory correctly describes light propagation in biological media.1, 2 Most often, the Monte Carlo method is used to solve the radiative transfer equation numerically. To obtain faster solutions describing the light propagation in biological tissue, the diffusion equation, which is an approximation to the radiative transfer theory, is applied.1
Several solutions of the layered diffusion equation have been reported in literature. Dayan applied a Fourier formalism to solve the two-layered diffusion equation. However, to be able to obtain analytical solutions several approximations were introduced in the inverse Fourier calculations.3 We derived exact solutions for the two-layered geometry in the steady-state, frequency, and time domains4, 5, 6 by numerically performing the inverse Fourier transform and introducing appropriate boundary conditions. Subsequently, different groups have used these solutions for further investigations.7, 8 Additionally, the solutions have been extended, for example, to a hybrid Monte Carlo/diffusion scheme,9 to a tilted interface between the layers,10 or to diffusing wave spectroscopy.11
Supplementally, different approaches to solve the layered diffusion equation have been made. Tualle applied the method of images to derive solutions of the two-layered geometry,12 whereas Martelli used the eigenfunction method together with the separation of variables technique to solve the time-dependent diffusion equation for two and three layers.13, 14, 15 Barnett solved the -layered diffusion equation via a finite-difference method.16 Recently, an approach based on the solution of a single layer to obtain the solutions for many layers using convolution was proposed.17
Building on our earlier work we solved the diffusion equation for -layered turbid media in the steady-state, frequency, and time domains and compared the results with Monte Carlo simulations. The solutions were derived for a semi-infinite layered and for a finite layered geometry, i.e., the bottom layer has an infinite and finite thickness, respectively. Recently, independently of us, solutions of the -layered diffusion equation have been published based on our earlier work for the special case of semi-infinite turbid media.18, 19 Whereas the extension of the Fourier formalism to layers is straightforward, the decisive point is the numerically correct and efficient 2-D inverse Fourier transform. To this end, we implemented four different methods to obtain the solutions in real space. Additionally, approximate solutions for the fluence rate in Fourier space were derived to avoid numerical errors due to the inverse Fourier transform.
In this paper, the light diffusion in turbid media is considered in the steady-state domain, whereas in the accompanying paper20 the light diffusion in the frequency and time domains is analyzed. Turbid media consisting of up to 20 layers having matched and mismatched refractive indices are studied, showing good agreement between the solutions of the diffusion theory and Monte Carlo simulations. The four methods for calculating the inverse Fourier transform are analyzed to optimize the performance of the solutions. Explicit solutions are given in the appendix for turbid media having two, three, and four layers for both the semi-infinite and the finite bottom layer. In addition, an executable program for calculation of the reflectance from an -layered turbid medium is provided on the Internet.21
In this section, the steady-state diffusion equation for an -layered turbid medium having laterally infinite extensions is derived using the Fourier transform formalism. A pencil light beam is incident in the direction perpendicular to the interface between the first layer and the nonscattering surrounding (see Fig. 1 ).
The origin of the coordinate system is the point of light impact on the turbid medium. As in our earlier paper, it is assumed that the light beam can be characterized by an isotropic point source at a depth in the first layer, , where and are the absorption and reduced scattering coefficients of the first layer, respectively.4 For this geometry, the diffusion equation becomesis the thickness of layer , is the fluence rate, is the diffusion coefficient, and and are the reduced scattering and absorption coefficients of layer , respectively.
We solved this system of equations using the Fourier transform approach and extrapolated boundary conditions. First, the equations are two-dimensionally Fourier transformed from the real space into the Fourier space with respect to the and coordinates. The boundary conditions depending only on the coordinate are applied in the Fourier space. The resulting system of equations is solved, yielding the solutions in the Fourier space. Finally, these solutions are two-dimensionally inverse Fourier transformed, delivering the solutions in the real space.
We start by applying the 2-D Fourier transform:1, 2, and obtain and . The fluence rate in the Fourier space in layer is denoted by . The following boundary conditions are used in the Fourier space for an -layered turbid medium having a finite thick ’th layer: and are the fluence rate in the first layer in the Fourier space above the source and below the delta source, respectively. The refractive index of layer is denoted by , and is the total thickness of the turbid medium. Equations 7, 8 describe the fluence rate and its derivative at the position of the source, and Eqs. 9, 10 at the interfaces of the different turbid layers. In Eqs. 6, 11, and are the positions of the extrapolated boundaries above the first and below the ’th layer, respectively (see Fig. 1), where the fluence rate is set to zero.22 These quantities are obtained using represents the fraction of photons that is internally diffusely reflected at the boundary to the nonscattering surrounding at the top or at the bottom layer. We calculated by solving the integrals derived by Haskell 22
In the case of an -layered turbid medium having an infinitely thick ’th layer, Eq. 11 becomes
We solved the resulting equation system for layers having a finite and semi-infinitely thick ’th layer by applying Cramer’s rule. The general solution for layers was obtained by induction using the solutions for two, three, and four layers. In the first case, we obtain for the fluence rate in the first layer (above the delta source)’th layer is and are obtained for an -layered turbid medium by iteratively using , , i.e., until and are obtained. This means that for , Eq. 17 can be directly applied without using Eq. 18. We give the fluence rate of the first and ’th layers, because we are interested in the reflectance from the first and the transmittance from the ’th layer, which are calculated from the fluence rates of these layers.
In the second case, the semi-infinite -layered turbid medium , the fluence rate is again calculated with Eq. 14, but the start terms for computing and are now
The recursion formula is the same as for the finite -layered turbid medium, see Eq. 18. For the semi-infinite -layered turbid medium we give only the fluence rate in the first layer, because the transmittance is zero. With these equations, the analytical formula of the fluence rate in the Fourier space is obtained for -layered turbid media for . In the appendix, we also give the solutions for two-layered turbid media for the finite and infinite cases. In addition, the explicit formulas for three and four layers are also presented for the finite and infinite thick bottom layer.
Next, the 2-D inverse Fourier transform of must be calculated to obtain the fluence rate in real space . We performed these transforms numerically as no general analytical solution was found up to now. For the accuracy of the calculation of the fluence rate in real space, this step is decisive. To minimize the numerical errors, we first used approximate solutions for large values (see Sec. 5.2). Second, we programmed and compared four different methods to accomplish the 2-D inverse Fourier transform:
In method A, Eq. 20 was solved numerically by evaluating it at discrete and values using a DFT (discrete Fourier transform). The FFT (fast Fourier transform) algorithm was applied to accelerate the calculations.23 Due to the rotational symmetry, it is sufficient to calculate the fluence rate in one direction.
In method B, we expressed in a 2-D Fourier series:
By accurate sampling of in the Fourier space and by regarding only the central solution of the resulting periodic function, is obtained without significant aliasing. Equation 21 can be simplified by using the rotational symmetry of the considered geometry:
Again, it is sufficient to evaluate in one direction (e.g., ). By simplifying the evaluation of the sums and assuming the same sampling rate in both directions we finally get
To obtain similar accuracy as for the other methods, we typically used 400 terms in the sums and a sampling rate in the range of . Note that the implementation of method B is the easiest of the four methods.
In the other two methods the rotational symmetry is used to obtain from Eq. 20 the 1-D inverse Hankel transform:denotes the Bessel function of first kind and zero order.
In method C, the integral is solved numerically by applying Gauss integration using typically 2 times 480 points.23 For the upper limit of the integral, we used values in the range of . In our original paper concerning the solution of the two-layered turbid medium, this method was used to calculate the inverse Fourier transform.4
In method D, Eq. 24 is transformed into a correlation integral, which is then represented by a discrete correlation.24 This can be evaluated efficiently with an FFT algorithm, which delivers the fluence rate at many (e.g., 512) distances. However, the distances are exponentially distributed. Alternatively, the discrete correlation can be evaluated at one distance without using the FFT algorithm.
By applying one or several of the four methods, the fluence rate is obtained. Then, the spatially resolved reflectance from the first layer is computed usingis the Fresnel reflection coefficient for a photon having an incident angle relative to the normal to the boundary. The use of the combined fluence rate and flux term in Eq. 25 is important to match better Monte Carlo simulations that use the Henyey-Greenstein function as a phase function with anisotropy factors typically found for biological media25 . For the finite, thick -layered tissue the transmittance is calculated with
To determine the transmittance, it is sufficient to use only the flux term to approximate the Monte Carlo simulations.
The solutions of the diffusion equations are compared with Monte Carlo simulations. The Monte Carlo method has been described in detail in the literature.26 The salient features of our code are as follows. A pencil beam is incident in direction perpendicular onto the turbid medium at the origin of the coordinate system. The Henyey-Greenstein function with an anisotropy factor of is applied as phase function. At the boundaries, the Fresnel reflection and transmission coefficients are used to consider the different refractive indices. Both the Monte Carlo simulations and the solutions of the diffusion equations were programmed in Delphi (Pascal). Extended variables (19 to 20 significant digits; ) are used to minimize numerical errors.
In this section, we first show the results for the solutions of the diffusion equation obtained with the different methods for computing the inverse Fourier transform. Second, the solutions of the diffusion equations are compared to Monte Carlo simulations.
Comparison of Different Methods for Calculation of the Inverse Fourier Transform
The spatially resolved reflectance from -layered turbid media was calculated using the four methods for computing the inverse Fourier transform, which are explained in the theory section. In general, we found that all four methods agreed excellently. Of course, the accuracy depends on the numerical efforts that have been applied. We chose the number of points used in the series or integral evaluations in the four methods such that the accuracy for calculation of was similar. As an example, we show the spatially resolved reflectance from a two-layered turbid medium calculated with method C (solid curve) and method D (circles) in Fig. 2 .
No differences can be seen in the figure between the results obtained from the two methods. To distinguish both curves we calculated the relative difference between them, which is shown in Fig. 3 .
The differences are for distances from the incident source. These are typical differences that we found by comparing all four methods. In most cases, the differences are smaller for small distances compared to large distances from the incident source. Note that a high accuracy for the calculation of is especially necessary for reconstruction problems, i.e., the determination of the optical coefficients with, e.g., nonlinear regression routines.
As the numerical accuracy of the four models is similar, it is interesting to compare their execution times. For a layered turbid medium having typical optical properties found in biological media we got following typical calculation times for the reflectance or transmittance in real space using a state-of-the-art computer. With methods A, B, C, and D computation times of (for distances; by using the FFT one obtains many distances in one calculation), (for one distance), (for one distance), and (for one distance) were needed, respectively. When the FFT was used to evaluate the correlation in method D (see Sec. 2) even distances could be evaluated in . The disadvantage of using the FFT in methods A and D is that the fluence rate cannot be obtained for arbitrarily chosen distances. In method A, the distances are equidistant, and in method C, exponentially distributed.
Comparison of the Solutions of the Diffusion Equation with Monte Carlo Simulations
Figure 4 shows the influence of different refractive indices of the first layer on the spatially resolved reflectance. A two-layered semi-infinite turbid medium was considered having (1) a refractive index of 1.4 for both layers [green (lighter) curve and circles] and (2) a refractive index of 1.8 for the first and 1.4 for the second layer (black curve and circles). For the nonscattering surrounding, a refractive index of 1.0 was assumed. A good agreement between the Monte Carlo simulations (circles) and diffusion theory (solid curves) can be observed.
Figure 5 shows a comparison of diffusion theory and Monte Carlo simulations for reflectance from a turbid medium having 20 layers. The optical properties of the layers were randomly chosen in the range between and , but were the same for the solution of the diffusion equation and Monte Carlo simulation. The thickness of each layer was . Thus, the thickness of the whole turbid medium was . The refractive index of each layer was 1.4, that of the nonscattering surrounding media was 1.0. Again, a good agreement between the two solutions can be stated except for small distances where the diffusion approximation is not valid.
Figure 6 gives the relative difference of both curves to better visualize the performance of the diffusion solution compared to the Monte Carlo simulation.
At small distances the differences are large, for intermediate distances the differences are about 5%, and for large distances the differences are smaller than 1% when the noise in the Monte Carlo simulations is accounted for. In general, the differences are even smaller for turbid media having larger thicknesses, because the diffusion approximation is better fulfilled.
Finally, Fig. 7 shows the transmittance from a four-layered turbid medium having a finite fourth layer. The solution of the diffusion equation (solid line) calculated with Eqs. 26, 24, 15, 16, 17, 18 was compared to Monte Carlo simulations (circles). The total thickness of the four layers was .
Figure 7 shows good agreement between the two calculations for all distances because in transmission, the diffusion approximation is valid also for small values due to the large value of the product , i.e., the transmitted photons have experienced a high number of scattering interactions in the turbid medium. The relative differences between the Monte Carlo simulation and the solutions from the diffusion equation are shown in Fig. 8 . The differences are largely caused by the statistics of the Monte Carlo simulations. The systematic differences due to the diffusion approximations are smaller than about 1% for all distances.
The solutions of the -layered diffusion equation were calculated in the steady-state domain for turbid media having a finite and an infinitely thick ’th layer. Analytical solutions were derived for the fluence rate in the Fourier space. The decisive point for obtaining precise values in real space is the accurate evaluation of the 2-D inverse Fourier transform. To this end, we first derived approximations of the analytical solutions in the Fourier space for large values (see Sec. 5.2). Second, we investigated four different methods for calculation of the inverse Fourier transform. Subsequently, we could estimate the errors of the obtained values for the reflectance and transmittance from the -layered turbid media by comparing the results of the four methods. We found relative differences that were typically smaller than . These values are a compromise between accuracy and speed of the solutions. Although we emphasized the accuracy in the calculations, we could achieve short calculation times in the range of even for a large number of reflectance or transmittance values at different distances from the source.
We compared the solutions of the -layered diffusion equation having an infinite and finite thick ’th layer with results from Monte Carlo simulations. The relative differences between both methods were in the same range as those found for homogeneous turbid media having a semi-infinite or slab geometry.25
Note that the derived solutions can be easily altered to obtain solutions of the heat conduction equation in an -layered material. Also, different boundary conditions, for example, the partial boundary condition, can be implemented with our method.22
In future work, the derived solutions will be used to study the inverse problem. For example, the determination of the optical properties of the brain using a layered model for the head will be investigated. We study whether a priori information of, for example, the optical properties of some of the involved tissue layers, is required to derive accurately the optical properties of the brain.
We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).
Appendix A. Fluence Rate in Fourier Space
In this subsection we present the analytical solutions of the fluence rate in the Fourier space for turbid media having two, three, and four layers . For turbid media having a finite thick ’th layer, the fluence rate in the first and ’th layer is given, whereas only the fluence rate in the first layer is shown for the case of an infinite thick ’th layer.
First, for the fluence rate in the first layer (above the source) of the turbid medium having both a finite and infinite thick bottom layer we obtain
The quantities and are for a finite two-layered turbid medium:
Second, for the fluence rate in the ’th layer we obtain for the finite two-layered turbid medium:is given by Eq. 29. For the finite three-layered turbid medium: is given by Eq. 33. For the finite four-layered turbid medium: is given by Eq. 37.
Appendix B. Approximate Solutions of the Fluence Rate in Fourier Space
To avoid numerical errors in the computation of the inverse Fourier transform, we derived approximate solutions for the fluence rate of -layered turbid media in the Fourier space. Note that in an earlier paper we derived corresponding results for a two-layered medium.6 For we gotand are iteratively calculated using Eq. 18. The starting values are obtained from Eq. 17 for a -layered turbid medium having a finite thick ’th layer and from Eq. 19 for a turbid medium having an infinitely thick ’th layer. This approximation was typically used for values larger than .