|
1.IntroductionOptical tomography is an imaging modality that utilizes multiple light projections through the media to quantitatively retrieve optical properties of interest for functional or molecular imaging.1,2 Accuracy of optical property reconstruction inside the media (the so-called “inverse problem”) is dependent on the forward model employed.3 The most popular light propagation model employed in the field of diffuse optics is the diffusion equation (DE), which is an approximation to the radiative transfer equation (RTE). The DE is well established for applications in which photons experience numerous scattering events. However, it is not appropriate when photons experience no scattering or only a few scattering events. With the rise of applications targeting small volumes or shallow depths,4,5 as well as imaging approaches that leverage low scattered photons (such as early time resolved data sets),6,7 alternative forward models need to be considered. For such imaging scenarios, one can solve the RTE or use a stochastic method like Monte Carlo (MC).8 Even though both approaches can produce comparable results, they differ in their ease of implementation. The MC method is straightforward to implement but computationally more costly, whereas the RTE is computationally fast but difficult to implement. Thanks to massively parallel implementations (GPU or CPU),9,10,11 the MC method has recently become the light model of choice for regimes in which the DE is not appropriate.12,13 It can be used in optical tomography to compute the forward model employed in the inverse problem within reasonable time frames. In mesoscopic fluorescence molecular tomography, the optical forward model can be computed in less than 5 min using voxel-based MC and symmetry in the imaging space.14,15 However, when complex boundaries are concerned, voxel-based techniques are less attractive due to the need for accurate modeling of the sample boundaries. This leads to an increased computational burden. For instance, in time-resolved preclinical optical tomography, full computations based on mesh MC and adjoint formulation can be achieved in around 30 min.12 Additionally, MC methods are not inherently suitable for adaptive image space discretization such as mesh adaptation and nonlinear inverse formulations.16,17 Hence, there is still incentive to develop fast and efficient deterministic RTE-based forward models. There are multiple challenges in solving the RTE analytically or numerically. The RTE can only be solved analytically in very few situations,18,19 such as when semi-infinite media are used.20 In such instances, the forward model can be computed in a matter of seconds. However, for media with complex boundaries, it is extremely difficult or even impossible to derive an analytical solution. Under these conditions, the RTE needs to be solved using numerical methods. Prior to applying the standard numerical methods, it is necessary to first process the angular domain, since the RTE retains the anisotropic nature of light propagation. Over the years, a few methods have been proposed to perform this task. Tarvainen et al.1,21 formed a partition on the angular domain and employed local linear basis functions to approximate the angular components of the solution. Thus, when forming the weak form of RTE by the finite element method (FEM), the high-dimensional integral with respect to the space variables and angular variables had to be calculated. Although this calculation could be transformed into two independent low-dimensional integrals, the tensor product of these two integral-related matrices still needs to be performed.1 In addition to using linear basis functions, Surya Mohan et al.22 also investigated the use of spherical harmonic functions to approximate the angular components of the solution. However, this formulation still required computation of the high-dimensional integral and tensor product of the matrices (e.g., when considering the three-dimensional (3-D) space and two-dimensional unit sphere in angular space, a five-dimensional integral was calculated in general). Another method to handle the angular domain of the RTE is to employ the discrete ordinate method (DOM) or its modified version. The standard DOM23 consists of transferring the RTE to a system of coupled partial differential equations (PDEs) by using numerical quadrature to calculate the integral term (the term that contains the phase function). Klose et al.,24 Fujii et al.,25 and Ren et al.26 employed level symmetric quadrature to directly discretize RTE into a number of different PDEs. The required number of PDEs was determined by the number of discrete angles of quadrature. Employing the modified DOM, Asllanaj et al.27 and Wang et al.28 partitioned the angular domain into many finite subdomains and performed the integration with respect to angular variables on both sides of RTE with the assumption that the solution to the RTE (photon radiance) was constant inside each subdomain. This would lead to easily computed integrals. However, a four-dimensional (4-D) integral, which could not be transformed into lower dimensional integrals, still needs to be computed when dealing with the term of the RTE containing the phase function.28 In comparison to dealing with the angular domain as mentioned above, the DOM is relatively easy to implement, as it does not require calculation of the integral with respect to angular variables.29 Consequently, the dimension of the RTE is reduced, allowing for computation in 3-D or 4-D. The main disadvantage of DOM implementations is the so-called “ray-effect,” which occurs due to angular discretization and is prominent in optical imaging since the anisotropy factor is large ( to 0.9).30 After dealing with the angular domain, multiple methods such as the finite difference method (FDM), the finite volume method (FVM), or the FEM can be employed to numerically solve the RTE. Klose et al.24 and Wang et al.28 employed an upwind FDM to approximately calculate the partial derivatives with respect to spatial variables. The FDM is relatively easy to implement compared with other numerical methods. The main disadvantage of the FDM is its use of structural grids, which cannot easily and accurately approximate complex, curved boundaries (although it is possible to deal with this situation as proposed in Ref. 31). For such scenarios, the FVM can be employed since it can deal with arbitrary shapes.24,27,30 To take into account the directional propagation of light, Asllanaj et al.27 proposed the inclusion of the nodal values located upstream at each integration point. In order to reduce the number of unknowns and have high resolution,27 the control elements were normally rebuilt to surround the nodes. The main limitation of the FVM is that it can be very involved to achieve higher order accuracy on unstructured grids, as the numerical scheme of the higher order FVM requires the use of multiple nearby elements.32 In terms of FEM-based approaches, it is widely known that employing continuous Galerkin FEM (CG-FEM) directly to solve the RTE will generate nonphysical oscillations33 due to the transport effect. Tarvainen et al.1 and Surya Mohan et al.22 proposed solving the RTE by CG-FEM, and in their work, the standard streamline diffusion modification technique was employed. The discontinuous Galerkin FEM (DG-FEM) was also employed by Eichholz34 and Gao and Zhao23 to solve the RTE. Compared with the CG-FEM, the DG-FEM has additional unknowns associated with the relaxation of continuous constraints of functions at the boundary between adjacent elements. However, the DG-FEM allows for a straightforward approach to solving transport equations such as the RTE.32 Herein, based on considerations of the ease of implementation and computational efficiency, an “intermediate” method, founded on the modification of the DOM with the CG-FEM, is proposed. In our work, several improvements were made, and the performance of the methods was illustrated numerically. First, various numerical quadratures were implemented to translate the RTE to a system of coupled PDEs with the DOM. A comparison was performed for us to understand the effect of different numerical quadratures on the solution of the RTE with multiple combinations of optical properties. Second, a phase function normalization35,36 strategy was employed in order to efficiently lessen the instability (oscillations) or “ray-effect” commonly encountered when using the DOM.23 More precisely, during the discretization of the angular space, the conservative properties of the phase function will generally be lost. Without employing large numbers of discrete points in angular space, this will generate spurious oscillations, as well as the “ray-effect” mentioned above. The phase function normalization techniques aim to preserve the conservative properties even after discretization of the angular space so as to reduce the oscillations of solutions with a smaller number of discrete points. Once the angular space is discretized, a streamline diffusion modified CG-FEM was selected and implemented due to its smaller number of unknowns compared with the DG-FEM as well as its ability to robustly deal with the transport nature of the RTE, following the work of Brooks and Hughes.37 The paper is arranged as follows. In Sec. 2, the modified DOM with the CG-FEM is discussed in detail. In addition, the normalized phase function technique, which reduces the error in computing the integral term of the RTE and mitigates the instability of the DOM, is illustrated. Since we want to compare the RTE to the DE, a brief introduction to solving the DE with the FEM is also provided. In Sec. 3, two types of sources as well as 3-D domain simulations are introduced. The quantities used to appraise the accuracy of solutions are also presented in detail. In Sec. 4, solutions obtained using the RTE, DE, and MC methods with different types of sources are compared for several combinations of tissue optical properties. Moreover, the importance of the normalized phase function technique is demonstrated. The last section summarizes our work and discusses further improvements to this method. 2.Algorithms2.1.Radiative Transfer EquationBased 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 as23,28 where the photon radiance () is defined on the domain , with the spatial domain being the entire 3-D Euclidian space or its subspace, and being the unit sphere in (direction or angular domain). In particular, the angular variable is expressed as , where is the polar angle and is the azimuthal angle. represents the attenuation coefficient () of the media, and , where and are the absorption coefficient () and scattering coefficient (), respectively. is the phase function, and it describes the probability that the photon scatters from direction to . In biomedical optics, the phase function in 3-D is usually modeled by the Henyey–Greenstein (H–G) function,23,38 where characterizes the anisotropy of the media, usually in the range of to 0.9 for biological tissue, indicating the highly forward scattering properties of the biological tissue. models the source inside the object.The boundary condition of Eq. (1) is defined so that no photon will travel into the object through the boundary1,21 except the region with irradiation by external light, that is, In Eq. (2), is the region on the boundary where the external light source exists or the region that is irradiated by an external light source (diffuse source39), and is the outward unit normal vector at the boundary . 2.2.Discretization of Angular SpaceFollowing 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), represents the weight in a discrete direction . The selection of discrete points on the unit sphere determines the accuracy of the solution to the RTE. Various numerical quadratures have been proposed34,40–42 in order to estimate the integral in Eq. (3) accurately. Some adaptive methods have also been developed,43 such as the ordinate splitting technique.41 In this paper, several widely employed quadratures, including level symmetric quadrature,40 product Gaussian quadrature,34,44 Legendre equal-weight quadrature,41 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,47 it still has not been widely used in optical imaging. Figure 1 shows the distribution of various quadrature points in the first octant. Fig. 1Distribution of discrete points on unit sphere of four numerical quadratures (first octant): (a) level symmetric quadrature (); (b) product Gaussian quadrature (); (c) Legendre equal-weight quadrature (); (d) Lebedev quadrature (). ![]() All the numerical quadratures mentioned above work well when the integrands are continuous and smooth. However, the H–G function becomes singular when is around 1 and the anisotropy factor is close to 1. Note that most biological tissues exhibit large . 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 TechniqueThe phase function normalization technique35,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,35 while the second one can be viewed as the definition of anisotropy . Generally, for any numerical quadrature discussed above, especially with large , 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 and weights of a chosen numerical quadrature; refer to Refs. 34, 40, 41, and 45 for how to calculate and . Let (based on the original discrete points) represent the angle between the th direction and the th direction , and let be the value of the phase function at . Namely, with the specific form of the H–G phase function. Then the normalized phase function, with its value at the same point denoted as , satisfies the following discrete form of Eq. (4): Here, the values of the normalized and the original phase function at are connected as Plugging Eq. (6) into Eq. (5), a system of linear equations will be obtained, with being the unknowns. Notice that there are conditions in Eq. (5) while the number of unknowns is . If it is further required that the normalized phase function be symmetric such that , the number of unknowns will be reduced to . In general, there are more unknowns than equations (when ), 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.48 Through normalization, the newly calculated weights are expressed as . 2.4.Discrete Ordinate Method with Continuous Galerkin Finite Element MethodAfter 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 , a finite element space is defined as Here, is the set of linear polynomials defined in , and denotes the space of the continuous functions in the closure of . Following the idea of the streamline diffusion method with a modified test function , the following CG-FEM method will be formed: look for , such that for any and any , Note that approximates the solution and . The parameter is related to the absorption and scattering coefficients,33 and it is chosen to be a piecewise constant, defined as , with being the diameter of a mesh element . (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 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 denote the ’th vertex of the mesh. Considering the Lagrange nodal basis of associated with vertices, satisfying , in which is the Dirac-Delta function. Each basis function is nonzero only within those elements sharing the vertex . The numerical solution and the test function can now be expanded in terms of basis, namely, In Eq. (10), represents the number of nodes (basis function) in discrete direction . Since the linear polynomial space is employed as discrete space, the gradient of the test function 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 . In the linear system, matrix contains the terms related with the stiffness matrix, the mass matrix and so on, . relates with the source term, which can be acquired either by employing the boundary condition in Eq. (2) or directly processing the source term25 . The linear equation was then solved using the generalized minimum residual method (GMRES),22 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.23 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 EquationBecause 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 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, approximation1 (in which case, one uses spherical harmonic functions whose index is 0 and 1). The quasisteady DE can be expressed as39 where is photon density () and is related to the photon radiance in Eq. (1) as follows:Here , where is the absorption coefficient (), while is the reduced scattering coefficient () and defined as . represents the light source. In this work, the Robin boundary condition is applied, that is, where is defined as , with expressed as and represents the refractive index of the object when the ambient medium is air ().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 MethodMC 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.52 Mesh-based MC methods have also been proposed12 to deal with complex shapes of media. In our simulation, mesh-based MC was employed to allow comparisons of identical meshes12,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.Settings for Numerical Simulations3.1.Ideal Pencil Beam SimulationIn our simulations, a 3-D rectangular computational domain is considered first, referred to as Rect. 1, which is similar to the setting54 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 with the direction orthogonal to the plane ; see Fig. 2(b). This type of source was modeled by the Dirac-Delta function . Fig. 2The 3-D rectangle used in simulations: (a) 3-D mesh of rectangle; (b) the slice at , where the arrow indicates the position in which the photon is launched perpendicularly to the plane . ![]() The optical properties of the 3-D rectangle were selected according to the most common values of human tissue55 (see Tables 3 and 4). For instance, the values for the healthy human brain are less than for the absorption coefficient and around for the scattering coefficient56 at 674 nm. The anisotropy of the tissue was set to and the refractive index 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 method53 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 value28 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 ( 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 . 3.2.Gaussian Shape Beam SimulationMany 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 in which represents the beam width. Two different 3-D rectangles were employed, Rect. 1 in Sec. 3.1 and Rect. 2, given as . Note that Rect. 2 is smaller than Rect. 1. Only the solutions to the RTE and MC were presented and compared, as the DE is not suitable for predicting photon propagation when the region of interest (ROI) is near the source. When the light source is changed from a point source to a Gaussian shape, solving the RTE changes only the boundary conditions presented in (2) or the source term. As a consequence, the matrix assembling time and computation time are almost unchanged. Similarly, extended illumination sources can be efficiently modeled using MC with overhead less than 5% of the total computational time.11,58 Here, photons were employed for MC simulation for each single point of Gaussian beam for both rectangles, Rect. 1 and Rect. 2.The anisotropy and refractive index were set to be the same values as in Sec. 3.1. The parameter was set to 0.5. In the first example (Rect. 1), the absorption coefficient was set to and the scattering coefficient was set to . 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 for absorption coefficient and for scattering coefficient, respectively. Fig. 33-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 where represents the number of points in comparison. Note that what RMSE measures is an absolute error. The relative error at each point between the solutions was defined asBased on this, the mean relative error (MRE, defined as ) 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.Results4.1.Three-Dimensional Rectangle Simulations with Ideal Pencil Beam4.1.1.Comparison between different numerical quadraturesIn order to evaluate the effect of the choices of numerical quadrature on the solution to the RTE, all the numerical quadratures reviewed in Sec. 2 were tested first for all sets of optical properties listed in Table 3. Tables 1 and 2 report only the errors for the optical settings with absorption coefficient and scattering coefficient . Table 1 lists the 3-D RMSE (, with defined in the corresponding tables) in different regions with different numerical quadratures. Note that the number in parentheses in the first column in Table 1 is the number of the discrete angles used for each numerical quadrature. Table 1RMSE with different numerical quadrature (μa=0.02 mm−1, μs=5 mm−1, g=0.9).
Table 2MRE with different numerical quadrature (μa=0.02 mm−1, μs=5 mm−1, g=0.9).
From Table 1, it is clear that when the ROI is far from the source and boundary, the error between the solutions to the RTE and from the MC simulations is smaller for all numerical quadratures. However, it can be seen from Table 1 that the error for Legendre equal-weight quadrature is larger in deeper regions compared with other quadratures. In order to find the quadrature with the least error in the simulations, the 3-D MRE was calculated and is shown in Table 2. It can be concluded that the MRE exhibits the same trend as the RMSE. From Table 2, one can also see that product Gaussian quadrature has a larger error () than Lebedev quadrature and level symmetric quadrature although it employs fewer discrete angles. Considering that the next product Gaussian quadrature will have 98 discrete angles and this will increase the computation cost, product Gaussian quadrature will not be our first choice. For Lebedev and level symmetric quadratures, the error for the whole rectangle is around 10%, while for deeper regions it is around 5%. Hence we can infer that level symmetric quadrature and Lebedev quadrature are the best among all tested quadratures, at least for the optical properties examined in this paper. It was observed that level symmetric quadrature can generate negative weights, which is nonphysical when the number of discrete angles is larger than 528.40 On the other hand, both numerical quadratures provide similar relative error estimations, with the difference less than 0.5%, at deeper depths. Thus, in this paper, Lebedev quadrature will be further investigated in the following simulations. 4.1.2.Ideal pencil beam simulationsFigure 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 and . With this combination of optical properties in tissue, the DE is not an accurate model to predict photon propagation.28 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 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. 4The 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 ; (b) the contours within the plane ; the value of the outermost curve is ; (c) the contours within the plane . ![]() Fig. 5Comparison between the solutions to RTE, DE, and MC method at different depths, along the line : (a) ; (b) ; (c) ; (d) ; (e) ; (f) . ![]() Table 3RMSE for different optical properties at different depth along the line x=0 mm (g=0.9).
Table 4MRE for different optical properties at different depth along line x=0 mm (g=0.9).
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 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, and , the errors are larger as ROI gets closer to the light source located at (0,0,0). At the extreme shallow region, , 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 are shown in Tables 3 and 4. For all examples tested and for depths smaller than , the MRE is less than or close to 5%. At a depth of , the MRE is less than 10%. Tables 3 and 4 show only the results along the line at different depths. For a comprehensive comparison of the solutions to the RTE and results from the MC, 3-D RMSE and MRE () were also calculated. The results are shown in Tables 5 and 6. Table 5RMSE for different optical properties to 3-D rectangle (g=0.9).
Table 6MRE for different optical properties to 3-D rectangle (g=0.9).
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 BeamThe 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. 6The 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) ; (b) , the value of the outermost contour is ; (c) . ![]() Fig. 7The 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) , the value of outermost curve is ; (b) , the value of outermost curve is ; (c) . ![]() Table 73-D RMSE and MRE for two rectangles with Gaussian shape beam (g=0.9).
4.3.Importance of the Normalized Phase Function TechniqueTo illustrate the importance of the normalized phase function technique to the stability of the proposed method, especially for large anisotropy with 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. 8Comparison 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. 5.Discussion and Future WorkIn this paper, the RTE was solved using our proposed algorithm with a modified DOM in angular variables and a streamline diffusion CG-FEM in space variables. The quadrature with the highest computational efficiency, Lebedev quadrature, was chosen to calculate the integral term, combined with the phase function normalization technique. 3-D rectangles with various common optical properties of human tissue were employed to simulate photon propagation. The solutions to the RTE and DE were compared to those found by MC methods. In summary, for a 3-D rectangle irradiated by an ideal pencil beam, the MRE for a region farther than 2 mm is around 5%. For a light source with Gaussian irradiation, even with a smaller computational domain, the worst MRE is still less than 10% for deeper regions. One advantage of our method is that, from a mathematical perspective, it reduces the dimension of the original RTE using the DOM. Furthermore, a streamline diffusion technique was employed to deal with the transport term in the RTE within a CG-FEM framework. The difference between the test function of the standard CG-FEM and streamline diffusion test function is the addition of a term that is the gradient of the original test function. Herein, it is quite easy to implement the modification. In our simulations of highly scattering media (a few tenths per mm of scattering coefficient), the solution was occasionally negative when the streamline diffusion technique was not used. In practice, with highly scattering tissue, there is no need to solve the RTE to obtain photon distribution, as DE provides a computationally more efficient model. The other advantage of using the streamline diffusion technique, at least with the optical properties examined in this work, is that fewer iterations were needed when the GMRES was used to solve the resulting linear equations. In our work, it is important to use the function normalization technique to calculate the integral term accurately. This technique mitigates the instability of DOM dramatically and further reduces the nonphysical oscillations of the solution. The oscillation of the solution to the RTE with DOM is more prominent with ideal pencil beam. Various techniques were proposed in literature for reducing the oscillation. For instance, in Ref. 23, the weight of numerical quadrature is calculated based on the integration of linear interpolation of the phase function. In comparison to the method discussed in Ref. 23, the phase function normalization technique is more direct and easier to apply. We want to point out the DG-FEM is known to be a good choice for solving problems with transport effects, like the RTE. It is also noted that the implementation of the DG-FEM can make the scheme parameters free and succinct. Although the DG-FEM involves more unknowns in general, it is still worth investigating in future work. Finally, we comment on the computational time of the proposed methods solving the RTE directly and the MC simulations. In the case of a single-point source such as pencil beam, solving the RTE directly is not advantageous in computational efficiency. When Gaussian beam or sources of arbitrary shape are considered, however, solving the RTE directly is more cost effective than conventional MC simulations. This is illustrated by Table 8, where computational times are reported for solving the RTE directly with a Gaussian beam and 50 discrete angles (including time to assemble matrices and solve linear equation), and MC simulations with or photons employed for each single point of Gaussian beam. If using MC simulations presented in Ref. 58, solving the RTE by the proposed algorithm will have similar computational efficiency when a large number of photons () is employed in MC simulations. We also want to point out that the use of parallel computation can reduce the time required to solve the forward photon propagation problem. In our code, CPU-based parallel computing techniques had already been implemented. In simulations, the time would be reduced by at least 40% if parallel computing was used (four-core CPU) in general. Still, using GPU computing can further decrease the time required for FEM implementation. Table 8Comparison of time consumption between solving RTE and MC simulations with Gaussian beam.
6.ConclusionA modified finite element method based on the streamline diffusion idea is applied to simulate the RTE with high accuracy. It is demonstrated numerically that the method provides an accurate and robust approach to solving the forward problem in the bio-optical regime. Future investigation includes evaluation of the performance of this algorithm on objects exhibiting a wider range of optical properties or objects with very general shapes. The output flux of photons calculated by this algorithm should also be verified by experiments. AcknowledgmentsThis work was funded by NSF CMMI award #1200270. F. Li was partially supported by NSF grant DMS-1318409. ReferencesT. Tarvainen, M. Vauhkonen and S. R. Arridge,
“Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation,”
J. Quantum Spectrosc. Radiat. Transfer, 109
(17–18), 2767
–2778
(2008). http://dx.doi.org/10.1016/j.jqsrt.2008.08.006 Google Scholar
A. P. Gibson, J. C. Hebden and S. R. Arridge,
“Recent advances in diffuse optical imaging,”
Phys. Med. Biol., 50
(4), R1
–R43
(2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 PHMBA7 0031-9155 Google Scholar
S. R. Arridge and J. C. Schotland,
“Optical tomography: forward and inverse problems,”
Inverse Prob., 25
(12), 123010
(2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar
V. Ntziachristos,
“Going deeper than microscopy: the optical imaging frontier in biology,”
Nat. Methods, 7
(8), 603
–614
(2010). http://dx.doi.org/10.1038/nmeth.1483 1548-7091 Google Scholar
M. Ozturk et al.,
“Mesoscopic fluorescence molecular tomography for evaluating engineered tissues,”
Ann. Biomed. Eng., 1
–13
(2015). http://dx.doi.org/10.1007/s10439-015-1511-4 ABPBBK 0084-6589 Google Scholar
M. J. Niedre et al.,
“Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice in vivo,”
Proc. Natl. Acad. Sci. U.S.A., 105
(49), 19126
–19131
(2008). http://dx.doi.org/10.1073/pnas.0804798105 Google Scholar
L. Zhao et al.,
“Lp regularization for early gate fluorescence molecular tomography,”
Opt. Lett., 39
(14), 4156
–4159
(2014). http://dx.doi.org/10.1364/OL.39.004156 OPLEDP 0146-9592 Google Scholar
C. Zhu and Q. Liu,
“Review of Monte Carlo modeling of light transport in tissues,”
J. Biomed. Opt., 18
(5), 050902
(2013). http://dx.doi.org/10.1117/1.JBO.18.5.050902 JBOPFO 1083-3668 Google Scholar
E. Alerstam, T. Svensson and S. Andersson-Engels,
“Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,”
J. Biomed. Opt., 13
(6), 060504
(2008). http://dx.doi.org/10.1117/1.3041496 JBOPFO 1083-3668 Google Scholar
Q. Fang and D. A. Boas,
“Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,”
Opt. Express, 17
(22), 20178
(2009). http://dx.doi.org/10.1364/OE.17.020178 OPEXFF 1094-4087 Google Scholar
Q. Fang and D. R. Kaeli,
“Accelerating mesh-based Monte Carlo method on modern CPU architectures,”
Biomed. Opt. Express, 3
(12), 3223
–3230
(2012). http://dx.doi.org/10.1364/BOE.3.003223 BOEICL 2156-7085 Google Scholar
J. Chen, Q. Fang and X. Intes,
“Mesh-based Monte Carlo method in time-domain widefield fluorescence molecular tomography,”
J. Biomed. Opt., 17
(10), 106009
(2012). http://dx.doi.org/10.1117/1.JBO.17.10.106009 JBOPFO 1083-3668 Google Scholar
M. S. Ozturk et al.,
“Mesoscopic fluorescence tomography of a photosensitizer (HPPH) 3D biodistribution in skin cancer,”
Acad. Radiol., 21
(2), 271
–280
(2014). http://dx.doi.org/10.1016/j.acra.2013.11.009 Google Scholar
L. Zhao et al.,
“The integration of 3-D cell printing and mesoscopic fluorescence molecular tomography of vascular constructs within thick hydrogel scaffolds,”
Biomaterials, 33
(21), 5325
–5332
(2012). http://dx.doi.org/10.1016/j.biomaterials.2012.04.004 BIMADU 0142-9612 Google Scholar
M. S. Ozturk et al.,
“Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue,”
J. Biomed. Opt., 18
(10), 100501
(2013). http://dx.doi.org/10.1117/1.JBO.18.10.100501 JBOPFO 1083-3668 Google Scholar
W. Bangerth,
“A framework for the adaptive finite element solution of large-scale inverse problems,”
SIAM J. Sci. Comput., 30
(6), 2965
–2989
(2008). http://dx.doi.org/10.1137/070690560 SJOCE3 1064-8275 Google Scholar
A. Edmans and X. Intes,
“Mesh optimization for Monte Carlo-based optical tomography,”
Photonics, 2
(2), 375
–391
(2015). http://dx.doi.org/10.3390/photonics2020375 Google Scholar
W. Cai, M. Xu and R. R. Alfano,
“Analytical form of the particle distribution based on the cumulant solution of the elastic Boltzmann transport equation,”
Phys. Rev. E, 71
(4), 041202
(2005). http://dx.doi.org/10.1103/PhysRevE.71.041202 Google Scholar
A. Kim and M. Moscoso,
“Chebyshev spectral methods for radiative transfer,”
SIAM J. Sci. Comput., 23
(6), 2074
–2094
(2002). http://dx.doi.org/10.1137/S1064827500382312 SJOCE3 1064-8275 Google Scholar
A. Liemert and A. Kienle,
“Exact and efficient solution of the radiative transport equation for the semi-infinite medium,”
Sci. Rep., 3
(2013). http://dx.doi.org/10.1038/srep02018 SRCEC3 2045-2322 Google Scholar
T. Tarvainen et al.,
“Finite element model for the coupled radiative transfer equation and diffusion approximation,”
Int. J. Numer. Methods Eng., 65
(3), 383
–405
(2006). http://dx.doi.org/10.1002/nme.1451 IJNMBH 0029-5981 Google Scholar
P. Surya Mohan et al.,
“Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,”
J. Comput. Phys., 230
(19), 7364
–7383
(2011). http://dx.doi.org/10.1016/j.jcp.2011.06.004 JCTPAH 0021-9991 Google Scholar
H. Gao and H. Zhao,
“A fast-forward solver of radiative transfer equation,”
Transport Theory Stat. Phys., 38
(3), 149
–192
(2009). http://dx.doi.org/10.1080/00411450903187722 TTSPB4 0041-1450 Google Scholar
A. D. Klose, V. Ntziachristos and A. H. Hielscher,
“The inverse source problem based on the radiative transfer equation in optical molecular imaging,”
J. Comput. Phys., 202
(1), 323
–345
(2005). http://dx.doi.org/10.1016/j.jcp.2004.07.008 JCTPAH 0021-9991 Google Scholar
H. Fujii et al.,
“Numerical modeling of photon propagation in biological tissue based on the radiative transfer equation,”
in AIP Conference Proceedings,
579
–585
(2013). http://dx.doi.org/10.1063/1.4794636 Google Scholar
K. Ren et al.,
“Algorithm for solving the equation of radiative transferin the frequency domain,”
Opt. Lett., 29
(6), 578
–580
(2004). http://dx.doi.org/10.1364/OL.29.000578 OPLEDP 0146-9592 Google Scholar
F. Asllanaj et al.,
“Radiative transfer equation for predicting light propagation in biological media: comparison of a modified finite volume method, the Monte Carlo technique, and an exact analytical solution,”
J. Biomed. Opt., 19
(1), 015002
(2014). http://dx.doi.org/10.1117/1.JBO.19.1.015002 JBOPFO 1083-3668 Google Scholar
Y. Wang et al.,
“Diffuse fluorescence tomography based on the radiative transfer equation for small animal imaging,”
Proc. SPIE, 8937 89370W
(2014). http://dx.doi.org/10.1117/12.2039017 Google Scholar
H.-S. Li, G. Flamant and J.-D. Lu,
“Mitigation of ray effects in the discrete ordinates method,”
Numer. Heat Transf. Part B Fundam., 43
(5), 445
–466
(2003). http://dx.doi.org/10.1080/713836241 Google Scholar
S. L. Jacques,
“Optical properties of biological tissues: a review,”
Phys. Med. Biol., 58
(11), R37
(2013). http://dx.doi.org/10.1088/0031-9155/58/11/R37 PHMBA7 0031-9155 Google Scholar
A. D. Klose,
“Transport-theory-based stochastic image reconstruction of bioluminescent sources,”
J. Opt. Soc. Am. A Opt. Image Sci. Vis., 24
(6), 1601
–1608
(2007). http://dx.doi.org/10.1364/JOSAA.24.001601 Google Scholar
J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods, Springer, New York, NY
(2008). Google Scholar
G. Kanschat,
“A robust finite element discretization for radiative transfer problems with scattering,”
East West J. Numer. Math., 6 265
–272
(1998). Google Scholar
J. A. Eichholz,
“Discontinuous Galerkin methods for the radiative transfer equation and its approximations,”
University of Iowa,
(2011). Google Scholar
B. Hunter and Z. Guo,
“A new and simple technique to normalize the hg phase function for conserving scattered energy and asymmetry factor,”
Numer. Heat Transf. Part B Fundam., 65
(3), 195
–217
(2014). http://dx.doi.org/10.1080/10407790.2013.849992 Google Scholar
B. Hunter and Z. Guo,
“Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media,”
Int. J. Heat Mass Transf., 55
(5–6), 1544
–1552
(2012). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2011.11.009 Google Scholar
A. N. Brooks and T. J. R. Hughes,
“Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations,”
Comput. Methods Appl. Mech. Eng., 32
(1–3), 199
–259
(1982). http://dx.doi.org/10.1016/0045-7825(82)90071-8 CMMECC 0045-7825 Google Scholar
A. Kienle, F. K. Forster and R. Hibst,
“Influence of the phase function on determination of the optical properties of biological tissue by spatially resolved reflectance,”
Opt. Lett., 26
(20), 1571
–1573
(2001). http://dx.doi.org/10.1364/OL.26.001571 OPLEDP 0146-9592 Google Scholar
M. Schweiger et al.,
“The finite element method for the propagation of light in scattering media: boundary and source conditions,”
Med. Phys., 22
(11), 1779
–1792
(1995). http://dx.doi.org/10.1118/1.597634 MPHYA6 0094-2405 Google Scholar
E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport, Wiley, New York
(1984). Google Scholar
G. Longoni and A. Haghighat,
“Development of new quadrature sets with the Ordinate Splitting technique,”
in Proc. of Int. Meeting on Mathematical Methods For Nuclear Application,
(2001). Google Scholar
V. Lebedev and D. Laikov,
“A quadrature formula for the sphere of the 131st algebraic order of accuracy,”
Dokl. Math., 59
(3), 477
–481
(1999). Google Scholar
J. C. Stone, Adaptive Discrete-Ordinates Algorithms and Strategies,
(2007). Google Scholar
K. Atkinson,
“Numerical integration on the sphere,”
ANZIAM J., 23
(3), 332
–347
(1982). http://dx.doi.org/10.1017/S0334270000000278 Google Scholar
V. I. Lebedev,
“Values of the nodes and weights of ninth to seventeenth order gauss-markov quadrature formulae invariant under the octahedron group with inversion,”
USSR Comput. Math. Math. Phys., 15
(1), 44
–51
(1975). http://dx.doi.org/10.1016/0041-5553(75)90133-0 CMMPA9 0965-5425 Google Scholar
V. I. Lebedev,
“Spherical quadrature formulas exact to orders 25–29,”
Siberian Math. J., 18
(1), 99
–107
(1977). http://dx.doi.org/10.1007/BF00966954 SMTJAW 0037-4466 Google Scholar
B. A. Gregersen and D. M. York,
“High-order discretization schemes for biochemical applications of boundary element solvation and variational electrostatic projection methods,”
J. Chem. Phys., 122
(19), 194110
(2005). http://dx.doi.org/10.1063/1.1899146 JCPSA6 0021-9606 Google Scholar
B. Hunter and Z. Guo,
“Comparison of quadrature schemes in DOM for anisotropic scattering radiative transfer analysis,”
Numer. Heat Transf. Part B Fundam., 63
(6), 485
–507
(2013). http://dx.doi.org/10.1080/10407790.2013.777644 Google Scholar
H. Dehghani et al.,
“Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction,”
Commun. Numer. Methods Eng., 25
(6), 711
–732
(2009). http://dx.doi.org/10.1002/cnm.1162 CANMER 0748-8025 Google Scholar
C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Courier Dover Publications, New York
(2012). Google Scholar
S. C. Brenner and R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York
(2008). Google Scholar
J. Chen and X. Intes,
“Comparison of Monte Carlo methods for fluorescence molecular tomography—computational efficiency,”
Med. Phys., 38
(10), 5788
–5798
(2011). http://dx.doi.org/10.1118/1.3641827 MPHYA6 0094-2405 Google Scholar
Q. Fang and D. A. Boas,
“Tetrahedral mesh generation from volumetric binary and gray-scale images,”
in Proc. of the Sixth IEEE Int. Conf. on Symp. on Biomedical Imaging: From Nano to Macro,
1142
–1145
(2009). http://dx.doi.org/10.1109/ISBI.2009.5193259 Google Scholar
F. Long et al.,
“Dental imaging using mesoscopic fluorescence molecular tomography: an ex vivo feasibility study,”
Photonics, 1
(4), 488
–502
(2014). http://dx.doi.org/10.3390/photonics1040488 Google Scholar
J. L. Sandell and T. C. Zhu,
“A review of in-vivo optical properties of human tissues and its impact on PDT,”
J. Biophotonics, 4
(11–12), 773
–787
(2011). http://dx.doi.org/10.1002/jbio.201100062 Google Scholar
F. P. Bevilacqua et al.,
“In-vivo local determination of tissue optical properties,”
Proc. SPIE, 3194 262
–268
(1998). http://dx.doi.org/10.1117/12.301063 PSISDG 0277-786X Google Scholar
K. Ren, G. Bal and A. Hielscher,
“Frequency domain optical tomography based on the equation of radiative transfer,”
SIAM J. Sci. Comput., 28
(4), 1463
–1489
(2006). http://dx.doi.org/10.1137/040619193 SJOCE3 1064-8275 Google Scholar
R. Yao, X. Intes and Q. Fang,
“Generalized mesh-based Monte Carlo for wide-field illumination and detection via mesh retessellation,”
Biomed. Opt. Express, 7
(1), 171
(2016). http://dx.doi.org/10.1364/BOE.7.000171 BOEICL 2156-7085 Google Scholar
BiographyFeixiao Long is a graduate student at Rensselaer Polytechnic Institute. He received his BE degree in biomedical engineering from Beijing University of Technology in 2006, and he received his ME degree in biomedical engineering from Tsinghua University in 2009. Fengyan Li is an associate professor in the Department of Mathematical Sciences at Rensselaer Polytechnic Institute. She received her PhD in applied mathematics from Brown University, and both her BSc and MSc degrees from Peking University, China. She was named a 2008 Alfred P. Sloan research fellow and a recipient of NSF CAREER award in 2009. She is currently serving on the editorial board of SIAM Journal of Numerical Computing. Xavier Intes is an associate professor in the Department of Biomedical Engineering at Rensselaer Polytechnic Institute. He obtained his PhD from Université de Bretagne Occidentale and postdoctoral training at the University of Pennsylvania. He acted as the chief scientist of Advanced Research Technologies Inc., from 2003 to 2006. His research interests are on the application of diffuse functional and molecular time-resolved optical techniques for biomedical imaging in preclinical and clinical settings. |