Hybrid diffusion-P3 equation in N-layered turbid media: steady-state domain

Zhenzhi Shi; Huijuan Zhao; Kexin Xu

doi:10.1117/1.3640810

1 October 2011 Hybrid diffusion-P3 equation in N-layered turbid media: steady-state domain

Zhenzhi Shi, Huijuan Zhao, Kexin Xu

Author Affiliations +

Journal of Biomedical Optics, Vol. 16, Issue 10, 105002 (October 2011). https://doi.org/10.1117/1.3640810

Abstract

This paper discusses light propagation in N-layered turbid media. The hybrid diffusion-P3 equation is solved for an N-layered finite or infinite turbid medium in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in Fourier space. Two inverse Fourier transform methods are developed to calculate the fluence rate in real space. In addition, the solutions of the hybrid diffusion-P3 equation are compared to the solutions of the diffusion equation and the Monte Carlo simulation. For the case of small absorption coefficients, the solutions of the N-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent and are in agreement with the Monte Carlo simulation. For the case of large absorption coefficients, the model of the hybrid diffusion-P3 equation is more precise than that of the diffusion equation. In conclusion, the model of the hybrid diffusion-P3 equation can replace the diffusion equation for modeling light propagation in the N-layered turbid media for a wide range of absorption coefficients.

1. Introduction

The in vivo determination of tissue optical properties and the study of light propagation in biological tissues are very important in a variety of biomedical fields that can be used to obtain knowledge of the physiological state of tissues.^{1, 2, 3, 4} Diffuse reflectance spectroscopy has emerged as one of the noninvasive optical techniques and it has been researched to obtain quantitative tissue characterization and optical properties.^{4, 5, 6, 7, 8, 9} In recent years, the study of light propagation in layered turbid media has gathered much attention because many parts of the body such as the skin, stomach, and head have layered tissue structures that require corresponding forward solutions.

Usually, the radiative transfer equation (RTE) is considered to be the most precise equation for describing light propagation in biological tissues and several solutions for the RTE have been reported in previous literature. However, the direct analytical solutions of the RTE cannot be easily obtained. Instead, the discrete ordinates (S _N) and the spherical harmonics equations (P _N) have been established to solve the RTE. The diffusion equation, which is the P ₁ approximation of RTE, has been widely used to obtain the solution for light propagation in layered turbid media and several solutions have been reported in literature. Dayan used the Fourier and Laplace transforms to acquire the solutions for the two-layered diffusion equation,¹⁰ while Kienle acquired the solutions for the two-layered diffusion equation in the steady-state, frequency, and time domains using the Fourier transforms.^{11, 12} Martelli presented the solutions for the case of two and three layers using the Eigenfunction method and the perturbation model.^{13, 14, 15, 16} For the case of the N-layered diffusion equation, several groups have also reported their solutions.^{17, 18, 19, 20, 21, 22} Liemert presented the solutions for N-layered infinite or finite turbid media in the steady-state, frequency, and time domains using the two-dimensional (2-D) Fourier transform formalism^{19, 20} and also presented the solutions for a cylindrical geometry.^{21, 22} Liemert's group gave us a complete set of solutions for the N-layered diffusion equation, and they also studied the correctly and efficiently numerical 2-D inverse Fourier transform and provided the executable program on the Internet using the Delphi (Pascal) language.

However, the diffusion equation has several limitations. First, the diffusion equation requires the optical absorption coefficient μ_a to be much smaller than the scattering coefficient μ_s. A typical criterion for the applicable region of the diffusion equation is that the reduced albedo a′ must be larger than 0.98.²³ Second, the diffusion equation is valid only when the radial distance is larger than 10 transport mean free paths.²⁴ Therefore, the diffusion equation can be satisfied in the therapeutic window (λ = 650 − 950 nm) of most tissues with small absorption coefficients. For example, the absorption coefficient of gray and white matter in the adult head at a wavelength of 800 nm is 0.025 and 0.005 mm⁻¹ and the reduced scattering coefficient is 2.5 and 6 mm⁻¹, respectively.²⁵ The absorption coefficient of human skin in vitro in the therapeutic window is less than 0.05 mm⁻¹ and the reduced scattering coefficient is almost 1.5 mm⁻¹.²⁶ However, in the near-IR (λ > 1000 nm) region of the spectrum, the absorption coefficients of tissues become large and the reduced albedo can be 0.5 or smaller. For example, the absorption coefficient of the dermis of the skin when λ > 1400 nm, is larger than 0.5 mm⁻¹ and the reduced scattering coefficient is less than 1.5 mm⁻¹.²⁷ Evidently, there is a significant need for methods that accurately describe light propagation in media with large absorption coefficients, especially in the application of noninvasive blood glucose monitoring using NIR light in the region of 1000 to 2000 nm because of the high absorption of water and skin tissues.²⁸ Hull and Foster²⁹ derived Green's function of the steady-state RTE in the P ₃ approximation for the infinite turbid media and demonstrated that the P ₃ approximation models the radiance in highly absorbing media or in the region close to the source with an accuracy that is superior to that of the diffusion equation. In order to develop the diffusion equation and to simplify the P ₃ approximation to be used for highly absorbing media, Hull and Foster investigated a modified diffusion approximation model called the hybrid diffusion-P3 equation.²⁹ Klose also developed a method of the simplified spherical harmonics equations (SP _N) to approximate the more complicated equation of radiative transfer for modeling light propagation in biological tissues and found that the SP _N significantly improves the diffusion solution in transport-like domains with high absorption and small geometries.³⁰

In this study, we present the solution for the hybrid diffusion-P3 equation for N-layered turbid media by combining the derivation of the N-layered diffusion equation in Liemert and Kienle's paper²⁰ and the modified diffusion approximation model investigated by Hull and Foster²⁹ in order to find a solution for light propagation in N-layered turbid media for the case of large absorption coefficients. The hybrid diffusion-P3 equation is solved for an N-layered finite or infinite turbid medium in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in the Fourier space. The 2-D inverse Fourier transform is numerically calculated. In addition, the solutions of hybrid diffusion-P3 equation are compared to those of the diffusion equation and Monte Carlo simulations.

2. Theory

In this section, we use the Fourier transform formalism to derive the steady-state hybrid diffusion-P3 equation for an N-layered turbid medium having finite or infinite extensions.

2.1.

N-layered Diffusion Equation

First, the N-layered diffusion equation can be described as follows,²⁰ when the turbid media in every layer is homogeneous:

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} D_1 \Delta \Phi _1 (x,y,z) - \mu _{a1} \Phi _1 (x,y,z) = - q(x,y,z),\;0 \le z < l_1,\hspace*{-6pt}\nonumber\\ \end{eqnarray}\end{document}

\begin{matrix} D_{1} Δ Φ_{1} (x, y, z) - μ_{a 1} Φ_{1} (x, y, z) = - q (x, y, z), 0 \leq z < l_{1}, \end{matrix}

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} D_k \Delta \Phi _k (x,y,z) - \mu _{ak} \Phi _k (x,y,z) &=& 0,\;\sum\limits_{j = 1}^{k - 1} {l_j } < z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document}

\begin{matrix} D_{k} Δ Φ_{k} (x, y, z) - μ_{a k} Φ_{k} (x, y, z) & = & 0, \sum_{j = 1}^{k - 1} l_{j} < z \leq \sum_{j = 1}^{k} l_{j}, \\ k & = & 2, 3, ..., N, \end{matrix}

where Φ_k is the fluence rate of the layer k, l _k is the thickness, and [TeX:] $\mu _{sk}^\prime $

μ_{s k}^{'}

and μ_ak are the reduced scattering and absorption coefficients, respectively. [TeX:] $D_k = 1/[3(\mu _{sk}^\prime + \mu _{ak})]$

D_{k} = 1 / [3 (μ_{s k}^{'} + μ_{a k})]

is the diffusion coefficient. q is the source term.

For a semi-infinite situation, the derivation of the N-layered diffusion equation involves characterizing the source term and satisfying the appropriate boundary condition. For the light propagation in biological tissues, a pencil beam is usually modeled as an infinite line of isotropic sources on a semi-infinite scattering medium. According to the formalism developed by Farrell,³¹ we can obtain the distribution equation of the source

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q(z) = \frac{a^{\prime}\mu_t^{\prime}}{4\pi }\exp (- \mu_t^{\prime}z), \end{equation}\end{document}

q (z) = \frac{a^{'} μ_{t}^{'}}{4 π} \exp (- μ_{t}^{'} z),

where [TeX:] $\mu _t^\prime = \mu _a + \mu _s^\prime $

μ_{t}^{'} = μ_{a} + μ_{s}^{'}

, and [TeX:] $a^\prime = \mu _s^\prime /(\mu _s^\prime + \mu _a)$

a^{'} = μ_{s}^{'} / (μ_{s}^{'} + μ_{a})

is the reduced albedo.

In order to represent the pencil beam in terms of simpler source distributions, isotropic point sources are adopted to describe the distributions of the sources. The distributions of isotropic point sources have a dipole moment with respect to an origin at the air–tissue interface as the distribution in Eq. 3.²⁹ To satisfy the dipole moment, a single point source is needed

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \int_0^\infty {za'\mu _t^{\prime} } \exp (- \mu _t^{\prime} z)dz = \int_0^\infty {za'} \delta (z - z_0)dz. \end{equation}\end{document}

\int_{0}^{\infty} z a^{'} μ_{t}^{'} \exp (- μ_{t}^{'} z) d z = \int_{0}^{\infty} z a^{'} δ (z - z_{0}) d z .

The right side of Eq. 4 indicates that an infinite line of isotropic sources can be modeled as a point source at an effective source depth [TeX:] $z_0 = 1/(\mu _s^\prime + \mu _a)$

z_{0} = 1 / (μ_{s}^{'} + μ_{a})

. Thus, the source item q(x, y, z)in Eq. 1 can be expressed as

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} q(x,y,z) = a_1 ^{\prime} \delta (x,y,z - z_0), \end{equation}\end{document}

q (x, y, z) = a_{1}^{'} δ (x, y, z - z_{0}),

where [TeX:] $z_0 = 1/(\mu _{a1} + \mu _{s1}^\prime)$

z_{0} = 1 / (μ_{a 1} + μ_{s 1}^{'})

must locate in the first layer.

To obtain the solutions of the radiance emitted from a semi-infinite medium, appropriate boundary conditions must be prescribed at the interface between the surrounding media and the biological tissue. The extrapolated boundary condition is one of the boundary conditions used for a semi-infinite scattering medium,^{32, 33} shown in Fig. 1, where Φ(x, y, z = −z _b) = 0 and z _b is the position of the extrapolated boundary.

Fig. 1

Scheme of the N-layered turbid media of one point source and the extrapolated boundary condition.

We then used the Fourier transform approach to solve the N-layered diffusion equation based on Liemert and Kienle's paper.²⁰ We apply the 2-D Fourier transform

Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _k (z,s_1,s_2) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\Phi _k (z,x,y)} } e^{i(s_1 x + s_2 y)} dxdy \end{equation}\end{document}

Φ_{k} (z, s_{1}, s_{2}) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} Φ_{k} (z, x, y) e^{i (s_{1} x + s_{2} y)} d x d y

to Eqs. 1, 2, and obtained the N-layered diffusion equation in Fourier space

Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _1 (z,s) - a_1^2 \Phi _1 (z,s) = - \frac{1}{{D_1 }}a_1 ^{\prime} \delta (z - z_0),\quad 0 \le z < l_1, \end{equation}\end{document}

\frac{\partial^{2}}{\partial z^{2}} Φ_{1} (z, s) - a_{1}^{2} Φ_{1} (z, s) = - \frac{1}{D_{1}} a_{1}^{'} δ (z - z_{0}), 0 \leq z < l_{1},

Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _k (z,s) - a_k^2 \Phi _k (z,s) &=& 0,\quad \sum\limits_{j = 1}^{k - 1} {l_j } < z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document}

\begin{matrix} \frac{\partial^{2}}{\partial z^{2}} Φ_{k} (z, s) - a_{k}^{2} Φ_{k} (z, s) & = & 0, \sum_{j = 1}^{k - 1} l_{j} < z \leq \sum_{j = 1}^{k} l_{j}, \\ k & = & 2, 3, ..., N, \end{matrix}

with [TeX:] $s = \sqrt {s_1^2 + s_2^2 } $

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

and [TeX:] $a_k^2 = (\mu _{ak} + D_k s^2)/D_k $

a_{k}^{2} = (μ_{a k} + D_{k} s^{2}) / D_{k}

.

2.2.

N-layered Hybrid Diffusion-P3 Equation

The N-layered hybrid diffusion-P3 equation is the method based on N-layered diffusion equation, combined with the coefficients of the P ₃ approximation, which is the high-order approximation of the solution of RTE. It has been demonstrated that the P ₃ approximation for modeling the radiance in highly absorbing media or in the region close to the source is more accurate than the diffusion equation.²⁹ v ⁻ is the high-order coefficient of the P3-Green function, which is named as the asymptotic attenuation coefficient and approximately equal to the effective attenuation coefficient [TeX:] $\mu _{\hbox{\scriptsize\textit{eff}}} $ $μ_{eff}$ of the diffusion-Green function for large albedos and increasingly deviates from [TeX:] $\mu _{\hbox{\scriptsize\textit{eff}}} $ $μ_{eff}$ as the albedo decreases.²⁹ In addition, with the use of [TeX:] $D = \mu _a /(\mu _{\hbox{\scriptsize\textit{eff}}})^2 $ $D = μ_{a} / {(μ_{eff})}^{2}$ , it has also been demonstrated that the hybrid diffusion-P3 approximation for semi-infinite media is more accurate than diffusion approximation for the case of large absorption coefficients.³⁴ So we can obtain that v ⁻ associates with [TeX:] $\mu _{\hbox{\scriptsize\textit{eff}}} $ $μ_{eff}$ of diffusion-Green function and [TeX:] ${{\mu _{ak} } / {D_k }}$ $μ_{a k} / D_{k}$ and D _k of diffusion equation associate with [TeX:] $(v_{^k }^ -)^2 $ ${(v_{k}^{-})}^{2}$ and [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ of P ₃ approximation, respectively. [TeX:] $(v_{^k }^ -)^2 $ ${(v_{k}^{-})}^{2}$ and [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ are the high-order coefficients compared to [TeX:] ${{\mu _{ak} } / {D_k }}$ $μ_{a k} / D_{k}$ and D _k of the diffusion equation.

Thus, the N-layered hybrid diffusion-P3 equation is obtained based on the N-layered diffusion equation by replacing [TeX:] ${{\mu _{ak} } / {D_k }}$ $μ_{a k} / D_{k}$ and D _k with [TeX:] $(v_{^k }^ -)^2 $ ${(v_{k}^{-})}^{2}$ and [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ , respectively.

Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &\displaystyle v_{_k }^ - = \frac{1}{{\sqrt {18} }}\left(\beta _k - \sqrt {\beta _k ^2 - \gamma _{ak} }\right)^{\frac{1}{2}} \nonumber\\ &\displaystyle \beta _k \equiv 27\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime) + 28\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \delta)\nonumber\\ &\displaystyle \hspace*{22pt}+\, 35\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \gamma _k)\mu _{ak} (\mu _{ak} + \mu _{sk} ^\prime \delta _k),\\ &\displaystyle \gamma _{ak} \equiv 3780(\mu _{ak} + \mu _{sk} ^\prime)(\mu _{ak} + \mu _{sk} ^\prime \gamma _k)(\mu _{ak} + \mu _{sk} ^\prime \delta _k), \nonumber\\ &\displaystyle \gamma _k = (1 - g_{2k})/(1 - g_{1k}),\nonumber\\ &\displaystyle\delta _k = (1 - g_{3k})/(1 - g_{1k}), \nonumber\\ &\displaystyle D_{asym\_k} = \mu _{ak} /(v_{_k }^ -)^2, \nonumber \end{eqnarray}\end{document}

\begin{matrix} v_{k}^{-} = \frac{1}{\sqrt{18}} {(β_{k} - \sqrt{β_{k}^{2} - γ_{a k}})}^{\frac{1}{2}} \\ β_{k} \equiv 27 μ_{a k} (μ_{a k} + μ_{s k}^{'}) + 28 μ_{a k} (μ_{a k} + μ_{s k}^{'} δ) \\ + 35 μ_{a k} (μ_{a k} + μ_{s k}^{'} γ_{k}) μ_{a k} (μ_{a k} + μ_{s k}^{'} δ_{k}), \\ γ_{a k} \equiv 3780 (μ_{a k} + μ_{s k}^{'}) (μ_{a k} + μ_{s k}^{'} γ_{k}) (μ_{a k} + μ_{s k}^{'} δ_{k}), \\ γ_{k} = (1 - g_{2 k}) / (1 - g_{1 k}), \\ δ_{k} = (1 - g_{3 k}) / (1 - g_{1 k}), \\ D_{a s y m_k} = μ_{a k} / {(v_{k}^{-})}^{2}, \end{matrix}

and g _1k, g _2k, and g _3k are the first-, second-, and third-moments of the phase function in the layer k, respectively. For the Henyey–Greenstein phase function, γ_k = 1 + g _1k and [TeX:] $\delta _k = 1 + g_{1k}\break + g_{1k} ^2 $

δ_{k} = 1 + g_{1 k} + g_{1 k}^{2}

.

Thus, we get the expression of the hybrid diffusion-P3 equation for N-layered turbid media in Fourier space

Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _{_1 }^H (z,s) - \big(a_1^H\big)^2 \Phi _{_1 }^H (z,s) &=& - \frac{1}{{D_{asym\_1} }}a_1 ^{\prime} \delta (z - z_0),\nonumber\\ && 0 \le z \le l_1, \end{eqnarray}\end{document}

\begin{matrix} \frac{\partial^{2}}{\partial z^{2}} Φ_{1}^{H} (z, s) - {(a_{1}^{H})}^{2} Φ_{1}^{H} (z, s) & = & - \frac{1}{D_{a s y m_1}} a_{1}^{'} δ (z - z_{0}), \\ 0 \leq z \leq l_{1}, \end{matrix}

Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \frac{{\partial ^2 }}{{\partial z^2 }}\Phi _{_k }^H (z,s) - \big(a_k^H\big)^2 \Phi _{_k }^H (z,s) &=& 0,\quad \sum\limits_{j = 1}^{k - 1} {l_j } \le z \le \sum\limits_{j = 1}^k {l_j },\nonumber\\ k &=& 2,3,\ldots,N, \end{eqnarray}\end{document}

\begin{matrix} \frac{\partial^{2}}{\partial z^{2}} Φ_{k}^{H} (z, s) - {(a_{k}^{H})}^{2} Φ_{k}^{H} (z, s) & = & 0, \sum_{j = 1}^{k - 1} l_{j} \leq z \leq \sum_{j = 1}^{k} l_{j}, \\ k & = & 2, 3, ..., N, \end{matrix}

with [TeX:] $s = \sqrt {s_1^2 + s_2^2 } $

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

and [TeX:] $(a_k^H)^2 = (v_{^k }^ -)^2 + s^2 $

{(a_{k}^{H})}^{2} = {(v_{k}^{-})}^{2} + s^{2}

. [TeX:] $\Phi _{_k }^H (z,s)$

Φ_{k}^{H} (z, s)

is the fluence rate in the Fourier space in the kth layer.

The following boundary conditions are used in the Fourier space for an finite N-layered turbid medium:²⁰

Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _1^H (- z_{b1},s) = 0, \end{equation}\end{document}

Φ_{1}^{H} (- z_{b 1}, s) = 0,

Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi_{1}^{{H}^{\prime}} (z_0,s) = \Phi_1^H (z_0,s),\;\;z_0 < l_1, \end{equation}\end{document}

Φ_{1}^{H^{'}} (z_{0}, s) = Φ_{1}^{H} (z_{0}, s), z_{0} < l_{1},

Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left. {\frac{{\partial \Phi _1^H (z,s)}}{{\partial z}}} \right|_{z = z_0 } - \left. {\frac{{\partial \Phi_{1}^{{H}^{\prime}} (z,s)}}{{\partial z}}} \right|_{z = z_0 } = \frac{{a_1 ^{\prime} }}{{D_{asym\_1} }}, \end{equation}\end{document}

{\frac{\partial Φ_{1}^{H} (z, s)}{\partial z}|}_{z = z_{0}} - {\frac{\partial Φ_{1}^{H^{'}} (z, s)}{\partial z}|}_{z = z_{0}} = \frac{a_{1}^{'}}{D_{a s y m_1}},

Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{{\Phi _k^H (L_k,s)}}{{\Phi _{k + 1}^H (L_k,s)}} = \left(\frac{{n_k }}{{n_{k + 1} }}\right)^2,\;L_k = \sum\limits_{j = 1}^k {l_j },\;1 \le k \le N - 1, \end{equation}\end{document}

\frac{Φ_{k}^{H} (L_{k}, s)}{Φ_{k + 1}^{H} (L_{k}, s)} = {(\frac{n_{k}}{n_{k + 1}})}^{2}, L_{k} = \sum_{j = 1}^{k} l_{j}, 1 \leq k \leq N - 1,

Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \left. {D_{asym\_k} \frac{{\partial \Phi _k^H (z,s)}}{{\partial z}}} \right|_{z = L_k } &=& \left. {D_{asym\_k + 1} \frac{{\partial \Phi _{k + 1}^H (z,s)}}{{\partial z}}} \right|_{z = L_k },\nonumber\\ L_k &=& \sum\limits_{j = 1}^k {l_j },\quad 1 \le k \le N - 1, \end{eqnarray}\end{document}

\begin{matrix} {D_{a s y m_k} \frac{\partial Φ_{k}^{H} (z, s)}{\partial z}|}_{z = L_{k}} & = & {D_{a s y m_k + 1} \frac{\partial Φ_{k + 1}^{H} (z, s)}{\partial z}|}_{z = L_{k}}, \\ L_{k} & = & \sum_{j = 1}^{k} l_{j}, 1 \leq k \leq N - 1, \end{matrix}

Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _N^H (L_N + z_{b2},s) = 0, \end{equation}\end{document}

Φ_{N}^{H} (L_{N} + z_{b 2}, s) = 0,

where z _b1 and z _b2 are the positions of the extrapolated boundaries above the first layer and below the last layer, respectively.

Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} z_{b1} = \frac{{1 + R_{{\hbox{\scriptsize\textit{eff}}}\_1} }}{{1 - R_{{\hbox{\scriptsize\textit{eff}}}\_1} }}2D_{asym\_1},\;z_{b2} = \frac{{1 + R_{{\hbox{\scriptsize\textit{eff}}}\_N} }}{{1 - R_{{\hbox{\scriptsize\textit{eff}}}\_N} }}2D_{asym\_N}, \end{equation}\end{document}

z_{b 1} = \frac{1 + R_{eff_1}}{1 - R_{eff_1}} 2 D_{a s y m_1}, z_{b 2} = \frac{1 + R_{eff_N}}{1 - R_{eff_N}} 2 D_{a s y m_N},

where [TeX:] $R_{{\hbox{\scriptsize\textit{eff}}}\_1} $

R_{eff_1}

and [TeX:] $R_{{\hbox{\scriptsize\textit{eff}}}\_N} $

R_{eff_N}

are the effective reflection coefficients at the top and bottom layers, respectively. We calculated [TeX:] $R_{{\hbox{\scriptsize\textit{eff}}}\_1} $

R_{eff_1}

and [TeX:] $R_{{\hbox{\scriptsize\textit{eff}}}\_N} $

R_{eff_N}

with the formula derived by Haskell ³³

2.3.

Solution of N-layered Hybrid Diffusion-P3 Equation

The solution of the N-layered hybrid diffusion-P3 equation [Eqs. 10, 11, 12, 13, 14, 15, 16, 17] for the fluence rate in the Fourier space is solved based on the derivation results from Liemert and Kienle's paper and by applying Cramer's rule.²⁰ The solution for the fluence rate in the first layer (above the single point source) is given by

Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _1^H (z,s) &=& a_1 ^{\prime} \left\{ \frac{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 - z)\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 - z)\big]}}{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 + z_{b1})\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 + z_{b1})\big]}}\right. \nonumber\\ && \left.\times\, \frac{{\sinh \big[a_1^H (z_0 + z_{b1})\big]}}{{a_1^H D_{asym\_1} }} - \frac{{\sinh \big[a_1^H (z_0 - z)\big]}}{{a_1^H D_{asym\_1} }}\right\}, \end{eqnarray}\end{document}

\begin{matrix} Φ_{1}^{H} (z, s) & = & a_{1}^{'} \{\frac{a_{1}^{H} β_{3} D_{a s y m_1} \cosh [a_{1}^{H} (l_{1} - z)] + a_{2}^{H} γ_{3} D_{a s y m_2} (n_{2}^{2} / n_{1}^{2}) \sinh [a_{1}^{H} (l_{1} - z)]}{a_{1}^{H} β_{3} D_{a s y m_1} \cosh [a_{1}^{H} (l_{1} + z_{b 1})] + a_{2}^{H} γ_{3} D_{a s y m_2} (n_{2}^{2} / n_{1}^{2}) \sinh [a_{1}^{H} (l_{1} + z_{b 1})]} \\ \times \frac{\sinh [a_{1}^{H} (z_{0} + z_{b 1})]}{a_{1}^{H} D_{a s y m_1}} - \frac{\sinh [a_{1}^{H} (z_{0} - z)]}{a_{1}^{H} D_{a s y m_1}}\}, \end{matrix}

and the fluence rate in the Nth layer (L _{N − 1} ⩽ z ⩽ L _N) can be given by

Eq. 20

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _N^H (z,s) = a_1 ^{\prime} \left\{ \frac{{\displaystyle\prod\nolimits_{k = 2}^{N - 1} {a_k^H D_{asym\_k} (n_N /n_1)^2 \sinh \big[a_1^H (z_0 + z_{b1})\big]\sinh \big[a_N^H (L_N + z_{b2} - z)\big]} }}{{a_1^H \beta _3 D_{asym\_1} \cosh \big[a_1^H (l_1 + z_{b1})\big] + a_2^H \gamma _3 D_{asym\_2} \big(n_2^2 /n_1^2\big)\sinh \big[a_1^H (l_1 + z_{b1})\big]}}\right\}. \end{equation}\end{document}

Φ_{N}^{H} (z, s) = a_{1}^{'} \{\frac{\prod_{k = 2}^{N - 1} a_{k}^{H} D_{a s y m_k} {(n_{N} / n_{1})}^{2} \sinh [a_{1}^{H} (z_{0} + z_{b 1})] \sinh [a_{N}^{H} (L_{N} + z_{b 2} - z)]}{a_{1}^{H} β_{3} D_{a s y m_1} \cosh [a_{1}^{H} (l_{1} + z_{b 1})] + a_{2}^{H} γ_{3} D_{a s y m_2} (n_{2}^{2} / n_{1}^{2}) \sinh [a_{1}^{H} (l_{1} + z_{b 1})]}\} .

In general, in the case of the N ⩾ 3 and finite turbid media, the quantities β₃ and γ₃ are obtained by using recursion formulas. The values of β₃ and γ₃ for finite N-layered turbid media are initiated by

Eq. 21

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \begin{array}{l} \left({\begin{array}{c} {\beta _N } \\ {\gamma _N } \\ \end{array}} \right) = a_{N - 1}^H D_{asym\_N - 1}^H \sinh \big[a_N^H (l_N + z_{b2})\big]\left[ {\begin{array}{c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right] \\[15pt] \qquad\qquad\; +\, a_N^H D_{asym\_N}^H \left(\displaystyle\frac{{n_N^2 }}{{n_{N - 1}^2 }}\right)\cosh \big[a_N^H (l_N + z_{b2})\big]\left[ {\begin{array}{c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right], \\ \end{array} \end{eqnarray}\end{document}

\begin{matrix} \begin{matrix} (\begin{matrix} β_{N} \\ γ_{N} \end{matrix}) = a_{N - 1}^{H} D_{a s y m_N - 1}^{H} \sinh [a_{N}^{H} (l_{N} + z_{b 2})] [\begin{matrix} \cosh (a_{N - 1}^{H} l_{N - 1}) \\ \sinh (a_{N - 1}^{H} l_{N - 1}) \end{matrix}] \\ + a_{N}^{H} D_{a s y m_N}^{H} (\frac{n_{N}^{2}}{n_{N - 1}^{2}}) \cosh [a_{N}^{H} (l_{N} + z_{b 2})] [\begin{matrix} \cosh (a_{N - 1}^{H} l_{N - 1}) \\ \sinh (a_{N - 1}^{H} l_{N - 1}) \end{matrix}], \end{matrix} \end{matrix}

and the recurrence relations are given by

Eq. 22

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left({\begin{array}{*{20}c} {\beta _{k - 1} } \\ {\gamma _{k - 1} } \\ \end{array}} \right) = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {a_{k - 2}^H D_{asym\_k - 2}^H \cosh \big(a_{k - 2}^H l_{k - 2}\big)} \\[5pt] {a_{k - 2}^H D_{asym\_k - 2}^H \sinh \big(a_{k - 2}^H l_{k - 2}\big)} \\[5pt] \end{array}} & {\begin{array}{*{20}c} {a_{k - 1}^H D_{asym\_k - 1}^H \big(n_{k - 1}^2 /n_{k - 2}^2\big)\sinh \big(a_{k - 2}^H l_{k - 2}\big)} \\[3pt] {a_{k - 1}^H D_{asym\_k - 1}^H \big(n_{k - 1}^2 /n_{k - 2}^2\big)\cosh \big(a_{k - 2}^H l_{k - 2}\big)} \\[3pt] \end{array}} \\ \end{array}} \right]\left({\begin{array}{*{20}c} {\beta _k } \\ {\gamma _k } \\ \end{array}} \right). \end{equation}\end{document}

(\begin{matrix} β_{k - 1} \\ γ_{k - 1} \end{matrix}) = [\begin{matrix} \begin{matrix} a_{k - 2}^{H} D_{a s y m_k - 2}^{H} \cosh (a_{k - 2}^{H} l_{k - 2}) \\ a_{k - 2}^{H} D_{a s y m_k - 2}^{H} \sinh (a_{k - 2}^{H} l_{k - 2}) \end{matrix} & \begin{matrix} a_{k - 1}^{H} D_{a s y m_k - 1}^{H} (n_{k - 1}^{2} / n_{k - 2}^{2}) \sinh (a_{k - 2}^{H} l_{k - 2}) \\ a_{k - 1}^{H} D_{a s y m_k - 1}^{H} (n_{k - 1}^{2} / n_{k - 2}^{2}) \cosh (a_{k - 2}^{H} l_{k - 2}) \end{matrix} \end{matrix}] (\begin{matrix} β_{k} \\ γ_{k} \end{matrix}) .

For the case of N ⩾ 3 and the semi-infinite N-layered turbid medium (L _N → ∞), the start term for calculating β₃ and γ₃ becomes

Eq. 23

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \left({\begin{array}{*{20}c} {\beta _N } \\ {\gamma _N } \\ \end{array}} \right) &=& a_{N - 1}^H D_{asym\_N - 1}^H \left[ {\begin{array}{*{20}c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right]\nonumber\\ && +\, a_N^H D_{asym\_N}^H \left(\frac{{n_N^2 }}{{n_{N - 1}^2 }}\right)\left[ {\begin{array}{*{20}c} {\cosh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] {\sinh \big(a_{N - 1}^H l_{N - 1}\big)} \\[3pt] \end{array}} \right],\quad \end{eqnarray}\end{document}

\begin{matrix} (\begin{matrix} β_{N} \\ γ_{N} \end{matrix}) & = & a_{N - 1}^{H} D_{a s y m_N - 1}^{H} [\begin{matrix} \cosh (a_{N - 1}^{H} l_{N - 1}) \\ \sinh (a_{N - 1}^{H} l_{N - 1}) \end{matrix}] \\ + a_{N}^{H} D_{a s y m_N}^{H} (\frac{n_{N}^{2}}{n_{N - 1}^{2}}) [\begin{matrix} \cosh (a_{N - 1}^{H} l_{N - 1}) \\ \sinh (a_{N - 1}^{H} l_{N - 1}) \end{matrix}], \end{matrix}

and the recursion formula is the same as Eq. 22.

In the case of the finite two-layered turbid medium (N = 2), the values of β₃ and γ₃ are

Eq. 24

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left({\begin{array}{*{20}c} {\beta _3 } \\ {\gamma _3 } \\ \end{array}} \right) = \left({\begin{array}{*{20}c} {\sinh \big[a_2^H (l_2 + z_{b2})\big]} \\[3pt] {\cosh \big[a_2^H (l_2 + z_{b2})\big]} \\[3pt] \end{array}} \right). \end{equation}\end{document}

(\begin{matrix} β_{3} \\ γ_{3} \end{matrix}) = (\begin{matrix} \sinh [a_{2}^{H} (l_{2} + z_{b 2})] \\ \cosh [a_{2}^{H} (l_{2} + z_{b 2})] \end{matrix}) .

For the two-layered semi-infinite N-layered turbid medium (L _N → ∞), β₃ = γ₃ = 1.

Next, we used the 2-D inverse Fourier transform to obtain the fluence rate Φ(x, y, z) in the real space. In this step, the numerical calculation is performed to obtain the solutions for the fluence rate in real space due to the difficulties of getting the analytical solutions. The expression for the 2-D inverse Fourier transform is given by

Eq. 25

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,y,z) &=& \frac{1}{{(2\pi)^2 }}\int_{ - \infty }^\infty \int_{ - \infty }^\infty \Phi _k^H (z,s)\nonumber\\ &&\times\exp [ - i(s_1 x + s_2 y)]{\mathop{\rm d}\nolimits} s_1 {\mathop{\rm d}\nolimits} s_2. \end{eqnarray}\end{document}

\begin{matrix} Φ_{k}^{H} (x, y, z) & = & \frac{1}{{(2 π)}^{2}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} Φ_{k}^{H} (z, s) \\ \times \exp [- i (s_{1} x + s_{2} y)] d s_{1} d s_{2} . \end{matrix}

The implementation of the 2-D inverse Fourier transform was completed in two methods using MATLAB language. In the first method (Inverse Fourier Transform, IFT1), [TeX:] $\Phi _k^H (x,y,z)$

Φ_{k}^{H} (x, y, z)

is expressed in a 2-D Fourier series:

Eq. 26

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} \Phi _k^H (x,y,z) = \displaystyle\frac{{\Delta s_1 \Delta s_2 }}{{(2\pi)^2 }}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\Phi _k^H (z,m\Delta s_1,n\Delta s_2)} } \\ \quad \quad \quad \quad \quad \; \times \exp [ - i(m\Delta s_1 x + n\Delta s_2 y)]. \\ \end{array} \end{equation}\end{document}

\begin{matrix} Φ_{k}^{H} (x, y, z) = \frac{Δ s_{1} Δ s_{2}}{{(2 π)}^{2}} \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} Φ_{k}^{H} (z, m Δ s_{1}, n Δ s_{2}) \\ \times \exp [- i (m Δ s_{1} x + n Δ s_{2} y)] . \end{matrix}

Equation 26 can be simplified by using rotational symmetry²⁰

Eq. 27

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,y,z) &=& \frac{{\Delta s_1 \Delta s_2 }}{{(2\pi)^2 }}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\Phi _k^H (z,m\Delta s_1,n\Delta s_2)} } \nonumber\\ && \times\, \cos (m\Delta s_1 x)\cos (n\Delta s_2 y). \end{eqnarray}\end{document}

\begin{matrix} Φ_{k}^{H} (x, y, z) & = & \frac{Δ s_{1} Δ s_{2}}{{(2 π)}^{2}} \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} Φ_{k}^{H} (z, m Δ s_{1}, n Δ s_{2}) \\ \times \cos (m Δ s_{1} x) \cos (n Δ s_{2} y) . \end{matrix}

Then, by simplifying the calculation of the sums, calculating in one direction (setting y = 0), and assuming the same sampling rate in both directions (Δs ₁ = Δs ₂ = Δs), it is sufficient to get the simpler expression of Eq. 27

Eq. 28

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \Phi _k^H (x,0,z) &=& \left(\frac{{\Delta s}}{{2\pi }}\right)^2 \left\{ \Phi _k^H (z,0,0)\right.\nonumber\\ &&+\, 2\sum\limits_{m = 1}^\infty {\Phi _k^H (z,m\Delta s,0)} [1 + \cos (m\Delta sx)]\nonumber\\ && +\left. 4\sum\limits_{m = 1}^\infty \sum\limits_{n = 1}^\infty \Phi _k^H (z,m\Delta s,n\Delta s) \cos (m\Delta sx)}.\nonumber\\ \end{eqnarray}\end{document}

\begin{matrix} Φ_{k}^{H} (x, 0, z) & = & {(\frac{Δ s}{2 π})}^{2} \{Φ_{k}^{H} (z, 0, 0) \\ + 2 \sum_{m = 1}^{\infty} Φ_{k}^{H} (z, m Δ s, 0) [1 + \cos (m Δ s x)] \\ + 4 \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} Φ_{k}^{H} (z, m Δ s, n Δ s) \cos (m Δ s x)\} . \end{matrix}

To obtain a better result, we typically use 6000 terms in the infinite sums and a sampling rate of [TeX:] $\Delta s = \mu _s^\prime /400$

Δ s = μ_{s}^{'} / 400

.

In the second method (IFT2), the one-dimensional inverse Hankel transform is obtained from Eq. 25

Eq. 29

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Phi _k^H (\rho,z) = \frac{1}{{2\pi }}\int_0^\infty {\Phi _k^H (z,s)sJ_0 (s\rho)ds}, \end{equation}\end{document}

Φ_{k}^{H} (ρ, z) = \frac{1}{2 π} \int_{0}^{\infty} Φ_{k}^{H} (z, s) s J_{0} (s ρ) d s,

where J ₀ is the Bessel function of the first kind and zero order. The discrete integral method is used to calculate the fluence rate. We typically use 2 × 10⁶ terms in the infinite sums to obtain high accuracy. The upper limit of the integral is set to [TeX:] $30\mu _s^\prime $

30 μ_{s}^{'}

. Additionally, the adaptive Gauss-Kronrod quadrature in the MATLAB Mathematics Toolbox³⁵ can be used to do numerical integration.

Using one of the two methods, the fluence rate in real space is obtained. Thereafter, the spatially resolved reflectance R ^H(ρ) of the first layer is given by³²

Eq. 30

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} R^H (\rho) &=& \int_{2\pi } {d\Omega [1 - R_{fres} (\theta)]} \frac{1}{{4\pi }}\big[\Phi _1^H (\rho,z = 0) \nonumber\\ && +\, 3D_{asym\_1}^H \frac{{\partial \Phi _1^H (\rho,z)}}{{\partial z}}|_{z = 0} \cos \theta\big]\cos \theta, \end{eqnarray}\end{document}

\begin{matrix} R^{H} (ρ) & = & \int_{2 π} d Ω [1 - R_{f r e s} (θ)] \frac{1}{4 π} [Φ_{1}^{H} (ρ, z = 0) \\ + 3 D_{a s y m_1}^{H} \frac{\partial Φ_{1}^{H} (ρ, z)}{\partial z} |_{z = 0} \cos θ] \cos θ, \end{matrix}

where R _fres(θ) is the Fresnel reflection coefficient for a photon with an incident angle θ relative to the boundary.

3. Results

In this section, we first compare the different methods for computing the inverse Fourier transform. Then, we compare the solutions of the steady-state hybrid diffusion-P3 equation and diffusion equation for a semi-infinite N-layered turbid medium with the Monte Carlo simulation. The principles of the Monte Carlo simulation of photon transport have been thoroughly described.³⁶ We use the MCML program to execute the Monte Carlo simulation developed by Lihong Wang.³⁶ A pencil photon beam is normally incident upon the semi-infinite turbid medium. The Henyey–Greenstein phase function is assumed for the calculation of the scattering angle. The spatial resolution of the steady-state Monte Carlo simulations is 0.1 mm and the number of initial photons is 10⁸.

3.1.

Comparison of the Different Methods for Computing the Inverse Fourier Transform

The spatially resolved reflectance from N-layered turbid media was calculated using the two methods (IFT1 and IFT2). Figure 2 compares the reflectance from a semi-infinite three-layered turbid medium calculated with IFT1 and IFT2. It can be seen that the differences between the two methods are less than 10⁻³ using MATLAB language (15 to 16 significant digits; 8 bytes). The results are ascribed to the fact that the Fourier and Henkel transforms are mathematically equivalent in circumstances of circular symmetry.

Fig. 2

Comparison of the reflectance from a semi-infinite three-layered turbid medium calculated with IFT1 and IFT2. The medium is three-layered with optical properties of μ_a1 = 0.1 mm⁻¹, μ_a2 = 0.01 mm⁻¹, μ_a3 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.3 mm⁻¹, n ₁ = 1.4, n ₂ = 1.7, and n ₃ = 1.4, respectively. The thicknesses of the three layers are l ₁ = 1 mm, l ₂ = 5 mm, and l ₃ = ∞, respectively.

Although these two methods are the mathematically simple methods to compute the inverse Fourier transform, they require lengthy calculation times. If a short calculation time is preferred, one can use the C or the Pascal language and the existing algorithms for a fast Fourier transform, which results in calculation times less than 10 ms.²⁰

3.2.

Results for the Case of Small Absorption Coefficients

Figure 3 shows a comparison of the semi-infinite three-layered diffusion equation and hybrid diffusion-P3 equation with the Monte Carlo simulations in the case of low absorption coefficients. The absorption coefficients of the three layers are μ_a1 = 0.01 mm⁻¹, μ_a2 = 0.1 mm⁻¹, and μ_a3 = 0.001 mm⁻¹, respectively. The reduced scattering coefficients of the three layers are μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, and μ_s3 ^′ = 1.3 mm⁻¹, respectively. The thicknesses of the three layers are l ₁ = 5 mm, l ₂ = 5 mm, and l ₃ = ∞, respectively. A good agreement among the diffusion equation, the hybrid diffusion-P3 equation, and the Monte Carlo simulations can be observed. Figure 4 gives the relative differences of the diffusion equation and the hybrid diffusion-P3 equation to the Monte Carlo simulations. At the radial distances of 0.5 mm < ρ < 12 mm, the differences are smaller than 5% with the noises of the Monte Carlo simulation and the numerical calculation accounted for. For large radial distances (ρ > 12 mm), the differences are smaller than 10% with the larger noises.

Fig. 3

Comparison of the solution of the semi-infinite three-layered diffusion equation (solid curve) and hybrid diffusion-P3 equation (dashed curve) with Monte Carlo simulations (circles). The medium is three-layered with optical properties of μ_a1 = 0.01 mm⁻¹, μ_a2 = 0.1 mm⁻¹, μ_a3 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.3 mm⁻¹, n ₁ = 1.4, n ₂ = 1.7, n ₃ = 1.4, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, respectively. The thicknesses of the three layers are l ₁ = 5 mm, l ₂ = 5 mm, l ₃ = ∞, respectively.

Fig. 4

Relative difference of the curves shown in Fig. 3.

We can observe from Figs. 3 and 4 that the solutions between the diffusion equation and the hybrid diffusion-P3 equation are nearly the same when the absorption coefficient is small. As can be seen in Sec. 2.2, the hybrid diffusion-P3 equation is obtained by replacing D _k with [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ . The main difference between the diffusion equation and the hybrid diffusion-P3 equation is the diffusion coefficient D of the diffusion equation and the asymptotic diffusion coefficient D _asym of the hybrid diffusion-P3 equation. Figure 5 gives the comparison between the diffusion coefficient D and the asymptotic diffusion coefficient D _asym in the range of μ_a. It can be observed from the results that as the absorption coefficient decreases, the values of the two coefficients (D _k and [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ ) become identical. Thus, it can be concluded that as the absorption coefficient decreases, the solutions for the N-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent. The conclusion of this section accords with that from Liemert and Kienle's paper.²⁰

Fig. 5

Comparison of the diffusion coefficient D and asymptotic diffusion coefficient D _asym as a function of the absorption coefficient assuming μ_s1 ^′ = 1.2 mm⁻¹, g = 0.8.

3.3.

Results for the Case of Large Absorption Coefficients

Figure 6 shows a comparison of the semi-infinite three-layered diffusion equation and hybrid diffusion-P3 equation with the Monte Carlo simulations in the case of large absorption coefficients. The absorption coefficients of the three layers are μ_a1 = 1 mm⁻¹, μ_a2 = 0.5 mm⁻¹, and μ_a3 = 0.001 mm⁻¹, respectively. The reduced scattering coefficients of the three layers are μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, and μ_s3 ^′ = 1.3 mm⁻¹, respectively. The thicknesses of the three layers are l ₁ = 1 mm, l ₂ = 5 mm, and l ₃ = ∞. The results in Fig. 6 show that in the case of large absorption coefficients, the hybrid diffusion-P3 equation agrees well with the Monte Carlo simulation results, while R(ρ) calculated with the N-layered diffusion equation is not close to that of the Monte Carlo simulations. Figure 7 gives the relative differences between the N-layered diffusion equation and the N-layered hybrid diffusion equation to the Monte Carlo simulations. At small distances (ρ < 2 mm), both the diffusion equation and the hybrid diffusion-P3 equation greatly deviate from the Monte Carlo simulations. At the radial distances of 2 mm < ρ < 4.7 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations are less than 20%. For the large distances of 4.7 mm < ρ < 6 mm, the differences are less than 8%. However, the differences between the diffusion equation and the Monte Carlo simulations are larger than 20%. The results suggest that the diffusion equation is not valid for the case of large absorption coefficients. From the results in Figs. 6 and 7, it can be concluded that the model of the hybrid diffusion-P3 equation is more precise than the diffusion equation for the case of a medium with a large absorption coefficient.

Fig. 6

Comparison of the solution of the semi-infinite three-layered diffusion equation (solid curve) and hybrid diffusion-P3 equation (dashed curve) with Monte Carlo simulations (circles). The medium is three-layered with optical properties of μ_a1 = 1 mm⁻¹, μ_a2 = 0.5 mm⁻¹, μ_a3 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.3 mm⁻¹, n ₁ = 1.4, n ₂ = 1.7, n ₃ = 1.4, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, respectively. The thicknesses of the three layers are l ₁ = 1 mm, l ₂ = 5 mm, l ₃ = ∞, respectively.

Fig. 7

Relative difference of the curves shown in Fig. 6.

Figure 8 shows another example of the semi-infinite five-layered turbid media for large absorption coefficients. The absorption coefficients of the five layers are μ_a1 = 0.5 mm⁻¹, μ_a2 = 0.2 mm⁻¹, μ_a3 = 0.1 mm⁻¹, μ_a4 = 0.3 mm⁻¹, and μ_a5 = 0.001 mm⁻¹, respectively. The thicknesses of the five layers are l ₁ = 1 mm, l ₂ = 1 mm, l ₃ = 1 mm, l ₄ = 2 mm, and l ₅ = ∞, respectively. The results show that the hybrid diffusion-P3 equation is more similar to the Monte Carlo simulations than the diffusion equation. Figure 9 gives the relative differences between the N-layered diffusion equation and the N-layered hybrid diffusion-P3 equation to the Monte Carlo simulations. At the radial distances of 1.6 mm < ρ < 2.5 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations is from 0 to 10%, and for the large distances of 2.5 mm < ρ < 6 mm, the differences between the hybrid diffusion-P3 equation and the Monte Carlo simulations are less than 20%. At small distances (ρ < 1.6 mm), the differences are largely affected by the light source. However, the differences between the diffusion equation and the Monte Carlo simulations are 0 to 45% for radial distances greater than 1.6 mm.

Fig. 8

Comparison of the solution of the semi-infinite five-layered diffusion equation (solid curve) and hybrid diffusion-P3 equation (dashed curve) with Monte Carlo simulations (circles). The medium is five-layered with optical properties of μ_a1 = 0.5 mm⁻¹, μ_a2 = 0.2 mm⁻¹, μ_a3 = 0.1 mm⁻¹, μ_a4 = 0.3 mm⁻¹, μ_a5 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.4 mm⁻¹, μ_s4 ^′ = 1.3 mm⁻¹, μ_s5 ^′ = 1.1 mm⁻¹, n ₁ = 1.4, n ₂ = 1.5, n ₃ = 1.6, n ₄ = 1.7, n ₅ = 1.4, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, g ₄ = 0.8, g ₅ = 0.8, respectively. The thicknesses of the five layers are l ₁ = 1 mm, l ₂ = 1 mm, l ₃ = 1 mm, l ₄ = 2 mm, l ₅ = ∞, respectively.

Fig. 9

Relative difference of the curves shown in Fig. 8.

From Figs. 6 and 8, it can be seen that the hybrid diffusion-P3 equation is more precise than the diffusion equation for modeling light propagation in the N-layered semi-infinite media with large absorption coefficients. In order to investigate the versatility of this conclusion, the medium with different optical properties and geometries have been investigated. Concerning the article length, only a few results are illustrated in Figs.10, 11, 12. The results also demonstrate the accuracy of the hybrid diffusion-P3 equation for modeling light propagation in the N-layered semi-infinite media with large absorption coefficients.

Fig. 10

(a) Comparisons among the solution of the diffusion equation (solid curve), the hybrid diffusion-P3 equation (dashed curve), and the Monte Carlo simulations (circles) for semi-infinite three-layered medium. The medium is three-layered with optical properties of μ_a1 = 0.6 mm⁻¹, μ_a2 = 0.8 mm⁻¹, μ_a3 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.3 mm⁻¹, n ₁ = 1.4, n ₂ = 1.7, n ₃ = 1.4, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, respectively. The thicknesses of the three layers are l ₁ = 2 mm, l ₂ = 3 mm, l ₃ = ∞, respectively. (b) Relative difference between the curves and circles shown in (a).

Fig. 11

(a) Comparisons among the solution of the diffusion equation (solid curve), the hybrid diffusion-P3 equation (dashed curve), and the Monte Carlo simulations (circles) for semi-infinite five-layered medium. The medium is five-layered with optical properties of μ_a1 = 1 mm⁻¹, μ_a2 = 0.4 mm⁻¹, μ_a3 = 0.2 mm⁻¹, μ_a4 = 0.1 mm⁻¹, μ_a5 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.4 mm⁻¹, μ_s4 ^′ = 1.3 mm⁻¹, μ_s5 ^′ = 1.1 mm⁻¹, n ₁ = 1.4, n ₂ = 1.5, n ₃ = 1.6, n ₄ = 1.7, n ₅ = 1.4, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, g ₄ = 0.8, g ₅ = 0.8, respectively. The thicknesses of the three layers are l ₁ = 3 mm, l ₂ = 1 mm, l ₃ = 1 mm, l ₄ = 1 mm, l ₅ = ∞, respectively. (b) Relative difference between the curves and circles shown in (a).

Fig. 12

(a) Comparisons among the solution of the diffusion equation (solid curve), the hybrid diffusion-P3 equation (dashed curve), and the Monte Carlo simulations (circles) for semi-infinite seven-layered medium. The medium is seven-layered with optical properties of μ_a1 = 0.8 mm⁻¹, μ_a2 = 0.5 mm⁻¹, μ_a3 = 0.3 mm⁻¹, μ_a4 = 0.4 mm⁻¹, μ_a5 = 0.2 mm⁻¹, μ_a6 = 0.1 mm⁻¹, μ_a7 = 0.001 mm⁻¹, μ_s1 ^′ = 1.2 mm⁻¹, μ_s2 ^′ = 1.1 mm⁻¹, μ_s3 ^′ = 1.3 mm⁻¹, μ_s4 ^′ = 1.2 mm⁻¹, μ_s5 ^′ = 1.1 mm⁻¹, μ_s6 ^′ = 1.3 mm⁻¹, μ_s7 ^′ = 1.2 mm⁻¹, n ₁ = 1.4, n ₂ = 1.5, n ₃ = 1.6, n ₄ = 1.4, n ₅ = 1.7, n ₆ = 1.6, n ₇ = 1.5, g ₁ = 0.8, g ₂ = 0.8, g ₃ = 0.8, g ₄ = 0.8, g ₅ = 0.8, g ₆ = 0.8, g ₇ = 0.8, respectively. The thicknesses of the three layers are l ₁ = 1 mm, l ₂ = 1 mm, l ₃ = 1 mm, l ₄ = 1 mm, l ₅ = 1 mm, l ₆ = 1 mm, l ₇ = ∞, respectively. (b) Relative difference between the curves and circles shown in (a).

4. Discussion

This paper discusses light propagation in N-layered turbid media. The solution of the hybrid diffusion-P3 equation is derived for N-layered finite or infinite turbid media. The solution is calculated in the steady-state domain for one point source using the extrapolated boundary condition. The Fourier transform formalism is applied to derive the analytical solutions of the fluence rate in the Fourier space based on the derivation of the N-layered diffusion equation from Liemert and Kienle's paper. The numerical calculation is performed to obtain the solutions for the fluence rate in real space because of the difficulty of obtaining analytical solutions. Two inverse Fourier transform methods are developed to calculate the fluence rate in real space.

In addition, the solution of the hybrid diffusion-P3 equation having an infinite or finite thick Nth layer is compared to that of the diffusion equation and Monte Carlo simulations. The main difference between the diffusion equation and the hybrid diffusion-P3 equation is the diffusion coefficient D of the diffusion equation and the asymptotic diffusion coefficient D _asym of the hybrid diffusion-P3 equation. The values of two coefficients (D _k and [TeX:] $D_{asym\_k} $ $D_{a s y m_k}$ ) become identical as the absorption coefficient decreases, thus it can be concluded that as the absorption coefficient decreases, the solutions of the N-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent. Simulation results show that, in the case of small absorption coefficients, the solutions of the N-layered diffusion equation and hybrid diffusion-P3 equation are almost equivalent and are in agreement with the Monte Carlo simulations. Additionally, we discussed the situation of the semi-infinite three-layered, five-layered, and seven-layered turbid media for the case of large absorption coefficients. It can be observed that the model of the hybrid diffusion-P3 equation is closer to the Monte Carlo simulation than the diffusion equation. Finally, we can make the conclusion that the model of the hybrid diffusion-P3 equation can replace the diffusion equation for light propagation in the turbid media for a wide range of absorption coefficients. The hybrid diffusion-P3 equation also has greater potential applications in the field of biomedical photonics than the diffusion equation.

Note that the derived solutions must satisfy the condition of one point source in the first layer l ₁ > z ₀. As a result, the model of the N-layered hybrid diffusion-P3 equation cannot be applied to situations in which the first layer is very thin. In the future, we will be committed to solve this problem and study other boundary conditions and source terms to obtain more accurate solutions at small radial distances (ρ < 2 mm).

Acknowledgments

The authors acknowledge the fruitful discussion with Dr. Liemert from Institut für Lasertechnologien in der Medizin und Meßtechnik at Universitat Ulm regarding the solution of the diffusion equation for N-layered turbid media. The authors would like to thank Anna H. Xue from the University of Washington for her help with revisions. The authors also acknowledge the funding supports from the National Natural Science Foundation of China (30870657, 60938002), and the Tianjin Municipal Government of China (09JCZDJC18200).

References

1.

A. Sassaroli, F. Martelli, and S. Fantini, “Perturbation theory for the diffusion equation by use of the moments of the generalized temporal point-spread function. III. Frequency-domain and time-domain results,” J. Opt. Soc. Am., 27 (7), 1723 –1742 (2010). https://doi.org/10.1364/JOSAA.27.001723 Google Scholar

2.

A. Sassaroli, F. Martelli, and S. Fantini, “Higher-order perturbation theory for the diffusion equation in heterogeneous media: application to layered and slab geometries,” Appl. Opt., 48 (10), D62 –D73 (2009). https://doi.org/10.1364/AO.48.000D62 Google Scholar

3.

M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “The Green's function for the radiative transport equation in the slab geometry,” J. Phys. A: Math. Theor., 43 065402 (2010). https://doi.org/10.1088/1751-8113/43/6/065402 Google Scholar

4.

I. Barman, N. C. Dingari, N. Rajaram, J. W. Tunnell, R. R. Dasari, and M. S. Feld, “Rapid and accurate determination of tissue optical properties using least-squares support vector machines,” Biomed. Opt. Express, 2 (3), 592 –599 (2011). https://doi.org/10.1364/BOE.2.000592 Google Scholar

5.

T. M. Bydlon, S. A. Kennedy, L. M. Richards, J. Q. Brown, B. Yu, M. K. Junker, J. Gallagher, J. Geradts, L. G. Wilke, and N. Ramanujam, “Performance metrics of an optical spectral imaging system for intra-operative assessment of breast tumor margins,” Opt. Express, 18 (8), 8058 –8076 (2010). https://doi.org/10.1364/OE.18.008058 Google Scholar

6.

G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt., 45 (5), 1062 –1071 (2006). https://doi.org/10.1364/AO.45.001062 Google Scholar

7.

R. Reif, O. A’Amar, and I. J. Bigio, “Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media,” Appl. Opt., 46 (29), 7317 –7328 (2007). https://doi.org/10.1364/AO.46.007317 Google Scholar

8.

K. Vishwanath, K. Chang, D. Klein, Y. F. Deng, V. Chang, J. E. Phelps, and N. Ramanujam, “Portable, fiber-based, diffuse reflection spectroscopy (DRS) systems for estimating tissue optical properties,” Appl. Spectrosc., 65 (2), 206 –215 (2011). https://doi.org/10.1366/10-06052 Google Scholar

9.

Q. Z. Wang, K. Shastri, and T. J. Pfefer, “Experimental and theoretical evaluation of a fiber-optic approach for optical property measurement in layered epithelial tissue,” Appl. Opt., 49 (28), 5309 –5320 (2010). https://doi.org/10.1364/AO.49.005309 Google Scholar

10.

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

11.

A. Kienle, T. Glanzmann, G. Wagnieres, and H. van den Bergh, “Investigation of two-layered turbid media with time-resolved reflectance,” Appl. Opt., 37 (28), 6852 –6862 (1998). https://doi.org/10.1364/AO.37.006852 Google Scholar

12.

A. Kienle, M. S. Patterson, N. Dognitz, R. Bays, G. Wagnieres, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt., 37 (4), 779 –791 (1998). https://doi.org/10.1364/AO.37.000779 Google Scholar

13.

F. Martelli, S. del Bianco, A. Sassaroli, and G. Zaccantia, “Retrieval of the optical properties of a layered medium based on an exact analytical solution of the time dependent diffusion equation,” Proc. SPIE, 4955 530 –535 (2003). https://doi.org/10.1117/12.478216 Google Scholar

14.

F. Martelli, S. Del Bianco, and G. Zaccanti, “Perturbation model for light propagation through diffusive layered media,” Phys. Med. Biol., 50 (9), 2159 –2166 (2005). https://doi.org/10.1088/0031-9155/50/9/016 Google Scholar

15.

F. Martelli, A. Sassaroli, S. Del Bianco, Y. Yamada, and G. Zaccanti, “Solution of the time-dependent diffusion equation for layered diffusive media by the eigenfunction method,” Phys. Rev. E, 67 (5), 056623 (2003). https://doi.org/10.1103/PhysRevE.67.056623 Google Scholar

16.

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 (10), 2827 –2843 (2007). https://doi.org/10.1088/0031-9155/52/10/013 Google Scholar

17.

A. H. Barnett, “A fast numerical method for time-resolved photon diffusion in general stratified turbid media,” J. Comput. Phys., 201 (2), 771 –797 (2004). https://doi.org/10.1016/j.jcp.2004.06.017 Google Scholar

18.

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 (6), 1382 –1390 (2006). https://doi.org/10.1364/JOSAA.23.001382 Google Scholar

19.

A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: frequency and time domains,” J. Biomed. Opt., 15 (2), 025002 (2010). https://doi.org/10.1117/1.3368682 Google Scholar

20.

A. Liemert and A. Kienle, “Light diffusion in N-layered turbid media: steady-state domain,” J. Biomed. Opt., 15 (2), 025003 (2010). https://doi.org/10.1117/1.3368685 Google Scholar

21.

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. II. Layered case,” Opt. Express, 18 (9), 9266 –9279 (2010). https://doi.org/10.1364/OE.18.009266 Google Scholar

22.

A. Liemert and A. Kienle, “Light diffusion in a turbid cylinder. I. Homogeneous case,” Opt. Express, 18 (9), 9456 –9473 (2010). https://doi.org/10.1364/OE.18.009456 Google Scholar

23.

J. B. Fishkin, S. Fantini, M. J. VandeVen, and E. Gratton, “Gigahertz photon density waves in a turbid medium: theory and experiments,” Phys. Rev. E, 53 (3), 2307 –2391 (1996). https://doi.org/10.1103/PhysRevE.53.2307 Google Scholar

24.

M. G. Nichols, E. L. Hull, and T. H. Foster, “Design and testing of a white-light, steady-state diffuse reflectance spectrometer for determination of optical properties of highly scattering systems,” Appl. Opt., 36 (1), 93 –104 (1997). https://doi.org/10.1364/AO.36.000093 Google Scholar

25.

V. V. Tuchin, H. Podbielska, and C. K. Hitzenberger, “Special section on coherence domain optical methods in biomedical science and clinics,” J. Biomed. Opt., 4 (1), 94 (1999). https://doi.org/10.1117/1.429950 Google Scholar

26.

A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys., 38 (15), 2543 –2555 (2005). https://doi.org/10.1088/0022-3727/38/15/004 Google Scholar

27.

K. Maruo, M. Tsurugi, J. Chin, T. Ota, H. Arimoto, Y. Yamada, M. Tamura, M. Ishii, and Y. Ozaki, “Noninvasive blood glucose assay using a newly developed near-infrared system,” IEEE J. Sel. Top. Quantum Electron., 9 (2), 322 –330 (2003). https://doi.org/10.1109/JSTQE.2003.811283 Google Scholar

28.

H. Cui, L. An, W. Chen, and K. Xu, “Quantitative effect of temperature to the absorbance of aqueous glucose in wavelength range from 1200 nm to 1700 nm,” Opt. Express, 13 (18), 6887 –6891 (2005). https://doi.org/10.1364/OPEX.13.006887 Google Scholar

29.

E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the P3 approximation,” J. Opt. Soc. Am. A, 18 (3), 584 –599 (2001). https://doi.org/10.1364/JOSAA.18.000584 Google Scholar

30.

A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys., 220 441 –470 (2006). https://doi.org/10.1016/j.jcp.2006.07.007 Google Scholar

31.

T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys., 19 (4), 879 –888 (1992). https://doi.org/10.1118/1.596777 Google Scholar

32.

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 (1), 246 –254 (1997). https://doi.org/10.1364/JOSAA.14.000246 Google Scholar

33.

R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A, 11 (10), 2727 –2741 (1994). https://doi.org/10.1364/JOSAA.11.002727 Google Scholar

34.

H. J. Tian, Y. Liu, L. J. Wang, Y. H. Zhang, and L. F. Xiao, “Hybrid diffusion approximation in highly absorbing media and its effects of source approximation,” Chin. Opt. Lett., 7 515 –518 (2009). https://doi.org/10.3788/COL20090706.0515 Google Scholar

35.

L. F. Shampine, “Vectorized adaptive quadrature in MATLAB,” J. Comput. Appl. Math., 211 (2), 131 –140 (2008). https://doi.org/10.1016/j.cam.2006.11.021 Google Scholar

36.

L. Wang, S. L. Jacques, and L. Zheng, “MCML–Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed., 47 (2), 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F Google Scholar

Citation Download Citation

Zhenzhi Shi, Huijuan Zhao, and Kexin Xu "Hybrid diffusion-P3 equation in N-layered turbid media: steady-state domain," Journal of Biomedical Optics 16(10), 105002 (1 October 2011). https://doi.org/10.1117/1.3640810

Published: 1 October 2011

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 5 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Diffusion

Absorption

Monte Carlo methods

Fourier transforms

Tissues

Optical properties

Scattering

1.

Introduction

2.

Theory