2.1.

Radiative Transfer Equation

Based on the setting of our mesoscopic imaging system (continuous wave),^13–15 in this paper, we investigate the quasisteady state of the RTE. Here, the photon propagation in the media is described through the RTE as^23,28

The boundary condition of Eq. (1) is defined so that no photon will travel into the object through the boundary^1,21 except the region with irradiation by external light, that is,

In Eq. (2), $S_0$ is the region on the boundary $\partial X$ where the external light source exists or the region that is irradiated by an external light source (diffuse source³⁹), and $n$ is the outward unit normal vector at the boundary $\partial X$.

2.2.

Discretization of Angular Space

Following the DOM to discretize the angular space, a numerical quadrature is applied first to compute the integral in (1), and this gives

In Eq. (3), $w_i$ represents the weight in a discrete direction $s_i$. The selection of discrete points on the unit sphere determines the accuracy of the solution to the RTE. Various numerical quadratures have been proposed^34,40–42 in order to estimate the integral in Eq. (3) accurately. Some adaptive methods have also been developed,⁴³ such as the ordinate splitting technique.⁴¹ In this paper, several widely employed quadratures, including level symmetric quadrature,⁴⁰ product Gaussian quadrature,^34,44 Legendre equal-weight quadrature,⁴¹ and Lebedev quadrature,^45,46 were investigated numerically. The reader can refer to the corresponding papers for details on the construction of these quadratures, advantages, and drawbacks. Here we want to emphasize that even though Lebedev quadrature has been reported to have the highest degree of accuracy and efficiency,⁴⁷ it still has not been widely used in optical imaging. Figure 1 shows the distribution of various quadrature points in the first octant.

Fig. 1

Distribution of discrete points on unit sphere $\mathrm{\Omega }$ of four numerical quadratures (first octant): (a) level symmetric quadrature ($N=80$); (b) product Gaussian quadrature ($N=72$); (c) Legendre equal-weight quadrature ($N=80$); (d) Lebedev quadrature ($N=86$).

All the numerical quadratures mentioned above work well when the integrands are continuous and smooth. However, the H–G function becomes singular when $s·s^{\prime }$ is around 1 and the anisotropy factor $g$ is close to 1. Note that most biological tissues exhibit large $g$. In such cases, a small number of quadrature points can cause a large error during the calculation of the integral at the right-hand side of Eq. (1). To reduce such an error, one can increase the number of quadrature points, at the expense of using more computational resources such as storage space and time, hence decreasing the efficiency of the algorithm. On the other hand, none of the numerical quadratures reviewed above preserve the key properties [see Eq. (4)] of the phase function, and this can lead to large errors in numerical solutions, as demonstrated in Sec. 4.3. All these motivate us to examine the phase function normalization technique, which is to further modify the numerical quadrature of our choice and will be discussed next.

2.3.

Phase Function Normalization Technique

The phase function normalization technique^35,36 is a method that preserves the conservative properties or laws of the phase function during the evaluation of the integral term in Eq. (1) with numerical quadrature by adjusting the quadrature weights. The following properties hold for the H–G phase function:

The first property is from the requirement that the scattered energy be conserved,³⁵ while the second one can be viewed as the definition of anisotropy $g$. Generally, for any numerical quadrature discussed above, especially with large $g$, the discretized form of the left-hand side of Eq. (4) is not equal to the right. A small difference between the two sides of Eq. (4) can generate a large error in the solution to (1); see numerical experiments in Sec. 4.3. The situation can worsen with large absorption and scattering coefficients, since in these cases numerical methods may even generate negative solutions that are nonphysical. Readers can refer to Refs. 35, 48 for more details on this technique.

Briefly, to apply the phase function normalization technique, we start with the original discrete points $s_i$ and weights $w_i$ of a chosen numerical quadrature; refer to Refs. 34, 40, 41, and 45 for how to calculate $s_i$ and $w_i$. Let ${\theta }^{ij}$ (based on the original discrete points) represent the angle between the $i^{\prime }$th direction $s_i$ and the $j^{\prime }$th direction $s_j$, and let $f(\cos {\theta }^{ij})$ be the value of the phase function at ${\theta }^{ij}$. Namely, $f(\cos {\theta }^{ij})=f(s_i,s_j)$ with the specific form of the H–G phase function. Then the normalized phase function, with its value at the same point ${\theta }^{ij}$ denoted as $f_{\text{norm}}^{ij}$, satisfies the following discrete form of Eq. (4):

Here, the values of the normalized and the original phase function at ${\theta }^{ij}$ are connected as

Plugging Eq. (6) into Eq. (5), a system of linear equations will be obtained, with ${\xi }^{ij}$ being the unknowns. Notice that there are $2N$ conditions in Eq. (5) while the number of unknowns is $N^2$. If it is further required that the normalized phase function be symmetric such that ${\xi }^{ij}={\xi }^{ji}$, the number of unknowns will be reduced to $(N^2+N)/2$. In general, there are more unknowns than equations (when $N>3$), hence there will be an infinite number of solutions. In order to deal with such an equation, a least squares solution is utilized, which minimizes the two-norm of the solution vector. A MATLAB^® built-in function lsqr.m is employed to solve Eq. (5). The sparse matrix technique is also used here for better efficiency in computation and storage.⁴⁸ Through normalization, the newly calculated weights are expressed as ${\overline{w}}_{ij}=w_i(1+{\xi }^{ij})/\sum _i{w}_i(1+{\xi }^{ij})$.

2.4.

Discrete Ordinate Method with Continuous Galerkin Finite Element Method

After discretizing the integral term in Eq. (1) using a modified numerical quadrature based on the phase function normalization technique, the RTE is turned into a system of coupled first-order PDEs. That is,

To solve Eq. (7), the CG-FEM with streamline diffusion modification of the test function is implemented.^33,37 For a given tetrahedron-based mesh $\mathcal{T}$, a finite element space is defined as

Here, $P^1(K)$ is the set of linear polynomials defined in $K$, and $C(X)$ denotes the space of the continuous functions in the closure of $X$.

Following the idea of the streamline diffusion method with a modified test function $v(r)+\delta s_j·\nabla v(r)$, the following CG-FEM method will be formed: look for ${\phi }_j(r)\in V_{\mathrm{h}}$, such that for any $v(r)$ and any $j=1,2,\dots ,N$,

Note that ${\phi }_j(r)\in V_{\mathrm{h}}$ approximates the solution $\phi (r,s_j)$ and $q_j(r)=q(r,s_j)$. The parameter $\delta$ is related to the absorption and scattering coefficients,³³ and it is chosen to be a piecewise constant, defined as $\delta {|}_K=C{h}_K/(1+h_K{\mu }_{\mathrm{t}})$, with $h_K$ being the diameter of a mesh element $K\in \mathcal{T}$. $C$ (0.3 to 0.6) is a constant, and the computed solution is not sensitive to its value for the range of absorption and scattering coefficients examined in this work on reasonably refined meshes. Our formula for $\delta$ satisfies the guiding principle for choosing this parameter, as outlined in Theorem 1 from Ref. 33. Performing integration by parts to the first term of Eq. (8), we obtain

Our final scheme is obtained with Eq. (9) plugged into Eq. (8).

Let $r_k$ denote the $k$’th vertex of the mesh. Considering the Lagrange nodal basis ${\{{\psi }_l^j(r)\}}_l$ of $V_{\mathrm{h}}$ associated with vertices, satisfying ${\psi }_l^j(r_k)={\delta }_{lk}$, in which ${\delta }_{lk}$ is the Dirac-Delta function. Each basis function ${\psi }_l^j(r)$ is nonzero only within those elements sharing the vertex $r_l$. The numerical solution ${\phi }_j(r)$ and the test function $v(r)$ can now be expanded in terms of basis, namely,

In Eq. (10), $M_j$ represents the number of nodes (basis function) in discrete direction $s_j$. Since the linear polynomial space is employed as discrete space, the gradient of the test function $v(r)$ will be constant in each element and can be calculated analytically. With these assumptions, the numerical method can be converted into a linear system, namely $Ax=q$. In the linear system, matrix $A$ contains the terms related with the stiffness matrix, the mass matrix and so on, $x={(a_1^1,a_2^1,\dots ,a_{M_1}^1,a_1^2,\dots ,a_{M_2}^2,\dots ,a_{M_N}^N)}^T$. $q$ relates with the source term, which can be acquired either by employing the boundary condition in Eq. (2) or directly processing the source term²⁵ $q(r,s)$. The linear equation was then solved using the generalized minimum residual method (GMRES),²² which is performed in the simulations with the MATLAB^® built-in function gmres.m. We want to point out that the source iteration method or improved source iteration method can also be employed.²³

When the object and ambient media have different refractive indices, the reflection and transmittance are modeled by Fresnel’s law.^23,25 In the simulations, the mismatched boundary condition will be applied when the refractive indices of the object and the ambient media are different.

2.5.

Diffusion Equation

Because of the computational difficulty of solving Eq. (1), an approximation model of the RTE, the DE, has been widely employed to simulate photon propagation in tissue.^39,49 Briefly, photon radiance $\phi (r,s)$ is expanded in terms of spherical harmonic functions. Then, the truncated spherical harmonic series, with a different number of terms, is used to approximate the RTE. For example, the DE can be viewed as the first order, $P_1$ approximation¹ (in which case, one uses spherical harmonic functions whose index is 0 and 1). The quasisteady DE can be expressed as³⁹

Here $\kappa (r)=1/[3({\mu }_{\mathrm{a}}+{\mu }_{\mathrm{s}}^{\prime })]$, where ${\mu }_{\mathrm{a}}$ is the absorption coefficient (${\mathrm{mm}}^{-1}$), while ${\mu }_{\mathrm{s}}^{\prime }$ is the reduced scattering coefficient (${\mathrm{mm}}^{-1}$) and defined as ${\mu }_{\mathrm{s}}^{\prime }=(1-g){\mu }_{\mathrm{s}}$. $q(\mathbf{r})$ represents the light source.

In this work, the Robin boundary condition is applied, that is,

Equation (11) is a second-order elliptic PDE with the Robin boundary condition. The CG-FEM is a well-established approach for solving such equations, and a standard finite element procedure can be applied. In our simulation, the continuous linear finite element space is employed and the corresponding matrices (stiffness matrix, mass matrix, surface integral related matrix, and load vector) are formed as in Refs. 49, 50, 51.

2.6.

Monte–Carlo Method

MC is a stochastic method in which a large number of photons are transported to obtain the photon distribution (photon density) in medium. This procedure utilizes a probability function to describe the stochastic events of photons, such as energy dissipation during collisions with media particles and the new direction of propagation when a photon hits the particles in media. It involves the absorption coefficient, scattering coefficient, and phase function for computation. Some authors have investigated the MC method for higher computational efficiency.⁵² Mesh-based MC methods have also been proposed¹² to deal with complex shapes of media. In our simulation, mesh-based MC was employed to allow comparisons of identical meshes^12,53 with solutions to RTE. Results obtained using the MC method are viewed as the ground truth standard to evaluate the accuracy of solutions to the RTE and DE.

3.1.

Ideal Pencil Beam Simulation

In our simulations, a 3-D rectangular computational domain $[-2.5,2.5]\times [-7.0,7.0]\times [-5.0,0]{\mathrm{mm}}^3$ is considered first, referred to as Rect. 1, which is similar to the setting⁵⁴ we previously used to simulate the photon propagation in media (teeth) with the MC method. A computational mesh was generated, as shown in Fig. 2(a), using the precompiled Computational Geometry Algorithms Library. The 3-D domain has 4858 nodes with 27,456 elements. The position where the photon was launched was $r_{\mathrm{s}}=(\mathrm{0,0},0)$ with the direction $s_0=(\mathrm{0,0},-1)$ orthogonal to the plane $z=0\mathrm{mm}$; see Fig. 2(b). This type of source was modeled by the Dirac-Delta function $\delta (r-r_s)\delta (s-s_0)$.

Fig. 2

The 3-D rectangle used in simulations: (a) 3-D mesh of rectangle; (b) the slice at $x=0\mathrm{mm}$, where the arrow indicates the position in which the photon is launched perpendicularly to the plane $z=0\mathrm{mm}$.

The optical properties of the 3-D rectangle were selected according to the most common values of human tissue⁵⁵ (see Tables 3 and 4). For instance, the values for the healthy human brain are less than $0.1{\mathrm{mm}}^{-1}$ for the absorption coefficient and around $10{\mathrm{mm}}^{-1}$ for the scattering coefficient⁵⁶ at 674 nm. The anisotropy of the tissue was set to $g=0.9$ and the refractive index $n$ was 1.5 for all simulations (pencil beam and Gaussian shape beam). Various combinations of optical properties (absorption and scattering coefficients) were tested. As mentioned above, results from the RTE were compared with the solution to the DE and MC simulations under the same settings. A tetrahedron mesh-based MC method⁵³ was employed. Since using different numbers of photons to simulate photon propagation in tissue by MC led to different absolute intensities, results were normalized with respect to a global maximum value²⁸ in order to compare the solutions obtained by the RTE, DE, and MC methods. Equation (12) was employed to transform the photon radiance to photon density by numerical integration to compare MC simulations and solutions to the DE. In order to easily process the data, a structural grid ($49\times 139\times 50$ with a resolution of 0.1 mm in each dimension) was generated. For the structural points inside one element, a linear interpolant through the values of four nodes of the element was employed.

Solutions to the RTE obtained by different numerical quadratures were compared to MC simulations in order to appraise their influence on the results. Since it was difficult to use the same number of discrete angles of quadratures with different constructions, comparable numbers of discrete angles were employed. Subsequently, solutions obtained with the most accurate quadrature (from the discussion in Sec. 4.1.1, Lebedev quadrature was mainly employed) for the RTE were compared with the results from the DE and MC simulations. Due to the limit of the storage space (RAM) of the computer used for simulation, no more than 100 angles of the DOM were employed. It is noted that the number of angles used in the simulation was less than those in Refs. 26 and 57 when the FVM was used with similar accuracy. The simulations were run on a computer with Intel Xeon 2.8GHz CPU, 24 GB RAM, and Windows 7 Professional 64-bit operating system. The number of photons in these MC simulations was set to ${10}^6$.

3.2.

Gaussian Shape Beam Simulation

Many light sources can be modeled by Gaussian distribution. In our next group of tests, a Gaussian-shaped beam, in which the intensity of the points at the cross-section of the beam obeys the Gaussian distribution, was employed as the light source. The beam intensity was modeled by

The anisotropy and refractive index were set to be the same values as in Sec. 3.1. The parameter $d$ was set to 0.5. In the first example (Rect. 1), the absorption coefficient was set to $0.08{\mathrm{mm}}^{-1}$ and the scattering coefficient was set to $5{\mathrm{mm}}^{-1}$. Since the light source was modeled by a Gaussian function, it was much smoother than the ideal pencil beam, thus fewer discrete angles were needed to obtain solutions with comparable accuracy to those gathered for the ideal pencil beam. For Rect. 1, a mesh with denser nodes was generated compared with the one in Sec. 3.1, and in total there are 7204 nodes with 41,559 elements.

For Rect. 2 [Fig. 3(a)], 10,996 nodes with 60,887 elements were generated. Figure 3(b) shows the Gaussian beam represented by a collection of single point sources for MC in the computation. The optical properties of the simulation were set to be $0.1{\mathrm{mm}}^{-1}$ for absorption coefficient and $5{\mathrm{mm}}^{-1}$ for scattering coefficient, respectively.

Fig. 3

3-D rectangle used in simulation 3.2: (a) mesh of 3-D rectangle; (b) irradiated light intensity modeled by Gaussian shape function.

Two quantities were used to assess the difference between the solution to the RTE and MC simulation, as shown below. One is the root mean square error between the solutions, defined as

Based on this, the mean relative error (MRE, defined as $\sum _i{e}_i/N$) and maximum relative error can be easily calculated. Note that both quantities can be calculated for any geometry, along any line in a specific plane, or in the entire 3-D rectangle.

4.1.

Three-Dimensional Rectangle Simulations with Ideal Pencil Beam

4.1.2.

Ideal pencil beam simulations

Figure 4 shows the contours of the logarithm of solutions to the RTE by the streamline diffusion modified CG-FEM and from MC simulations within three planes. The optical properties for these simulations were set as ${\mu }_{\mathrm{a}}=0.02{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}=5{\mathrm{mm}}^{-1}$. With this combination of optical properties in tissue, the DE is not an accurate model to predict photon propagation.²⁸ Eighty-six angles were employed in Lebedev quadrature to obtain the results for the RTE. It is clear that at deeper depths, the difference between the solutions to the RTE and MC simulations is small. The photon densities along the line $x=0$ mm with different depths are shown in Fig. 5. The RMSE and MRE along each line at different depths are shown in Tables 3 and 4.

Fig. 4

The contours of logarithm of solutions to RTE and MC results, solid curve for the solution to RTE and dashed curve for MC results: (a) the contours within the plane $x=0\mathrm{mm}$; (b) the contours within the plane $y=0\mathrm{mm}$; the value of the outermost curve is $-1.5$; (c) the contours within the plane $z=-3.0\mathrm{mm}$.

Fig. 5

Comparison between the solutions to RTE, DE, and MC method at different depths, along the line $x=0\mathrm{mm}$: (a) $z=-0.4\mathrm{mm}$; (b) $z=-1.0\mathrm{mm}$; (c) $z=-2.0\mathrm{mm}$; (d) $z=-3.0\mathrm{mm}$; (e) $z=-4.0\mathrm{mm}$; (f) $z=-5.0\mathrm{mm}$.

Table 3

RMSE for different optical properties at different depth along the line x=0  mm (g=0.9).

(μa,μs)  mm−1	−0.4  mm	−1.0  mm	−2.0  mm	−3.0  mm	−4.0  mm	−5.0  mm
(0.02, 5)	$6.84\times {10}^{-2}$	$9.07\times {10}^{-2}$	$1.03\times {10}^{-2}$	$1.01\times {10}^{-3}$	$4.19\times {10}^{-4}$	$2.89\times {10}^{-4}$
(0.02, 8)	$8.88\times {10}^{-2}$	$2.58\times {10}^{-2}$	$3.25\times {10}^{-3}$	$1.02\times {10}^{-3}$	$5.69\times {10}^{-4}$	$2.23\times {10}^{-4}$
(0.08, 5)	$7.08\times {10}^{-2}$	$9.16\times {10}^{-2}$	$1.02\times {10}^{-2}$	$9.49\times {10}^{-4}$	$2.80\times {10}^{-4}$	$1.97\times {10}^{-4}$
(0.01, 10)	$8.92\times {10}^{-2}$	$2.33\times {10}^{-2}$	$2.90\times {10}^{-3}$	$7.40\times {10}^{-4}$	$3.24\times {10}^{-4}$	$1.94\times {10}^{-4}$

Table 4

MRE for different optical properties at different depth along line x=0  mm (g=0.9).

(μa,μs)  mm−1	−0.4  mm	−1.0  mm	−2.0  mm	−3.0  mm	−4.0  mm	−5.0  mm
(0.02, 5)	17.76%	17.24%	8.25%	4.73%	3.84%	3.59%
(0.02, 8)	11.94%	7.31%	4.10%	3.19%	3.46%	2.71%
(0.08, 5)	19.51%	18.98%	9.22%	5.05%	3.72%	3.44%
(0.01, 10)	13.80%	7.99%	5.54%	3.54%	3.01%	3.09%

In comparing MC to the DE, it is quite clear that the discrepancy between the MC solution (blue-dashed line in Fig. 5) and the DE solution (green line in Fig. 5) is large. In comparing MC to the RTE, the solutions to the RTE and MC coincide quite well except in shallow regions [Figs. 5(a) and 5(b)]. The MRE for $z=-3.0,z=-4.0,z=-5.0\mathrm{mm}$ is less than 5% (Table 4). At these depths, the maximum relative error for all points is 10.04%, 12.03%, and 9.06%, respectively. In the shallow region, for instance, $z=-1.0\mathrm{mm}$ and $z=-2.0\mathrm{mm}$, the errors are larger as ROI gets closer to the light source located at (0,0,0). At the extreme shallow region, $z=-0.4\mathrm{mm}$, the difference between the solution to the RTE and the results of the MC method is large, which leads to large RMSE and MRE.

To further demonstrate the performance of the algorithms, three more combinations of optical properties were tested. The RMSE and MRE between the solutions to the RTE and the results of the MC along line $x=0\mathrm{mm}$ are shown in Tables 3 and 4. For all examples tested and for depths smaller than $-3.0\mathrm{mm}$, the MRE is less than or close to 5%. At a depth of $z=-2.0\mathrm{mm}$, the MRE is less than 10%.

Tables 3 and 4 show only the results along the line $x=0\mathrm{mm}$ at different depths. For a comprehensive comparison of the solutions to the RTE and results from the MC, 3-D RMSE and MRE ($x\in [-2.5,2.5]\mathrm{mm},y\in [-7.0,7.0]\mathrm{mm}$) were also calculated. The results are shown in Tables 5 and 6.

Table 5

RMSE for different optical properties to 3-D rectangle (g=0.9).

(μa,μs)  mm−1	−5.0≤z<0.0  mm	−5.0≤z<−1.0  mm	−5.0≤z<−2.0  mm
(0.02, 5)	$2.81\times {10}^{-2}$	$1.07\times {10}^{-2}$	$2.16\times {10}^{-3}$
(0.02, 8)	$2.54\times {10}^{-2}$	$6.39\times {10}^{-3}$	$1.81\times {10}^{-3}$
(0.08, 5)	$2.88\times {10}^{-2}$	$1.06\times {10}^{-2}$	$2.02\times {10}^{-3}$
(0.01, 10)	$2.60\times {10}^{-2}$	$5.77\times {10}^{-3}$	$1.58\times {10}^{-3}$

Table 6

MRE for different optical properties to 3-D rectangle (g=0.9).

(μa,μs)  mm−1	−5.0≤z<0.0  mm	−5.0≤z<−1.0  mm	−5.0≤z<−2.0  mm
(0.02, 5)	10.05%	6.82%	5.45%
(0.02, 8)	8.33%	5.51%	4.78%
(0.08, 5)	10.90%	7.42%	5.82%
(0.01, 10)	8.87%	5.76%	4.92%

From Tables 5 and 6, one can see that for each set of optical properties (corresponding to each row in the table), the 3-D RMSE and MRE decrease when the regions where the errors are measured are reduced. In our simulations, the largest error occurs near the source. For the optical properties listed in Tables 5 and 6, the MRE for deeper regions is around 5%.

4.2.

Three-Dimensional Rectangle Simulations with Gaussian Modeled Intensity Beam

The contours of the logarithm of photon densities within three planes of Rect. 1, obtained by solving the RTE and using MC simulations, are shown in Fig. 6. The contours of the logarithm of photon densities of Rect. 2 are shown in Fig. 7. For both Rect. 1 and Rect. 2, Lebedev quadrature with 50 angles was employed. 3-D RMSE and MRE for different regions were also calculated and are shown in Table 7. When the volume of the ROI is reduced, there is a subsequent decrease in both the RMSE and MRE. Because the Gaussian-shaped intensity beam is smoother than the point source, the solution to the RTE contains fewer oscillations near the light source. This can be easily seen through comparison of the ROI near the source in Figs. 4 and 6.

Fig. 6

The contours of logarithm of photon densities of Rect. 1 at three planes, in which the solid curves represent the solution to RTE while the dashed curves represent MC results: (a) $x=0\mathrm{mm}$; (b) $y=0\mathrm{mm}$, the value of the outermost contour is $-1.5$; (c) $z=-3\mathrm{mm}$.

Fig. 7

The contours of logarithm of photon densities of Rect. 2 at three planes, in which the solid curves represent the solution to RTE while the dashed curves represent MC results: (a) $x=0\mathrm{mm}$, the value of outermost curve is $-2.5$; (b) $y=0\mathrm{mm}$, the value of outermost curve is $-2.5$; (c) $z=-3\mathrm{mm}$.

Table 7

3-D RMSE and MRE for two rectangles with Gaussian shape beam (g=0.9).

4.3.

Importance of the Normalized Phase Function Technique

To illustrate the importance of the normalized phase function technique to the stability of the proposed method, especially for large anisotropy with $g$ around 0.8 to 0.9, we compare in Fig. 8 the numerical solutions of the methods with or without using the normalization. Rect. 1 with Gaussian beam was considered, with all the optical properties the same as in Sec. 3.2. The comparison is shown in Fig. 8.

Fig. 8

Comparison between the logarithm of contours with and without phase function normalization: (a) with phase function normalization and (b) without phase function normalization.

When processing the CG-FEM without the phase function normalization, it is observed from Fig. 8(b) that the algorithm is less stable and the computed solution oscillates between positive and negative values. In the plot, the negative value is shown as the machine epsilon in MATLAB^® for logarithm calculation.