1.

Introduction

The use of visible or near-infrared light provides many possibilities, particularly for applications in biomedical treatment and diagnosis, such as optical tomography, laser surgery, and photodynamic therapy. To develop an efficient tool, it is essential to have high-order accuracy models for predicting light propagation within inhomogeneous, absorbing, and highly forward-scattering media, typical for biological tissues. Modeling of light propagation through biological tissues has traditionally been done with the diffusion approximation (DA)^1,2 or the statistical Monte Carlo (MC) method.^3–11 The diffusion theory presumes that the scattering is predominate and that the medium is optically diffuse so that the angle-dependent radiant intensity is replaced by an angle-independent photon flux and the radiative transfer equation (RTE) is simplified as a diffusion equation. However, the DA is hardly applicable to heterogeneous biological tissues with nonscattering or low-scattering regions; moreover, experiments² have shown that it fails to match experimental data when the tissue sample is not optically diffuse. The MC technique was broadly used because of its simple algorithm and flexibility to include real physical conditions. Nethertheless, this method requires a large number of photons to obtain smooth accurate solutions, which is highly time consuming. Alternatively, different (deterministic) numerical methods have been developed to solve the RTE in the spatial domain, for example, the finite difference method,^12–15 the finite element method (FEM),^16,17 or the finite volume method (FVM).¹⁸ By using an angular discretization of the RTE, the cited methods are considered as a discrete ordinates method (DOM).¹⁹ It is known that one of the more serious shortcomings of the DOM is false scattering, which is a consequence of spatial discretization errors. Another serious drawback of the method is the so-called ray effect, which is a consequence of angular discretization (see, for example, chapter 16 in Ref. 19). In our earlier work, a modified FVM (MFVM) for solving the two-dimensional (2-D) RTE has been developed.²⁰ Recently, two authors of the present paper succeeded in deriving exact analytical solutions to the RTE for a variety of different geometries.^21–23 It has to be emphasized that these exact analytical solutions are very useful for validation of the numerical methods.

This paper presents a comparison between three different solutions of the RTE based on the MFVM,²⁰ the analytical approach,^21,22 and the MC technique⁶ for predicting light propagation within absorbing and highly forward-scattering media, subjected to a collimated light beam. Few works can be found in the literature on the same topic.^24–28 In Ref. 24, calculations based on the DA, the DOM for solving the RTE in the time domain, and the MC technique were conducted and compared for predicting light transport in multidimensional biological tissues. In that case, the DA were solved with the FEM using a commercial package FEMLAB. Comparisons among the three methods were performed over a broad range of parameters, such as the scattering and absorption coefficients, the heterogeneity of the tissues, the CPU time, etc. It was concluded that the DOM was found to closely match the MC simulation. Da Silva et al.²⁵ compared the auxiliary function method (giving an exact solution) and the DOM for a heterogeneous absorbing and scattering (one-dimensional) slab composed of a single flat layer or a multilayer in the steady-state domain. The comparison was applied to two different media presenting two typical and extreme scattering behaviors: Rayleigh and Mie scattering with smooth or very anisotropic phase functions, respectively. A good agreement between the two methods was observed in both cases. Mishra et al.²⁶ studied numerically a square short-pulse laser having pulse-width of the order of a femtosecond on transmittance and reflectance signals in case of an absorbing and scattering planar layer. A comparison of the performance of the three numerical methods for solving the RTE, namely, the discrete transfer method, the DOM, and the FVM, was performed. Effects of the optical properties (extinction coefficient, scattering albedo, anisotropy factor) and the laser properties, such as the pulse-width and the angle of incidence, on the transmittance and the reflectance signals were compared. In all the cases, results of the three methods were found to compare very well with each other. Computationally, the DOM was found to be the most efficient. In Ref. 27, a comparison of light transport models in view of optical tomography applications was performed. These models were based on two FEMs associated with the DOM (for solving the 2-D RTE in the steady-state domain), namely, the least square and the discontinuous Galerkin finite element (DGFE) formulations. The comparison of both methods using the $S_n$ and the $T_n$ angular quadratures shows that the DGFE gives more accurate solutions for problems with strong discontinuities, but may exhibit some oscillations due to the Galerkin procedure. A test on a collimated irradiation showed that both methods give the same accuracy due to the separation of the radiance into collimated and diffuse components, which removes the discontinuities in the implementation of the boundary conditions. Recently, the RTE in the time domain was solved for an axisymmetric cylindrical medium using both the DOM and the FVM.²⁸ Steady and transient flux profiles were determined for absorbing and scattering media. Results for each solution method were compared and shown for various grid numbers, scattering albedos, and optical thicknesses. A comparison of CPU time and memory usage among the methods was presented. It was found that the FVM uses more memory and has a longer convergence time than the DOM for all cases due to the difference in angular treatment.

The remainder of this paper is organized as follows. Section 2 deals with the problem statement for 2-D geometries and gives details of the MFVM. In our previous publication,²⁰ only homogeneous biological media have been considered and the boundaries were assumed to be transparent. Also, some extensions of the numerical method, regarding the space-varying optical properties and the interface between two media (refractive index mismatch with Fresnel reflection at the interface), have been developed for the present work. Moreover, more complex collimated light beams and a new time semi-discretization of the RTE are employed here. The third section analyzes and discusses the comparison results among the different methods cited above. The last section reports our conclusions and suggests further developments for this work.

2.

Numerical Treatment of the Radiative Transfer Equation

2.1.

Modified Finite Volume Method

2.1.1.

Problem statement

A 2-D absorbing and highly forward-scattering medium, such as biological tissue exposed to arbitrary collimated external radiation, is analyzed. Except for the refractive index, which is constant, the medium has nonuniform optical properties to model, in the present study, a homogeneous or a two-layered medium and, in the future, the complex heterogeneities inside biological tissues. The external radiation that penetrates from the outside into the participating medium is brought by a perpendicular collimated light beam. One part of it propagates in the medium without being deviated, while the other part is scattered in all directions of space. Then, it is convenient to split the radiance into two components denoted as ${\psi }_c(s,\mathbf{\Omega },t)$, where ${\psi }_c$ is the collimated radiance ($\mathrm{W}{\mathrm{mm}}^{-2}{\mathrm{sr}}^{-1}$), $s$ is the spatial position (mm), $\mathbf{\Omega }$ is the direction vector, and $t$ is the time (ps), and $\psi (s,\mathbf{\Omega },t)$, which are, respectively, the collimated and diffuse components of light.¹⁹ The intensity ($\mathrm{W}{\mathrm{mm}}^{-2}$) of the collimated light beam is $\mathrm{\Upsilon }(s_w,t)$, given at any location point $s_w$ on the bounding surface that is illuminated. Then, the ${\psi }_c(s,\mathbf{\Omega },t)$ collimated radiance is governed by Bouguer-Beer-Lambert attenuation law.

Eq. (1)

$${\psi }_c(s,\mathbf{\Omega },t)=\mathrm{\Upsilon }(s_w,t-\frac{n{\mathrm{\Delta }}_s}{c})\exp [-{\int }_{s_w}^s{\mu }_t(u)\mathrm{d}u]\delta (\mathbf{\Omega }-{\mathbf{\Omega }}_c)\mathrm{for}t\ge \frac{n{\mathrm{\Delta }}_s}{c}\text{and}{\psi }_c(s,\mathbf{\Omega },t)=0\text{otherwise},$$

where $n$ is the refractive index, $c$ is the speed of light in vacuum (0.299793 $\mathrm{mm}{\mathrm{ps}}^{-1}$), ${\mu }_t$ is the attenuation coefficient (${\mathrm{mm}}^{-1}$), $\delta$ is the Dirac-delta function, ${\mathbf{\Omega }}_c$ is the direction of the collimated light beam, and ${\mathrm{\Delta }}_s=|s-s_w|$. The diffuse radiance $\psi (s,\mathbf{\Omega },t)$ is solution of the RTE with an additional radiation source term ($S$) ($\mathrm{W}{\mathrm{mm}}^{-2}$) due to the scattered part of the collimated light beam.^19,20

Eq. (2)

$$S_c(s,\mathbf{\Omega },t)={\mu }_s(s)p({\mathbf{\Omega }}_c\to \mathbf{\Omega })\mathrm{\Upsilon }(s_w,t-\frac{n{\mathrm{\Delta }}_s}{c})\exp [-{\int }_{s_w}^s{\mu }_t(u)\mathrm{d}u],$$

where ${\mu }_s$ is the scattering coefficient (${\mathrm{mm}}^{-1}$) and $p$ is the scattering phase function (${\mathrm{sr}}^{-1}$). The scattering phase function $p({\mathbf{\Omega }}^{\prime }\to \mathbf{\Omega })$ describes the probability that a light beam from one direction, ${\mathbf{\Omega }}^{\prime }$, will be scattered into a certain other direction $\mathbf{\Omega }$. The fluence $\varphi (s,t)$ ($\mathrm{W}{\mathrm{mm}}^{-2}$) and the outgoing flux at the (transparent) bounding surface of the medium $Q_{\text{out}}(s,t)$ can be calculated, respectively, as

Eq. (3)

$$\phi (s,t)={\int }_{2\pi }\psi (s,\mathbf{\Omega },t)+{\psi }_c(s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }={\int }_{2\pi }\psi (s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }+\mathrm{\Upsilon }(s_w,t-\frac{n{\mathrm{\Delta }}_s}{c})\exp [-{\int }_{s_w}^s{\mu }_t(u)\mathrm{d}u]Q_{\text{out}}(s,t)={\int }_{\mathbf{\Omega }·n_{\text{out}}>0}\psi (s,\mathbf{\Omega },t)(\mathbf{\Omega }·{\mathit{n}}_{\text{out}})\mathrm{d}\mathbf{\Omega }\{+\pi \mathrm{\Upsilon }(s_w,t-\frac{n{\mathrm{\Delta }}_s}{c})\exp [-{\int }_{s_w}^s{\mu }_t(u)\mathrm{d}u]\text{if}{\mathbf{\Omega }}_c·{\mathit{n}}_{\text{out}}>0\},$$

where the contribution of $\mathrm{\Upsilon }$ is valid only for $t\ge \frac{n{\mathrm{\Delta }}_s}{c}$. $\mathit{n}$ is the local unit surface normal pointing into the medium. It should be noted that the reflectance (respectively transmittance) corresponds to the outgoing flux at the bounding surface (resp. opposite bounding surface) that is illuminated.

2.1.2.

Discretization of the time, angular, and spatial domains

The time domain was discretized with a $\mathrm{\Delta }t$ constant step such that $t_k=k\mathrm{\Delta }t$, where $k$ was the number of iteration. The angular space ($2\mathrm{\Pi }\mathrm{Sr}$) was uniformly subdivided into ${\mathbf{\Omega }}^m$ discrete directions, $m\in \{1,\dots ,N_{\theta }\}$. The computational spatial domain was divided into three-node triangular elements using unstructured (or structured) meshes. As in Ref. 20, a cell-vertex formulation²⁹ was adopted in this work to have a high resolution in space of the RTE. It consists of building control volumes around each node of the mesh and computing the solution at the nodes of the mesh (nodes of triangles). The polygonal control volumes connected to each node were built by joining the centroids of the elements to the midpoints of the corresponding sides (Fig. 1). Control volumes surrounding the nodes at the boundaries of the medium were built as illustrated in Ref. 30. The surface of the control volume related to node $P$ was subdivided into $N_f$ surface elements. For a surface element $f$, $i_f$ was the integration point located at the center of the surface element, $A_f$ was its surface area (${\mathrm{mm}}^2$) and, $n_f$ was the local unit outward normal vector. It should be noted that the $i_f$ integration points of panels $f(f=1,2,\dots ,12)$ are defined only from the coordinates of the vertices of triangles of the mesh (Fig. 1).

Fig. 1

Control volume related to an interior node.

2.1.3.

Full discretization of the RTE

FVM applied to the RTE for the spatial domain

The integration of the RTE over a $V_P$ control volume (${\mathrm{mm}}^3$) (Fig. 1) and into a $\mathrm{\Delta }{\mathbf{\Omega }}^m$ discrete solid angle centered around an ${\mathbf{\Omega }}^m$ discrete direction yields^19,20

Eq. (4)

$$\begin{gathered} \frac{n}{c}{\int }_{V_P}{\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\frac{\partial \psi (s,\mathbf{\Omega },t)}{\partial t}\mathrm{d}\mathbf{\Omega }\mathrm{d}V+{\int }_{V_P}{\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\mathbf{\Omega }·\nabla \psi (s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }\mathrm{d}V+{\int }_{V_P}{\mu }_t(s){\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\psi (s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }\mathrm{d}V \\ ={\int }_{V_P}{\mu }_s(s){\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}{\int }_{{\mathbf{\Omega }}^{\prime }=2\pi }p({\mathbf{\Omega }}^{\prime }\to \mathbf{\Omega })\psi (s,{\mathbf{\Omega }}^{\prime },t)\mathrm{d}{\mathbf{\Omega }}^{\prime }\mathrm{d}\mathbf{\Omega }\mathrm{d}V+{\int }_{V_P}{\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}{S}_c(s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }\mathrm{d}V, \end{gathered}$$

where $\psi (s,\mathbf{\Omega },t)$ is the radiance defined at time $t$ and at position $s$ traveling in the direction of unit vector $\mathbf{\Omega }$ [located here in the plane $(Oxy)$] where light propagates at the velocity $v=c/n$ in the medium. $S_c$ is the radiation source term given by Eq. (2). By applying the divergence theorem to the second term of the left member of Eq. (4), it follows that

Eq. (5)

$${\int }_{V_P}{\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\mathbf{\Omega }·\nabla \psi (s,\mathbf{\Omega },t)\mathrm{d}\mathbf{\Omega }\mathrm{d}V={\int }_{{\mathrm{\Gamma }}_P}{\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\psi (s,\mathbf{\Omega },t)(\mathbf{\Omega }·{\mathit{n}}_{\text{out}})\mathrm{d}\mathbf{\Omega }\mathrm{d}S.$$

Let ${\psi }_{i_f}^{k,m}$ be an approximation of the radiance at time $t_k$, and at the $i_f$ integration position traveling in the ${\mathbf{\Omega }}^m$ discrete direction. Then, Eq. (5) is approximated by

Eq. (6)

$$\sum _{f=1}^{N_f}{\psi }_{i_f}^{k,m}{A}_f{\mathrm{\Delta }}_f^m\text{with}{\mathrm{\Delta }}_f^m={\int }_{\mathrm{\Delta }{\mathbf{\Omega }}^m}\mathbf{\Omega }·n_f\mathrm{d}\mathbf{\Omega }.$$

It should be noted that ${\mathrm{\Delta }}_f^m$ is an integral that depends only on the orientation of surface element $f$ for the direction considered. This quantity is calculated in an exact way.³⁰ We assume that ($\psi$, $S$) and optical properties (${\mu }_a$, ${\mu }_s$), where ${\mu }_a$ is the absorption coefficient (${\mathrm{mm}}^{-1}$), are constant inside a sufficiently small $V_P$ control volume (taking only one value at node $P$) and inside $\mathrm{\Delta }{\mathbf{\Omega }}^m$. The full discretization of the steady-state Eq. (4) yields the following algebraic equation:

Eq. (7)

$$\sum _{f=1}^{N_f}{\psi }_{i_f}^m{A}_f{\mathrm{\Delta }}_f^m=\{-{\mu }_{tP}{\psi }_P^m+{\mu }_{sP}\sum _{m^{\prime }=1}^{N_{\theta }}{\tilde{p}}^{m^{\prime }\to m}{\psi }_P^{m^{\prime }}\mathrm{\Delta }{\mathbf{\Omega }}^{m^{\prime }}+S_{c,P}^m\}\mathrm{\Delta }{\mathbf{\Omega }}^m{V}_P,$$

where

Eq. (8)

$$S_{c,P}^m={\mu }_{sP}{\tilde{p}}^{m_c\to m}{\psi }_c(P,{\mathbf{\Omega }}_c),$$

and ${\tilde{p}}^{m^{\prime }\to m}$ represents the average part of the scattered energy from the control solid angle $\mathrm{\Delta }{\mathbf{\Omega }}^{m^{\prime }}$ toward the control solid angle $\mathrm{\Delta }{\mathbf{\Omega }}^m$.²⁰

Spatial discretization of the transport term of the radiance

To solve Eq. (7), closure relations between the ${\psi }_{i_f}^m$ integration-point values and the nodal values of the radiance are required (Fig. 1). The directional nature of radiative transfer needed to be taken into account in order to establish the closure relations. Thus, for a specific direction of propagation of light, only the nodal values located upstream from the integration point had to be considered. In the same way, as presented in Ref. 20, a locally rigorous one-dimensional integration of the RTE along the optical path $(u_f,i_f)$ has to be used.

Eq. (9)

$$\psi (i_f,\mathbf{\Omega })=\psi (u_f,\mathbf{\Omega })\exp [-{\int }_{u_f}^{i_f}{\mu }_t(s)\mathrm{d}s]+{\int }_{u_f}^{i_f}\{S_c(s,\mathbf{\Omega })+{\mu }_s(s){\int }_{{\mathbf{\Omega }}^{\prime }=2\pi }p({\mathbf{\Omega }}^{\prime }\to \mathbf{\Omega })\psi (s,{\mathbf{\Omega }}^{\prime })\mathrm{d}{\mathbf{\Omega }}^{\prime }\}\exp [-{\int }_s^{i_f}{\mu }_t(u)\mathrm{d}u]\mathrm{d}s.$$

Here, $u_f$ and $i_f$ are assumed to be on the same optical path of $\mathbf{\Omega }$ direction with $u_f$ located upstream from $i_f$ and ${\mathrm{\Delta }}_{s_f}=|u_f-i_f|$. The following approximation can be used, $\psi (s,{\mathbf{\Omega }}^{\prime })\approx \psi (u_f,{\mathbf{\Omega }}^{\prime })$ $s\in (u_f,i_f)$.

To compact the equations, the following notations are introduced:

Eq. (10)

$$D_{u_f}^m=A_f{\mathrm{\Delta }}_f^m\exp [-{\int }_{u_f}^{i_f}{\mu }_t(s)\mathrm{d}s];E_{u_f}^m=A_f{\mathrm{\Delta }}_f^m{\int }_{u_f}^{i_f}{\mu }_s(s)\exp [-{\int }_s^{i_f}{\mu }_t(u)\mathrm{d}u]\mathrm{d}s{C}_{u_f}^m=A_f{\mathrm{\Delta }}_f^m{\int }_{u_f}^{i_f}{S}_c^m(s)\exp [-{\int }_s^{i_f}{\mu }_t(u)\mathrm{d}u]\mathrm{d}s.$$

From Eqs. (2) and (8), an approximate value of $C_{u_f}^m$ can be deduced.

Eq. (11)

$$C_{u_f}^m\approx A_f{\mathrm{\Delta }}_f^m\exp [-{\int }_{s_w}^{i_f}{\mu }_t(u)du]{\tilde{p}}^{m_c\to m}|i_f-u_f|\frac{1}{2}\{{\mu }_{s_{i_f}}\mathrm{\Upsilon }(s_w)+{\mu }_{s_{u_f}}\mathrm{\Upsilon }(s_w)\}.$$

The computation of integrals in Eq. (10) is detailed in Appendix B. Combining Eqs. (9) and (10) in Eq. (6) gives

Eq. (12)

$$\sum _{f=1}^{N_f}{\psi }_{i_f}^m{A}_f{\mathrm{\Delta }}_f^m=\sum _{f=1}^{N_f}\{{\psi }_{u_f}^m{D}_{u_f}^m+C_{u_f}^m+\left[\sum _{m^{\prime }=1}^{N_{\theta }}{\tilde{p}}^{m^{\prime }\to m}{\psi }_{u_f}^{m^{\prime }}\mathrm{\Delta }{\mathbf{\Omega }}^{m^{\prime }}\right]E_{u_f}^m\},$$

where $\theta$ is the azimuthal angle (rad). The projections and linear interpolations presented in Ref. 30 were expected to improve the closure relations and the accuracy of the results. Also, they were used and generalized in this study to link ${\psi }_{u_f}^m$ defined at the $u_f$ point with the nodal values of the radiance (see Appendix A).

Full discretization of the RTE

The formula for the fourth-order Runge-Kutta (RK4)³¹ was used for the time semi-discretization of Eq. (4) combined to discrete Eq. (7). We also tested the formula for the Runge-Kutta-Fehlberg method (RK45), which is known to be more precise than RK4, but the results were similar.

An iterative procedure is required to solve equations in the time domain (and also in the steady state). Inspection of these equations shows that they can only be solved for one prescribed $m$ discrete direction. To evaluate the sums over $m^{\prime }$ in Eqs. (7) and (12), it is necessary to know, as preliminary at time $t_k$ (respectively at iteration $k$ in the steady state), the specific radiances ${\psi }_P^{k,m^{\prime }}$ and ${\psi }_{u_f}^{k,m^{\prime }}$ in all discrete directions $m^{\prime }\in \{1,\dots ,N_{\theta }\}$. In our case, these specific radiances were changed by the same one computed at the previous iteration $k-1$ since the values are known at time $t_{k-1}$ (respectively at iteration $k-1$ in the steady state). The convergence criterion was related to a relative difference dealing with the values of the outgoing flux at the bounding surface of the medium. In all results presented further, the iterations have been achieved with a relative difference that does not exceed a prescribed tolerance set to ${10}^{-5}$. Moreover, our numerical method is designed to take into account mismatched interfaces as described in Appendix C. The MFVM described above was coded using C language. A particular effort has been done over the optimization of data structures and computations in the sequential code. As all real numbers manipulated in the code are coded in double precision to ensure a high accuracy of the results, a trade-off has been required between memory consumption and computational time in the context of sequential execution on a single computer. However, this constraint may be at least partially relaxed in a context of cluster computing, implying the use of several computers and thus, a larger aggregated memory. Moreover, even in the context of a single computer, there is a great potential of computational time reduction by using the multiple cores present in current computers and/or GPUs.

3.

Results and Discussion

In this section, the obtained solutions with the MFVM are validated against the analytical approach and the MC technique. The Henyey-Greenstein phase function has been used to account for the anisotropic nature of scattering by the medium. This function, which depends only on the scalar product between the incident direction ${\mathbf{\Omega }}^{\prime }$ and that scattered $\mathbf{\Omega }$ and the anisotropy factor $g$, is expressed for 2-D media as³²

In all the cases presented further, $g$ is set to 0.9, which corresponds to highly forward-scattering media, such as biological tissue. The refractive index of the outside medium (air) is equal to 1. In the first three cases, the steady state is considered and the bounding surface of the medium at $x=0\mathrm{mm}$ was illuminated with a perpendicular collimated Gaussian beam having a spatial intensity of

3.1.

Homogenous Slab (Case 1)

The first test case, taken from Ref. 22, deals with a homogenous slab. The boundaries of the medium are assumed to be transparent. The optical properties of the medium [Fig. 2(a)] are typical for biological tissue in the near-infrared spectral range, where the absorption coefficient is particularly low. It should be noted that the medium is optically thick ($\tau =$10.01), and thus the regime of multiple scattering has to be considered [$l_s\ll L<l_a$, where $L$ is the thickness of the medium (mm), $l_s$ is the scattering mean free path ($=1/{\mu }_s$), and $l_a$ is the absorption mean free path ($=1/{\mu }_a$)]. To correctly take into account the first scattering events in the medium ($l_s=0.1\mathrm{mm}$), the first node of the grid [according to the $(Ox)$ axis] from the bounding surface of the medium at $x=0\mathrm{mm}$ was located at 0.033 mm for mesh 1, 0.025 mm for mesh 2, and 0.016 mm for meshes 3 and 4.

Fig. 2

Physical model.

The steady-state spatially resolved reflectance (at $x=0\mathrm{mm}$) and transmittance (at $x=1\mathrm{mm}$) are shown in Fig. 3(a). Due to the relatively small thickness, the magnitudes for reflectance and transmittance decrease strongly when spatial position $y$ increases. It can be seen that the numerical solution given by the MFVM is in good agreement with the exact analytical solution.²² False scattering and ray effect cause mainly the relative differences between the two solutions, and a potential way to reduce it is to use a spatial mesh refinement [as illustrated in Fig. 3(b)] accompanied by a finer angular quadrature. The differences are higher for transmittance. The relative differences for reflectance are $<2\%$ (meshes 1 and 2) and 1% (mesh 3), respectively. For transmittance, they are $<\sim 4.3\%$ (mesh 1), 3.4% (mesh 2), and 2.4% (mesh 3), respectively.

Fig. 3

Steady-state spatially resolved reflectance and transmittance (case 1).

The diffuse fluence inside the medium is depicted in Fig. 4, and we can seen a symmetry plane $(y=0)$ mm and forward light scattering, which is in agreement with the value of $g(=0.9)$. The diffuse fluence is concentrated in the region between $y=-1\mathrm{mm}$ and $y=1\mathrm{mm}$. These imply that the relative difference for transmittance has peaked regions [Fig. 3(b)] between $y=0\mathrm{mm}$ and $y=1\mathrm{mm}$. A uniform angular quadrature was used in this work. A potential way to reduce the relative difference for transmittance is to use a nonuniform angular quadrature and to refine it near the zone of the collimated direction.

Fig. 4

Steady-state diffuse fluence computed with the modified finite volume method (case 1).

3.2.

Refractive Index Mismatch with Fresnel Reflection at the Interface (Case 2)

The medium studied in this section has a refractive index equal to 1.4 and the boundaries of the medium are assumed to be semitransparent (see Appendix C). The other parameters are the same as in case 1. The steady-state spatially resolved reflectance (at $x=0\mathrm{mm}$) and transmittance (at $x=1\mathrm{mm}$) are depicted in Fig. 5(a). Compared to case 1, the magnitudes for reflectance and transmittance are lower for $y$ near to 0. It is due to the reflected light at the boundaries of the medium directed back into the medium. At the same time, it contributes to give higher magnitudes for reflectance and transmittance when $y$ increases. The numerical solution given by the MFVM compares well with the MC solution⁶ taken as the reference solution in this case. Also, the relative differences between the two solutions [Fig. 5(b)] decrease with spatial mesh refinement, and the differences are higher for transmittance. The relative differences for reflectance are almost always $<1\%$ (meshes 1, 2, 3) except for $y>2.4\mathrm{mm}$, but this is due to the statistical errors of the MC method, which increase with larger $y$ values. For transmittance, they are $<\sim 4.1\%$ (mesh 1), 3% (mesh 2), and 2.2% (mesh 3), respectively.

Fig. 5

Steady-state spatially resolved reflectance and transmittance (case 2).

The number of simulated photons with MC was 3 ${10}^8$ and the computational time was 112 min for this case. The computational times required by the MFVM were, respectively, 31 min (mesh 1), 80 min (mesh 2), and 131 min (mesh 3).

3.3.

Two-Layered Medium (Case 3)

A two-layered medium having a relatively thin first layer of thickness 0.1 mm and a second layer of thickness 0.9 mm is studied. The boundaries of the medium are assumed to be transparent. The values of optical properties of the medium are shown in Figs. 2(b) and 6. In Fig. 6, the steady-state spatially resolved reflectance (at $x=0\mathrm{mm}$) are presented for three different absorption coefficients of the first layer (${\mu }_{a1}=0.01$, 1, and $10{\mathrm{mm}}^{-1}$). As expected, the reflectance decreases when the absorption coefficient of the first layer increases. The results show a good agreement between the numerical solution given by the MFVM and the MC solution.⁶ The relative differences between the two solutions are higher for ${\mu }_{a1}=10{\mathrm{mm}}^{-1}$ and lie between 0 and 5% [Fig. 6(b)]. For the two other absorption coefficients of the first layer, the relative differences are globally $<2\%$.

Fig. 6

Steady-state spatially resolved reflectance (case 3).

3.4.

Results in the Time Domain (Case 4)

The last test case was chosen to evaluate the performance in the time domain of the numerical scheme of the RTE. The optical properties of the homogeneous medium and the boundaries of the medium are the same as in case 1. The Courant-Friedrich-Levy (CFL) condition needs to be satisfied, namely, $\mathrm{\Delta }t\le (n/c)\mathrm{\Delta }{s}_{\text{min}}$, where $\mathrm{\Delta }{s}_{\text{min}}$ is the size of the smallest cell of the spatial mesh.¹³ In other words, since a light beam always travels with a velocity $c/n$ corresponding to the speed of light in the medium, the traveling distance between two consecutive time steps should not exceed the size of the smallest cell. The CFL condition was satisfied in our simulations by setting $\mathrm{\Delta }t=0.01\mathrm{ps}$ for mesh 3. As initial condition ($t=0$), there is no light (or photons) in the medium before light impingement. At subsequent times, the bounding surface of the medium at $x=0\mathrm{mm}$ was illuminated with a perpendicular collimated Gaussian beam in space and in time.

Eq. (15)

$$\mathrm{\Upsilon }(y,t)=\tilde{\mathrm{\Upsilon }}(y)\widehat{\mathrm{\Upsilon }}(t)\text{for}t\ge 0,$$

where $\tilde{\mathrm{\Upsilon }}$ is given by Eq. (14), and

Eq. (16)

$$\widehat{\mathrm{\Upsilon }}(t)=\frac{1}{{ϵ}_t\sqrt{2\pi }}\exp [-\frac{{(t-t_c)}^2}{2{ϵ}_t^2}]\text{for}t\ge 0,$$

where $t_c$ is a constant positive time value. For the simulations presented further, ${ϵ}_y=0.3\mathrm{mm}$, ${ϵ}_t=5\mathrm{ps}$, and $t_c=25\mathrm{ps}$. The curves of time-resolved reflectance (at $x=0\mathrm{mm}$) are presented for different spatial positions, namely, $y=0.5,1.5,2.5,3.5,4.5\mathrm{mm}$ (Fig. 7). A close shape of the Gaussian pulse in time [Eq. (16)] is found. The curves of reflectance decrease with the increase of $y$-coordinate due to the light propagation. An agreement between the MFVM-RK4 and MC solutions can be seen in Fig. 7(a). However, the relative differences in the time domain are higher than in the steady state. The relative differences increase with the $y$-coordinate and for long times. They are related to the time discretization error of the RTE but also to the statisical error in the MC simulation.

Fig. 7

Time and spatially resolved reflectance (case 4).

The number of simulated photons with MC was $4.1·{10}^7$ and the computational time was 148 min. The computational time required by the MFVM-RK4 to carry out 7000 iterations (direct solution in the time domain) was approximately five days.

4.

Conclusion

An MFVM for predicting light propagation within 2-D absorbing and highly forward-scattering medium subjected to a collimated Gaussian light beam has been presented. The accuracy of the MFVM has been examined by comparing the results with an exact analytical solution of the RTE and the MC technique. Numerical results for the spatially resolved reflectance and transmittance over different situations (refractive index mismatch with Fresnel reflection at the interface, homogeneous and two-layered media) have shown that the MFVM yields accurate results. In the steady state, relative differences of the different solutions were found to decrease with spatial mesh refinement obtaining $<2.4\%$ compared to analytical and MC solutions. In the time domain, the RK4 has been used for the time semi-discretization of the RTE. An agreement between the MFVM-RK4 and MC solutions has been shown, but the relative differences were found to be higher than in the steady state.

In our future work, we are interested to implement another technique for the time semi-discretization of the RTE as recently proposed in Ref. 33. The technique is based on solving the RTE in the frequency domain at multiple modulation frequencies. It has been reported by the authors that using the Fourier-series representation of the radiance, the solution can be obtained in the time domain with good accuracy and with significantly fewer computational resources than are needed in the direct solution in the time domain. Moreover, the MFVM presented in this paper will be extended to 3-D geometry for biomedical applications, such as diffuse optical tomography for tumor diagnosis. Several aspects of parallel computing for 2-D/3-D MFVM code will be explored to speed up the computation, such as the multithreaded computation (the use of multiple cores present in the same computer), the potential use of GPUs, and the larger scale of parallelism making use of multiple computers in order to manage larger and more complex problems. Our preliminary studies on the parallel computing of the MFVM over direction (in our simulations, typically 128 discrete directions were used) show that we can expect significant gains in computing speed.