## 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] and in frequency domain it is given by

$$\mathbf{J}\left(\mathbf{r}\right)=-D\left(\mathbf{r}\right)\nabla \varphi \left(\mathbf{r}\right),$$

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–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 *ϕ*
^{(i)} 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̂**
_{1} is a unit surface normal vector directed from the medium 1 to the medium 2 at the interface and vice versa for **n̂**
_{2}.

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.

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

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**)=${\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̂**
_{1} and **n̂**
_{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 **n̂**
_{1}·**J**=ζ′^{(1)}
*ϕ*
^{(1)}-ζϕ^{(2)}ϕ^{(2)} and **n̂**
_{2}·**J**=ζ′^{(2)}ϕ^{(2)}-ζϕ^{(1)}ϕ^{(1)} 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*
_{0} is the light speed in the vacuum.

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*^{n}
*$\tilde{\varphi}$*(**r**,*t*)d*t* 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*)d*t* 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

with

and

then **Φ**
_{0} satisfies

The equation **Φ**
_{1} 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**)=[-*∂ϕ*
_{0}(**r**)/*∂µ*_{a,i}
]/*ϕ*
_{0}(**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 -*∂ϕ*
_{0}/*∂µ*_{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 **Φ**
_{1} 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}$/**Φ**
_{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),

with the delta function Robin source,

located at the detection position m on the measurement boundary Γ [8, 10] with 2*A*=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*^{n}
*$\tilde{\psi}$*(**r**,*t*)d*t* 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*)d*t*, 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

where the vector **q** has the entry *q*_{i}
=*ζφ*_{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:

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 (**Φ**
_{0},**Φ**
_{1}) and (**Ψ**
_{0},**Ψ**
_{1}) with the finite element basis functions.

## 5. Results

#### 5. 1. 2D Results

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×10^{8} photons and the isotropic scattering was assumed throughout the domain (i.e., *g*=0).

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*
_{1}≠*n*
_{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.

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*
_{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] 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.

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*
_{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.

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

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*
_{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|*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}≡1/ζ′^{(1)} and *J*_{n}
≡**n̂**
_{1}·**J**, where *ϕ*
_{1} and *ϕ*
_{1} are the fluence rate in the layer 1 and layer 2, respectively, at the interface. In principle, any set of (*ϕ*
_{1}, *ϕ*
_{2}, *J*_{n}
) is allowed under the constraint imposed by Eq. (31). We have collected (*ϕ*
_{1}, *ϕ*
_{2}, *J*_{n}
) values at the nodes along the interface for the 2D mesh shown in Fig. 3 and plotted *ϕ*
_{1}/3*J*_{n}
with respect to *ϕ*
_{2}/3*J*_{n}
in Fig. 12. As shown in Fig. 12, (*ϕ*
_{1}/3*J*_{n}
, *ϕ*
_{2}/3*J*_{n}
) pairs are located at diverse positions on the locus defined by Eq. (31) with *C*
_{1}=0.431 and 0.608 for (*n*
_{1}, *n*
_{2})=(1.33, 1.58) and (1.58, 1.33), respectively. For (*n*
_{1}, *n*
_{2})=(1.33, 1.58), The minimum value of *ϕ*
_{1}/3|*J*_{n}
| and *ϕ*
_{2}/3|*J*_{n}
| is 2.463 and 3.275, respectively, and such minima occur at x=0. When (*ϕ*
_{1}/*ϕ*
_{2})/(*n*
_{1}/*n*
_{2})^{2} becomes close to 1, *ϕ*
_{1}/3|*J*_{n}
| and *ϕ*
_{2}/3|*J*_{n}
| reaches the maximum value over hundreds. At a large x, *ϕ*
_{1}/3|*J*_{n}
| and *ϕ*
_{2}/3|*J*_{n}
| tend to saturate to roughly 7.5 and 10.8, respectively. For (*n*
_{1}, *n*
_{2})=(1.58, 1.33), the mimimum value of *ϕ*
_{1}/3|*J*_{n}
| and *ϕ*
_{2}/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, (*ϕ*
_{1}/*ϕ*
_{2})/(*n*
_{1}/*n*
_{2})^{2} 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.

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.

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