Research Article  Open Access
Extended Finite Element Method with Simplified Spherical Harmonics Approximation for the Forward Model of Optical Molecular Imaging
Abstract
An extended finite element method (XFEM) for the forward model of 3D optical molecular imaging is developed with simplified spherical harmonics approximation (). In XFEM scheme of equations, the signed distance function is employed to accurately represent the internal tissue boundary, and then it is used to construct the enriched basis function of the finite element scheme. Therefore, the finite element calculation can be carried out without the timeconsuming internal boundary mesh generation. Moreover, the required overly fine mesh conforming to the complex tissue boundary which leads to excess time cost can be avoided. XFEM conveniences its application to tissues with complex internal structure and improves the computational efficiency. Phantom and digital mouse experiments were carried out to validate the efficiency of the proposed method. Compared with standard finite element method and classical Monte Carlo (MC) method, the validation results show the merits and potential of the XFEM for optical imaging.
1. Introduction
Light propagation model in biological tissue is the foundation of optical imaging. An accurate forward model is important for location and quantification of target distribution in the fields of optical imaging modalities, such as diffusion optical tomography (DOT), fluorescence molecular tomography (FMT), bioluminescence tomography (BLT), and Cerenkov luminescence tomography (CLT) [1–5]. The propagation of the emission photons in tissue can be accurately represented by the radiative transfer equation (RTE) or Monte Carlo (MC) models, but they are extremely computationally expensive. Therefore, the commonly used mathematical model in optical imaging field is the diffusion approximation (DA) to RTE. However, the DA model can be used only in the highly scattering property of the biological tissue, and is not suitable for the real mouse with complex internal tissues. To reach a compromise between the accuracy and efficiency, simplified spherical harmonics approximation () to RTE is employed due to its capacity in improving the solution in transportlike domains with high absorption and small geometries [6, 7].
Owing to the complex and curvilinear geometries associated with the biological tissues, the classical finite element methods (FEM) with approximation become necessary for optical imaging, especially for heterogeneous tissues [8–10]. In the FEM scheme, the region of heterogeneous tissue is divided into small tetrahedron elements. The linear functions of the tetrahedron element are employed in the standard finite element basis function which requires the homogeneity of tissue in one element. The fine triangle mesh between the two tissues is required to conform to the complex internal boundary for ensuring the calculation accuracy of FEM. However, the generation of the internal boundary mesh is a hard work and time consuming task. Moreover, the fine mesh yields the excess cost in finite element computation, especially with approximation in three dimensions. Fortunately, the extended finite element method inherits all the advantages of standard FEM and exempts the internal boundary mesh generation. Consequently, the extendfinite element method (XFEM) may deal with the problems of fine triangle mesh generation perfectly.
XFEM was first introduced in the literature [11], and it has many applications in the area of mechanics [11, 12]. To the best of our knowledge, it has not been used with approximation for optical imaging. In this study, we establish the mathematical framework of XFEM with approximation as the optical forward model. Specifically, a mesh without the internal boundary conformation is employed in the XFEM scheme. The distance relationship between the tetrahedron vertex and the real tissue boundary is reflected by a signed distance function which is used to construct the enriched basis function. Then the enriched function is added to the standard finite element basis function, therefore a standard approximation of FEM is thus “enriched” in the discretized region (usually cut by the internal tissue boundary) of interest. In the weak form of equations, Gaussian quadrature is employed to calculate the integrals and the linear system equations are established. As a result, the calculation of XFEM can be carried out without the timeconsuming internal boundary mesh generation. Moreover, the required overly fine mesh conforming to the complex tissue boundary which leads to excess time cost can be avoided. XFEM conveniences its application to tissues with complex internal structure and improves the computational efficiency. Numerical experiments with a phantom and a digital heterogeneous mouse were carried out to evaluate the performance of the proposed method. The results were compared with that of standard FEM and MC method to demonstrate the efficiency of XFEM.
The paper is organized as follows. In Section 2, the detailed procedure of using XFEM for solving equations is derived. In Section 3, we evaluate the performance of the proposed method by comparing with the standard FEM and MC method in the experiments, and demonstrate the efficiency of the method. Conclusions and discussions are given in the last section.
2. Method
2.1. Approximations
The general form of the equations and its boundary conditions for optical imaging in three dimensions are [13]: where is the composite moments of the radiances in RTE, and , , , , , and denote the coefficients which are related to , , , and for equations at each point in the region or boundary . These can be calculated by the absorption coefficient , scattering coefficient , anisotropy parameter , and refractive index . is the internal source, and it represents the bioluminescence source in BLT or the fluorophore in and FMT. is the external source which represent the external laser source in DOT and FMT. represents the unit normal vector outward to the boundary.
The exiting partial current can be obtained from detector readings at the tissue boundary which can be acquired by CCD camera in practical application, and we have the general formulation: where and are the coefficients which can be calculated for equations. The detailed derivation of (1)–(2) refers to the literature [13].
2.2. Extended Finite Element Discretization
The equations and its boundary conditions can be solved using Galerkin finite element scheme, the weak form of equations can be written as follows: where and are the coefficient matrices with respect to and . is the test function, and it is the same as the standard linear basis function in elements. For finite element analysis in three dimensions, tetrahedral elements have become popular in numerical computation, because of their ability to describe complex geometries such as heterogeneous tissues. Thus the volumetric domain is discretized into small tetrahedral elements, and the elements need to conform to the internal tissue boundary as shown in Figures 1(a) and 1(b). When the internal boundary is not smooth, the size of the element must be small enough to ensure the accuracy and this may cause difficulty in mesh generation and lead to huge computation burden. In XFEM framework, mesh is generated as a region consists many different tissues, and the elements “cut” by the actual tissue boundary are enriched as shown in Figures 1(a) and 1(c). Using the enriched basis functions constituted by the signed distance functions, the boundary can be determined.
2.3. Enriched Strategy
In XFEM scheme, the signed distance function is adopted to depict the internal boundary of different tissues. The definition of the signed distance function is where is the point at the internal boundary , is the signed function, and denotes the unit normal vector outward to the internal boundary as shown in Figure 1(a). Herein, the position relation between the discrete point and the continuous boundary is completely reflected by the function .
The enriched basis function composed of the signed distance function of the element at tetrahedron node can be derived, and the continuous is discretized with its value at the node. The enrichment functions introduce a discontinuity in the gradient of the radiances field or the distribution of optical parameters to , thus the following integral of enriched function is more accurate than that of linear basis function: where is the linear basis function of the tetrahedron element. Using the enriched basis function , the enriched approximation can be written as the following form: where is the conventional approximation and is the number of nodes. is the extended approximation, is the enriched degrees of freedom, and is the number of the enriched nodes. It is clear that includes as a special case, and each enriched element has eight degrees of freedom. Rewrite the extended linear basis functions and the discrete point value of in matrix form, we obtain where and have components for each variable of equations. This form can directly use to assemble the finite element matrix.
2.4. Integrals and System Matrix
For finite element discretization, the matrix form of the equations (3) can be written as follows: where denotes all the elements in the region, for the variables and , each of which has components in total according to (7). in the left coefficient matrix is the small element matrix and can be written in details: where and are the corresponding matrix elements of in (7) in the finite element . are the number marks of the four points in element . Similarly, in the right is
In finite element framework, for ease of calculation, the integrand is assumed continuous, and the Gaussian quadrature can be used. In this paper, on the one hand, exact integrals can be obtained for the standard linear elements. On the other hand, the second order Gaussian quadrature with four quadrature points is adopted for the remaining enriched elements as showed in Figure 2. The Gaussian quadrature formula in the standard element is where is the general integrand, is the coefficient and all is 0.25, is the Gaussian quadrature point as showed in Figure 2.
For the integrand , the integrals of enriched function can be calculated by the linear basis function , and at the Gaussian points is a constant as follows:
Incorporating (7) and (8), and assembling all the element matrixes (8), after using the Gaussian quadrature to calculate the enriched elements of (9), the linear system equations of equations is established and can be rewritten in the matrix form: where is a matrix including rows and columns, is the unknowns including components, and is the source term including components. Obtain the of (13) and instead of in (2). The linear relationship between the surface detector readings and the source or is established.
The goal of forward problem is to get the relationship between the surface detector readings and the internal sources. To have the general and particular comparison, the correlation coefficient and the well known mean relative numerical error of are both defined to quantitatively evaluate the performance of XFEM or FEM on the coarse mesh with respect to the FEM on the fine mesh: where is the compute by XFEM or FEM, is compute by FEM on the fine mesh. is the number of the sampling point, and is mean value of in (14). = 1 illustrates the two data is identical after normalization and can be used to assess the degree of closeness between the two data. Similarly, by substituting the for the in (14) and (15), the mean relative numerical error and correlation coefficient of XFEM and FEM with respect to MC method can be obtained for evaluation.
3. Results and Discussion
3.1. Regular Phantom Experiment Compared with Standard FEM
The validation studies were performed using a cylindrical phantom of 30 mm height and 10 mm radius to model a mouse. It consisted of ellipsoids or cylinders to represent the tissues of mouse as shown in Figure 3(a). A solid sphere source of 1 mm radius and 0.238 nW/mm^{3} power density was centered at (3, 5, 0) inside the right lung. The relevant optical parameters from literature [14] are listed in Table 1. Numerical simulations are carried out to compare the standard FEM and XFEM for equations. The coarse mesh containing 3459 nodes without internal mesh boundary generation was used for XFEM and another coarse mesh contained 3573 nodes for FEM as shown in Figures 3(b) and 3(c). Because the precise analytic solutions for heterogamous phantom is difficult to obtained, the FEM on the fine mesh can get relative accurate solution according to the classical finite element analysis: where the numerical solution of FEM converges to the exact solution as the mesh size decreases. We choose the result of the FEM on the fine mesh containing 12312 nodes as the standard for comparison. The program of FEM and XFEM is coded in MATLAB on the desktop computer (Intel(R) Xeon(R) 2 CPU E5430 @ 2.66 GHz, and 8 G RAM).

(a)
(b)
(c)
345 interpolate points whose value is nonzero are uniformly sampled around the phantom surface. Choosing the FEM on the fine mesh as the standard, the absolute value of exiting partial current on the sampling points are arranged in ascending order. Then the results of XFEM, and FEM on the comparable coarse mesh at these sampling points can be obtained by interpolation. The comparison results with Diffusion, , approximations are showed in Figures 4(a), 4(b), and 4(c), respectively. Although the curve has fluctuation caused by the discretization error, the results have similar tendency. It is clear that the blue curve of XFEM on coarse mesh is closer to the green curve of FEM on fine mesh than that of FEM on the coarse mesh. Compared with the FEM, the of XFEM is closer to 1 for DA, , and with 0.94, 0.95, and 0.90. The two methods is carried out on comparative mesh even the XFEM on the slightly less coarse mesh (3459 nodes). Thus the solution using XFEM is nearer to the true solution. The of XFEM is 14%, 15%, and 16%, which is much smaller than FEM, whose is 31%, 22%, and 22% for DA, , and equations respectively. The results listed in Table 2 using the two evaluation indexes demonstrate the validity of the proposed method for equations. Moreover, it indicates that the XFEM is superior to the standard FEM and the solution of XFEM is more closer to the true solution for equations.

(a)
(b)
(c)
(d)
3.2. Digital Mouse Experiment Compared with Monte Carlo Method
In this experiment, a digital heterogeneous mouse from CT and cryosection data consisting of several organs is adopted to evaluate the performance of the XFEM. The mouse is shown in Figure 5(a) [15]. A cylinder source of 0.8 mm radius, 1.6 mm height, and 0.311 nW/mm^{3} power density was centered at (12, 7.5, 50.5) inside the liver. The optical parameters of the organs at 670 nm wavelength computed by the literature [14] are employed and shown in Table 3. The experiments were conducted using XFEM on the coarse mesh, FEM on the fine mesh, and classical MC method. The surface detector readings using MC method was achieved from MOSE [16] using 5 million photons. Since there are still errors between the equation and MC to depict the light propagation, the XFEM is compared with the FEM and MC simultaneously.

(a)
(b)
(c)
As shown in Figure 4(d), it is clear that the blue and green curves agree well. Considering the difference of the results between and equations are very small, for simplicity, only the results are presented in Figure 6. The absolute values on the surface using the FEM, XFEM and MC methods agree well and have the similar tendency in general. For detailed comparison, 521 points with nonzero values on the mouse surface are sampled, and then the surface value of these points using the three methods are obtained by interpolation and shown in Figure 7. Choosing the MC method as standard, the mean relative numerical error and correlation coefficient of XFEM and FEM with respect to MC method are obtained and shown in Table 4. The and of FEM is 0.86 and 44%, while those of XFEM is 0.93 and 45%. The computational time of several modules that perform specific computational tasks, mesh generation, matrix assembly and solver. Since the XFEM does not have the complex mesh generation as FEM, the total time cost of matrix assembly and solver is considered for comparison. It is clear that the time cost of XFEM is only 367 seconds on the coarser mesh while that of FEM is 2675 seconds with the similar accuracy results. All the time cost of two methods are far less than that of MC method. The XFEM has a distinct strength on timeefficiency and this makes it more practical in imaging process.

(a)
(b)
(c)
4. Conclusion
We have derived the extended finite element method with approximations for the forward model of the three dimensional optical imaging. Considering the complex geometric object, it is necessary to have the fined mesh to conform to the internal boundary. And the mesh conformation is a difficult issue in the pretreatment in the FEM, moreover the standard FEM on the fine mesh for approximations cost too much especially for the high order approximation. Fortunately, the XFEM can deal with this problem. The XFEM includes the standard FEM as a special case, which does not require a geometric representation of the interface or any boundary mesh generation. Use the signed distance functions add to the standard basis functions, the number of nodes can keep unchanged. The interface can be well depicted by the enriched functions, also the solution of equation is more accurate.
The XFEM is validated through the numerical experiments. Phantom experiment was conducted and its results agreed well with the well known classical FEM. Moreover, compared with the FEM, XFEM can get more accurate result even on the slightly coarse mesh, which is closer to the FEM on the fine mesh. Digital mouse experiments further indicate that XFEM is superior to the standard FEM. The MC method was employed to evaluate the performance the proposed method. Although the relative errors and correlation coefficient using the XFEM with respect to MC method is comparable to that of standard FEM, the time cost of XFEM is greatly decreased due to the adoption of coarser element mesh. All these indicate that the XFEM is more suitable for equations especially in complex heterogeneous tissues.
Adaptive FEM method is often adopted for its high efficiency [17, 18]. But nearly all the adaptive method is based on the mesh refinement which is extremely complicate especially in three dimensions. XFEM can also be seen as an adaptive FEM method because it increases the degrees of freedom near the internal boundary while using fixed mesh, therefore the mesh refinement can be avoided. Thus the process in each adaptive level is simplified and the calculation cost is decreased. Moreover, the discretization error of the forward model is improved, which is important for the quantification reconstruction. XFEM may provide a potential tool for reconstruction algorithm.
Acknowledgments
This work was supported by the Program of the National Basic Research and Development Program of China (973) under Grant no. 2011CB707702, the National Natural Science Foundation of China under Grant nos. 81090272, 81101083, 81101084, 81101100, 81000632, and 30900334, the National Key Technology Support Program under Grant no. 2012BAI23B06, the Natural Science Basic Research Plan in Shaanxi Province of China under Grant no. 2012JQ4015, and the Fundamental Research Funds for the Central Universities.
References
 V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of wholebody photonic imaging,” Nature Biotechnology, vol. 23, no. 3, pp. 313–320, 2005. View at: Publisher Site  Google Scholar
 M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Physics in Medicine and Biology, vol. 50, no. 12, pp. 2837–2858, 2005. View at: Publisher Site  Google Scholar
 A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Physics in Medicine and Biology, vol. 50, no. 4, pp. R1–R43, 2005. View at: Publisher Site  Google Scholar
 W. Cong, G. Wang, D. Kumar et al., “Practical reconstruction method for bioluminescence tomography,” Optics Express, vol. 13, no. 18, pp. 6756–6771, 2005. View at: Publisher Site  Google Scholar
 A. E. Spinelli, C. Kuo, B. W. Rice et al., “Multispectral Cerenkov luminescence tomography for small animal optical imaging,” Optics Express, vol. 19, no. 13, pp. 12605–12618, 2011. View at: Publisher Site  Google Scholar
 A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” Journal of Computational Physics, vol. 220, no. 1, pp. 441–470, 2006. View at: Publisher Site  Google Scholar
 M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using threedimensional frequencydomain simplified spherical harmonics equations,” Physics in Medicine and Biology, vol. 54, no. 8, pp. 2493–2509, 2009. View at: Publisher Site  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,” Journal of Computational Physics, vol. 202, no. 1, pp. 323–345, 2005. View at: Publisher Site  Google Scholar
 M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Medical Physics, vol. 22, no. 11 I, pp. 1779–1792, 1995. View at: Publisher Site  Google Scholar
 H. Jiang, “Frequencydomain fluorescent diffusion tomography: a finiteelementbased algorithm and simulations,” Applied Optics, vol. 37, no. 22, pp. 5337–5343, 1998. View at: Google Scholar
 N. Moës, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999. View at: Google Scholar
 J. Chessa, P. Smolinski, and T. Belytschko, “The extended finite element method (XFEM) for solidification problems,” International Journal for Numerical Methods in Engineering, vol. 53, no. 8, pp. 1959–1977, 2002. View at: Publisher Site  Google Scholar
 Y. Lu, B. Zhu, H. Shen, J. C. Rasmussen, G. Wang, and E. M. SevickMuraca, “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Physics in Medicine and Biology, vol. 55, no. 16, pp. 4625–4645, 2010. View at: Publisher Site  Google Scholar
 G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined opticalPET (OPET) system: a computer simulation feasibility study,” Physics in Medicine and Biology, vol. 50, no. 17, pp. 4225–4241, 2005. View at: Publisher Site  Google Scholar
 B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Physics in Medicine and Biology, vol. 52, no. 3, pp. 577–587, 2007. View at: Publisher Site  Google Scholar
 MOSE, http://www.mosetm.net/.
 Y. Lv, J. Tian, W. Cong et al., “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Optics Express, vol. 14, no. 18, pp. 8211–8223, 2006. View at: Publisher Site  Google Scholar
 A. Joshi, W. Bangerth, and E. M. SevickMuraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Optics Express, vol. 12, no. 22, pp. 5402–5417, 2004. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2012 Wei Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.