1. Introduction

It is well established that near-infrared (NIR) light propagations in tissue can be described by diffusion approximation (DA) [1] and in frequency domain it is given by

where ϕ is the fluence rate, J is the diffuse power flux, D(r)=[3(µ_a (r)+(1-g)µ_s (r))]^-1 is the diffusion coefficient, µ_a and µ_s are the absorption and the scattering coefficients, and q ₀ is an isotropic source distribution. The anisotropy factor g is typically of the order of 0.9 for biological tissue, indicating strongly forward biased scattering. The term µ′_s =(1-g)µ_s is called the reduced scattering coefficient that incorporates the anisotropic scattering effect into the isotropic diffusion equation. For brevity, ϕ(r,ω) has been written to ϕ(r) in Eq. (1).

Within DA, the boundary condition at air/tissue boundary was derived [1–5] and can be written to

where n̂ is a unit outward surface normal vector at the tissue surface. ζ is a coefficient that accounts for the Fresnel reflections at the boundary and is given by

where R_ϕ and R_J are diffuse reflectances for the fluence rate ϕ and the normal flux J_n =n̂·J, respectively. R_ϕ and R_J are given by

θ is the angle of incidence, and R_p is the Fresnel reflection coefficient for unpolarized light. Most works on NIR imaging of tissue have assumed that the refractive index is constant throughout the medium. In that case, both ϕ and J_n are continuous for the overall domain. When there is an index mismatched internal boundary inside the medium, ϕ is discontinuous across the boundary although J_n is still continuous across it [4–6]. At the refractive index mismatched interface between the medium 1 and the medium 2, the boundary condition at the interface in Ref. 4–6 can be written to

with

for (i, j)=(1, 2) and (2, 1) where ϕ ⁽ⁱ⁾ and ($R_{\varphi }^{\left(i\right)}$ ,$R_J^{\left(i\right)}$ ) are the fluence rate and the diffuse reflectances, respectively, in the medium i. n̂ ₁ is a unit surface normal vector directed from the medium 1 to the medium 2 at the interface and vice versa for n̂ ₂.

The numerical implementation of DA incorporating the interfacial index mismatch was presented first by using the boundary integral equation [5]. The finite element method was also introduced by ignoring the term J_n in Eq. (6) and the amplitude and the phase of the diffuse photon density wave (DPDW) were numerically investigated by FEM and by MC simulations for a two-layer slab model [7]. In this work, we introduce the exact FEM implementation of Eq. (6) without any approximations and present the results of its application to the slab model.

2. Finite element method

The Galerkin weighted residual method was used to build the FEM system matrices and source vectors. We derive the FEM formulation only for the case in which only one interfacial index mismatch exists in the medium, but the method can be easily extended to the case in which there are several interfacial index mismatches.

 

Fig. 1. Splitting of the domain Ω into Ω₁ and Ω₂ and its external boundary Γ into Γ₁ and Γ₂ at the index-mismatched interface S by duplications of nodes and facets in 3D (or edges in 2D) at S. As a result, Ω _L has the boundary given by Γ _L ∪ S for L=1, 2. The domains Ω₁ and Ω₂ are coupled by the discontinuity condition for the fluence rate on S.

Download Full Size | PDF

In the preprocessing step, we split the domain Ω into Ω₁ and Ω₂ at the interface S with refractive index mismatch. The exterior boundary Γ of Ω is also split into Γ₁ and Γ₂ as shown in Fig. 1. After splitting the domain, we have

in each domain Ω _L for L=1, 2. To meet the discontinuity condition given by Eq. (6), we should duplicate each node on S. If the original node is chosen to be possessed by Ω₁, the duplicated node must be assigned to Ω₂ and vice versa. We expanded the fluence rate ϕ ^(L)(r) by ϕ ^(L)(r)=${\sum }_{i=1}^{N\left[L\right]}$${\mathrm{\Phi }}_i^{\left(L\right)}$ ${\phi }_i^{\left(L\right)}$ (r) using the piece-wise linear nodal basis functions {${\phi }_i^{\left(L\right)}$ (r)} for Ω _L , where N[L] is the total number of nodes assigned to Ω _L . The FEM discretization of Eq. (8) gives

for L=1 and 2 where Λ^(L) is given by

n̂ ₁ and n̂ ₂ are the unit vectors directed toward Ω₂ and Ω₁, respectively, on S as illustrated in Fig. 1 and ζ^(L) is ζ at Γ _L determined by Eq. (3). From Eq. (6), we have n̂ ₁·J=ζ′⁽¹⁾ ϕ ⁽¹⁾-ζϕ⁽²⁾ϕ⁽²⁾ and n̂ ₂·J=ζ′⁽²⁾ϕ⁽²⁾-ζϕ⁽¹⁾ϕ⁽¹⁾ and their substitutions to Eq. (10) result in the following form:

with F ^(L)=K ^(L)+C ^(L)-iω B ^(L)+ζ^(L) A ^(L) for L=1, 2. The entries of the matrices are given by

and the source vector has terms

for L=1, 2 and M=1, 2. In Eq. (14), n ^(L) is the refractive index of Ω _L and c ₀ is the light speed in the vacuum.

 

Fig. 2. Reproduction of boundary integral results [5] by FEM. The data show the normalized |R-R _s |ϕ where the detector location is R=(0, 0, z), the source lacation is R _s =(0, 0, 10) in mm, and ϕ is the fluence rate. S is the index-mismatched interface between medium 1 with refractive index n ₁ and medium 2 with n ₂. The scattering coefficient µ_s , the absorption coefficient µ_a , and the anisotropy factor g are the same as those given in [5]: µ_{s ,1}=7.5 mm^-1, µ_{a ,1}=0.0035 mm^-1, g ₁=0.8 for medium 1 and µ_{s ,2}=5.0 mm^-1, µ_{a ,2}=0.024 mm^-1, g ₂=0.8 for medium 2. In the figure, (a) is the result for n ₂=1.0 and (b) is the result for n ₂=2.0 with n ₁=1.333 in both cases.

Download Full Size | PDF

To ensure that our FEM formulation is correct, we tested it for two different media with a simple planar interface [5]. We found that our FEM formulation reproduces the results computed previously by the boundary integral method [5]. The FEM results for the normalized products of the source-detector separation and the fluence rate are shown in Fig. 2 and they are the same as Fig. 2 and Fig. 3 in Ref. 5. The details of the optical parameters were taken from Ref. 5.

3. Extraction of time domain measurement data

The features of the frequency domain data using DPDW mode were previously investigated for the two-layer slab model [7]. We report the features available in the time domain data for the model. The (integrated) intensity and the mean time of flight can be derived from the frequency domain DA at the zero frequency limit with the aid of the temporal moments [8]. The n-th temporal moment is defined by ϕ_n (r)≡${\int }_{-\infty }^{\infty }$ tⁿ $\tilde{\varphi }$(r,t)dt where $\tilde{\varphi }$(r,t) is the fluence rate in the time domain. Because ϕ(r,ω)=${\int }_{-\infty }^{\infty }$ $\tilde{\varphi }$(r,t)exp(iωt)dt we have ϕ_n (r)=(-i∂/∂ω) ⁿ ϕ(r,ω)|_ω=0 which is the n-th order derivative of ϕ(r,ω) with respect to ω at ω=0. ϕ ₀(r) is nothing else but the DC fluence rate. It is well established that the integrated intensity E(r) and the mean time of flight <t(r)> are given by E(r)=ζϕ ₀(r) and <t(r)>=ϕ ₁(r)/ϕ ₀(r), respectively. In the finite element domain, we have E=ζΦ ₀ and <t>=Φ ₁/Φ ₀ where Φ ₀ and Φ ₁ are the DC and the first temporal moment solution vectors for the fluence rate, respectively. If we rewrite Eq. (11) to

with

and

then Φ ₀ satisfies

The equation Φ ₁ is derived from the differentiation of Eq. (18) with respect to ω at ω=0 and it becomes

From Ref. [9], the partial mean optical path length in the i-th layer can be written to c_i <t_i >(r)=[-∂ϕ ₀(r)/∂µ_a,i ]/ϕ ₀(r) where c_i is the speed of light, µ_a,i is the absorption coefficient, and <t_i (r)> is the mean time of flight, respectively, in the i-th layer. When we denote the term -∂ϕ ₀/∂µ_a,i as ${\varphi }_1^{\mathit{\mu }\mathit{,}i}$, the finite element equation for its vector ${\mathbf{\Phi }}_1^{\mathit{\mu }\mathit{,}i}$ is derived in the same fashion as Φ ₁ by differentiating Eq. (21) with respect to µ_a,i . If one finds the solution vector, ${\mathbf{\Phi }}_1^{\mathit{\mu }\mathit{,}i}$, the partial mean time vector <t _i > and the partial mean optical path length vector c_i <t _i > for the i-th layer are determined from c_i <t _i >=${\mathbf{\Phi }}_1^{\mathit{\mu }\mathit{,}i}$/Φ ₀.

4. Extraction of absorption sensitivity functions

To investigate the effects of the interfacial index mismatch on the tomographic reconstruction of the absorption profiles of the medium, the absorption sensitivities in intensity and mean time measurements need to be computed and these can be computed also by FEM. The sensitivity functions are derived from the solution of the adjoint problem of Eq. (1),

with the delta function Robin source,

located at the detection position m on the measurement boundary Γ [8, 10] with 2A=1/ζ. Like Eq. (1), ψ(r,ω) has been written to ψ(r) for brevity. When the interfacial index mismatch exists, ψ(r) must also satisfy the saltus condition given by Eq. (6). The n-th temporal moment of the adjoint solution is given by ψ_n (r)≡${\int }_{-\infty }^{\infty }$ tⁿ $\tilde{\psi }$(r,t)dt where $\tilde{\psi }$(r, t) is the adjoint solution in the time domain. Since ψ(r,ω)=${\int }_{-\infty }^{\infty }$ $\tilde{\psi }$(r,t)exp(-iωt)dt, we have ψ_n (r)=(i∂/∂ω) ⁿ ψ(r,ω)|_ω=0. In the finite element domain, the adjoint problem given by Eq. (23), Eq. (24), and Eq. (6) with the replacement of ϕ(r) by ψ(r) becomes

where the vector q has the entry q_i =ζφ_i (m) and Ψ is the adjoint solution vector. Just as we did for Φ ₀ and Φ ₁, the DC adjoint solution vector Ψ ₀ and the first temporal moment vector Ψ ₁ are obtained from the solutions of the following equations:

According to Ref. 8 and 10, the absorption sensitivity function J ^(E) in the intensity measurement and the absorption sensitivity functions J ^(<t>) in the mean time measurement are given by [8, 10]

and

where J ^(T) is the Jacobian of the first temporal moment and is given by

In Eqs. (28)–(30), all the values of the spatially varying functions are obtained from (Φ ₀,Φ ₁) and (Ψ ₀,Ψ ₁) with the finite element basis functions.

5. Results

5. 1. 2D Results

 

Fig. 3. A planar two-layer slab model [7] and the generated FEM mesh for the truncated domain. In (b), the truncated domain has 160 mm width and 80 mm height. 22474 triangle elements and 11496 nodes were initially generated over the domain. For the index-mismatched case, 161 nodes on the interface were duplicated and total 11657 nodes were used. The source was applied to the centre of the surface of the layer 1. n ₁ and n ₂ are indices of refraction of the layer 1 and 2, respectively.

Download Full Size | PDF

The two-layer slab model and the FEM mesh are shown in Fig. 3 and we kept the optical parameters used in Ref. 7: µ′_s =1.0-mm and µ_a =0.01mm^-1 for the whole computing domains. Four different sets of refractive indices (n ₁, n ₂) were generated by permutation of 1.33 and 1.58 like Ref. 7. We have assumed the collimated pencil beam incidence and created a unit delta function source at one reduced scattering length 1/µ′_s below the irradiated surface as is usually done. The upper layer has 5 mm thickness and infinite width. The lower layer has infinite thickness and width in our model. The computing domain was truncated into a rectangle region of 160 mm width and 80 mm height. 22474 triangle elements and 11496 nodes were initially generated for the domain and they were used as the input mesh for the index-matched case. For the index-mismatched case, 161 nodes on the interface were duplicated and total 11657 nodes were used for FEM computations.

The infinite regions outside the rectangle region were transformed into conforming shells with finite thickness and attached to the truncated boundary. The fluence rate at the exterior boundary of the shells was set to zero because such outermost boundary maps to the infinity. By the exterior domain mapping technique used in this work, we can obtain the correct solutions of the fluence rate in the finite regions of interest when the medium is infinite or semi-infinite. To the top surface of the upper layer, Eq. (2) was applied. MC simulations were performed by launching 1.5×10⁸ photons and the isotropic scattering was assumed throughout the domain (i.e., g=0).

 

Fig. 4. Plots of log10 [fluence rate (W/mm²)] for (a) n ₁=1.33, n ₂=1.33, (b) n ₁=1.33, n ₂=1.58, (c) n ₁=1.58, n ₂=1.33, and (d) n ₁=1.58, n ₂=1.58. n ₁ and n ₂ are indices of refraction of the layer 1 and 2, respectively. Maximum value of the contour is 0.0 and contour step is 0.5 in logarithmic scale.

Download Full Size | PDF

The fluence rates for the four different (n ₁, n ₂) sets are plotted in Fig. 4. One can observe that the correct fluence rates in the semi-infinite homogeneous domain are obtained for n ₁=n ₂. When the refractive index is mismatched (i.e., n ₁≠n ₂) at the interface between the layers, the discontinuity of the fluence rate arises at the interface. The fluence rate jumps for n ₁<n ₂ and it drops for n ₁>n ₂ when it goes across the interface from the medium 1 to the medium 2, which is a feature consistent with that shown in Fig. 2.

 

Fig. 5. Intensity and mean time (of flight) data computed by FEM and MC. n ₁ and n ₂ are refractive indices of the layer 1 and the layer 2 respectively.

Download Full Size | PDF

 

Fig. 6. Partial mean optical pathlength in the layer 1 and layer 2.

Download Full Size | PDF

The (integrated) intensity and the mean time (of flight) of detected photons are shown in Fig. 5. In the intensity, the data computed by FEM are the raw data and they are less than the MC data. However, with a proper normalization [7] or shifting upward in a logarithmic scale, they show good matching with MC data. The data show that the intensity is insensitive to the refractive index mismatch at the interface. On the other hand, in mean time, a significant difference between the data with and without the index mismatch is observed. In comparison with the case of n ₁=n ₂, the mean time of flight increases for n ₁<n ₂ and it decreases for n ₁>n ₂. We found the excellent agreements of the mean time between the FEM data and the MC data. The partial mean optical path length [9] in each layer is shown in Fig. 6. Near the source, the paths in the layer 1 are dominantly taken by photons while the paths in the layer 2 are significantly suppressed. For large source-detector separations, the majority paths change to the paths in the layer 2. We have again the good agreements between the FEM data and the MC data for the partial mean optical path length.

 

Fig. 7. Absorption sensitivity in intensity measurements for (a) n ₁=1.33, n ₂=1.33, (b) n ₁=1.33, n ₂=1.58, (c) n ₁=1.58, n ₂=1.33, and (d) n ₁=1.58, n ₂=1.58. Source-detector separations are 10, 30, and 50 mm from the left column to the right.

Download Full Size | PDF

The absorption sensitivity profiles for intensity and for mean time are shown in Fig. 7 and Fig. 8, respectively. The sensitivity functions have been computed at each node on the finite element mesh and have been plotted using the finite element interpolation functions (i.e. basis functions) throughout the domain. The sensitivity functions in intensity measurements shown in Fig. 7 are more localized than those in mean time measurements shown in Fig. 8, which is a general feature of the sensitivity functions [8]. For the index-matched case of n ₁=n ₂, the sensitivity distributions are essentially identical for the index of 1.33 and 1.58 and only the absolute magnitudes differ between the two. For the index-mismatched case, the discontinuities of the sensitivity functions occur for both the intensity and the mean time of flight measurements. In comparison with the case of n ₁=n ₂, the sensitivity in the layer 2 increases for n ₁<n ₂ and it decreases for n ₁>n ₂ below the interface. This indicates that the interfacial index mismatch will have significant influences on the quality of the tomographic reconstruction of the absorption coefficients in the layered structure often encountered in real tissues. The suppression of the sensitivity in the deeper tissue layer (i.e. layer 2) for n ₁>n ₂ implies that the measured intensity and mean time data will be less altered by the absorption changes there. The enhancement of sensitivity in the deeper tissue layer for n ₁<n ₂ implies that the measured data will be greatly altered by the absorption changes there in comparion with the index-matched case. Therefore, the absorption inhomogeneities in the deeper layer will be easy or hard to reconstruct depending upon the refractive index mismatch at the interface.

 

Fig. 8. Absorption sensitivity in mean time measurements for (a) n ₁=1.33, n ₂=1.33, (b) n ₁=1.33, n ₂=1.58, (c) n ₁=1.58, n ₂=1.33, and (d) n ₁=1.58, n ₂=1.58. Source-detector separations are 10, 30, and 50 mm from the left column to the right.

Download Full Size | PDF

5.2. 3D Results

We have also computed the intensity, the mean time of flight, and the mean optical path length for the 3D two-layer slab model. The computing domain was truncated into a cylindrical region of 50 mm radius and 50 mm height. Just as we did in 2D, the infinite regions outside the truncated domain were transformed into conforming shells of the finite thickness attached to the artificial boundary of the truncated domain. 291762 tetrahedral elements and 55106 nodes were initially generated for the domain and they were used as the input mesh for the index-matched case. For the index-mismatched case, 1870 nodes on the interface were duplicated and total 56976 nodes were used for the FEM computations.

 

Fig. 9. Cut-through view of the tetrahedral mesh and log10[fluence rate (W/mm²)] computed by FEM for the 3D two-layer slab model: (a) tetrahedral mesh, (b) fluence rate for (n ₁, n ₂)=(1.33, 1.58), and (c) fluence rate for (n ₁, n ₂)=(1.58, 1.33). n ₁ and n ₂ are refractive indices of the layer 1 and the layer 2, respectively. Maximum value of the contour is 0.0 and contour step is 0.5 in logarithmic scale.

Download Full Size | PDF

 

Fig. 10. Intensity, mean time (of flight), and partial mean optical path length data computed by FEM and MC in 3D. Symbols are results from MC simulations and solid lines are results from FEM. n ₁ and n ₂ are refractive indices of the layer 1 and the layer 2, respectively.

Download Full Size | PDF

 

Fig. 11. Absorption sensitivity in mean time measurements at the detection position located at 30 mm radial distance computed by 3D FEM for (a) (n ₁, n ₂)=(1.33, 1.58) and for (b) (n ₁, n ₂)=(1.58, 1.33). n ₁ and n ₂ are refractive indices of the layer 1 and the layer 2, respectively.

Download Full Size | PDF

The results are shown in Fig. 9–11. The mesh and the fluence rates shown in Fig. 9 and sensitivity functions in Fig. 11 are plotted as cut-through views with respect to a plane including the axis of symmetry of the cylinder. In Fig. 9 and Fig. 11, we have plotted only the index-mismatched cases because the fluence rate and the sensitivity functions in 3D for the index-matched case are well understood by the photon migration community. The intensity, the mean time, and the partial mean optical path length data obtained from the FEM and the MC simulations are shown in Fig. 10 and the good agreement between the results from the FEM and the MC simulations are found again in 3D. In the qualitative aspects, however, there are no essential differences between the results from 2D and from 3D calculations.

6. Discussions

The mean time varies significantly by the internal refractive index mismatch. Its counterpart in frequency domain may be the phase of the DPDW, which also exhibits the influence of the internal refractive index mismatch [7]. When n ₂ varies, the mean time increases for n ₂>n ₁ and decreases for n ₂<n ₁ in comparison with the case for n ₂=n ₁. This phenomenon comes partly from the variation of the light speed. But the investigation of the optical path length shows that the distance the average photon travels in each layer is significantly altered by the Fresnel reflections at the internal boundary. For n ₂>n ₁, if photons entered the layer 2 from the layer 1 they have difficulties in coming back to the layer 1 due to the total internal reflections and tend to reside in the layer 2; as a result, the absorption sensitivity in the layer 2 increases. For n ₂<n ₁, the transmission of light into the layer 2 from the layer 1 tends to be prohibited again due to the total internal reflections and photons tend to reside in the layer 1; in this case, the absorption sensitivity in the layer 2 decreases. For n ₂=n ₁, there are no restrictions for the photons to move across the interface. The unbalance of the incoming and the outgoing photons across the interface due to the index mismatch is responsible for the path length variations, for the mean time variations and also for the measurement sensitivity variations.

Because we have still used the diffusion approximation (DA) in the presence of the refractive index mismatch at the internal boundary, the validity of DA may be addressed at the interface. The validity of DA near the boundary has often been stated in terms of the ratio of the fluence rate to the normal flux at the boundary ϕ/3|J_n |. In order for the DA to be valid, ϕ/3|J_n | should be large compared to 1. The air/tissue boundary condition given by Eq. (2) is often called the partial current boundary condition (PCBC) and the existence of such boundary strains the DA [3]. For PCBC, the value of ϕ/3|J_n | is deterministic by the boundary and depends only on the refractive index of the underlying tissue. For the tissue/tissue boundary, however, ϕ/3|J_n | is not deterministic. It is generally dependent upon the positions along the interface and the geometry of the interface. Rearranging the terms in Eq. (6) gives rise to

with C ₁≡1/ζ′⁽¹⁾ and J_n ≡n̂ ₁·J, where ϕ ₁ and ϕ ₁ are the fluence rate in the layer 1 and layer 2, respectively, at the interface. In principle, any set of (ϕ ₁, ϕ ₂, J_n ) is allowed under the constraint imposed by Eq. (31). We have collected (ϕ ₁, ϕ ₂, J_n ) values at the nodes along the interface for the 2D mesh shown in Fig. 3 and plotted ϕ ₁/3J_n with respect to ϕ ₂/3J_n in Fig. 12. As shown in Fig. 12, (ϕ ₁/3J_n , ϕ ₂/3J_n ) pairs are located at diverse positions on the locus defined by Eq. (31) with C ₁=0.431 and 0.608 for (n ₁, n ₂)=(1.33, 1.58) and (1.58, 1.33), respectively. For (n ₁, n ₂)=(1.33, 1.58), The minimum value of ϕ ₁/3|J_n | and ϕ ₂/3|J_n | is 2.463 and 3.275, respectively, and such minima occur at x=0. When (ϕ ₁/ϕ ₂)/(n ₁/n ₂)² becomes close to 1, ϕ ₁/3|J_n | and ϕ ₂/3|J_n | reaches the maximum value over hundreds. At a large x, ϕ ₁/3|J_n | and ϕ ₂/3|J_n | tend to saturate to roughly 7.5 and 10.8, respectively. For (n ₁, n ₂)=(1.58, 1.33), the mimimum value of ϕ ₁/3|J_n | and ϕ ₂/3|J_n | is 4.915 and 3.337, respectively, and their overall variation patterns with respect to x are qualitatively identical to the previous case. Even the index-matched case shows the similar patterns. As shown in Fig. 12, (ϕ ₁/ϕ ₂)/(n ₁/n ₂)² is close to 1 in the vicinity of x=11 mm, where the sign of J_n changes and the maximum value of ϕ _1,2/3|J_n | appears.

 

Fig. 12. Scattered positions of (ϕ ₁/3J_n , ϕ ₂/3J_n ) on the locus imposed by the constraint and the spatial variation of the ratio ϕ ₁/ϕ ₂ along the interface. In the left figure, the pairs outside the range have not been shown. In the right figure, x is the distance from the centre of the interface.

Download Full Size | PDF

With this point in mind, one should not be confused by the issue on the strain of the DA for the saltus conditions given by Eq. (6) or Eq. (31). When the interface is in the vicinity of the source, the minimum value of ϕ _1,2/3|J_n | may reduce and possibly the DA is strained at the positions around which the minimum occurs. But the strain of the DA in this case occurs not because of the presence of the interface but because of the location of the interface near the source.

7. Conclusions

We have formulated the FEM implementation of DA incorporating the effect of the refractive index mismatch at the internal boundary of diffusive medium without any approximations. The reliability of the proposed FEM formulation has been demonstrated by reproducing the results generated previously by the boundary integral equation technique and by comparison with the MC results for the two-layer slab model. By investigating the light propagations in the two-layer slab model with the FEM and MC simulations, we have shown that the mean time and the absorption sensitivity functions vary significantly as the index mismatch develops at the internal boundary.