1.

Introduction

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.¹

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.³ We derived exact solutions for the two-layered geometry in the steady-state, frequency, and time domains^{4, 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,⁹ to a tilted interface between the layers,¹⁰ or to diffusing wave spectroscopy.¹¹

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,¹² 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 $N$ -layered diffusion equation via a finite-difference method.¹⁶ Recently, an approach based on the solution of a single layer to obtain the solutions for many layers using convolution was proposed.¹⁷

Building on our earlier work we solved the diffusion equation for $N$ -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 $N$ -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 $N$ 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 paper²⁰ 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 $N$ -layered turbid medium is provided on the Internet.²¹

2.

Theory

In this section, the steady-state diffusion equation for an $N$ -layered turbid medium having laterally infinite extensions is derived using the Fourier transform formalism. A pencil light beam is incident in the $z$ direction perpendicular to the interface between the first layer and the nonscattering surrounding (see Fig. 1 ).

Fig. 1

Scheme of the $N$ -layered turbid medium including the position of the source and the extrapolated boundaries.

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 $z=z_0=1∕({\mu }_{a1}+{\mu }_{s1}^{\prime })$ in the first layer, $\delta (x,y,z-z_0)$ , where ${\mu }_{a1}$ and ${\mu }_{s1}^{\prime }$ are the absorption and reduced scattering coefficients of the first layer, respectively.⁴ For this geometry, the diffusion equation becomes

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 $x$ and $y$ coordinates. The boundary conditions depending only on the $z$ 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:

In the case of an $N$ -layered turbid medium having an infinitely thick $N$ ’th layer, Eq. 11 becomes

We solved the resulting equation system for $N$ layers having a finite and semi-infinitely thick $N$ ’th layer by applying Cramer’s rule. The general solution for $N$ 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)

The quantities ${\beta }_3$ and ${\gamma }_3$ in Eqs. 14, 16 are obtained by using recursion formulas. By applying ${\beta }_N$ and ${\gamma }_N$ ,

In the second case, the semi-infinite $N$ -layered turbid medium $(L_N\to \infty )$ , the fluence rate ${\Phi }_1(z,s)$ is again calculated with Eq. 14, but the start terms for computing ${\beta }_3$ and ${\gamma }_3$ are now

The recursion formula is the same as for the finite $N$ -layered turbid medium, see Eq. 18. For the semi-infinite $N$ -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 $N$ -layered turbid media for $N⩾3$ . 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 $\Phi (z,s)$ must be calculated to obtain the fluence rate in real space $\Phi (x,y,z)$ . 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 $s$ 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 $s_1$ and $s_2$ values using a DFT (discrete Fourier transform). The FFT (fast Fourier transform) algorithm was applied to accelerate the calculations.²³ Due to the rotational symmetry, it is sufficient to calculate the fluence rate in one direction.

In method B, we expressed ${\Phi }_k(x,y,z)$ in a 2-D Fourier series:

By accurate sampling of ${\Phi }_k$ in the Fourier space and by regarding only the central solution of the resulting periodic function, ${\Phi }_k(x,y,z)$ 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 ${\Phi }_k$ in one direction (e.g., $y=0$ ). By simplifying the evaluation of the sums and assuming the same sampling rate in both directions $(\Delta s_1=\Delta s_2=\Delta s)$ 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 $\Delta s\approx 20{\mu }_s^{\prime }∕400$ . 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:

In method C, the integral is solved numerically by applying Gauss integration using typically 2 times 480 points.²³ For the upper limit of the integral, we used values in the range of $30\times {\mu }_s^{\prime }$ . In our original paper concerning the solution of the two-layered turbid medium, this method was used to calculate the inverse Fourier transform.⁴

In method D, Eq. 24 is transformed into a correlation integral, which is then represented by a discrete correlation.²⁴ 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 ${\Phi }_1(\rho ,z)$ is obtained. Then, the spatially resolved reflectance $R\left(\rho \right)$ from the first layer is computed using

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.²⁶ The salient features of our code are as follows. A pencil beam is incident in $z$ direction perpendicular onto the turbid medium at the origin of the coordinate system. The Henyey-Greenstein function with an anisotropy factor of $g=0.8$ 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; $10\text{bytes}$ ) are used to minimize numerical errors.

3.

Results

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.

3.1.

Comparison of Different Methods for Calculation of the Inverse Fourier Transform

The spatially resolved reflectance from $N$ -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 $\Phi (\rho ,z)$ 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 .

Fig. 2

Comparison of the reflectance from a semi-infinite two-layered turbid medium calculated with method C (solid curve) and method D (circles). The parameters are ${\mu }_{s1}^{\prime }=1.0{\mathrm{mm}}^{-1}$ , ${\mu }_{a1}=0.02{\mathrm{mm}}^{-1}$ , ${\mu }_{s2}^{\prime }=1.5{\mathrm{mm}}^{-1}$ , ${\mu }_{a2}=0.005{\mathrm{mm}}^{-1}$ , $l_1=5\mathrm{mm}$ , $l_2=\infty \mathrm{mm}$ , and $n_1=n_2=1.4$ .

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 .

Fig. 3

Relative difference of the curves shown in Fig. 2 [(|method $\mathrm{C}-\text{method}$ D|)/method C].

The differences are $<{10}^{-7}$ for distances $\lesssim 40\mathrm{mm}$ 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 $\Phi (\rho ,z)$ 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 $\approx 3\mathrm{s}$ (for $\approx 1000$ distances; by using the FFT one obtains many distances in one calculation), $\approx 100\mathrm{ms}$ (for one distance), $\approx 35\mathrm{ms}$ (for one distance), and $\approx 10\mathrm{ms}$ (for one distance) were needed, respectively. When the FFT was used to evaluate the correlation in method D (see Sec. 2) even $\approx 1000$ distances could be evaluated in $\approx 10\mathrm{ms}$ . 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.

3.2.

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.

Fig. 4

Comparison of the solution of the semi-infinite two-layered diffusion equation and Monte Carlo simulations. The parameters are ${\mu }_{s1}^{\prime }=1.3{\mathrm{mm}}^{-1}$ , ${\mu }_{a1}=0.04{\mathrm{mm}}^{-1}$ , ${\mu }_{s2}^{\prime }=0.8{\mathrm{mm}}^{-1}$ , ${\mu }_{a2}=0.005{\mathrm{mm}}^{-1}$ , $l_1=3\mathrm{mm}$ , and $l_2=\infty \mathrm{mm}$ .

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 $0.5{\mathrm{mm}}^{-1}<{\mu }_s^{\prime }<1.5{\mathrm{mm}}^{-1}$ and $0.0{\mathrm{mm}}^{-1}<{\mu }_a<0.04{\mathrm{mm}}^{-1}$ , but were the same for the solution of the diffusion equation and Monte Carlo simulation. The thickness of each layer was $1\mathrm{mm}$ . Thus, the thickness of the whole turbid medium was $20\mathrm{mm}$ . 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.

Fig. 5

Comparison of the solution of a 20-layered turbid medium calculated with the diffusion equation (solid curve) with Monte Carlo simulations (circles). Optical coefficients were randomly chosen but were the same for both methods. The refractive index of all layers was 1.4.

Figure 6 gives the relative difference of both curves to better visualize the performance of the diffusion solution compared to the Monte Carlo simulation.

Fig. 6

Relative difference of the curves shown in Fig. 5 [(Monte $\text{Carlo}-\text{diffusion}$ )/Monte Carlo].

At small distances $(\rho <1\mathrm{mm})$ the differences are large, for intermediate distances the differences are about 5%, and for large distances $(\rho >8\mathrm{mm})$ 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 $24\mathrm{mm}$ .

Fig. 7

Comparison of the solution of the semi-infinite four-layered diffusion equation (solid curve) with Monte Carlo simulations (circles). The parameters are ${\mu }_{s1}^{\prime }=1.6{\mathrm{mm}}^{-1}$ , ${\mu }_{a1}=0.005{\mathrm{mm}}^{-1}$ , ${\mu }_{s2}^{\prime }=0.9{\mathrm{mm}}^{-1}$ , ${\mu }_{a2}=0.015{\mathrm{mm}}^{-1}$ , ${\mu }_{s3}^{\prime }=0.7{\mathrm{mm}}^{-1}$ , ${\mu }_{a3}=0.01{\mathrm{mm}}^{-1}$ , ${\mu }_{s4}^{\prime }=1.2{\mathrm{mm}}^{-1}$ , ${\mu }_{a4}=0.02{\mathrm{mm}}^{-1}$ , $n_1=1.4$ , $n_2=1.3$ , $n_3=1.5$ , $n_4=1.4$ , $l_1=4\mathrm{mm}$ , $l_2=5\mathrm{mm}$ , $l_3=3\mathrm{mm}$ , and $l_4=12\mathrm{mm}$ .

Figure 7 shows good agreement between the two calculations for all distances because in transmission, the diffusion approximation is valid also for small $\rho$ values due to the large value of the product ${\mu }_s^{\prime }{L}_4$ , 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.

Fig. 8

Relative difference [(Monte $\text{Carlo}-\text{diffusion}$ )/Monte Carlo] of the curves shown in Fig. 7.

4.

Discussion

The solutions of the $N$ -layered diffusion equation were calculated in the steady-state domain for turbid media having a finite and an infinitely thick $N$ ’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 $s$ 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 $N$ -layered turbid media by comparing the results of the four methods. We found relative differences that were typically smaller than ${10}^{-7}$ . 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 $10\mathrm{ms}$ even for a large number of reflectance or transmittance values at different distances from the source.

We compared the solutions of the $N$ -layered diffusion equation having an infinite and finite thick $N$ ’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.²⁵

Note that the derived solutions can be easily altered to obtain solutions of the heat conduction equation in an $N$ -layered material. Also, different boundary conditions, for example, the partial boundary condition, can be implemented with our method.²²

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.