Light diffusion in N-layered turbid media: frequency and time domains

André Liemert; Alwin Kienle

doi:10.1117/1.3368682

1 March 2010 Light diffusion in N-layered turbid media: frequency and time domains

André Liemert, Alwin Kienle

Author Affiliations +

Journal of Biomedical Optics, Vol. 15, Issue 2, 025002 (March 2010). https://doi.org/10.1117/1.3368682

Abstract

We deal with light diffusion in mismatched N-layered turbid media having a finite or an infinitely thick Nth layer. We focus on time-resolved light propagation in both the frequency and time domains. Based on our results for the steady-state domain, solutions of the N-layered diffusion equations in the frequency and time domains are obtained by applying the Fourier transform technique. Different methods for calculation of the inverse Fourier transform are studied to validate the solutions, showing relative differences typically smaller than 10^-6. The solutions are compared to Monte Carlo simulations, revealing good agreement. Finally, by applying the Laplace and Fourier transforms we derive a fast (≈ 1 ms) and accurate analytical solution for the time domain reflectance from a two-layered turbid medium having equal reduced scattering coefficients and refractive indices in both layers.

1. Introduction

In the first of our two papers dealing with light diffusion in $N$ -layered turbid media we studied the light propagation in the steady-state domain.¹ 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.² Contrarily, time-resolved measurements are more complex and expensive, but they contain more information on the investigated tissues.³ 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 $100 MHz$ and the demodulation and phase shift of the remitted or transmitted light is detected.⁶

In the literature, experiments in all three measurement regimes were reported assuming that the investigated turbid media can be considered homogeneous or heterogeneous.⁷ 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.¹

In this paper, we present the derived solutions of the time-resolved reflectance and transmittance from an $N$ -layered turbid medium having a finite thick $N$ ’th layer, and the time-resolved reflectance from an $N$ -layered turbid medium having an infinitely thick $N$ ’th layer. The solutions in the frequency domain are obtained by using the Fourier transform technique, as reported in the accompanying paper.¹ 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.

2. Theory

In this section, the derivation of the solutions of the $N$ -layered diffusion equation in the frequency and time domains is presented. We start with the diffusion equation in the time domain:

Eq. 1

(\frac{1}{c} \frac{\partial}{\partial t} + μ_{a} - D Δ) Φ (r, t) = S (r, t),

where

Φ (r, t)

is the fluence rate in real space,

μ_{a}

is the absorption coefficient,

D = 1 ∕ (3 μ_{s}^{'})

the diffusion coefficient, and

μ_{s}^{'}

is the reduced scattering coefficient of the turbid medium. By solving Eq. 1 for a

δ

source in time and space

[S (r, t) = δ (r) δ (t)]

Green’s function

G (r, t)

of the problem is obtained with

Eq. 2

(\frac{1}{c} \frac{\partial}{\partial t} + μ_{a} - D Δ) G (r, t) = δ (r) δ (t) .

For an arbitrary source in time and space $S (r, t)$ , the fluence rate can be calculated with the knowledge of Green’s function by applying

Eq. 3

Φ (r, t) = \int_{- \infty}^{\infty} \int_{V} G (r^{'}, t^{'}) S (r - r^{'}, t - t^{'}) d^{3} r^{'} d t^{'} .

We propose two different methods to calculate $G (r, t)$ . 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 $t$ coordinate and the 2-D Fourier transform relative to the $x$ and $y$ coordinates to obtain an ordinary differential equation depending on the $z$ coordinate in the inverse Laplace and Fourier space. Then, the boundary conditions are applied relative to the $z$ 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 $t$ 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.⁸ 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.¹ 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:

Eq. 4

S (r, t) = S (r) e^{i ω t},

with

i = \sqrt{- 1}

and

ω = 2 π f

, where

f

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

Eq. 5

Φ (r, t) = \int_{- \infty}^{\infty} \int_{V} G (r^{'}, t^{'}) S (r - r^{'}) \exp [i ω (t - t^{'})] d^{3} r^{'} d t^{'} .

Subsequently, the 3-D integration in the spatial domain is performed. By assuming a $δ$ source in space, $S (r - r^{'}) = δ (r - r^{'})$ , we obtain

Eq. 6

Φ (r, t) = e^{i ω t} \int_{- \infty}^{\infty} G (r, t^{'}) \exp (- i ω t^{'}) d t^{'} = G (r, ω) e^{i ω t},

where

G (r, ω)

is Green’s function in the frequency domain. We can derive

G (r, ω)

by solving the following differential equation, which is obtained from the Fourier transform relative to the

t

coordinate of Eq. 2:

Eq. 7

(Δ - α^{2}) G (r, ω) = - \frac{1}{D} δ (r), α^{2} = \frac{μ_{a}}{D} + i \frac{ω}{D c},

where

Δ

is the Laplace operator, and

c

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:

Eq. 8

G (r, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} G (r, ω) e^{i ω t} d ω .

For a vanishing intensity oscillation $(ω = 0 Hz)$ we get from Eq. 7 the steady-state diffusion equation, which was solved for $N$ -layered turbid media in the accompanying paper.¹ Thus, to obtain the solution of the $N$ -layered diffusion equation in the frequency domain assuming a $δ$ -source at $x = y = 0$ and $z = z_{0}$ :

Eq. 9

Δ Φ_{1} (x, y, z, ω) - (\frac{μ_{a 1}}{D_{1}} + i \frac{ω}{c_{1} D_{1}}) Φ_{1} (x, y, z, ω) = - \frac{1}{D_{1}} δ (x, y, z - z_{0}), 0 ⩽ z < l_{1},

Eq. 10

Δ Φ_{k} (x, y, z, ω) - (\frac{μ_{a k}}{D_{k}} + i \frac{ω}{c_{k} D_{k}}) Φ_{k} (x, y, z, ω) = 0, \sum_{j = 1}^{k - 1} l_{j} < z ⩽ \sum_{j = 1}^{k} l_{j}, k = 2, 3, \dots, N,

all the solutions of the

N

-layered steady-state diffusion equation derived in the accompanying study¹ can be used simply by replacing

μ_{a k} ∕ D_{k}

with

μ_{a k} ∕ D_{k} + i ω ∕ (c_{k} D_{k})

, where

c_{k}

is the velocity of light in layer

k

. In addition, the four methods¹ 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 $N$ -layered turbid medium in the frequency domain is computed from the fluence rate by

Eq. 11

R (ρ, ω) = \int_{2 π} d Ω [1 - R_{f} (θ)] \frac{1}{4 π} [Φ_{1} (ρ, z = 0, ω) + 3 D_{1} {\frac{\partial Φ_{1} (ρ, z, ω)}{\partial z} |}_{z = 0} \cos (θ)] \cos (θ),

where

R_{f} (θ)

is¹ the Fresnel reflection coefficient and

ρ = {(x^{2} + y^{2})}^{1 ∕ 2}

. For large distances from the source

(ρ ⪢ 1 ∕ μ_{s}^{'})

it is sufficient to use the Fick’s law for calculation of the reflectance⁹

Eq. 12

R (ρ, ω) = D_{1} {\frac{\partial Φ_{1} (ρ, z, ω)}{\partial z} |}_{z = 0} .

In the frequency domain, the intensity-modulated source causes that the reflectance measured at the detector position $ρ$ is also modulated by the same frequency $f$ 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 $M$ :

Eq. 13

θ (ρ, ω) = \arctan {\frac{Im [R (ρ, ω)]}{Re [R (ρ, ω)]}},

Eq. 14

M (ρ, ω) = {\frac{Re {[R (ρ, ω)]}^{2} + Im {[R (ρ, ω)]}^{2}}{Re {[R (ρ, ω = 0)]}^{2}}}^{1 ∕ 2} .

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 paper¹ are used by replacing $μ_{a k} ∕ D_{k}$ with $μ_{a k} + i ω ∕ (c_{k} D_{k})$ . The solutions in the time domain are calculated by computing the solution in the frequency domain for many frequencies (e.g. 512) and by using¹⁰ an FFT. Note that because $Φ (r, t)$ 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 $Φ (r, ω)$ . Then, $Φ (r, t)$ is obtained by an FFT of two times the real part or two times the imaginary part of $Φ (r, ω)$ , accelerating the calculations by a factor of 2.

The solutions of the $N$ -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,¹ 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.¹¹ The results in the frequency domain are obtained by Fourier transforming (via an FFT) the Monte Carlo data calculated in the time domain.

3. Results

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.

3.1.

Comparison of Different Methods for Solving the Diffusion Equation

In the accompanying paper,¹ 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 $ρ = 22 mm$ .

Fig. 1

Comparison of the reflectance from a semi-infinite two-layered turbid medium calculated with two different methods for computing the inverse Fourier transform. The parameters are $μ_{s 1}^{'} = 1.1 {mm}^{- 1}$ , $μ_{a 1} = 0.009 {mm}^{- 1}$ , $μ_{s 2}^{'} = 1.3 {mm}^{- 1}$ , $μ_{a 2} = 0.03 {mm}^{- 1}$ , $l_{1} = 5 mm$ , $l_{2} = \infty mm$ , $n_{1} = 1.4$ , and $n_{2} = 1.6$ .

Both curves are indistinguishable in the figure. Thus, Fig. 2 shows the relative difference between both curves. The differences are smaller than $10^{- 6}$ 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.

Fig. 2

Relative difference of the data [(|method $C - method$ D|)/method C] shown in Fig. 1.

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 $ρ = 30 mm$ .

Fig. 3

Comparison of the reflectance from a semi-infinite two-layered turbid medium calculated with an analytical solution (circles) [see Eq. 41] and an numerical solution (method D, solid curve). The parameters are $μ_{s 1}^{'} = 1.2 {mm}^{- 1}$ , $μ_{a 1} = 0.005 {mm}^{- 1}$ , $μ_{s 2}^{'} = 1.2 {mm}^{- 1}$ , $μ_{a 2} = 0.05 {mm}^{- 1}$ , $l_{1} = 3 mm$ , $l_{2} = \infty mm$ , $n_{1} = 1.4$ , and $n_{2} = 1.4$ .

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 $10^{- 6}$ 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.

Fig. 4

Relative difference of the data shown in Fig. 3 [(|method $C - method$ D|)/method C].

3.2.

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 $20.5 mm$ 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).

Fig. 5

Comparison of the solution of the four-layered diffusion equation (solid curve) with Monte Carlo simulations (circles) in the time domain. The parameters are $μ_{s 1}^{'} = 1.5 {mm}^{- 1}$ , $μ_{a 1} = 0.04 {mm}^{- 1}$ , $μ_{s 2}^{'} = 0.8 {mm}^{- 1}$ , $μ_{a 2} = 0.02 {mm}^{- 1}$ , $μ_{s 3}^{'} = 1.2 {mm}^{- 1}$ , $μ_{a 3} = 0.01 {mm}^{- 1}$ , $μ_{s 4}^{'} = 1.2 {mm}^{- 1}$ , $μ_{a 4} = 0.015 {mm}^{- 1}$ , $l_{1} = 2 mm$ , $l_{2} = 3 mm$ , $l_{3} = 5 mm$ , $l_{4} = \infty mm$ , and $n_{k} = 1.4$ .

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 $0.5 ns$ , 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%.

Fig. 6

Relative differences between the Monte Carlo simulation and the solution of the diffusion equation [(Monte $Carlo - diffusion$ )/Monte Carlo] for the data shown in Fig. 5.

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 $20.5 mm$ 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 $f$ . 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.

Fig. 7

Comparison of the solution of the four-layered diffusion equation with Monte Carlo simulations in the frequency domain for the same turbid medium that is discussed in Fig. 5. The phase values calculated by the Monte Carlo method (circles) and those calculated by the solution of the diffusion equation (solid curve) are shown.

Fig. 8

Modulation calculated by the Monte Carlo method (circles) and those calculated by the solution of the diffusion equation (solid curve) are shown for the same turbid medium considered in Fig. 7.

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.

4. Discussion

Solutions of the $N$ -layered diffusion equation were derived for turbid media having a finite and infinite thick $N$ ’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 paper¹ 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 $1 ms$ 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 $10^{- 6}$ 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 $N$ -layered diffusion equation in the time and steady-state domains are available on the Internet.¹²

Acknowledgment

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076).

References

1.

A. Liemert and A. Kienle, “Light diffusion in

N

-layered turbid media: steady state domain,” J. Biomed. Opt., 15 (2), 025003 (2010). 1083-3668 Google Scholar

2.

M. S. Patterson, B. C. Wilson, and D. R. Wyman, “The propagation of optical radiation in tissue. II: Optical properties of tissues and resulting fluence distributions,” Lasers Med. Sci., 6 379 –389 (1991). https://doi.org/10.1007/BF02042460 0268-8921 Google Scholar

3.

S. R. Arridge, “Optical tomography in medical imaging inverse problems,” Inverse Probl., 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar

4.

T. Svensson, J. Swartling, P. Taroni, A. Torricelli, P. Lindblom, C. Ingvar, and S. Andersson-Engles, “Characterization of normal breast tissue heterogeneity using time-resolved near-infrared spectroscopy,” Phys. Med. Biol., 50 2559 –2571 (2005). https://doi.org/10.1088/0031-9155/50/11/008 0031-9155 Google Scholar

5.

A. Pifferi, A. Torricelli, L. Spinelli, D. Contini, R. Cubeddu, F. Martelli, G. Zaccanti, A. Tosi, A. Dalla Mora, F. Zappa, and S. Cova, “Time-resolved diffuse reflectance using small source-detector separation and fast single-photon gating,” Phys. Rev. Lett., 100 138101 (2008). https://doi.org/10.1103/PhysRevLett.100.138101 0031-9007 Google Scholar

6.

A. Cerussi, N. Shah, D. Hsiang, A. Durkin, J. Butler, and B. J. Tromberg, “In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy,” J. Biomed. Opt., 11 044005 (2006). https://doi.org/10.1117/1.2337546 1083-3668 Google Scholar

7.

S. L. Jacques and B. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt., 13 041302 (2008). https://doi.org/10.1117/1.2967535 1083-3668 Google Scholar

8.

I. Dayan, S. Havlin, and G. H. Weiss, “Photon migration in a two-layer turbid medium. A diffusion analysis,” J. Mod. Opt., 39 1567 –1582 (1992). https://doi.org/10.1080/09500349214551581 0950-0340 Google Scholar

9.

A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A, 14 246 –254 (1997). https://doi.org/10.1364/JOSAA.14.000246 0740-3232 Google Scholar

10.

W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Pascal, Cambridge University, Cambridge, UK (1990). Google Scholar

11.

A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol., 41 2221 –2227 (1996). https://doi.org/10.1088/0031-9155/41/10/026 0031-9155 Google Scholar

12.

. http://www.uni-ulm.de/ilm/index.php?id=10020200 Google Scholar

13.

M. S. Patterson, B. Chance, and B. C. Wilson, “Time-resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 0003-6935 Google Scholar

14.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York (1971). Google Scholar

Appendices

Appendix

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,

Eq. 15

L {G (x, y, z, t)} = G (x, y, z, p) = \int_{0}^{\infty} G (x, y, z, t) e^{- p t} d t,

we obtain

Eq. 16

\frac{p}{c} G (x, y, z, p) - D \nabla^{2} G (x, y, z, p) + μ_{a} G (x, y, z, p) = δ (x, y, z) .

Then, the 2-D Fourier transform with respect to

x

and

y

is employed:

Eq. 17

G (s_{1}, s_{2}, z, p) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} G (x, y, z, p) \exp [i (s_{1} x + s_{2} y)] d x d y,

delivering

Eq. 18

\frac{\partial^{2} G (s_{1}, s_{2}, z, p)}{\partial z^{2}} - α^{2} G (s_{1}, s_{2}, z, p) = - \frac{1}{D} δ (z),

where

G (s_{1}, s_{2}, z, p)

is the fluence rate in the Laplace-Fourier space, and

α^{2} = s_{1}^{2} + s_{2}^{2} + p ∕ D c + μ_{a} ∕ D

. Next, the boundary conditions are introduced. First,

G (s_{1}, s_{2}, z, p)

is zero infinitely far away from the source:

Eq. 19

G (s_{1}, s_{2}, z = \infty, p) = G (s_{1}, s_{2}, z = - \infty, p) = 0 .

Second, the source term at

z = 0

must be considered similar as in the accompanying paper:¹

Eq. 20

G (s_{1}, s_{2}, z = 0^{-}, p) = G (s_{1}, s_{2}, z = 0^{+}, p),

Eq. 21

\frac{\partial G (s_{1}, s_{2}, z = 0^{-}, p)}{\partial z} - \frac{\partial G (s_{1}, s_{2}, z = 0^{+}, p)}{\partial z} = \frac{1}{D} .

The ordinary differential Eq. 18 can be solved using

Eq. 22

G (s_{1}, s_{2}, z, p) = c_{1} e^{α z} + c_{2} e^{- α z} z < 0,

Eq. 23

G (s_{1}, s_{2}, z, p) = c_{3} e^{α z} + c_{4} e^{- α z} z ⩾ 0 .

With the preceding four boundary conditions, the unknown constants are determined resulting in the following solutions:

Eq. 24

G (s_{1}, s_{2}, z, p) = \frac{1}{2 α D} e^{α z} z < 0,

Eq. 25

G (s_{1}, s_{2}, z, p) = \frac{1}{2 α D} e^{- α z} z ⩾ 0 .

Inserting $α$ in Eq. 24 yields

Eq. 26

G (s_{1}, s_{2}, z, p) = {(\frac{c}{4 D})}^{1 ∕ 2} \frac{e^{(- z ∕ \sqrt{D c}) {[p + c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})]}^{1 ∕ 2}}}{{[p + c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})]}^{1 ∕ 2}} z < 0,

Eq. 27

G (s_{1}, s_{2}, z, p) = {(\frac{c}{4 D})}^{1 ∕ 2} \frac{e^{(- z ∕ \sqrt{D c}) {[p + c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})]}^{1 ∕ 2}}}{{[p + c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})]}^{1 ∕ 2}} z ⩾ 0 .

Subsequently, these equations are two-dimensionally inverse Fourier and Laplace transformed to obtain $G (x, y, z, t)$ . 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:

Eq. 28

L^{- 1} [\frac{\exp (- k \sqrt{p})}{\sqrt{p}}] = \frac{1}{\sqrt{π t}} \exp (- k^{2} ∕ 4 t),

Eq. 29

L [e^{- a t} G (t)] = G (p + a), G (p) = L {G (t)},

where

a = c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})

, one gets

Eq. 30

G (s_{1}, s_{2}, z, t) = {(\frac{c}{4 D π t})}^{1 ∕ 2} \exp {- [c μ_{a} + D c (s_{1}^{2} + s_{2}^{2})] t} \exp (- \frac{z^{2}}{4 D c t}) .

Finally, the 2-D inverse Fourier transform is employed using the formula (and its corresponding equation with respect to $y$ and $s_{2}$ ):

Eq. 31

F [\exp (- \frac{b x^{2}}{2})] = \frac{\sqrt{2 π}}{\sqrt{b}} \exp (- \frac{s^{2}}{2 b}),

which delivers the well-known Green’s function for the time-resolved diffusion equation:¹³

Eq. 32

G (x, y, z, t) = \frac{c}{{(4 π D c t)}^{3 ∕ 2}} \exp (- μ_{a} c t) \exp (- \frac{x^{2} + y^{2} + z^{2}}{4 D c t}) .

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 $z$ . Then the boundary conditions for the two-layered turbid medium were applied. The resulting equation system was solved, as presented in the accompanying paper.¹ 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 $D = D_{1} = D_{2}$ and $n_{1} = n_{2}$ ,

Eq. 33

G (s, z, p) = \frac{\sinh [α_{1} (z_{0} + z_{b})]}{D α_{1}} \frac{α_{1} \cosh [α_{1} (l_{1} - z)] + α_{2} \sinh [α_{1} (l_{1} - z)]}{α_{1} \cosh [α_{1} (l_{1} + z_{b})] + α_{2} \sinh [α_{1} (l_{1} + z_{b})]} - \frac{\sinh [α_{1} (z_{0} - z)]}{D α_{1}},

with

α_{k}^{2} = s^{2} + μ_{a k} ∕ D + p ∕ D c

,

s^{2} = s_{1}^{2} + s_{2}^{2}

, and

c

is the velocity of light in both layers.

We make use of the following auxiliary equation

Eq. 34

\frac{α_{1} \cosh [α_{1} (l_{1} - z)] + α_{2} \sinh [α_{1} (l_{1} - z)]}{α_{1} \cosh [α_{1} (l_{1} + z_{b})] + α_{2} \sinh [α_{1} (l_{1} + z_{b})]} = \frac{(α_{1} + α_{2}) \exp [α_{1} (l_{1} - z)] + (α_{1} - α_{2}) \exp - α_{1} (l_{1} - z)}{(α_{1} + α_{2}) \exp [α_{1} (l_{1} + z_{b})] + (α_{1} - α_{2}) \exp [- α_{1} (l_{1} + z_{b})]} = \frac{\exp [α_{1} (l_{1} - z)] + (α_{1} - α_{2} ∕ α_{1} + α_{2}) \exp [- α_{1} (l_{1} - z)]}{\exp [α_{1} (l_{1} + z_{b})] + (α_{1} - α_{2} ∕ α_{1} + α_{2}) \exp [- α_{1} (l_{1} + z_{b})]} = \exp [- α_{1} (z + z_{b})] \frac{1 + D_{a} \exp [- 2 α_{1} (l_{1} - z)]}{1 + D_{a} \exp [- 2 α_{1} (l_{1} + z_{b})]} = \exp [- α_{1} (z + z_{b})] {1 + D_{a} \exp [- 2 α_{1} (l_{1} - z)]} \sum_{m = 0}^{\infty} {(- D_{a})}^{m} \exp [- 2 m α_{1} (l_{1} + z_{b})],

using

D_{a} = (α_{1} - α_{2}) ∕ (α_{1} + α_{2})

and the geometric series in the last step. Inserting Eq. 34 into Eq. 33 yields

Eq. 35

G (s, z, p) = \frac{\exp [α_{1} (z_{0} - z)] - \exp [- α_{1} (z_{0} + z + 2 z_{b})]}{2 D α_{1}} (1 + {1 - \exp [2 α_{1} (z + z_{b})]} \sum_{m = 1}^{\infty} {(- D_{a})}^{m} \exp [- 2 m α_{1} (l_{1} + z_{b})]) - \frac{\exp [α_{1} (z_{0} - z)] - \exp [- α_{1} (z_{0} - z)]}{2 D α_{1}} .

By rearranging and simplifying Eq. 35 it follows that

Eq. 36

G (s, z, p) = \frac{\exp [α_{1} (z_{0} - z)] - \exp [- α_{1} (z_{0} + z + 2 z_{b})]}{2 D α_{1}} + \frac{1}{2 D α_{1}} \sum_{m = 1}^{\infty} [(\exp {- α_{1} [z_{0} - z + 2 m (l + z_{b})]} - \exp {- α_{1} [z_{0} + z + 2 z_{b} + 2 m (l + z_{b})]} + \exp {- α_{1} [- z_{0} + z + 2 m (l + z_{b})]} - \exp {- α_{1} [- z_{0} - z - 2 z_{b} + 2 m (l + z_{b})]}) {(- D_{a})}^{m}] .

Subsequently, we perform the inverse transforms. Employing the Laplace transform formulas¹⁴

Eq. 37

\frac{\exp [- a {(p + b)}^{1 ∕ 2}]}{{(p + b)}^{1 ∕ 2}}, a > 0 \Leftrightarrow \frac{1}{\sqrt{π t}} e^{- b t} \exp (- \frac{a^{2}}{4 t}) ε (t),

Eq. 38

{[\frac{{(p + a)}^{1 ∕ 2} - {(p + b)}^{1 ∕ 2}}{{(p + a)}^{1 ∕ 2} + {(p + b)}^{1 ∕ 2}}]}^{m}, m > 0 \Leftrightarrow \frac{m}{t} \exp [- \frac{1}{2} (a + b) t] I_{m} (\frac{a - b}{2} t) ε (t),

Eq. 39

X (p) \cdot Y (p) \Leftrightarrow X (t) * Y (t) = \int_{0}^{t} X (τ) Y (t - τ) d τ,

where

ε (t)

is the Heaviside step function,

I_{m}

is the modified Bessel function of the

m

’th order, and the inverse Fourier transform with respect to the

x

coordinate (and accordingly to the

y

coordinate)

Eq. 40

\exp (- D c t \cdot s_{1}^{2}) \Leftrightarrow \frac{1}{{(4 π D c t)}^{1 ∕ 2}} \exp (- \frac{x^{2}}{4 D c t}),

we, finally, obtain Green’s function in the real space

Eq. 41

G (r, t) = \frac{c}{{(4 π D c t)}^{3 ∕ 2}} \exp (- μ_{a 1} c t) \exp (- \frac{ρ^{2}}{4 D c t}) {\exp [- \frac{{(z - z_{0})}^{2}}{4 D c t}] - \exp [- \frac{{(z + z_{0} + 2 z_{b})}^{2}}{4 D c t}]} + \frac{1}{4 π D c t} \exp (- \frac{ρ^{2}}{4 D c t}) \int_{0}^{t} \frac{1}{2} {[\frac{c}{D π (t - τ)}]}^{1 ∕ 2} \exp [- μ_{a 1} c (t - τ)] \exp [- \frac{1}{2} (μ_{a 1} + μ_{a 2}) c τ] \times \sum_{m = - \infty}^{\infty} (\exp {- \frac{{[z + z_{0} + 2 z_{b} + 2 m (l_{1} + z b)]}^{2}}{4 D c (t - τ)}} - \exp {- \frac{{[z_{0} - z + 2 m (l_{1} + z b)]}^{2}}{4 D c (t - τ)}}) \frac{| m |}{τ} I_{| m |} (\frac{μ_{a 2} - μ_{a 1}}{2} c τ) d τ .

Equation 41 was evaluated by Gauss integration (see Fig. 3). The first term of Eq. 41 is identical to the solution of the diffusion equation for a semi-infinite homogeneous geometry.

To accelerate the calculations we derive an alternative formula for $G (r, t)$ . 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 $G (r, t)$ given in Eq. 41. The second part of Eq. 36 is transformed as follows:

Eq. 42

\frac{1}{2 D α_{1}} {\exp [- α_{1} (z_{0} - z)] - \exp [- α_{1} (z_{0} + z + 2 z_{b})] + \exp [α_{1} (z_{0} - z)] - \exp [α_{1} (z_{0} + z + 2 z_{b})]} \sum_{m = 1}^{\infty} {- D_{a} \exp [- 2 α_{1} (l + z_{b})]}^{m} = \frac{1}{2 D α_{1}} {\cosh [α_{1} (z_{0} - z)] - \cosh [α_{1} (z_{0} + z + 2 z_{b})]} {\frac{1}{1 + D_{a} \exp [- 2 α_{1} (l + z_{b})]} - 1} = \frac{1}{2 D α_{1}} {\cosh [α_{1} (z_{0} + z + 2 z_{b})] - \cosh [α_{1} (z_{0} - z)]} \frac{D_{a}}{D_{a} + \exp [2 α_{1} (l + z_{b})]} .

Then, Eq. 42 is inverse Fourier transformed with respect to the $x$ and $y$ 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 $a = D c s^{2}$ . In the inverse Laplace transform with respect to the time variable we used $p = i ω$ , because the imaginary axis is within the region of convergence. Thus, we obtain

Eq. 43

G (r, t) = \frac{c}{{(4 π D c t)}^{3 ∕ 2}} \exp (- μ_{a 1} c t) \exp (- \frac{ρ^{2}}{4 D c t}) {\exp [- \frac{{(z - z_{0})}^{2}}{4 D c t}] - \exp [- \frac{{(z + z_{0} + 2 z_{b})}^{2}}{4 D c t}]} + \frac{1}{4 π D c t} \exp (- \frac{ρ^{2}}{4 D c t}) \frac{1}{2 π} \int_{- \infty}^{\infty} \frac{{\cosh [α_{1} (z_{0} + z + 2 z_{b})] - \cosh [α_{1} (z_{0} - z)]} D_{a}}{D α_{1} {D_{a} + \exp [2 α_{1} (l + z_{b})]}} e^{i ω t} d ω,

where

Eq. 44

α_{1} = {(\frac{μ_{a 1}}{D} + i \frac{ω}{D c})}^{1 ∕ 2} .

The integral in Eq. 43 is identical to the Fourier transform and can be evaluated efficiently with the FFT method, which results in calculation times of the order of

1 ms

.

Citation Download Citation

André Liemert and Alwin Kienle "Light diffusion in N-layered turbid media: frequency and time domains," Journal of Biomedical Optics 15(2), 025002 (1 March 2010). https://doi.org/10.1117/1.3368682

Published: 1 March 2010

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 61 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Diffusion

Fourier transforms

Monte Carlo methods

Reflectivity

Modulation

Refractive index

Scattering

1.

Introduction

2.

Theory

Eq. 1

Eq. 2

Eq. 3

Eq. 4

Eq. 5

Eq. 6

Eq. 7

Eq. 8

Eq. 9

Eq. 10

Eq. 11

Eq. 12

Eq. 13

Eq. 14

3.

Results

3.1.

Comparison of Different Methods for Solving the Diffusion Equation

Fig. 1

Fig. 2

Fig. 3

Fig. 4

3.2.

Comparison of the Solutions of the Diffusion Equation with Monte Carlo Simulations

Fig. 5

Fig. 6

Fig. 7

Fig. 8

4.

Discussion

Acknowledgment

References

Appendices

Appendix

Appendix A: Infinite Turbid Medium

Eq. 15

Eq. 16

Eq. 17

Eq. 18

Eq. 19

Eq. 20

Eq. 21

Eq. 22

Eq. 23

Eq. 24

Eq. 25

Eq. 26

Eq. 27

Eq. 28

Eq. 29

Eq. 30

Eq. 31

Eq. 32

Appendix B: Two-Layered Turbid Medium

Eq. 33

Eq. 34

Eq. 35

Eq. 36

Eq. 37

Eq. 38

Eq. 39

Eq. 40

Eq. 41

Eq. 42

Eq. 43

Eq. 44

Show All Keywords

Keywords/Phrases

Search In:

Publication Years