Modeling boundary measurements of scattered light using the corrected diffusion approximation.

We study the modeling and simulation of steady-state measurements of light scattered by a turbid medium taken at the boundary. In particular, we implement the recently introduced corrected diffusion approximation in two spatial dimensions to model these boundary measurements. This implementation uses expansions in plane wave solutions to compute boundary conditions and the additive boundary layer correction, and a finite element method to solve the diffusion equation. We show that this corrected diffusion approximation models boundary measurements substantially better than the standard diffusion approximation in comparison to numerical solutions of the radiative transport equation.

1. Introduction

Non-invasive boundary measurements of light scattered by tissues are important for biomedical applications [1, 2]. By extracting information from these measurements regarding the optical properties of tissues (e.g. absorption and scattering), one may gain valuable insight into tissue health. For example, in diffuse optical tomography (DOT) and fluorescence diffuse optical tomography (fDOT), one seeks to reconstruct the optical properties of tissues from measurements of light at the boundary of the domain. The applications include, for example, detection and classification of breast cancer, monitoring of infant brain tissue oxygenation level and functional brain activation studies, for reviews see e.g. [3-5]. In quantitative photoacoustic tomography (QPAT) one seeks to estimate concentration of chromophores, such as hemoglobin and melanin, inside tissues by combining the optical contrast and ultrasound propagation [6-8].

Image reconstruction problems in DOT, fDOT and QPAT are non-linear ill-posed inverse problems. There are no direct methods for the solution of these problems, and thus they are typically stated as minimization problems such as regularized output least squares. The iterative solution of this problem requires repetitive solutions of the forward model. Therefore, it is essential to have a computationally feasible forward model that describes light propagation in tissues accurately.

The theory of radiative transport governs light propagation in tissues [1, 2]. This theory takes into account absorption and scattering due to inhomogeneities in the medium. The major challenge in using radiative transport theory to study light propagation in tissues is that it is mathematically complicated due mostly to the large number of variables in the radiative transport equation (RTE). The large dimensionality of the RTE makes even computational methods challenging.

Tissues typically scatter light strongly and absorb light weakly. For that reason, one often replaces the RTE by the diffusion approximation (DA) [1-4]. In the DA, one assumes that the light becomes nearly isotropic due to strong multiple scattering. The DA is much simpler to solve than the RTE. However, applying the DA to model boundary measurements is problematic. It is well known that the DA is not valid near sources or boundaries. This is because the assumption that light is nearly isotropic is too restrictive to take into account sources and boundary conditions. Nonetheless, the DA has been used to model boundary measurements with some success despite these limitations. Regardless, there still exists a need for more accurate models of boundary measurements. Consequently, the prescription of “correct” boundary conditions [9-13] and source terms [14-16] for the DA has been a long-standing issue.

There have been some works that have taken into account sources and boundaries correctly by combining the solutions of the RTE and the DA to form a hybrid method. Wang and Jacques [17] used Monte Carlo simulations for the RTE in combination with the DA. Tarvainen et al [18] developed a coupled method combining solutions of the RTE and DA both within a finite element framework for both space and angle. This coupled method can take into account correctly boundaries, sources as well as low-scattering regions in the interior of the domain. Recently, Gao and Zhao [19] developed a sophisticated numerical method to solve the RTE. This method employs multigrid methods in both space and angle where the coarsest grid level is consistent with the DA when scattering is strong and absorption is weak.

Recently, Kim [20] presented an asymptotic analysis of the RTE leading to the so-called corrected diffusion approximation (cDA). The additive correction to the DA is given by a boundary layer solution. This boundary layer solution corrects for the error made by the DA near the boundary. It vanishes rapidly away from the boundary (on the order of one scattering mean free path) so that the standard DA approximates the solution deep within the interior of the domain. The asymptotic analysis used to derive this boundary layer solution was established in the early 1970’s in the neutron transport community [21, 22]. In fact, Pomraning and Ganapol [23] used this asymptotic analysis to study boundary conditions for the DA in detail.

In this paper, we use the cDA to model boundary measurements. Using a finite element method (FEM) to solve the DA, we are able to consider general spatial domains. By computing the boundary layer solution only for boundary points where we are modeling measurements, we show that the cDA provides a superior approximation to the solution of the RTE requiring only a small amount of more work than solving the DA itself. In particular, we consider a spatial domain Ω ⊂ ℝ2 with boundary ∂Ω. Although the simulations presented here are limited to two-dimensional (2D) case, the theory derived in [20] and reviewed here is represented in dimensional independent form allowing for the realization of the method both in two and three dimensions (3D). In 3D, the cDA consist of solving 3D diffusion equation and one-dimensional radiative transport equation for each of the measurements locations. Thus, the extension of method to 3D is straightforward.

The remainder of the paper is as follows. In Section 2, we give an overview of the asymptotic analysis of the RTE leading to the cDA. In Section 3 we give details of the calculations needed to compute the boundary condition coefficients for the diffusion equation and the boundary layer solution. In addition, we give details of the numerical method used to solve the diffusion equation. In Section 4 we show results from our computational simulations. Section 5 is the conclusions.

2. Asymptotic analysis of the radiative transport equation

The steady-state radiative transport equation governs continuous light propagation in an absorbing and scattering medium. The radiance ϕ := ϕ(r,ŝ) gives the power at position r ∈ Ω ⊂ ℝ flowing in direction ŝ ∈ S−1 with n denoting the number of spatial dimensions and S−1 denoting the unit sphere. The absorption and scattering coefficients are denoted by μ := μ(r) and μ := μ(r), respectively. The scattering operator 𝒧 is defined as The scattering phase function Θ gives the fraction of light scattered in direction ŝ due to light incident in direction ŝ′. We assume that scattering is rotationally invariant so that the scattering phase function Θ depends only ŝ · ŝ′.

To solve Eq. (1) in Ω × S−1, we must supplement boundary conditions of the form where ϕ0 is the source, ∂Ω the boundary of domain Ω and n̂ denotes the unit outward normal on ∂Ω. Here, ℛϕ denotes the reflection of light due to a mismatch in the refractive index at the boundary (see Appendix B, Eq. (B.4)). Boundary condition Eq. (3) prescribes the radiance over only the directions pointing into the domain.

Boundary measurements are given by the exitance Γ(r) defined as where T = 1 − R is the Fresnel transmission coefficient in the case of mismatched refractive indices at the boundary and mapping K implements the Snell’s law in a vector form (see Appendix B, Eq. (B.15)). Furthermore, other quantity of interest is the fluence rate ϒ(r), defined as an integral of the radiance over angular directions

We consider the case in which scattering is strong and absorption is weak. To make this assumption explicit, we introduce a small, dimensionless parameter 0 < ε ≪ 1 according to By substituting Eq. (6) into Eq. (1), we obtain We seek an asymptotic solution of Eq. (7) in the limit as ε → 0+ of the form Here, Φ denotes the interior solution and Ψ denotes the boundary layer solution. This asymptotic analysis was done recently in [20]. In what follows, we summarize the results from that work.

Seeking the interior solution Φ in the form we find that Φ0 = Φ0(r) and Φ1 = −nκŝ · ∇Φ0 with Φ0 satisfying the diffusion equation The diffusion coefficient κ is defined as with g denoting the anisotropy factor defined as the mean of the cosine of the scattering phase function The leading order behavior of the interior solution Φ ∼ Φ0(r) − εnκŝ · ∇Φ0(r) is a weakly linear function of ŝ. In general, a weakly linear function of ŝ is not sufficient to satisfy boundary condition Eq. (3). For that reason, we add a boundary layer solution Ψ to correct the interior solution near the boundary. This boundary layer solution decays rapidly away from the boundary on a length scale that is O(ε). Thus, ϕ ∼ Φ0 − εnκŝ · ∇Φ0 deep in the interior of Ω far away from the boundary ∂Ω.

To compute Ψ near a particular boundary point r ∈ ∂Ω, consider a coordinate system (ρ,z) where ρ is a vector parallel to the tangent plane at r and z is the coordinate along −n̂. Let z = εζ, μ = ŝ · (−n̂) and ŝ⊥ = ŝ + μn̂. Then, each of the terms in the boundary layer solution Ψ(ρ,ζ, ŝ) = Ψ0(ρ,ζ, ŝ) + εΨ1(ρ,ζ,ŝ) + O(ε2) satisfies boundary value problems for the one-dimensional radiative transport equations of the form and with σ̄ = σ(r) and ∇⊥ denoting the gradient with respect to ρ. Both Ψ0 and Ψ1 must satisfy the asymptotic matching condition

To ensure asymptotic matching condition Eq. (15) is satisfied, one must set Here, 𝓟 is the operator that projects boundary sources for the one-dimensional half space problem onto the mode of the solution that does not decay as ζ → 0. We give more details about 𝓟 in Section 3.

For the special case in which the boundary source is axisymmetric about n̂ so that ϕ0 = ϕ0(r,μ), the right-hand side of Eq. (14a) vanishes identically and Eq. (16) reduces to the familiar Robin boundary condition with The function R(μ) is the Fresnel reflection coefficient defined with respect to μ which is defined on the local coordinate system with respect to r. Equation (10) together with the boundary condition Eq. (17) is known as the diffusion approximation to the RTE. The standard diffusion approximation, which is derived from the spherical harmonics expansion of the RTE, is shown in Appendix A.

Upon solution of Eq. (10) subject to boundary condition Eq. (17), the boundary layer solution in axisymmetric case Ψ(ζ,μ) = Ψ0 + εΨ1 satisfies

Upon solution of boundary value problem Eq. (21), the corrected diffusion approximation evaluated at r ∈ ∂Ω is then given by Equation (22) gives the asymptotic solution up to O(ε2).

3. Numerical implementations

In this section, we give a method used to compute the boundary condition coefficients given by Eqs. (18) – (20) and the solution of the boundary value problem Eq. (21) given the diffusion approximation. Furthermore, we describe briefly a finite element method to solve Eq. (10) subject to boundary condition Eq. (17).

3.1. Computing boundary condition coefficients and boundary layer solutions

To compute the coefficients of the boundary conditions for the diffusion approximation and the boundary layer solution for Ω ⊂ ℝ2, we need to study the canonical half space problem: Here, θ parameterizes S1 and cosθ = ŝ · (−n̂). The boundary source ψ is assumed to be even so that ψ(θ) = ψ(−θ) which corresponds to a boundary source that is axisymmetric with respect to n̂. We use the Henyey-Greenstein scattering phase function [24, 25]

Since boundary condition Eq. (23b) prescribes an even function of θ at ζ = 0 and because the scattering phase function given in Eq. (24) is rotationally invariant, the solution is an even function of θ: For that case, we have Let μ = cosθ. Then, we can rewrite boundary value problem Eq. (23) as Here, the redistribution function h is defined as

3.1.1. Plane wave solutions

We use plane wave solutions to solve boundary value problem Eq. (27). Plane wave solutions are special solutions of Eq. (27a) of the form Ψ = e(μ). Substituting this ansatz into Eq. (27a), we obtain the eigenvalue problem There are several properties of plane wave solutions that are useful for computing solutions of the radiative transport equation [26, 27].

To calculate plane wave solutions numerically, we use the discrete ordinate method. In particular, we use the Gauss-Chebyshev quadrature rule with Replacing the integral operation in Eq. (27a) with the Gauss-Chebyshev quadrature rule and evaluating that result at μ, we obtain Equation (32) is a discrete eigenvalue problem suitable for numerical computations. Notice that in Eq. (32) we have added a small regularization parameter 0 < δ ≪ 1 which is equivalent to adding a small amount of absorption. This ensures that the eigenvalues we calculate numerically are real and distinct. The numerical error incurred by introducing δ is exponentially small. For our numerical calculations, we typically have chosen δ = 10−8.

Suppose we solve numerically Eq. (32). We will obtain N eigenvalues λ and eigenvectors V(μ). For each pair [λ,V(μ)] satisfying Eq. (32), the pair [−λ,V(−μ)] satisfies Eq. (32) also. As a result, we order and index the eigenvalues according to Using this indexing the symmetry of the plane wave solutions corresponds to λ− = −λ and V−(μ) = V(−μ). The eigenvectors are orthogonal according to We normalize the eigenvectors according to

3.1.2. Boundary condition coefficients

In Eqs. (18) – (20), the coefficients a, b and f are defined in terms of a projection operator 𝓟. This operator is given as an expansion in plane wave solutions derived in [20]. Here, we state the result. Let p for i = N/2 + 1,···,N be defined as where y satisfies the N/2 × N/2 linear system of equations Then, we compute the boundary condition coefficients through evaluation of

3.1.3. Boundary layer solution

Now that we have computed the boundary condition coefficients in boundary condition Eq. (17), let us suppose that we have solved Eq. (10) subject to boundary condition Eq. (17). We show how we solve this boundary value problem numerically in Section 3.2. Then, the boundary layer solution Ψ satisfying boundary value problem Eq. (21) can be computed as an expansion in plane wave solutions. This expansion is derived in [20]. Here, we give the result for Ψ(0,μ) for −1 ≤ μ ≤ 1 which is what we need to correct the DA at the boundary. Let H be defined as Then Ψ(0,μ) for the points corresponding to −1 ≤ μ ≤ 1 is given by

3.2. Computing the diffusion approximation

In this work, the FEM is used to solve the DA, Eq. (10) subject to Eq. (17). We follow the same procedure as in the case of the standard DA, see e.g. [28-32]. The variational formulation of the DA is where v is a test function. By representing Φ0 in a finite dimensional basis the problem is discretized. The FE-approximation of the DA can be written in the form where where k, p = 1,···,N, ϑ(r) is the nodal basis function and N is the number of spatial nodes. The source vector is Vector c = (c1, ···, c) ∈ ℝ is the solution of the DA at the nodes of the spatial grid.

3.3. Summary of the algorithm

To summarize, the procedure for a numerical solution of the cDA is given as

For a given optical parameters μ and μ, compute the asymptotic parameter ε using Eq. (50) and scaled optical parameters α and σ using Eq. (6).

Solve the eigenvalue problem, Eq. (32), at the measurement points r.

Solve the matrix y using Eq. (37). Evaluate the discrete projection operator from the Eq. (36) and the coefficients of the Robin boundary condition of the DA from Eqs. (38)–(40).

Solve the DA Eq. (10) subject to the boundary condition Eq. (17) using the FEM, Eq. (44). Evaluate the solution and the gradient at the measurement points r.

Compute the boundary layer correction using Eq. (42).

Compute the approximation to the radiance ϕ at the measurement points r using the Eq. (22).

4. Numerical results

The performance of the cDA was tested with 2D simulations. Simulation domain Ω was a circle with radius of 20 mm centered at the origin. The source ϕ0(r,ŝ) was located at (x,y) = (−20,0) mm with cosine shape giving the largest value in inward direction and value zero in direction of tangent to the boundary

Three types of test cases were considered: a homogeneous medium with matched refractive indices inside and outside the domain for three different values of scattering coefficient μ, the same cases with mismatched refractive indices, and a heterogeneous medium.

The results of the cDA were compared with the results of the standard DA (see Appendix A) and the RTE. The cDA was solved as explained above in Section 3.3. The standard DA was solved with the FEM similarly as in [32]. The FE-approximation of the RTE was implemented similarly as in [32,33] in the case of matched refractive indices. The implementation is explained in more detail in Appendix B.

The FE-mesh for the spatial discretization of the domain contained approximately 4687 nodal points and 9196 triangular elements for the homogeneous test cases and 6114 nodal points and 11 988 elements for the heterogeneous test case. For the angular discretization of the RTE 64 equally spaced angular directions were used. Parameter ε was chosen as the ratio of the total mean free path l = (μ + μ)−1 to the size of the domain L [34]

4.1. Matched refractive indices

For the first case, we consider a medium with matched refractive indices at the boundary. The refractive indices inside and outside the medium were nin = 1 and nout = 1, respectively. The scattering coefficient was given three value μ = 50, 5 and 0.5 mm−1. The absorption coefficient and the anisotropy factor were constants, μ = 0.01 mm−1 and g = 0.8, respectively.

Radiances at the boundary point (x,y) = (0,−20) for μ = 5 mm−1 computed using the cDA, the DA and the RTE are shown in left image of Fig. 1. To compare the performance of the cDA against the DA, two different quantities were investigated. First, the fluence rate inside the domain was computed using Eq. (5), and secondly the exitance at the boundary was computed using Eq. (4). The exitances computed using the cDA, the DA and the RTE are shown in Fig. 2. In addition, the relative errors of fluence rates computed using the cDA (top row) and the DA (bottom row) against the RTE for different values of μ are shown in Fig. 3. Furthermore, the relative errors of exitances are shown in Fig. 4. To give a quantitative estimate of the errors, the means of the relative errors of fluence rates ΔΦ and the exitances ΔΓ were computed. The results are given in Table 1. We also recorded the computation times of the different models. These are given in Table 2.

Fig. 1

Radiance at the boundary point (x,y) = (0,−20) computed using the cDA, the DA and the RTE for μ = 5 mm−1 with matched (left) and mismatched refractive indices (right).

Fig. 2

Logarithm of exitance at the boundary of the domain computed using the cDA, the DA and the RTE for different values of μ with matched refractive indices at the boundary.

Fig. 3

Percent relative error of fluence rate computed using the cDA (top row) and the DA (bottom row) for different values of μ with matched refractive indices at the boundary. The values are cut at 10%.

Fig. 4

Percent relative error of exitance at the boundary computed using the cDA and the DA for different values of μ with matched refractive indices.

Table 1

The mean of the relative error of fluence rate ΔΦ (%) and exitance ΔΓ(%) computed using the cDA and the DA for different values of μ(mm−1) and asymptotic parameter ε with matched refractive indices.

		ΔΦ0		ΔΓ
μs	ε	cDA	DA	cDA	DA
50	0.001	2.10	5.21	2.11	4.21
5	0.01	1.43	4.19	4.86	3.04
0.5	0.1	17.32	14.10	42.31	42.54

Table 2

The computation times of the models for different values of μ.

	Computation time			Relative computation time
μs	cDA	DA	RTE	cDA/RTE	DA/RTE
50	8.9 s	3.3 s	2.4 min	6.2 %	2.3 %
5	8.8 s	3.2 s	2.2 min	6.5 %	2.4 %
0.5	8.6 s	3.3 s	2.3 min	6.4 %	2.5 %

As it can be seen from Fig. 1, the radiance computed using the cDA agree relatively well with the RTE and satisfy the zero boundary condition in inward direction. Note that the RTE may give small negative or positive values in inward direction due to numerical reasons. The DA gives negative radiance in inward direction which are unphysical. The relative error of fluence rate is approximately 2 % for the cDA and between 4–6 % for the standard DA when scattering coefficient is μ = 50 and μ = 5 mm−1 as it can be seen from Fig. 3 and Table 1. For μ = 0.5 mm−1 both the cDA and the standard DA give large errors since the assumptions of the DA are not valid anymore. Figure 4 shows that as scattering becomes large compared to absorption, the relative error of exitance decreases for the cDA just as the asymptotic theory predicts.

The computation times in Table 2 show that solving the cDA is feasible when compared with the standard DA. In addition, solving both the cDA and the DA are much faster than solving the RTE. Thus, the cDA models light propagation more accurately than the standard DA and requires only a small amount of more work.

4.2. Mismatched refractive indices

For the second case, we consider a medium with mismatched refractive indices at the boundary. The refractive indices inside and outside the medium were nin = 1.33 and nout = 1, respectively. Other optical parameters were the same as before. Radiances at the boundary point (x,y) = (0,−20) for μ = 5 mm−1 computed using the cDA, the DA and the RTE are shown in right image of Fig. 1. Again, the quantities of interest were the fluence rate inside the domain and the exitance at the boundary. The exitances computed using different models are shown in Fig. 5. In addition, the relative errors of fluence rates and exitances are shown in Figs. 6 and 7, respectively. As before, the means of the relative errors of fluence rates and exitances were computed and the results are given in Table 3.

Fig. 5

Logarithm of exitance at the boundary of the domain computed using the cDA, the DA and the RTE for different values of μ with mismatched refractive indices.

Fig. 6

Percent relative error of fluence rate computed using the cDA (top row) and the DA (bottom row) for different values of μ with mismatched refractive indices. The values are cut at 10%.

Fig. 7

Percent relative error of exitance computed using the cDA and the DA for different values of μ with mismatched refractive indices.

Table 3

The mean of the relative error of fluence rate ΔΦ (%) and exitance ΔΓ(%) computed using the cDA and the DA for different values of μ and asymptotic parameter ε with mismatched refractive indices.

		ΔΦ0		ΔΓ
μs	ε	cDA	DA	cDA	DA
50	0.001	2.63	3.73	1.96	22.14
5	0.01	1.59	6.31	4.38	13.33
0.5	0.1	11.18	11.74	39.20	26.26

As it can be seen Fig. 1, the approximation to the radiance given by the cDA agrees better with the RTE than the approximation given by the DA. Furthermore, the results in Fig. 6 and Table 3 show the relative error of fluence rate is between 2–3 % for the cDA when μ = 50 mm−1 and μ = 5 mm−1. For the DA the relative error is larger giving the largest errors at the boundary. For the case μ = 0.5 mm−1 neither of the approximations are valid. Figure 7 shows that the relative error of exitance decreases for the cDA when the ratio of absorption and scattering decreases due to the asymptotic theory behind the model. In contrast, the relative error of the standard DA may not have a global error bound over the whole domain. In that case, the relative error might be large close to the boundary while giving satisfactory results inside the domain. In addition, decrease in ratio of absorption and scattering does not ensure decrease in relative error for the DA. Thus, the cDA models light propagation more accurately than the standard DA when compared with the RTE.

4.3. Heterogeneous medium

For the third case, we consider a medium with heterogeneous optical properties. Simulated optical parameters of the medium are shown in Fig. 8. The scattering and absorption coefficients of the background were (μ,μ) = (10,0.01) mm−1. Furthermore, optical parameters of scattering and absorbing inclusions were (μ,μ) = (20,0.01) mm−1 and (μ,μ) = (10,0.02) mm−1, respectively. The anisotropy factor in the whole domain was g = 0.8. The refractive indices inside and outside the medium were nin = 1.33 and nout = 1, respectively. Asymptotic parameter ε = 0.005 was computed from Eq. (50) using the optical parameters of the background. The exitances are shown in Fig. 9. The relative errors of fluence rates computed using the cDA and the DA are shown in Fig. 10. Furthermore, the relative errors of exitances are shown in Fig. 11.

Fig. 8

Optical parameters for the heterogeneous test case.

Fig. 9

Logarithm of exitance at the boundary of the domain computed using the cDA, the DA and the RTE for the heterogeneous test case.

Fig. 10

Percent relative error of fluence rate computed using the cDA (left) and the DA (right) for the heterogeneous test case. The values are cut at 10%.

Fig. 11

Percent relative error of exitance for the heterogeneous test case.

The mean of the relative error of fluence rate is 1.69 % for the cDA over whole domain. In contrast, the mean of the relative error of fluence rate is 5.07 % for the DA. In addition, the DA gives larger errors close to the boundary than inside the domain as it can be seen from Fig. 10. Furthermore, the relative error of exitance is smaller for the cDA (mean 5.38 %) than that of for the DA (mean 18.73 %). Thus, the cDA gives more accurate results over whole domain than the standard DA due to a global error bound given by the asymptotic theory.

5. Conclusions

In this work, recently introduced corrected diffusion approximation was numerically implemented. In the cDA, an additive correction term is computed for the DA at the boundary based on asymptotic analysis of the RTE. The procedure for computing the cDA requires only small modifications to the existing solvers for the DA. In particular, one only needs to modify the spatial coefficients of the DA and the Robin boundary condition. In addition, an additive correction term, which satisfies a one-dimensional radiative transport equation, is readily solved using the plane wave solutions.

The performance of the cDA was tested with 2D simulations. The results were compared with the results of the DA and the RTE. The results show that the cDA models boundary measurements of scattered light more accurately than the standard DA with only a small increase in computation time.