Light diffusion in a turbid cylinder. II. Layered case

doi:10.1364/OE.18.009266

« Show journal navigation

Light diffusion in a turbid cylinder. II. Layered case

André Liemert and Alwin Kienle »View Author Affiliations

Optics Express, Vol. 18, Issue 9, pp. 9266-9279 (2010)
http://dx.doi.org/10.1364/OE.18.009266

View Full Text Article

Acrobat PDF (915 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Theory
3. Results
4. Discussion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (22)
Cited By
Figures (9)
Metrics

Abstract

This paper is the second of two dealing with light diffusion in a turbid cylinder. The diffusion equation was solved for an N-layered finite cylinder. Solutions are given in the steady-state, frequency, and time domains for a point beam incident at an arbitrary position of the first layer and for a circular flat beam incident at the middle of the cylinder top. For special cases the solutions were compared to other solutions of the diffusion equation showing excellent agreement. In addition, the derived solutions were validated by comparison with Monte Carlo simulations. In the time domain we also derived a fast solution (≈ 10ms) for the case of equal reduced scattering coefficients and refractive indices in all layers.

1. Introduction

In the accompanying paper we described the light diffusion in a homogeneous cylindrical turbid medium [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

]. In this study we present solutions of the diffusion equation for cylindrical turbid media having an arbitrary number of layers.

The light propagation in layered turbid media has been investigated for therapeutical or diagnostic applications because many parts of the body exhibit a layered structure, e.g. the arm, the leg, or the head. The transport theory and an approximation of it, the diffusion theory, have been normally applied to obtain quantitative data [2

A. Ishimaru, “Wave Propagation and Scattering in Random Media,” Academic Press, New York (1978).

]. The advantage of the diffusion equation is that it can be solved analytically for layered geometries yielding fast and accurate solutions compared to the normally applied numerical solutions.

For laterally infinitely extended turbid media Dayan et al. presented approximative solutions of the two-layered diffusion equation [3

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). [CrossRef]

]. We derived exact solutions of the two-layered diffusion equation in the steady-state, frequency, and time domains using the Fourier transform formalism [4

A. Kienle, M.S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

, 5

A. Kienle, T. Glanzmann, G. Wagnières, and H. van den Bergh, “Investigation of Two-Layered Turbid Media with Time-Resolved Reflectance,” Appl. Opt. 37, 6852–6862 (1998). [CrossRef]

, 6

A. Kienle and T. Glanzmann, “In Vivo Determination of the Optical Properties of Muscle Using a Layered-Model,” Phys. Med. Biol. 44, 2689–2702 (1999). [CrossRef] [PubMed]

]. This approach was, recently, extended by us and others for an infinite number of layers [8

X.C. Wang and S.M. Wang, “Light Transport Modell in a N-Layered Mismatched Tissue,” Waves Rand. Compl. Media 16, 121–135 (2006). [CrossRef]

, 7

J.M. Tualle, H.M. Nghiem, D. Ettori, R. Sablong, E. Tinet, and S. Avrillier, “Asymptotic Behavior and Inverse Problem in Layered Scattering Media,” J. Opt. Soc. Am. A 21, 24–34 (2004). [CrossRef]

, 9

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt. , accepted. [PubMed]

, 10

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. , accepted. [PubMed]

], whereas the diffusion equation for a finite number of layers was investigated also by several other groups [11

G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium,” Appl. Opt. 37, 7401–7409 (1998). [CrossRef]

, 12

S-H. Tseng, C. Hayakawa, B.J. Tromberg, J. Spanier, and A.J. Durkin, “Quantitative Spectroscopy of Superficial Turbid Media,” Opt. Lett. 23, 3165–3167 (2005). [CrossRef]

, 13

G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Monte Carlo Diffusion Hybrid Model for Photon Migration in a Two-Layer Turbid Medium in the Frequency Domain,” Appl. Opt. 39, 2235–2244 (2000). [CrossRef]

, 14

M. Das, C. Xu, and Q. Zhu, “Analytical Solution for Light Propagation in a Two-Layer Tissue Structure with a Tilted Interface for Breast Imaging,” Appl. Opt. 45, 5027–5036 (2006). [CrossRef] [PubMed]

, 15

A.H. Barnett, “A Fast Numerical Method for Time-Resolved Photon Diffusion in General Stratified Turbid Media,” J. Comp. Phys. 201, 771–797 (2004). [CrossRef]

, 16

C. Donner and H.W. Jensen, “Rapid Simulations of Steady-State Spatially Resolved Reflectance and Transmittance Profiles of Multilayered Turbid Materials,” J. Opt. Soc. Am. A 23, 1382–1390 (2006). [CrossRef]

A solution of the layered diffusion equation for a cylindrical geometry has been presented by Martelli et al. for the case of two and three layers [17

F. Martelli, A. Sassaroli, S. Del Bianco, Y. Yamada Y, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for Layered Diffusive Media by the Eigenfunction Method,” Phys. Rev. E 67, 056623 (2003). [CrossRef]

, 18

F. Martelli, A. Sassaroli, S. Del Bianco, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for a Three-Layer Medium: Application to Study Photon Migration through a Simplified Adult Head Model,” Phys. Med. Biol. 52, 2827–2843 (2007). [CrossRef] [PubMed]

, 19

F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation Model for Light Propagation through Diffusive Layered Media,” Phys. Med. Biol. 50, 2159–2166 (2005). [CrossRef] [PubMed]

]. They used the eigenfunction method together with the separation of variables to derive solutions in the time domain.

In this study we present solutions of the diffusion equation for an N-layered cylinder. The solutions are given in the steady-state, frequency, and time domains. For equal reduced scattering coefficients and refractive indices in all layers we also derived a simplified solution in the time domain. For special cases, for which solutions of the diffusion equation are described in literature, we compared these solutions to those derived in this study. In addition, we validated the derived solutions with Monte Carlo simulations.

Besides in vivo applications the derived solutions can be used especially for experiments on layered phantoms, because they often have a cylindrical geometry. Thus, the phantoms’ volume does not have to be so large that the assumption of a semi-infinite geometry is valid. In addition, the obtained solutions can be easily transferred to solutions of the heat conduction equation, where multi-layered cylindrical geometries are often encountered.

2. Theory

2.1. Diffusion Theory

Solutions of the diffusion theory for an N-layered cylinder are, first, derived in the frequency domain, which include for the special case of zero modulation frequency the solutions in the steady-state domain. Then, solutions in the time domain are presented. The solutions are derived for a point beam incident onto an arbitrary location of the first layer. Thus, the beam can be directed onto the cylinder top as well as onto the cylinder barrel of the first layer. Additionally, we give the solutions for a circular flat beam incident onto the middle of the cylinder top.

Fig. 1 shows a scheme of the considered N-layered cylinder. The thickness, the refractive index, the reduced scattering and absorption coefficients of layer k are denoted by l_k , n_k , μ′ _sk , and μ_ak , respectively. The radius of the cylinder is indicated with a. As usual, it is assumed that the incident beam can be represented by an isotropic source at a distance of 1/(μ′_s1 + μ _a1) from the location of incidence at the boundary of the cylinder. The point source is located in the first layer of the cylinder.

Fig. 1. Scheme of the N-layered cylinder used in the calculations. In the figure the beam is incident onto the cylinder top, but it can also be incident onto the cylinder barrel of the first layer.

As in the accompanying paper the extrapolated boundary condition is used for the top, bottom, and barrel boundaries [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

]. The extrapolation length is calculated with

z_{bk} = \frac{1 + R_{eff, k}}{1 - R_{eff, k}} 2 D_{k},

(1)

where R _eff,k is the angle-averaged probability for reflection at the boundary between layer k and the surrounding medium and D_k = 1/(3μ′_sk) is the diffusion coefficient of layer k. For the boundary at the top and the bottom of the cylinder k equals 1 and N, respectively. For the boundary at the cylinder barrel k is assumed to be the layer at which the remitted or transmitted light is calculated, because it is more involved to solve the diffusion equation for N layers using a different boundary lengths at the cylinder barrel for each layer.

2.1.1. Solution in the frequency domain

The diffusion equation in the frequency domain is given by

ΔΦ (\vec{r}, ω) - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ (\vec{r}, ω) = - \frac{1}{D} S (\vec{r}, ω),

(2)

where Φ, c, and ω denote the fluence rate, the speed of light, and the angular frequency of the intensity modulated light, respectively. The source function for a point source in cylindrical coordinates is

S (\vec{r}, ω) = \frac{1}{ρ} δ (ρ - ρ_{0}) δ (ϕ - ϕ_{0}) δ (z - z_{0})

. Eq. 2 can be rewritten by using cylindrical coordinates as follows:

\frac{\partial^{2}}{\partial ρ^{2}} Φ + \frac{1}{ρ} \frac{\partial}{\partial ρ} Φ + \frac{1}{ρ^{2}} \frac{\partial^{2}}{\partial ϕ^{2}} Φ + \frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + i \frac{ω}{Dc}) Φ = - \frac{1}{D_{ρ}} δ (ρ - ρ_{0}) δ (ϕ - ϕ_{0}) δ (z - z_{0}) .

(3)

In the accompanying paper [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

] we presented an integral transform, which can be used to simplify Eq. 3 to an ordinary differential equation. Due to the extrapolated boundary condition this transform is used in the following form

Φ (s_{n}, ϕ, m) = \int_{0}^{2 π} \int_{0}^{a'} ρΦ (ρ, ϕ') J_{m} (s_{n} ρ) \cos (m (ϕ - ϕ')) d ρ d ϕ',

(4)

where s_n are determined from the following equation

J_{m} (a' s_{n}) = 0, n = 1,2, \dots,

(5)

a′ = a + z_bk , and J_m is the first kind Bessel function of order m. By applying Eq. 4 to Eq. 3 the resulting ordinary differential equation for Φ = Φ(s_n , ν, m, z, ω) is

\frac{\partial^{2}}{\partial z^{2}} Φ - (\frac{μ_{a}}{D} + s_{n}^{2} + i \frac{ω}{Dc}) Φ = - \frac{1}{D} J_{m} (s_{n} ρ_{0}) \cos (m (ϕ - ϕ_{0})) δ (z - z_{0}) .

(6)

The solution of Eq. 6 is obtained by solving

\frac{\partial^{2}}{\partial z^{2}} G (s_{n}, z, ω) - α_{k}^{2} G (s_{n}, z, ω) = - \frac{1}{D} δ (z - z_{0}),

(7)

where

α = \sqrt{μ_{a} / D + s_{n}^{2} + iω / Dc}

, and by using

Φ (s_{n}, ϕ, m, z, ω) = G (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) \cos (m (ϕ - ϕ_{0})) .

(8)

For the first layer the solution of Eq. 7 is composed of the homogeneous and particular part, for the other layers only by the homogeneous part. For deriving the particular solution of the first layer G ₁ ^(p)(s_n , z, ω) we use the Fourier transform to find the one-dimensional Green’s function for an unbounded region in z-direction. By applying

G_{1}^{(p)} (s_{n}, k, ω) = \int_{- \infty}^{\infty} G_{1}^{(p)} (s_{n}, z, ω) e^{- ikz} d z

(9)

to Eq. 7 we get the Green’s function in the Fourier-space as

G_{1}^{(p)} (s_{n}, k, ω) = \frac{1}{D_{1}} \frac{1}{k^{2} + α_{1}^{2}} e^{- ik z_{0}} .

(10)

The inverse Fourier transform

G_{1}^{(p)} (s_{n}, z, ω) = \frac{1}{2 π} \int_{- \infty}^{\infty} G_{1}^{(p)} (s_{n}, k, ω) e^{ikz} d k

(11)

is used to obtain the particular solution of Eq. 7 for k = 1

G_{1}^{(p)} (s_{n}, z, ω) = \frac{e^{- α_{1} ∣ z - z_{0} ∣}}{2 D_{1} α_{1}} .

(12)

Thus, the complete solution of Eq. 7 is

G_{1} (s_{n}, z, ω) = A_{1} e^{α_{1} z} + B_{1} e^{- α_{1} z} + \frac{e^{- α_{1} ∣ z - z_{0} ∣}}{2 D_{1} α_{1}}, k = 1

G_{k} (s_{n}, z, ω) = A_{k} e^{α_{k} z} + B_{k} e^{- α_{k} z}, k = 2, \dots, N .

(13)

The constant coefficients A_k und B_k are determined by using the following boundary conditions in z-direction

G_{1} (s_{n}, z = - z_{b 1}, ω) = 0

n_{k + 1}^{2} G_{k} (s_{n}, z = L_{k}, ω) = n_{k}^{2} G_{k + 1} (s_{n}, z = L_{k}, ω)

D_{k} \partial G_{k} (s_{n}, z = L_{k}, ω) / \partial z = D_{k + 1} \partial G_{k + 1} (s_{n}, z = L_{k}, ω) / \partial z

G_{N} (s_{nN}, z = L_{N} + z_{bN}, ω) = 0,

(14)

where

L_{k} = \sum_{i = 1}^{k} l_{i}, k = 1, \dots, N

. The solution for the fluence rate in the first layer (0 ≤ z < l ₁) is given by

G_{1} (s_{n}, z, ω) = \frac{e^{- α_{1} ∣ z - z_{0} ∣} - e^{- α_{1} (z + z_{0} + 2 z b_{1})}}{2 D_{1} α_{1}} + \frac{\sinh [α_{1} (z_{0} + z_{b 1})] \sinh [α_{1} (z + z_{b 1})]}{D_{1} α_{1} e^{α_{1} (l_{1} + z_{b 1})}}

\times \frac{D_{1} α_{1} n_{1}^{2} β_{3} - D_{2} α_{2} n_{2}^{2} γ_{3}}{D_{1} α_{1} n_{1}^{2} β_{3} \cosh [α_{1} (l_{1} + z_{b 1})] + D_{2} α_{2} n_{2}^{2} γ_{3} \sinh [α_{1} (l_{1} + z_{b 1})]}

(15)

and for the fluence rate in the N^th layer (L _N-1 ≤ z < L_N ) by

G_{N} (s_{n}, z, ω) = \frac{n_{N}^{2} \prod_{i = 2}^{N - 1} (D_{i} α_{i} n_{i}^{2}) \sinh [α_{1} (z_{0} + z_{b 1})] \sinh [α_{N} (L_{N} + z_{bN} - z)]}{D_{1} α_{1} n_{1}^{2} β_{3} \cosh [α_{1} (l_{1} + z_{b 1})] + D_{2} α_{2} n_{2}^{2} γ_{3} \sinh [α_{1} (l_{1} + z_{b 1})]} .

(16)

In general, the quantities β ₃ and γ ₃ are obtained by recurrence relations. The start values are

β_{N} = D_{N - 1} α_{N - 1} n_{N - 1}^{2} \cosh (α_{N - 1} l_{N - 1}) \sinh [α_{N} (l_{N} + z_{b 2})]

+ D_{N} α_{N} n_{N}^{2} \sinh (α_{N - 1} l_{N - 1}) \cosh [α_{1} (l_{N} + z_{b 2})]

γ_{N} = D_{N - 1} α_{N - 1} n_{N - 1}^{2} \sinh (α_{N - 1} l_{N - 1}) \sinh [α_{N} (l_{N} + z_{b 2})]

+ D_{N} α_{N} n_{N}^{2} \cosh (α_{N - 1} l_{N - 1}) \cosh [α_{N} (l_{N} + z_{b 2})]

(17)

and the recurrence relations are given by

β_{k - 1} = D_{k - 2} α_{k - 2} n_{k - 2}^{2} \cosh (α_{k - 2} l_{k - 2}) β_{k} + D_{k - 1} α_{k - 1} n_{k - 1}^{2} \sinh (α_{k - 2} l_{k - 2}) γ_{k}

γ_{k - 1} = D_{k - 2} α_{k - 2} n_{k - 2}^{2} \sinh (α_{k - 2} l_{k - 2}) β_{k} + D_{k - 1} α_{k - 1} n_{k - 1}^{2} \cosh (α_{k - 2} l_{k - 2}) γ_{k} .

(18)

For a two-layered cylinder (N=2) Eqs. 17 and 18 are not needed, instead β ₃ = sinh[α ₂(l ₂+z _b2)] und γ ₃ = cosh[α ₂(l ₂+z _b2)] have to be applied. For a three-layered cylinder (N = 3) only Eq. 17 is required. We note that we derived the above presented solutions already for the case of a laterally infinitely extended medium [9

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt. , accepted. [PubMed]

], but now they are presented in a different way, so that no approximations are needed anymore for avoiding numerical errors as was necessary for the previous solutions. Thus, we recommend to use Eqs. 15–18 also for the laterally infinitely extended medium. The inverse relation of Eq. 4 is given by [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

]

Φ (ρ, ϕ) = \frac{1}{π {a'}^{2}} \sum_{m = - \infty}^{\infty} (2 - δ_{m, 0}) \sum_{n = 1}^{\infty} Φ (s_{n}, ϕ, m) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(19)

Finally, the solution of Eq. 3 in the frequency domain for layer k is obtained using Eq. 8 and the following recurrence relation

J_{m + 1} (a' s_{n}) = \frac{2 m}{a' s_{n}} J_{m} (a' s_{n}) - J_{m - 1} (a' s_{n}) = - J_{m - 1} (a' s_{n})

(20)

Φ_{k} (\vec{r}, ω) = \frac{1}{π {a'}^{2}} \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} G_{k} (s_{n}, z, ω) J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) k = 1,2 \dots N .

(21)

For a point source that is incident onto the center of the cylinder top Eq. 21 is simplified to

Φ_{k} (\vec{r}, ω) = \frac{1}{π {a'}^{2}} \sum_{n = 1}^{\infty} G_{k} (s_{n}, z, ω) J_{0} (s_{n} ρ) J_{1}^{- 2} (a' s_{n}) .

(22)

For the derivation of the solutions of the N-layered diffusion equation for a circular flat beam with radius ρ_w incident onto the middle of the cylinder top we refer to the accompanying paper [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

] where we gave the corresponding solution for a homogeneous cylinder. Accordingly, the solution of an N-layered cylinder is

Φ_{k} (\vec{r}, ω) = \frac{2}{π {a'}^{2} ρ_{w}} \sum_{n = 1}^{\infty} G_{k} (s_{n}, z, ω) \frac{J_{0} (s_{n} ρ) J_{1} (s_{n} ρ_{w})}{s_{n} J_{1}^{2} (s_{n} a')} .

(23)

We note that the solutions in the steady-state domain can be easily obtained from the above presented solutions in the frequency domain by setting ω= 0.

2.1.2. Solution in the time domain

In order to obtain the solutions of the diffusion equation in the time domain

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

(24)

we use the inverse Fourier transform

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

(25)

Thus, we can apply the solutions derived in the frequency domain to obtain the solutions in the time domain. For example, for a point source incident at an arbitrary position onto the first layer of the N-layered cylinder we obtain by inserting Eq. 21 into Eq. 25

Φ_{k} (\vec{r}, t) = \frac{1}{2 π^{2} {a'}^{2}} \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} \int_{- \infty}^{\infty} G_{k} (s_{n}, z, ω) e^{iωt} dω J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (a' s_{n}) .

(26)

The inverse Fourier integral is calculated by using the inverse discrete Fourier transform. For equal reduced scattering coefficients and refractive indices in all N layers faster solutions can be derived in the time domain. For this case, the complex numbers α_k in the frequency domain have the form

α_{k} = \sqrt{s_{n}^{2} + \frac{μ_{ak}}{D} + i \frac{ω}{Dc}} = \frac{\sqrt{μ_{ak} c + Dc s_{n}^{2} + iω}}{\sqrt{Dc}} .

(27)

The definition of the Laplace transform with the complex frequency s = σ+iω is

X (σ + iω) = \int_{- \infty}^{\infty} x (t) e^{- (σ + iω) t} dt, ROC : Re {s} > - σ,

(28)

where ROC is the region of convergence. By choosing σ=μ_akc+Dcs ² _n one obtains

X (μ_{ak} c + Dc s_{n}^{2} + iω) = \int_{- \infty}^{\infty} x (t) e^{- (μ_{ak} c + Dc s_{n}^{2} + iω) t} dt = \int_{- \infty}^{\infty} x (t) e^{- Dc s_{n}^{2} t} e^{- (μ_{ak} c + iω) t} dt .

(29)

Now it can be seen that Eq. 29 describes the Laplace transform of the function x(t)e ^{-Dcs_n

²
t}. The region of convergence for this case is Re{s} > -μ_akc. The imaginary axis is within the region of convergence, so that we can calculate the inverse transform of X(μ_ac+iω) by means of the Fourier transform, in particular the FFT algorithm. For this special case we obtain the solution as

Φ_{k} (\vec{r}, t) = \frac{1}{2 π^{2} {a'}^{2}} \int_{- \infty}^{\infty} G_{k} (z, ω) e^{iωt} dω \sum_{m = - \infty}^{\infty} \cos (mφ) \sum_{n = 1}^{\infty} e^{- Dc s_{n}^{2} t} J_{m} (s_{n} ρ_{0}) J_{m} (s_{n} ρ) J_{m + 1}^{- 2} (s_{n} a') .

(30)

Finally, we give the results for a point source incident onto the center of the cylinder top

Φ_{k} (\vec{r}, t) = \frac{1}{2 π^{2} {a'}^{2}} \int_{- \infty}^{\infty} G_{k} (z, ω) e^{iωt} dω \sum_{n = 1}^{\infty} e^{- Dc s_{n}^{2} t} J_{0} (s_{n} ρ) J_{1}^{- 2} (s_{n} a') .

(31)

We note that the solutions for equal reduced scattering coefficients and refractive indices in each layer can be quickly evaluated (≈ 10ms) using the above equations.

The derived solutions of the diffusion equation give the fluence rate in the steady-state, frequency, and time domains. In order to compare this to Monte Carlo simulations or experiments the reflectance and the transmittance have to be calculated. For the reflectance in the steady-state domain the fluence rate and its derivative in z-direction (flux-term) are used, whereas for all other solutions only the flux term is applied [9

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt. , accepted. [PubMed]

, 10

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. , accepted. [PubMed]

, 1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

2.2. Monte Carlo simulations

Monte Carlo simulations are used for the validation of the derived solutions. The basics of our Monte Carlo program were described in the accompanying paper [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

] for the calculation of the light propagation in a homogeneous cylinder. For the simulations of the layered cylinder we, additionally, consider the different optical properties of the involved layers using, e.g., Fresnel equations at the boundaries. 10⁷ photons were used in Monte Carlo simulations.

3. Results

The derived solutions of the diffusion equation for an N-layered cylinder are first compared to other available solutions of the diffusion equation that are valid for special cases. Then, the derived solutions are verified against results obtained by Monte Carlo simulations.

3.1. Comparison with other solutions of the diffusion theory

In literature so far a solution of the layered cylinder was derived for two and three layers [17

, 18

]. In Fig. 2 this solution is compared to Eq. 26 for the time resolved reflectance from a finite two-layered cylinder (solid curve). A point source is incident perpendicular onto the center of the cylinder top and the reflectance is calculated at a distance of ρ = 17mm.

Fig. 2. Comparison of the reflectance from a two-layered cylinder with radius a = 22mm calculated with Eq. 26 and by the eigenvalue method [17

,18

] (solid curves). The parameters are μ′_s1 = 1.1mm^-1, μ _a1 = 0.02mm^-1, μ′_s2 = 1.3mm^-1, μ _a2 = 0.008mm^-1, l ₁ = 7mm, l ₂ = 8mm, n ₁ = 1.4, and n ₂ = 1.4. For the non-scattering surrounding the refractive index is 1.0. The distance between the incident beam and detection is ρ = 17mm. The dashed curve shows the reflectance for a laterally infinitely extended medium.

The disagreement between the time resolved reflectance from the two-layered cylinder calculated with the two different methods is very small, so that we also plot the relative differences, see Fig. 3. These differences are in the range of 10^-5 and can normally be neglected for applications in the field of biophotonics. If necessary, they can be decreased by putting more effort into the numerical evaluation.

For comparison, also the reflectance from a laterally infinitely extended two-layered turbid medium having the same optical properties as the cylinder is shown in Fig. 2 (dashed curve). As expected, for small times the relative difference is small, whereas for larger times the differences increase due to the fact that the photons remitted at larger times have experienced larger volumes of the turbid media and, thus, have explored the differences of the two considered geometries.

In order to validate the derived solutions for an arbitrary number of layers solutions for different geometries have to be used, because no other solutions for an N-layered cylinder are known so far. Thus, we used the derived solutions of the N-layered turbid cylinder for an infinite large radius (a → ∞) and compared it to solutions derived for N-layered turbid media that are laterally infinitely extended [9

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt. , accepted. [PubMed]

, 10

A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. , accepted. [PubMed]

]. Fig. 4 shows the time resolved transmittance from a 7-layered turbid medium at three different detector distances from the middle of the cylinder bottom (ρ = 50, 60, 70mm). The point source is incident at the center of the cylinder top. The optical properties and thicknesses of the seven layers used for the calculations shown in Fig. 4 are listed in Table 1.

Tab. 1: Reduced scattering and absorption coefficients, refractive index, and thickness of the 7-layered turbid medium.

Fig. 3. Relative difference (|Martelli et al. - Eq. 26|/Martelli et al.) of the time resolved reflectance calculated with the two methods for the cylinder shown in Fig. 2.

Fig. 4. Comparison of the time resolved transmittance from a finite 7-layered turbid medium. The distances between the center of the cylinder bottom and the detection positions are ρ= 50mm (black curve), ρ= 60mm (green curve) and ρ= 70mm (red curve). The radius of the cylinder is a = 150mm.

layer	μ′_s[mm⁻¹]	μ_a [mm⁻¹]	n	l ₁[mm]
1	1.2	0.01	1.4	3
2	1.1	0.03	1.3	2
3	1.3	0.02	1.5	2
4	1.6	0.015	1.6	2
5	1.5	0.008	1.1	2
6	1.4	0.035	1.7	2
7	1.7	0.025	1.4	3

The agreement between the two solutions for all three distances is again excellent. Fig. 5 shows also the relative difference of the two solutions.

Fig. 5. Relative difference of the two solutions (|cylinder - laterally infinite|/cylinder) shown in Fig. 4.

For the most time values the relative differences are in the range of 10^-10 and 10^-8. Finally, we compare the derived solutions of the two-layered diffusion equation using the same optical properties in both layers to the solution of a homogeneous cylinder [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

]. The point source is incident at the center of the cylinder and time resolved transmittance from the cylinder barrel at a depth of z = 7mm is calculated, see Fig. 6. The results are given for two radii of the cylinder (a = 18mm (red curve), a = 40mm (green curve)). As in the other comparisons the relative difference between the curves calculated with the two models depend on the absolute reflectance data. In general, the larger the reflectance the smaller the relative differences.

The agreement is very good so that we also show the relative difference, see Fig. 7. The red curve shows the differences for a = 18mm and the green curve those for a = 40mm.

For comparison purpose the transmittance data obtained by a homogeneous parallelepiped model [21

A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]

] having the same optical properties as the cylinder and a length of 36mm or a = 80mm in x- and y-direction and a length of 10mm in z-direction are also shown in Fig. 6 (dashed curves). Similar as for the cylinder the parallelepiped is illuminated at the center and the transmittance is calculated at a depth of z = 7mm in the middle of a lateral side, so that the distance between the source and the detector is the same as for the cylinder. For small times the solutions are close, for larger times they again diverge due to the different geometries.

Fig. 6. Comparison of the time resolved transmission from the cylinder barrel of a homogeneous cylinder. The parameters are μ′ _s = 1.5mm^-1, μ_a = 0.02mm^-1, l ₁ = 4mm, l ₂ = 6mm, n = 1.4 and the two radii are a = 18mm (red curve) and a = 40mm (green curve). The dashed lines show the transmittance calculated with the solution for a homogeneous parallelepiped.

Fig. 7. Relative difference between time resolved transmittance from a two-layered and a homogeneous cylinder for a = 18mm (red curve) and a = 40mm (green curve) shown in Fig. 6.

3.2. Comparison with Monte Carlo simulations

In order to test the applicability of the derived solutions for evaluating experimental data we compare them to Monte Carlo simulations. First, the spatially resolved transmittance from a 7-layered cylinder which is illuminated at the center of the cylinder top is investigated. The optical and geometrical properties are those given in Tab. 1, but the refractive index is now n = 1.0 for all layers and the surrounding. The transmittance from the center of the cylinder bottom is calculated for three different radii a = 4mm (red curve), a = 5mm (green curve) and a = 10mm (blue curve). Fig. 8 shows the results obtained from the solution of the diffusion equation (solid curves) and those obtained from the Monte Carlo simulations (points).

Fig. 8. Comparison of the solution of the 7-layered diffusion equation (solid curves) with Monte Carlo simulations (circles) in the steady-state domain. The three radii are a = 4mm (red curve), a = 5mm (green curve), and a = 10mm (blue curve).

Good agreement is found for all three radii. This is expected because the distances between the incident source and the detector positions are large so that the diffusion approximation is valid.

Finally, the solution of the diffusion equation is compared to Monte Carlo simulations for a point beam which is non-centric incident onto the cylinder top at r⃗ ₀ = (x,y,z) = (4√2mm, 4√2mm, 1/(μ _a1 + μ′_s1)mm). Fig. 9 shows the spatially resolved reflectance along the x-axis (blue curve), along the negative y-axis (green curve), and along the bisecting line between the x-axis and y-axis (red line), see inset. The optical and geometrical properties are those used for the calculations shown in Fig. 8 besides the radius of the cylinder, which is changed to a = 12mm.

The spatially resolved reflectance along the bisecting line has a maximum at a distance of 8mm which corresponds to the location of the incident point source, where the diffusion equation fails to describe correctly the reflectance. At large distances, close to 12mm the reflectance from all cylinders rapidly decreases due to the influence of the cylinder border. In general, the agreement between diffusion theory and Monte Carlo simulations is in the expected range.

Fig. 9. Comparison of the reflectance from a 7-layered turbid medium calculated with the solution of the diffusion equation (solid curves) and with the Monte Carlo method (points).

4. Discussion

In this study the solutions of the diffusion equation for an N-layered cylinder were derived in the steady-state, frequency, and time domains applying formulae that we presented in the accompanying paper [1

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.

]. Solutions were obtained for a point source incident onto an arbitrary location of the first layer and a circular flat beam incident onto the middle of the cylinder top. In the time domain we derived fast solutions (≈ 10ms) for the special case of equal reduced scattering coefficient and refractive index in each layer. For the general case in the time domain the calculation time is about 1s when using a state-of-the-art computer and the Pascal (Delphi) programming language, whereas in the steady-state domain the calculation times are orders of magnitude faster.

A solution up to three layers was reported in the literature using a different approach [18

]. We compared this solution to the one derived in this paper for two-layers (and also three-layers, data not shown) and found excellent agreement. We also compared the derived solutions for the special cases of infinitely large radii of a 7-layered turbid medium and of homogeneous cylinders and observed very good accordance.

For the validation against transport theory we compared the derived solutions with Monte Carlo simulations. Calculations in transmission and remission were performed showing the expected small relative differences in the range of a few percent.

In general, although the derivation of the above presented solutions of the diffusion equation might be somewhat cumbersome, the obtained final results can be easily implemented. In addition, we provide executable files of the solutions in the steady-state and time domains on the internet [22

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

In the future, the solutions will be used to examine the inverse problem, the determination of the optical properties (and layer thickness) of (some of) the layers from measurements in the different domains.

Acknowledgement

We acknowledge the support by the European Union (nEUROPt, grant agreement no. 201076). We like to thank Fabrizio Martelli for providing us his code to solve the two-layered diffusion equation used for the comparison in Fig. 2.

References and links

1.	A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder. I Homogeneous Case,” submitted.
2.	A. Ishimaru, “Wave Propagation and Scattering in Random Media,” Academic Press, New York (1978).
3.	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). [CrossRef]
4.	A. Kienle, M.S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]
5.	A. Kienle, T. Glanzmann, G. Wagnières, and H. van den Bergh, “Investigation of Two-Layered Turbid Media with Time-Resolved Reflectance,” Appl. Opt. 37, 6852–6862 (1998). [CrossRef]
6.	A. Kienle and T. Glanzmann, “In Vivo Determination of the Optical Properties of Muscle Using a Layered-Model,” Phys. Med. Biol. 44, 2689–2702 (1999). [CrossRef] [PubMed]
7.	J.M. Tualle, H.M. Nghiem, D. Ettori, R. Sablong, E. Tinet, and S. Avrillier, “Asymptotic Behavior and Inverse Problem in Layered Scattering Media,” J. Opt. Soc. Am. A 21, 24–34 (2004). [CrossRef]
8.	X.C. Wang and S.M. Wang, “Light Transport Modell in a N-Layered Mismatched Tissue,” Waves Rand. Compl. Media 16, 121–135 (2006). [CrossRef]
9.	A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt. , accepted. [PubMed]
10.	A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt. , accepted. [PubMed]
11.	G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium,” Appl. Opt. 37, 7401–7409 (1998). [CrossRef]
12.	S-H. Tseng, C. Hayakawa, B.J. Tromberg, J. Spanier, and A.J. Durkin, “Quantitative Spectroscopy of Superficial Turbid Media,” Opt. Lett. 23, 3165–3167 (2005). [CrossRef]
13.	G. Alexandrakis, T.J. Farrell, and M.S. Patterson, “Monte Carlo Diffusion Hybrid Model for Photon Migration in a Two-Layer Turbid Medium in the Frequency Domain,” Appl. Opt. 39, 2235–2244 (2000). [CrossRef]
14.	M. Das, C. Xu, and Q. Zhu, “Analytical Solution for Light Propagation in a Two-Layer Tissue Structure with a Tilted Interface for Breast Imaging,” Appl. Opt. 45, 5027–5036 (2006). [CrossRef] [PubMed]
15.	A.H. Barnett, “A Fast Numerical Method for Time-Resolved Photon Diffusion in General Stratified Turbid Media,” J. Comp. Phys. 201, 771–797 (2004). [CrossRef]
16.	C. Donner and H.W. Jensen, “Rapid Simulations of Steady-State Spatially Resolved Reflectance and Transmittance Profiles of Multilayered Turbid Materials,” J. Opt. Soc. Am. A 23, 1382–1390 (2006). [CrossRef]
17.	F. Martelli, A. Sassaroli, S. Del Bianco, Y. Yamada Y, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for Layered Diffusive Media by the Eigenfunction Method,” Phys. Rev. E 67, 056623 (2003). [CrossRef]
18.	F. Martelli, A. Sassaroli, S. Del Bianco, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for a Three-Layer Medium: Application to Study Photon Migration through a Simplified Adult Head Model,” Phys. Med. Biol. 52, 2827–2843 (2007). [CrossRef] [PubMed]
19.	F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation Model for Light Propagation through Diffusive Layered Media,” Phys. Med. Biol. 50, 2159–2166 (2005). [CrossRef] [PubMed]
20.	F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, “Light Propagation through Biological Tissue and other Diffusive Media,” SPIE Press, Bellingham (2010).
21.	A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
22.	http://www.uni-ulm.de/ilm/index.php?id=10020200.

OCIS Codes
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.5280) Medical optics and biotechnology : Photon migration

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: February 12, 2010
Revised Manuscript: April 7, 2010
Manuscript Accepted: April 14, 2010
Published: April 19, 2010

Virtual Issues
Vol. 5, Iss. 9 Virtual Journal for Biomedical Optics

Citation
André Liemert and Alwin Kienle, "Light diffusion in a turbid cylinder. II. Layered case," Opt. Express 18, 9266-9279 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-9-9266

Sort: Author | Year | Journal | Reset

References

A. Liemert and A. Kienle, “Light Diffusion in a N-layered Turbid Cylinder.I Homogeneous Case,” submitted.
A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, New York, 1978).
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). [CrossRef]
A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagniàres, and H. van den Bergh, “Noninvasive Determination of the Optical Properties of Two-Layered Turbid Media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]
A. Kienle, T. Glanzmann, G. Wagnières, and H. van den Bergh, “Investigation of Two-Layered Turbid Media with Time-Resolved Reflectance,” Appl. Opt. 37, 6852–6862 (1998). [CrossRef]
A. Kienle and T. Glanzmann, “In Vivo Determination of the Optical Properties of Muscle Using a Layered-Model,” Phys. Med. Biol. 44, 2689–2702 (1999). [CrossRef] [PubMed]
J. M. Tualle, H. M. Nghiem, D. Ettori, R. Sablong, E. Tinet, and S. Avrillier, “Asymptotic Behavior and Inverse Problem in Layered Scattering Media,” J. Opt. Soc. Am. A 21, 24–34 (2004). [CrossRef]
X. C. Wang and S. M. Wang, “Light Transport Modell in a N-Layered Mismatched Tissue,”Waves Rand. Compl. Media 16, 121–135 (2006). [CrossRef]
A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Steady-State Domain,” J. Biomed. Opt.accepted. [PubMed]
A. Liemert and A. Kienle, “Light Diffusion in N-layered Turbid Media: Frequency and Time Domains,” J. Biomed. Opt.accepted. [PubMed]
G. Alexandrakis, T. J. Farrell, and M. S. Patterson, “Accuracy of the Diffusion Approximation in Determining the Optical Properties of a Two-Layer Turbid Medium,” Appl. Opt. 37, 7401–7409 (1998). [CrossRef]
S-H. Tseng, C. Hayakawa, B. J. Tromberg, J. Spanier, and A. J. Durkin, “Quantitative Spectroscopy of Superficial Turbid Media,” Opt. Lett. 23, 3165–3167 (2005). [CrossRef]
G. Alexandrakis, T. J. Farrell, and M. S. Patterson, “Monte Carlo Diffusion Hybrid Model for Photon Migration in a Two-Layer Turbid Medium in the Frequency Domain,” Appl. Opt. 39, 2235–2244 (2000). [CrossRef]
M. Das, C. Xu, and Q. Zhu, “Analytical Solution for Light Propagation in a Two-Layer Tissue Structure with a Tilted Interface for Breast Imaging,” Appl. Opt. 45, 5027–5036 (2006). [CrossRef] [PubMed]
A. H. Barnett, “A Fast Numerical Method for Time-Resolved Photon Diffusion in General Stratified Turbid Media,” J. Comp. Phys. 201, 771–797 (2004). [CrossRef]
C. Donner and H. W. Jensen, “Rapid Simulations of Steady-State Spatially Resolved Reflectance and Transmittance Profiles of Multilayered Turbid Materials,” J. Opt. Soc. Am. A 23, 1382–1390 (2006). [CrossRef]
F. Martelli, A. Sassaroli, S. Del Bianco, Y. Yamada Y, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for Layered Diffusive Media by the Eigenfunction Method,” Phys. Rev. E 67, 056623 (2003). [CrossRef]
F. Martelli, A. Sassaroli, S. Del Bianco, and G. Zaccanti, “Solution of the Time-Dependent Diffusion Equation for a Three-Layer Medium: Application to Study Photon Migration through a Simplified Adult Head Model,” Phys. Med. Biol. 52, 2827–2843 (2007). [CrossRef] [PubMed]
F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation Model for Light Propagation through Diffusive Layered Media,” Phys. Med. Biol. 50, 2159–2166 (2005). [CrossRef] [PubMed]
F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, Light Propagation through Biological Tissue and other Diffusive Media (SPIE Press, Bellingham, 2010).
A. Kienle, “Light Diffusion Through a Turbid Parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
http://www.uni-ulm.de/ilm/index.php?id=10020200.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper