Approximate P3 solution for the semi-infinite medium: steady state and time domain

Xichang Wang

doi:10.1117/1.JBO.22.9.095003

14 September 2017 Approximate P₃ solution for the semi-infinite medium: steady state and time domain

Xichang Wang

Author Affiliations +

Journal of Biomedical Optics, Vol. 22, Issue 9, 095003 (September 2017). https://doi.org/10.1117/1.JBO.22.9.095003

Abstract

The steady-state solution of the Green’s function obtained by the P₃ equation in a semi-infinite medium is presented, the proposed solution is a diffusion-based model. Two time-domain solutions are established: one is the solution under extrapolation boundary condition, which we call the optical parameter method, and the other corresponds to the diffusion equation, which we call the double-diffusion coefficient method. The spatial-resolved reflectance and the time-resolved reflectance are calculated. The Monte Carlo simulation is used to verify the P₃ equation. The results show that the P₃ steady-state equation and the two time-domain equations are in good agreement with the Monte Carlo simulation. In the steady state, when the distance between the detector and the light source is less than several free paths, the P₃ equation is more accurate than the diffusion equation. In other cases, the P₃ model and the diffusion model have similar results. However, when the absorption coefficient is large, P₃ is more accurate. In the time domain, the optical parameter method is more accurate, and the double-diffusion coefficient method is more consistent with the diffusion equation.

1. Introduction

Spectroscopy is a common and effective method for calculating the absorption and scattering coefficients of biological tissues and other turbid media. The method requires measuring the reflectance of the medium. An inversion algorithm is then used to reconstruct the optical properties of the medium, such as the absorption coefficient and the scattering coefficient. For this purpose, the diffusion equation is employed,¹^–⁴ which, however, has several limitations in the application. First, the diffusion theory is applied when the light absorption coefficient is smaller than the scattering coefficient. Second, the diffusion theory becomes invalid when the distance between the detector and the light source increases to several transport mean free paths.

Some researchers have published methods for extracting tissue optical properties in low-albedo biological media, or methods for small source–detector separations. Venugopalan et al.⁵ developed an equation estimating approximate increases in the diffusion of the radiation transport. Kienle and Patterson⁶ developed a method of fitting the optical properties of a reflectivity measurement of a low exposure medium close to the source by fitting the data to a series of reflection curves produced by a single Monte Carlo simulation. Jia et al.⁷ established the method of “virtual source” diffuse approximation. These methods, to some extent, improved the accuracy of the diffusion equations. However, researchers are still looking for new methods to alleviate some of their limitations further. For example, using the radiative transfer equation,⁵ called the $P_{N}$ approximation model, it is difficult to establish a mathematical expression. However, it can be simplified as a diffusion model ( $P_{1}$ approximation) and $P_{3}$ approximation. $P_{3}$ is considered more accurate than the diffuse approximation and, therefore, has recently become a major direction of research.

Boas et al.⁸ established a frequency-dependent solution for the $P_{3}$ equation in the media of infinite thickness. Hull and Foster⁹ obtained the Green’s function of the steady-state radiation transport equation in the $P_{3}$ equation. They proved that the $P_{3}$ model had better accuracy than the diffusion theory in the high absorption medium or near the source of the radiation. Dickey et al.¹⁰ pointed out that the radiance and the fluence expressions based on the $P_{3}$ equation were used to optically characterize an intralipid tissue phantom, varying in concentration of the scatterer and absorber. Chai et al.¹¹ established the $P_{3}$ equation of the spherical harmonics method. Gao et al.¹² derived the equation set of the $P_{3}$ equation and its Green’s function approximate solution. Further, a complete solution of the $P_{3}$ equation with the extrapolated boundary conditions was also obtained. However, this progression was the continuation of the theory first proposed by Boas et al.,⁸ and the solutions of the $P_{3}$ equation had already been derived by Liemert and Kienle¹³ in a semi-infinite medium illuminated by a collimated beam of light. The solution was evaluated with nearly the same computational effort as needed for solving the standard diffusion equation. For this purpose, the method used a combination of the classic spherical harmonics and rotated reference frames for solving the $P_{3}$ equations in closed form.

In this paper, a similar method as that of the diffusion equation is used. The steady-state solutions of the Green’s function by the $P_{3}$ equation in the semi-infinite media under extrapolation boundary condition are estimated, the proposed solution is a diffusion-based model. Further, two frequency-domain solutions are estimated and transformed into the time domain by Fourier transform. The spatially resolved reflectance and the time-resolved reflectance of the $P_{3}$ equation are calculated, and the Monte Carlo simulation is used to verify the $P_{3}$ equation. The $P_{3}$ equation and the diffusion equation are compared.

2. Theory

2.1.

Solution of the Green’s Function of P₃ Equation in the Steady State

We present the solution of the Green’s function of the $P_{3}$ approximate steady-state equation in the infinite medium (see Fig. 1).

Fig. 1

Schematic of the infinite medium. The parameter $z$ is taken in the vertical direction, and the interface is at $z = 0$ . The initial position of the photon is at the point ( $0,0, z_{0}$ ).

We assume that an infinitely thin beam is vertically incident on the medium. All the photons undergo first scatter at depth $z_{0}$ of one transport mean free path. An isotropic photon source is established at that point, and therefore the beam is scattered isotropically at a depth $z = z_{0} = 1 / (μ_{s}^{'} + u_{a})$ , where $μ_{s}^{'}$ and $u_{a}$ are the deduced scattering coefficients and the absorption coefficients, $μ_{s}^{'} = μ_{s} (1 - g)$ , $μ_{s}$ and $g$ are the scattering coefficients and anisotropy factor. The origin of the coordinate system is a point where the beam enters the turbid medium, and $z$ coordinate has the same direction as the incident beam. The $x$ and $y$ coordinates lie on the surface of the turbid sample and $ρ = {(x^{2} + y^{2})}^{1 / 2}$ .

If we use the Henyey–Greenstein phase function, all its higher-order moments $g_{l}$ can be determined.⁹ We assume

Eq. (1)

μ_{t}^{(l)} = μ_{a} + μ_{s} (1 - g_{l}), μ_{t}^{(0)} = μ_{a} .

Some scholars¹^,³^,¹⁴^,¹⁵ used inverse Fourier transform to solve the diffusion equation, that is, the two solutions $ϕ_{0}$ ( $ρ, z$ ) and $ϕ_{1}$ ( $ρ, z$ ) of the diffusion equation were obtained. On the basis of the diffuse model, we obtain the four solutions $ϕ_{0} (ρ, z)$ , $ϕ_{1} (ρ, z)$ , $ϕ_{2} (ρ, z)$ and $ϕ_{3} (ρ, z)$ of the $P_{3}$ equation by the inverse Fourier transform in the same way as the diffusion equation. The four solutions of the $P_{3}$ equation of Green’s functions are calculated by Eq. (2) as follows:

Eq. (2)

ϕ_{i} (ρ, z) = \frac{1}{2 π} \int_{0}^{\infty} ϕ_{i} (z, s) s J_{0} (s ρ) d s i = 0,1, 2,3 .

Equation (2) is the inverse Fourier transform.

ϕ_{i} (z, s)

is the result of the Fourier transform that is written in the form of exponentials when

z < z_{0}

. The specific expression is shown by

Eq. (3)

ϕ_{0} (z, s) = A \exp (- z v^{-}) + B \exp [(z - z_{0}) v^{-}] + C \exp (- z v^{+}) + E \exp [(z - z_{0}) v^{+}] ϕ_{1} (z, s) = A h_{1} (- v^{-}) \exp (- z v^{-}) + B h_{1} (v^{-}) \exp [(z - z_{0}) v^{-}] + C h_{1} (- v^{+}) \exp (- z v^{+}) + E h_{1} (v^{+}) \exp [(z - z_{0}) v^{+}] ϕ_{2} (z, s) = A h_{2} (- v^{-}) \exp (- z v^{-}) + B h_{2} (v^{-}) \exp [(z - z_{0}) v^{-}] + C h_{2} (- v^{+}) \exp (- z v^{+}) + E h_{2} (v^{+}) \exp [(z - z_{0}) v^{+}] ϕ_{3} (z, s) = A h_{3} (- v^{-}) \exp (- z v^{-}) + B h_{3} (v^{-}) \exp [(z - z_{0}) v^{-}] + C h_{3} (- v^{+}) \exp (- z v^{+}) + E h_{3} (v^{+}) \exp [(z - z_{0}) v^{+}] .

In Eq. (3),

ν^{+}

and

ν^{-}

are expressed as

Eq. (4)

v^{+} = \sqrt{\frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}}, v^{-} = \sqrt{\frac{- b - \sqrt{b^{2} - 4 a c}}{2 a}},

where

a, b

, and

c

are given by

Eq. (5)

a = 27 μ_{t}^{(1)}, b = - [27 μ_{t}^{(1)} (α_{0}^{2} + α_{3}^{2}) - 28 μ_{t}^{(3)} (α_{2}^{2} - α_{0}^{2})] c = 27 μ_{t}^{(1)} α_{0}^{2} α_{3}^{2} + 28 μ_{t}^{(3)} (α_{0}^{2} α_{1}^{2} - α_{2}^{2} α_{1}^{2}) .

In Eq. (5),

Eq. (6)

α_{0}^{2} = (μ_{a} + D s^{2}) / D, D = \frac{1}{3 μ_{t}^{(1)}} α_{1}^{2} = α_{2}^{2} = s^{2} a_{3}^{2} = [μ_{t}^{(3)} + 3 D_{2} s^{2}] / 3 D_{2}, D_{2} = \frac{3}{35 μ_{t}^{(2)}},

where the constant items

A, B, C, E

, and

h_{i} (ν)

in Eq. (2) can be expressed as follows:

Eq. (7)

h_{1} (v) = \frac{3 (v^{2} - α_{3}^{2})}{2 (v^{2} - α_{1}^{2})} \frac{3 v}{7 μ_{t}^{(3)}} \frac{(v^{2} - α_{0}^{2})}{2 (v^{2} - α_{2}^{2})} h_{2} (v) = - \frac{(v^{2} - α_{0}^{2})}{2 (v^{2} - α_{2}^{2})} h_{3} (v) = + 3 v / 7 μ_{t}^{(3)} \frac{(v^{2} - α_{0}^{2})}{2 (v^{2} - α_{2}^{2})},

Eq. (8)

A = - \frac{[- 1 / 2 h_{3} (v^{+}) - 1 / 3 h_{1} (v^{+})] \exp {- [(z_{0} + z_{b 1}) v^{-}]}}{[h_{3} (v^{+}) h_{1} (v^{-}) - h_{3} (v^{-}) h_{1} (v^{+})] \exp (z_{b 1} v^{-})} B = \frac{- 1 / 2 h_{3} (v^{+}) - 1 / 3 h_{1} (v^{+})}{h_{3} (v^{+}) h_{1} (v^{-}) - h_{3} (v^{-}) h_{1} (v^{+})} C = - \frac{[1 / 3 h_{1} (v^{-}) + 1 / 2 h_{3} (v^{-})] \exp {- [(z_{0} + z_{b 2}) v^{+}]}}{[h_{3} (v^{+}) h_{1} (v^{-}) - h_{3} (v^{-}) h_{1} (v^{+})] \exp (z_{b 2} v^{+})} E = \frac{1 / 3 h_{1} (v^{-}) + 1 / 2 h_{3} (v^{-})}{h_{3} (v^{+}) h_{1} (v^{-}) - h_{3} (v^{-}) h_{1} (v^{+})} .

In Eq. (4), if

s = 0

, Eq. (4) is expressed as

ν^{+} (0)

and

ν^{-} (0)

. Thus,

z_{b 1}

and

z_{b 2}

in Eq. (8) can be expressed as shown in

Eq. (9)

z_{b 1} = 2 \frac{1 + R_{eff}}{1 - R_{eff}} \frac{μ_{a}}{v_{-}^{2} (0)}, z_{b 2} = 2 \frac{1 + R_{eff}}{1 - R_{eff}} \frac{μ_{a}}{v_{+}^{2} (0)},

If the refractive index of the biological medium is 1.4,

R_{eff} = 0.493

. Whereas, if the refractive index is 1.33,

R_{eff} = 0.431

.

R_{eff}

was calculated according to Haskell et al.² We present four solutions of the

P_{3}

equation in the steady-state equation.

To calculate the $P_{3}$ equation, we present the expression of the reflectance of the $P_{3}$ equation, which is also a continuation of the diffusion equation. The solution of the $P_{3}$ equation steady-state spatially resolved reflectance $R (ρ)$ is shown by Eq. (10)

Eq. (10)

R (ρ) = k_{0} ϕ_{0} (ρ, z) + k_{1} ϕ_{1} (ρ, z) + k_{2} ϕ_{2} (ρ, z) + k_{3} ϕ_{3} (ρ, z) |_{z = 0},

where the coefficients

k_{1}, k_{2}, k_{3}

, and

k_{4}

are the constants associated with the refractive index. If the refractive index

n = 1.4

, then

k_{1} = 0.1177

,

k_{2} = - 0.3056

,

k_{3} = 0.4444

, and

k_{4} = - 0.5487

. If

n = 1.33

, then

k_{1} = 0.1320

,

k_{2} = - 0.3359

,

k_{3} = 0.4806

, and

k_{4} = - 0.5845

.

2.2.

Solution of the Green’s Function of P₃ Equation in the Time Domain

The basic idea of the time-domain method is to first establish the frequency-domain equation, through the Fourier transform, and then convert it to the time-domain equation. Here, we present two methods in the time domain.

The first method, we call the optical parameter method. If we use Eq. (11) instead of Eq. (1), we can get the equation in the frequency domain

Eq. (11)

μ_{t}^{(l)} = μ_{a} + μ_{s} (1 - g_{l}) + j \frac{ω}{c},

where

ω

is an angular frequency of the incident light and

c

is the speed at which light travels in the medium. Equations (2)–(10) can still be applied to the

P_{3}

frequency-domain equation.

The second method we call the double-diffusion coefficient method. In Eqs. (1)–(10), only Eq. (6) is modified by Eq. (12),

Eq. (12)

α_{0}^{2} = (μ_{a} + D s^{2} + j \frac{ω}{c}) / D, D = \frac{1}{3 μ_{t}^{(1)}} α_{1}^{2} = α_{2}^{2} = s^{2} a_{3}^{2} = [μ_{t}^{(3)} + 3 D_{2} s^{2} + + j \frac{ω}{c}] / 3 D_{2}, D_{2} = \frac{3}{35 μ_{t}^{(2)}} .

Equation (13) is used in the time-domain and frequency-domain methods of the diffusion equation

Eq. (13)

α_{0}^{2} = (μ_{a} + D s^{2} + j \frac{ω}{c}) / D, D = \frac{1}{3 μ_{t}^{(1)}},

where

D

is the coefficient of the diffusion equation. To correspond to the diffuse equation, we adopt

a_{3}^{2} = [μ_{t}^{(3)} + 3 D_{2} s^{2} + j \frac{ω}{c}] / 3 D_{2}, D_{2} = \frac{3}{35 μ_{t}^{(2)}}

. Since the new parameter

D_{2}

is similar to the diffusion coefficient, we name this approximation method a double-diffusion coefficient method.

We estimate different frequencies and then carry out Fourier transform for the solution in the time domain.

2.3.

Relationship Between P₃ Equation and Related Theory

We compared the relationship between the $P_{3}$ model and other related models. Kienle et al.¹⁴ used the Fourier approach to solve the diffusion equation for a two-layered turbid medium. They then used this method to publish a series of articles including the study,¹⁵ where Eq. (3) was widely used. We also used this method to solve the diffusion equation ³ and establish the $P_{3}$ equation. The spatially resolved reflectance as utilized in the diffusion equation is given by

Eq. (14)

R (ρ) = C_{1} ϕ (ρ, z = 0) + C_{2} D \frac{\partial ϕ (ρ, z = 0)}{\partial z} .

The refractive index of the biological tissue is given by 1.4 and constants

C_{1} = 0.1117

,

C_{2} = 0.3056

. From the diffusion equation, we have

ϕ_{1} = - D \partial ϕ (ρ, z = 0) / \partial z

. Using these values and expression, Eq. (14) takes the form,

R (ρ) = C_{1} ϕ_{0} (ρ, z = 0) + C_{2} ϕ_{1}

. This reflectance represents the first two terms of Eq. (10). Hull and Foster⁹ resolved a solution of the spatially resolved reflectance. However, the expression was more complex and completely different from the one we presented. Gao et al.¹² also employed the equation of reflectance given by Hull and Foster⁹

The literature⁹ presents yet another expression that is different from Eq. (10) in our study. The $k_{n}$ in the literature⁹ depends on the numerical aperture of the detector and the refractive-index mismatch. The relationship between Eq. (10) and the literature⁹ is expressed by

Eq. (15)

k_{i} = \frac{(2 i + 1)}{4 π} k_{i}^{'} (i = 0,1, 2, 3),

where

k_{i}

is the coefficient in Eq. (10) and

k_{i}^{'}

is the coefficient in Ref. 9.

In the diffuse equation, $D$ is the diffusion coefficient. Equation (6) not only has the same definition but also a further definition, in order to distinguish the diffusion coefficient, we named $D_{2}$ .

Hull and Foster⁹ presented the specific expression of $ν$ , which was also used by Gao et al.¹² Boas et al.⁸ established a specific expression of $ν$ in the frequency domain. In the limit, when the frequency domain is approaching 0, the result⁸ becomes similar to that of Hull and Foster⁹ This is expressed as follows:

Eq. (16)

v {}_{+}= \frac{1}{\sqrt{18}} {(β + \sqrt{β^{2} - γ_{a}})}^{\frac{1}{2}}, v {}_{-}= \frac{1}{\sqrt{18}} {(β - \sqrt{β^{2} - γ_{a}})}^{\frac{1}{2}} β = 27 μ_{a} μ_{t}^{(1)} + 28 μ_{a} μ_{t}^{(3)} + 35 μ_{t}^{(2)} μ_{t}^{(3)} γ_{a} = 3780 μ_{a} μ_{t}^{(1)} μ_{t}^{(2)} μ_{t}^{(3)} .

Equations (4) and (5), in the limit,

s = 0

, reduce to Eq. (16). This shows that Eq. (16) has a special form of our model.

Equation (9) is called the extrapolated boundary condition. The extrapolated boundary conditions for the $P_{3}$ equation in Refs. 9 and 12 can be modified to Eq. (17)

Eq. (17)

z_{b 1} = 2 A \frac{μ_{a}}{v_{-}^{2}}, z_{b 2} = 2 A \frac{μ_{a}}{v_{+}^{2}} .

This is consistent with our Eq. (12). Hull and Foster⁹ pointed out the drawback of extrapolated boundary conditions. The drawback is that it is does not rigorously satisfy boundary conditions for the higher-order terms. However, it can be demonstrated that this approach yields excellent agreement with experimental data, suggesting that the boundary condition is not especially sensitive to these terms.

The solution of the Green’s function by the $P_{3}$ equation is similar to the diffusion equation. The $P_{3}$ model is an extension of the diffusion equation. Our model has some similarities to the published $P_{3}$ equation,⁹^,¹² but they have essential differences. The same similarities are the extrapolation boundary condition [Eq. (17)], the interrelated coefficients [Eq. (15)].

The differences are that: (a) the solution of Hull and Foster⁹ is obtained by polynomial approximation, and our equation is obtained by strict Fourier transform and inverse transformation, (b) the solution of Hull is derived assuming the approximation of a spherically symmetric medium. Our solution is derived assuming the approximation of a nonspherically symmetric medium: (c) the calculated reflectance equation is completely different between our equation and the equation of Hull; (d) using our equations, we can derive the characteristic parameters of Hull and Foster; and (e) our equation is established on the basis of the diffusion equation, which is the continuation of the diffusion equation.

3. Results

3.1.

Comparison of P₃ Equation with Monte Carlo Simulation and Diffusion Equation in the Steady State

Monte Carlo simulation is used as a universal means of verification. The diffusion equation is often used in biomedical optics.³^,¹⁶ We use the Monte Carlo simulation¹⁷ to verify our $P_{3}$ model and compare the accuracy with the diffuse model. The Monte Carlo simulation uses a pencil photon that was normally incident onto the semi-infinite media, and use the Henyey–Greenstein phase function. The reflectance of Monte Carlo simulations depicts the results of Monte Carlo simulations of isotropic point sources in semi-infinite media.

First, we compare the $P_{3}$ model with the diffusion equation of the semi-infinite medium using the Monte Carlo simulation. We use the diffuse equation method of the paper.³^,¹⁶

If the detection distance is 0.5 to 30 mm, it is assumed that the scattering coefficient $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , and the absorption coefficients are $μ_{a} = 0.002 {mm}^{- 1}$ and $μ_{a} = 0.0005 {mm}^{- 1}$ . The results are shown in Fig. 2(a). The straight part is the $P_{3}$ equation result, and the square part is the Monte Carlo simulation result. From Fig. 2(a), we can see that the $P_{3}$ equation result is almost the same as the Monte Carlo simulation.

Fig. 2

Comparison of $P_{3}$ approximation model with Monte Carlo simulation and diffuse approximation. The optical parameters of (a) are $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , $μ_{a} = 0.002 {mm}^{- 1}$ , and $μ_{a} = 0.0005 {mm}^{- 1}$ . The optical parameter of (b) is $μ_{a} = 0.0005 {mm}^{- 1}$ . The optical parameter of (c) is $μ_{a} = 0.002 {mm}^{- 1}$ . The optical parameters of (d) are $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , $μ_{a} = 0.002 {mm}^{- 1}$ , and $μ_{a} = 0.0005 {mm}^{- 1}$ . The optical parameter of (e) is $μ_{a} = 0.0005 {mm}^{- 1}$ . The optical parameter of (f) is $μ_{a} = 0.002 {mm}^{- 1}$ .

We compare the relative error (unit: %) among the $P_{3}$ equation, the Monte Carlo simulation, and the diffuse equation. As shown in Figs. 2(b) and 2(c), the relative error is larger except for the front 1 mm. The relative error of the results is very small. In other words, the $P_{3}$ equation is very consistent with the Monte Carlo simulation results. At the same time, in addition to the near light source, the relative error between the $P_{3}$ equation and the diffuse approximation is close to zero, that is, the $P_{3}$ equation has the same result as the diffusion equation.

It is assumed that the scattering coefficient $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , and the absorption coefficients are $μ_{a} = 0.002 {mm}^{- 1}$ and $μ_{a} = 0.0005 {mm}^{- 1}$ . The results are shown in Fig. 2(d), and the relative errors are shown in Figs. 2(e) and 2(f). These confirm the previous conclusions.

If the detection distance is 0.05 to 3 mm, it is assumed that the scattering coefficient $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , and the absorption coefficient is $μ_{a} = 0.0005 {mm}^{- 1}$ . The results are shown in Fig. 3(a). The straight line is the $P_{3}$ equation result, the dotted line is the Monte Carlo simulation result. From Fig. 3(a), we can see that when the detection distance is greater than a certain value, the result of the $P_{3}$ equation is almost the same as the Monte Carlo simulation.

Fig. 3

Comparison of $P_{3}$ approximation model with Monte Carlo simulation and diffuse approximation. The optical parameters of (a) and (b) are $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , and $μ_{a} = 0.0005 {mm}^{- 1}$ . The optical parameters of (c) and (d) are $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , and $μ_{a} = 0.0005 {mm}^{- 1}$ .

We compare the relative error (unit: %) among the $P_{3}$ equation, the Monte Carlo simulation, and the diffuse equation. As shown in Fig. 3(b), when the detection distance is between 0.6 and 1.5 mm, the relative error of $P_{3}$ is small, and when the detection distance is greater than 1.5 mm, the $P_{3}$ equation has the same result as the diffuse equation. That is, when the detection distance is greater than 0.6 mm, the $P_{3}$ equation is more accurate.

It is assumed that the scattering coefficient $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , and the absorption coefficients are $μ_{a} = μ_{a} = 0.0005 {mm}^{- 1}$ . The results are shown in Fig. 3(c), and the relative error is shown in Fig. 3(d), confirming our conclusions.

The diffusion equation is usually used under the condition that the absorption coefficient is small. In the same parameters of Fig. 2(a), we increase the absorption coefficient. The specific parameters are $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , and $μ_{a} = 0.05 {mm}^{- 1}$ . We compare the $P_{3}$ equation, diffuse approximation, and Monte Carlo simulation when the distance of the detector is 0.05 to 3 mm.

The result is shown in Fig. 4(a). Figure 4(b) shows the relative error between Monte Carlo simulation, $P_{3}$ equation, and the diffusion approximation.

Fig. 4

(a) Comparison of $P_{3}$ approximation with the diffuse approximation and Monte Carlo simulation and (b) comparison of relative error of $P_{3}$ approximation with diffusion model and Monte Carlo simulation.

Figure 4(b) shows that when the detection distance is more than 2 mm, $P_{3}$ has a small error. When the detection distance is more than 2.5 mm, the relative error of $P_{3}$ is less than 4%, and the error of diffusion is greater than 8%. This shows that $P_{3}$ is more accurate.

3.2.

Comparison of the P₃ Equation with Monte Carlo Simulation in Time Domain

First, the double-diffusion coefficient method of the $P_{3}$ equation in the time domain is compared with Monte Carlo simulation. Since the double-diffusion coefficient corresponds to the diffusion equation, we add a $D_{2}$ based on the diffusion equation, so we first verify this method.

Assuming that the three semi-infinite dielectrics have a refractive index of 1.4 and their parameters are $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , and $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{a} = 0.005 {mm}^{- 1}$ , and $μ_{a} = 0.05 {mm}^{- 1}$ , we set two detection distances of 5.5 and 8.5 mm, respectively. The results are shown in Figs. 5(a) and 5(b), where the detection distances are 5.5 and 8.5 mm, respectively.

Fig. 5

Comparison of the time-resolved reflectance calculated by $P_{3}$ theory with Monte Carlo simulation (a) refractive index $n = 1.4$ , $ρ = 5.5 mm$ ; (b) $n = 1.4$ , $ρ = 8.5 mm$ , and (c) $n = 1.33$ , $ρ = 8.5 mm$ .

It can be seen from Figs. 5(a) and 5(b) that the time-resolved reflectance calculated by the double-diffusion coefficient method is consistent with the Monte Carlo simulation, and it is proved that the time-resolved reflectance calculated by the double-diffusion coefficient method is correct.

To further verify the correctness of the double-diffusion coefficient method, we compared the three media with a refractive index of 1.33. The parameters are $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , respectively. The absorption coefficients are $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{a} = 0.005 {mm}^{- 1}$ , and $μ_{a} = 0.05 {mm}^{- 1}$ , respectively. The detection distance is 8.5 mm. The results are shown in Fig. 5(c).

It can be seen from Fig. 5(c) that the medium with a refractive index of 1.33 is consistent with the Monte Carlo simulation, further proving that the time-domain method of the double-diffusion coefficient method is correct.

We have shown that the time-resolved reflectance calculated by the double-diffusive coefficient method is consistent with the Monte Carlo simulation. We have also verified that the time-resolved reflectance using the optical parametric method is correct. We could have first used the Monte Carlo simulation to verify it, but we found that the results are the same as in Fig. 5.

We further compare the two time-domain methods. The parameters in Fig. 6 are the same as the first two parameters in Fig. 5(c), using the same detection distance, we can see that the image is completely coincident, without any difference. This proves that the optical parameter method is accurate. Because the optical parameter method does not have any approximation, it is more accurate than the double-diffusion coefficient method which is the solution corresponding to the diffusion model.

Fig. 6

Comparison of two time-domain methods of $P_{3}$ approximation.

In order to compare the differences between the two methods, we calculate the relative error. In time-resolved reflectance, the relative error is usually calculated using the logarithmic relative error. The relative error is given as follows:

Eq. (18)

E = | \frac{\log R_{D} - \log R_{u}}{\log R_{D}} | \times 100 %,

where log

R_{D}

represents the time-resolved reflectance calculated by the double-diffusion coefficient method and log

R_{u}

represents the time-resolved reflectance calculated by the optical parametric method. The results are shown in Fig. 7(a). The results show that if the time exceeds 0.2 nm, the error is less than 0.3%.

Fig. 7

Comparison of relative errors of two time-domain methods in $P_{3}$ approximation (a) calculated by Eq. (18) and (b) calculated by Eq. (19).

Further, we use the relative error of the steady state as follows:

Eq. (19)

E = | \frac{R_{D} - R_{u}}{R_{D}} | \times 100 % .

The relative error calculated by Eq. (19), as shown in Fig. 7(b), comes out to be about 4% when the time is greater than 0.2 ns. Therefore, the two time-domain methods are correct. The optical parameter method is better in terms of more accurate solution whereas the double-diffusion coefficient method excels in its simplicity.

3.3.

Comparison of the P₃ Equation with the Diffuse Equation in Time Domain

We have not only verified the $P_{3}$ equation but also compared the time-domain equation of the $P_{3}$ equation and diffusion model. Since the two time-domain methods using the $P_{3}$ model and the time-domain method of the diffusion model are coincident, we simply compare the relative errors between the two $P_{3}$ models and the diffusion model.

In Figs. 8(a) and 8(b), the relative error equation in the steady state is used. The parameters of Fig. 8(a) are refractive index $n = 1.4$ , $μ_{s} = 10 {mm}^{- 1}$ , $g = 0.9$ , $μ_{a} = 0.005 {mm}^{- 1}$ . The parameters of Fig. 8(b) are refractive index $n = 1.33$ , $μ_{s} = 15 {mm}^{- 1}$ , $g = 0.9$ , $μ_{a} = 0.01 {mm}^{- 1}$ .

Fig. 8

Relative error of $P_{3}$ approximation and diffusion equation: (a) $n = 1.4$ and (b) $n = 1.33$ .

It can be seen from Fig. 8 that when the time exceeds 0.2 ns, the relative error between the optical parameter method and the diffuse equation is about 2%; the error between the double-diffusion coefficient method and the diffuse equation is very small. In fact, the double-diffusive coefficient method is almost identical to that estimated by the diffuse model. Since the optical parameter method is an accurate $P_{3}$ equation time-domain equation, the result is more accurate than the time-domain equation of the diffusion model.

If the logarithmic relative error equation is used, it is also shown that the diffusion equation is relatively consistent with the double-diffusion coefficient method. Therefore, the accuracy of the diffusion equation is not as good as that of $P_{3}$ approximating the time-domain equation, but the difference is very small. This is shown in Fig. 9. The parameters are the same as in Fig. 8(a).

Fig. 9

Logarithmic relative error of $P_{3}$ approximation and diffusion equation.

4. Conclusion

The diffusion equation becomes increasingly inaccurate with the increase of the absorption coefficient value. This justifies the need to establish a model such as the $P_{3}$ equation, which is more accurate than the diffusion equation. However, due to relative complexity, it is not easy to implement as the mathematical model. In this paper, the solution of the $P_{3}$ equation in the steady state under extrapolation boundary condition is presented in the form of Fourier transform and inverse transformation, the proposed solution is a diffusion-based model. The solution is consistent with the diffusion model.

Monte Carlo simulation is used to demonstrate the accuracy of the $P_{3}$ model in the steady state. The relative errors of the $P_{3}$ equation, diffuse equation, and Monte Carlo simulation are calculated, which indicates that $P_{3}$ is more accurate. In the steady state, the $P_{3}$ equation is similar to the diffuse equation, and it is shown that the diffusion equation and $P_{3}$ equation have almost the same value under the condition of satisfying the diffusion equation. Monte Carlo simulation shows that when there is a large absorption coefficient, the relative error of the diffusion equation is quite large, whereas the relative error of the $P_{3}$ equation is comparatively small.

Based on the steady-state equation of the $P_{3}$ equation, two kinds of frequency-domain equations are established: the optical parameter method and double-diffusion coefficient method, they are transformed into two time-domain equations by Fourier transform in the frequency domain.

By utilizing Monte Carlo simulation, we have verified two different time-domain methods with different refractive indices and different parameters. The results show that our two methods are correct.

The time-domain method of the diffusion equation is compared. The results show that the numerical and error values calculated by the double-diffusive coefficient method are consistent with the diffusion equation. Since the optical parameter method of the $P_{3}$ equation does not have any approximation under extrapolation boundary condition, its accuracy is higher than that of the double-diffusion coefficient method. Therefore, the accuracy of the optical parameter method is higher than that of the diffusion equation. However, the double-diffusion coefficient method is simpler, and the logarithmic relative error of the time domain is calculated at about 0.2%, so the error is almost negligible.

Our research method can lay the foundation for the solution of the $P_{5}$ model. We estimate that the $P_{5}$ model needs to establish the cubic equation of the $v$ in Eq. (4) and the reconstruction of Eq. (7). Based on this paper, we can study the $P_{3}$ approximate multilayer model, which can replace the diffusion equation. ³^,¹⁶^,¹⁸

Disclosures

The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.

References

1.

A. Kienle et al., “Investigation of two layered turbid medium with time-resolved reflectance,” Appl. Opt., 37 (28), 6852 –6861 (1998). http://dx.doi.org/10.1364/AO.37.006852 APOPAI 0003-6935 Google Scholar

2.

R. C. Haskell et al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A, 11 (10), 2727 –2741 (1994). http://dx.doi.org/10.1364/JOSAA.11.002727 Google Scholar

3.

X. Wang and S. M. Wang, “Light transport model in a n-layered mismatched tissue,” Waves Random Complex Medium, 16 (2), 121 –135 (2006). http://dx.doi.org/10.1080/17455030600683321 Google Scholar

4.

J. B. Fishkin et al., “Gigahertz photon density waves in a turbid medium: theory and experiments,” Phys. Rev. E, 53 (3), 2307 –2319 (1996). http://dx.doi.org/10.1103/PhysRevE.53.2307 Google Scholar

5.

V. Venugopalan, J. S. You and B. J. Tromberg, “Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source–detector separations,” Phys. Rev. E, 58 (2), 2395 –2407 (1998). http://dx.doi.org/10.1103/PhysRevE.58.2395 Google Scholar

6.

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

7.

M. Jia et al., “Virtual-source diffusion approximation for enhanced near-field modeling of photon-migration in low-albedo medium,” Opt. Express, 23 (2), 1337 –1353 (2015). http://dx.doi.org/10.1364/OE.23.001337 OPEXFF 1094-4087 Google Scholar

8.

D.A. Boas et al., “Photon migration within the

P_{3}

approximation,” Proc. SPIE, 2389 240 (1995). http://dx.doi.org/10.1117/12.209973 PSISDG 0277-786X Google Scholar

9.

E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the

P_{3}

approximation,” J. Opt. Soc. Am. A, 18 (3), 584 –599 (2001). http://dx.doi.org/10.1364/JOSAA.18.000584 JOAOD6 0740-3232 Google Scholar

10.

D. J. Dickey, R. B. Moore and J. Tulip, “Using radiance predicted by the

P_{3}

approximation in a spherical geometry to predict tissue optical properties,” Proc. SPIE, 4156 181 (2001). http://dx.doi.org/10.1117/12.413701 PSISDG 0277-786X Google Scholar

11.

C. Chai, Q. Luo and Q. Liu, “Analytical solution of

P_{3}

approximation to radiative transfer equation for an infinite homogenous media and its validity,” J. Mod. Opt., 55 (21), 3611 –3624 (2008). http://dx.doi.org/10.1080/09500340802372801 JMOPEW 0950-0340 Google Scholar

12.

Z. Gao et al., “Study of spatially resolved diffuse reflection based on P₃ approximation,” Acta Opt. Sin., 26 (8), 1220 –1225 (2006). GUXUDC 0253-2239 Google Scholar

13.

A. Liemert and A. Kienle, “Explicit solutions of the radiative transport equation in the

P_{3}

approximation,” Med. Phys., 41 (11), 111916 (2014). http://dx.doi.org/10.1118/1.4898097 MPHYA6 0094-2405 Google Scholar

14.

A. Kienle et al., “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt., 37 (4), 779 –791 (1998). http://dx.doi.org/10.1364/AO.37.000779 APOPAI 0003-6935 Google Scholar

15.

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

16.

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

17.

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

18.

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

Biography

Xichang Wang graduated from the Changchun Institute of Optics and Fine Mechanics, Chinese Academy of Sciences, in 1999. He is a professor of Yantai University. His main research interest is biomedical optics.

Citation Download Citation

Xichang Wang "Approximate P₃ solution for the semi-infinite medium: steady state and time domain," Journal of Biomedical Optics 22(9), 095003 (14 September 2017). https://doi.org/10.1117/1.JBO.22.9.095003

Received: 20 June 2017; Accepted: 21 August 2017; Published: 14 September 2017

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 4 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Diffusion

Monte Carlo methods

Reflectivity

Absorption

Refractive index

Fourier transforms

Scattering

1.

Introduction

2.

Theory

2.1.

Solution of the Green’s Function of P3 Equation in the Steady State

Fig. 1

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

2.2.

Solution of the Green’s Function of P3 Equation in the Time Domain

Eq. (11)

Eq. (12)

Eq. (13)

2.3.

Relationship Between P3 Equation and Related Theory

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

3.

Results

3.1.

Comparison of P3 Equation with Monte Carlo Simulation and Diffusion Equation in the Steady State

Fig. 2

Fig. 3

Fig. 4

3.2.

Comparison of the P3 Equation with Monte Carlo Simulation in Time Domain

Fig. 5

Fig. 6

Eq. (18)

Fig. 7

Eq. (19)

3.3.

Comparison of the P3 Equation with the Diffuse Equation in Time Domain

Fig. 8

Fig. 9

4.

Conclusion

Disclosures

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

Solution of the Green’s Function of P₃ Equation in the Steady State

Solution of the Green’s Function of P₃ Equation in the Time Domain

Relationship Between P₃ Equation and Related Theory

Comparison of P₃ Equation with Monte Carlo Simulation and Diffusion Equation in the Steady State

Comparison of the P₃ Equation with Monte Carlo Simulation in Time Domain

Comparison of the P₃ Equation with the Diffuse Equation in Time Domain