OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 12, Iss. 8 — Apr. 19, 2004
  • pp: 1727–1740
« Show journal navigation

Finite element method for diffusive light propagations in index-mismatched media

Jae Hoon Lee, Seunghwan Kim, and Youn Tae Kim  »View Author Affiliations


Optics Express, Vol. 12, Issue 8, pp. 1727-1740 (2004)
http://dx.doi.org/10.1364/OPEX.12.001727


View Full Text Article

Acrobat PDF (2576 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Near-infrared (NIR) light propagations in strongly scattering tissue have been studied in the past few decades and diffusion approximations (DA) have been extensively used under the assumption that the refractive index is constant throughout the medium. When the index is varying, the discontinuity of the fluence rate arises at the index-mismatched interface. We introduce the finite element method (FEM) incorporating the refractive index mismatch at the interface between the diffusive media without any approximations. Intensity, mean time, and mean optical path length were computed by FEM and by Monte Carlo (MC) simulations for a two-layer slab model and a good agreement between the data from FEM and from MC was found. The absorption sensitivity of intensity and mean time measurements was also analyzed by FEM. We have shown that mean time and absorption sensitivity functions vary significantly as the refractive index mismatch develops at the interface between the two layers.

© 2004 Optical Society of America

1. Introduction

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

1. A. Ishimaru, Wave propagation and scattering in random media (Academic, 1978), Chap 7–9.

] and in frequency domain it is given by

·D(r)ϕ(r)+[μa(r)+i(ωc)]ϕ(r)=q0(r),
J(r)=D(r)ϕ(r),
(1)

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 0 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

1. A. Ishimaru, Wave propagation and scattering in random media (Academic, 1978), Chap 7–9.

5

5. J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]

] and can be written to

n̂·J=ζϕ
(2)

where 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

ζ=12(1Rϕ1+RJ)
(3)

where Rϕ and RJ are diffuse reflectances for the fluence rate ϕ and the normal flux Jn =·J, respectively. Rϕ and RJ are given by

Rϕ=20π2Rp(cosθ)cosθsinθdθ,
(4)
RJ=30π2Rp(cosθ)cos2θsinθdθ,
(5)

θ is the angle of incidence, and Rp 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 Jn are continuous for the overall domain. When there is an index mismatched internal boundary inside the medium, ϕ is discontinuous across the boundary although Jn is still continuous across it [4

4. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]

6

6. G. W. Faris, “Diffusion equation boundary conditions for the interface between turbid media: a comment,” J. Opt. Soc. Am. A 19, 519–520 (2002) [CrossRef]

]. 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

n̂i·J=ζ(i)ϕ(i)ζ(j)ϕ(j)
(6)

with

ζ(i)=1Rϕ(i)2(RJ(i)+RJ(j))
(7)

for (i, j)=(1, 2) and (2, 1) where ϕ (i) and (Rϕ(i) ,RJ(i) ) are the fluence rate and the diffuse reflectances, respectively, in the medium i. 1 is a unit surface normal vector directed from the medium 1 to the medium 2 at the interface and vice versa for 2.

The numerical implementation of DA incorporating the interfacial index mismatch was presented first by using the boundary integral equation [5

5. J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]

]. The finite element method was also introduced by ignoring the term Jn 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

7. H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,” Phys. Med. Biol. 48, 2713–2727 (2003) [CrossRef] [PubMed]

]. 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 Ω1 and Ω2 and its external boundary Γ into Γ1 and Γ2 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 Γ LS for L=1, 2. The domains Ω1 and Ω2 are coupled by the discontinuity condition for the fluence rate on S.

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

·D(L)(r)ϕ(L)(r)+[μa(L)(r)+i(ωc(L))]ϕ(L)(r)=q0(L)(r)
(8)

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 Ω1, the duplicated node must be assigned to Ω2 and vice versa. We expanded the fluence rate ϕ (L)(r) by ϕ (L)(r)=i=1N[L]Φi(L) φi(L) (r) using the piece-wise linear nodal basis functions {φi(L) (r)} for Ω L , where N[L] is the total number of nodes assigned to Ω L . The FEM discretization of Eq. (8) gives

[K(L)+C(L)iωB(L)+ζ(L)A(L)]Φ(L)+Λ(L)[n̂L·J]=Q(L)
(9)

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

Λi(L)[n̂L·J]=sdSφi(L)(r)n̂L(r)·J(r).
(10)

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

[F(1)+ζ'(1)Σ(11)ζ'(2)Σ(12)ζ'(1)Σ(2)F(2)+ζ'(2)Σ(22)][Φ(1)Φ(2)]=[Q(1)Q(2)]
(11)

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

Kij(L)=ΩLdΩD(L)(r)φi(L)(r)φj(L)(r),
(12)
Cij(L)=ΩLdΩμa(L)(r)φi(L)(r)φj(L)(r),
(13)
Bij(L)=1c0ΩLdΩn(L)(r)φi(L)(r)φj(L)(r),
(14)
Aij(L)=ΓLdΓφi(L)(r)φj(L)(r),
(15)
ij(LM)=sdSφi(L)(r)φj(M)(r).
(16)

and the source vector has terms

Qi(L)=ΩLφi(L)(r)q0(L)(r)
(17)

for L=1, 2 and M=1, 2. In Eq. (14), n (L) is the refractive index of Ω L and c 0 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 1 and medium 2 with n 2. 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 1=0.8 for medium 1 and µs ,2=5.0 mm-1, µa ,2=0.024 mm-1, g 2=0.8 for medium 2. In the figure, (a) is the result for n 2=1.0 and (b) is the result for n 2=2.0 with n 1=1.333 in both cases.

To ensure that our FEM formulation is correct, we tested it for two different media with a simple planar interface [5

5. J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]

]. We found that our FEM formulation reproduces the results computed previously by the boundary integral method [5

5. J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]

]. 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

7. H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,” Phys. Med. Biol. 48, 2713–2727 (2003) [CrossRef] [PubMed]

]. 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

8. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed]

]. The n-th temporal moment is defined by ϕn (r)≡ tn ϕ˜(r,t)dt where ϕ˜(r,t) is the fluence rate in the time domain. Because ϕ(r,ω)= ϕ˜(r,t)exp(iωt)dt we have ϕn (r)=(-i∂/∂ω) n ϕ(r,ω)|ω=0 which is the n-th order derivative of ϕ(r,ω) with respect to ω at ω=0. ϕ 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)=ζϕ 0(r) and <t(r)>=ϕ 1(r)/ϕ 0(r), respectively. In the finite element domain, we have EΦ 0 and <t>=Φ 1/Φ 0 where Φ 0 and Φ 1 are the DC and the first temporal moment solution vectors for the fluence rate, respectively. If we rewrite Eq. (11) to

[K+CiωB+A+Σ]Φ=Q
(18)

with

A(ζ(1)A(1)00ζ(2)A(2))
(19)

and

Σ(ζ(1)Σ(11)ζ(2)Σ(12)ζ(1)Σ(21)ζ(2)Σ(22)),
(20)

then Φ 0 satisfies

[K+C+A+Σ]Φ0=Q.
(21)

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

[K+C+A+Σ]Φ1=BΦ0.
(22)

From Ref. [9

9. M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. Van Der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. 38, 1859–1876 (1993) [CrossRef] [PubMed]

], the partial mean optical path length in the i-th layer can be written to ci <ti >(r)=[-∂ϕ 0(r)/∂µa,i ]/ϕ 0(r) where ci is the speed of light, µa,i is the absorption coefficient, and <ti (r)> is the mean time of flight, respectively, in the i-th layer. When we denote the term -∂ϕ 0/∂µa,i as ϕ1µ,i, the finite element equation for its vector Φ1µ,i is derived in the same fashion as Φ 1 by differentiating Eq. (21) with respect to µa,i . If one finds the solution vector, Φ1µ,i, the partial mean time vector <t i > and the partial mean optical path length vector ci <t i > for the i-th layer are determined from ci <t i >=Φ1µ,i/Φ 0.

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),

·D(r)ψ(r)+[μa(r)i(ωc)]ψ(r)=0,
(23)

with the delta function Robin source,

ψ(r)+2AD(r)ψ(r)n̂=δ(rm),
(24)

located at the detection position m on the measurement boundary Γ [8

8. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed]

, 10

10. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999) [CrossRef]

] 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)≡ tn ψ˜(r,t)dt where ψ˜(r, t) is the adjoint solution in the time domain. Since ψ(r,ω)= ψ˜(r,t)exp(-iωt)dt, we have ψn (r)=(i∂/∂ω) n ψ(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

[K+C+iωB+A+Σ]Ψ=q
(25)

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

[K+C+A+Σ]Ψ0=q,
(26)
[K+C+A+Σ]Ψ1=BΨ0.
(27)

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

8. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed]

, 10

10. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999) [CrossRef]

]

J(E)(r)=ψ0(r)ϕ0(r)
(28)

and

J(<t>)(r)=[J(T)(r)<t(m)>J(E)(r)]E(m)
(29)

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

J(T)(r)=ψ1(r)ϕ0(r)+ψ0(r)ϕ1(r).
(30)

In Eqs. (28)(30), all the values of the spatially varying functions are obtained from (Φ 0,Φ 1) and (Ψ 0,Ψ 1) 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 1 and n 2 are indices of refraction of the layer 1 and 2, respectively.

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 1, n 2) 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×108 photons and the isotropic scattering was assumed throughout the domain (i.e., g=0).

Fig. 4. Plots of log10 [fluence rate (W/mm2)] for (a) n 1=1.33, n 2=1.33, (b) n 1=1.33, n 2=1.58, (c) n 1=1.58, n 2=1.33, and (d) n 1=1.58, n 2=1.58. n 1 and n 2 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.

The fluence rates for the four different (n 1, n 2) sets are plotted in Fig. 4. One can observe that the correct fluence rates in the semi-infinite homogeneous domain are obtained for n 1=n 2. When the refractive index is mismatched (i.e., n 1n 2) at the interface between the layers, the discontinuity of the fluence rate arises at the interface. The fluence rate jumps for n 1<n 2 and it drops for n 1>n 2 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 1 and n 2 are refractive indices of the layer 1 and the layer 2 respectively.
Fig. 6. Partial mean optical pathlength in the layer 1 and layer 2.

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

7. H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,” Phys. Med. Biol. 48, 2713–2727 (2003) [CrossRef] [PubMed]

] 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 1=n 2, the mean time of flight increases for n 1<n 2 and it decreases for n 1>n 2. We found the excellent agreements of the mean time between the FEM data and the MC data. The partial mean optical path length [9

9. M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. Van Der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. 38, 1859–1876 (1993) [CrossRef] [PubMed]

] 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=1.33, n 2=1.33, (b) n 1=1.33, n 2=1.58, (c) n 1=1.58, n 2=1.33, and (d) n 1=1.58, n 2=1.58. Source-detector separations are 10, 30, and 50 mm from the left column to the right.

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

8. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed]

]. For the index-matched case of n 1=n 2, 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 1=n 2, the sensitivity in the layer 2 increases for n 1<n 2 and it decreases for n 1>n 2 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 1>n 2 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 1<n 2 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=1.33, n 2=1.33, (b) n 1=1.33, n 2=1.58, (c) n 1=1.58, n 2=1.33, and (d) n 1=1.58, n 2=1.58. Source-detector separations are 10, 30, and 50 mm from the left column to the right.

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/mm2)] computed by FEM for the 3D two-layer slab model: (a) tetrahedral mesh, (b) fluence rate for (n 1, n 2)=(1.33, 1.58), and (c) fluence rate for (n 1, n 2)=(1.58, 1.33). n 1 and n 2 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.
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 1 and n 2 are refractive indices of the layer 1 and the layer 2, respectively.
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 1, n 2)=(1.33, 1.58) and for (b) (n 1, n 2)=(1.58, 1.33). n 1 and n 2 are refractive indices of the layer 1 and the layer 2, respectively.

The results are shown in Fig. 911. 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

7. H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,” Phys. Med. Biol. 48, 2713–2727 (2003) [CrossRef] [PubMed]

]. When n 2 varies, the mean time increases for n 2>n 1 and decreases for n 2<n 1 in comparison with the case for n 2=n 1. 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 2>n 1, 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 2<n 1, 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 2=n 1, 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|Jn |. In order for the DA to be valid, ϕ/3|Jn | 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

3. R. C. Haskell, L. O. Svaasand, T. -T. Tsay, T. -C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

]. For PCBC, the value of ϕ/3|Jn | is deterministic by the boundary and depends only on the refractive index of the underlying tissue. For the tissue/tissue boundary, however, ϕ/3|Jn | 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

ϕ1Jn=(n1n2)2(ϕ2Jn)+C1
(31)

Fig. 12. Scattered positions of (ϕ 1/3Jn , ϕ 2/3Jn ) on the locus imposed by the constraint and the spatial variation of the ratio ϕ 1/ϕ 2 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.

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|Jn | 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.

Acknowledgments

We thank the reviewers for helpful comments to revise our work. This work has been supported by the Ministry of Information and Communication of Korea.

References and links

1.

A. Ishimaru, Wave propagation and scattering in random media (Academic, 1978), Chap 7–9.

2.

M. Keizer, M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. 27, 1820–1824 (1988). [CrossRef]

3.

R. C. Haskell, L. O. Svaasand, T. -T. Tsay, T. -C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

4.

R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]

5.

J. Ripoll and M. Nieto-Vesperinas, “Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,” J. Opt. Soc. Am. A 16, 1947–1957 (1999). [CrossRef]

6.

G. W. Faris, “Diffusion equation boundary conditions for the interface between turbid media: a comment,” J. Opt. Soc. Am. A 19, 519–520 (2002) [CrossRef]

7.

H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,” Phys. Med. Biol. 48, 2713–2727 (2003) [CrossRef] [PubMed]

8.

S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite-element-method calculations,” Appl. Opt. 34, 8026–8037 (1995) [CrossRef] [PubMed]

9.

M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. Van Der Zee, and D. T. Delpy, “A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,” Phys. Med. Biol. 38, 1859–1876 (1993) [CrossRef] [PubMed]

10.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999) [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(170.3660) Medical optics and biotechnology : Light propagation in tissues
(170.7050) Medical optics and biotechnology : Turbid media
(290.1990) Scattering : Diffusion

ToC Category:
Research Papers

History
Original Manuscript: February 18, 2004
Revised Manuscript: April 2, 2004
Published: April 19, 2004

Citation
Jae Hoon Lee, Seunghwan Kim, and Youn Kim, "Finite element method for diffusive light propagations in index-mismatched media," Opt. Express 12, 1727-1740 (2004)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-12-8-1727


Sort:  Journal  |  Reset  

References

  1. A. Ishimaru, Wave propagation and scattering in random media (Academic, 1978), Chap 7-9.
  2. M. Keizer, M. Star, and P. R. M. Storchi, �??Optical diffusion in layered media,�?? Appl. Opt. 27, 1820-1824 (1988). [CrossRef]
  3. R. C. Haskell, L. O. Svaasand, T. �??T. Tsay, T. �??C. Feng, M. S. McAdams, and B. J. Tromberg, �??Boundary conditions for the diffusion equation in radiative transfer,�?? J. Opt. Soc. Am. A 11, 2727-2741 (1994). [CrossRef]
  4. R. Aronson, "Boundary conditions for diffusion of light,�?? J. Opt. Soc. Am. A 12, 2532-2539 (1995). [CrossRef]
  5. J. Ripoll and M. Nieto-Vesperinas, �??Index mismatch for diffuse photon density waves at both flat and rough diffuse-diffuse interfaces,�?? J. Opt. Soc. Am. A 16, 1947-1957 (1999). [CrossRef]
  6. G. W. Faris, �??Diffusion equation boundary conditions for the interface between turbid media: a comment,�?? J. Opt. Soc. Am. A 19, 519-520 (2002) [CrossRef]
  7. H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. Paulsen, �??The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach,�?? Phys. Med. Biol. 48, 2713-2727 (2003) [CrossRef] [PubMed]
  8. S. R. Arridge and M. Schweiger, �??Photon measurement density functions. Part 2: Finite-element-method calculations,�?? Appl. Opt. 34, 8026-8037 (1995) [CrossRef] [PubMed]
  9. M. Hiraoka, M. Firbank, M. Essenpreis, M. Cope, S. R. Arridge, P. Van Der Zee, and D. T. Delpy, �??A Monte Carlo investigation of optical pathlength in inhomogeneous tissue and its application to near-infrared spectroscopy,�?? Phys. Med. Biol. 38, 1859-1876 (1993) [CrossRef] [PubMed]
  10. S. R. Arridge, �??Optical tomography in medical imaging,�?? Inverse Problems 15, R41-R93 (1999) [CrossRef]

Cited By

Alert me when this paper is cited

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


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited