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 (SP_𝑁). In XFEM scheme of SP_𝑁 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 time-consuming 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 (SP_𝑁) to RTE is employed due to its capacity in improving the solution in transport-like 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 SP_𝑁 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 SP_𝑁 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 extend-finite 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 SP_𝑁 approximation for optical imaging. In this study, we establish the mathematical framework of XFEM with SP_𝑁 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 SP_𝑁 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 time-consuming 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 SP_𝑁 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. SP_𝑁 Approximations

The general form of the (𝑁+1)/2SP_𝑁 equations and its (𝑁+1)/2 boundary conditions for optical imaging in three dimensions are [13]: −∇⋅𝐶_{𝑖,∇𝜑𝑖}(𝑟)∇𝜑_𝑖(𝑟)+^(𝑁+1)/2_𝑗=1𝐶_{𝑖,𝜑𝑖}(𝑟)𝜑_𝑖(𝑟)=𝐶_𝑖,𝑄(𝑟)𝑄(𝑟)(),𝑟∈Ω,𝑖,𝑗∈1,𝑁+12_{(𝑁+1)/2𝑗=1}.𝐶^𝑏_{𝑖,∇𝜑𝑗}(𝑟)𝑛⋅∇𝜑_𝑗(𝑟)=^(𝑁+1)/2_𝑗=1𝐶^𝑏_{𝑖,𝜑𝑗}(𝑟)𝜑_𝑗(𝑟)+𝐶_𝑖,𝑆𝑆_𝑖(𝑟),()𝑟∈𝜕Ω,𝑖,𝑗∈1,𝑁+12, (1) where 𝜑_𝑖(𝑟) is the SP_𝑁 composite moments of the radiances in RTE, and 𝐶_{𝑖,∇𝜑𝑖}(𝑟), 𝐶_{𝑖,𝜑𝑖}(𝑟), 𝐶_𝑖,𝑄(𝑟), 𝐶^𝑏_{𝑖,∇𝜑𝑗}(𝑟), 𝐶^𝑏_{𝑖,𝜑𝑗}(𝑟), and 𝐶_𝑖,𝑆 denote the coefficients which are related to ∇𝜑_𝑖(𝑟), 𝜑_𝑖(𝑟), 𝑄, and 𝑆 for SP_𝑁 equations at each point 𝑟 in the region Ω or boundary 𝜕Ω. These can be calculated by the absorption coefficient 𝜇_𝑎[mm⁻¹], scattering coefficient 𝜇_𝑠[mm⁻¹], 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 𝐽⁺[nW/mm²] 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: 𝐽⁺=^(𝑁+1)/2_𝑗=1𝐶^𝐽_∇𝜑𝑗(𝑟)𝑛⋅∇𝜑_𝑗(𝑟)+^(𝑁+1)/2_𝑗=1𝐶^𝐽_𝜑𝑗(𝑟)𝜑_𝑗(𝑟), (2) where 𝐶^𝐽_∇𝜑𝑗(𝑟) and 𝐶^𝐽_𝜑𝑗(𝑟) are the coefficients which can be calculated for SP_𝑁 equations. The detailed derivation of (1)–(2) refers to the literature [13].

2.2. Extended Finite Element Discretization

The SP_𝑁 equations and its boundary conditions can be solved using Galerkin finite element scheme, the weak form of SP_𝑁 equations can be written as follows: _Ω𝐶_{𝑖,∇𝜑𝑖}(𝑟)⋅∇𝜑_𝑖(𝑟)(𝑟)+⋅∇𝑣𝑑𝑟_{(𝑁+1)/2𝑗=1Ω}𝐶_{𝑖,𝜑𝑗}(𝑟)𝜑_𝑗(𝑟)(𝑟)−⋅𝑣𝑑𝑟_{(𝑁+1)/2𝑗=1𝜕Ω}𝐶_{𝑖,∇𝑓𝜑𝑗}(𝑟)∇𝑓_{𝑛⋅𝜑𝑖}𝜑_𝑗(𝑟)(𝑟)=⋅𝑣𝑑𝑟_Ω𝐶_𝑖,𝑄(𝑟)𝑄(𝑟)⋅𝑣(𝑟)𝑑𝑟+_𝜕Ω𝐶_{𝑖,𝑓𝑆𝑖}(𝑟)𝑓_𝑆𝑖𝑆_𝑖(𝑟)⋅𝑣(𝑟)𝑑𝑟, (3) 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.

Figure 1 

(a) Heterogeneous tissues. (b) Standard finite element mesh boundary conforms to the interface. (c) Enriched elements and nodes of XFEM “cut” by the interface.

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 𝑁(𝑟)=min_{𝑟𝑏∈Γ}𝑟−𝑟^𝑏sign𝑛_𝑠⋅𝑟−𝑟^𝑏.sign(𝜉)=1if𝜉≥0,−1if𝜉<0, (4) where 𝑟^𝑏 is the point at the internal boundary Γ, sign 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: 𝜓_𝑗(𝑟)=𝑣_𝑗(𝑟)|𝑁(𝑟)|−𝑁_𝑗.𝑁(𝑟)=⁴_𝑘=1𝑣_𝑘(𝑟)𝑁_𝑘.𝑁_𝑗=𝑁𝑟_𝑗, (5) 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: 𝜑(𝑟)≈^𝑁c_𝑖=1𝜙_𝑖𝑣_𝑖(𝑟).𝜑(𝑟)≈𝜑(𝑟)+^𝑁e_𝑗=1𝑎_𝑗𝜓_𝑗(𝑟)≜𝑣(𝑟)𝜑, (6) where 𝜑(𝑟) is the conventional approximation and 𝑁_c is the number of nodes. 𝜑(𝑟) is the extended approximation, 𝑎_𝑗 is the enriched degrees of freedom, and 𝑁_e 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 𝑣(𝑟)=𝑣₁(𝑟),𝑣₂(𝑟)⋯𝑣_𝑁c+𝑁e(𝑟=𝑣₁)(𝑟),𝑣₂(𝑟)⋯𝑣_𝑁c(𝑟),𝜓₁(𝑟),𝜓₂(𝑟)⋯𝜓_𝑁e(𝑟).𝜑_𝑖=𝜙_𝑖,1,𝜙_𝑖,2⋯𝜙_𝑖,𝑁c,𝑎_𝑖,1,𝑎_𝑖,2⋯𝑎_𝑖,𝑁e^𝑇,()𝑖∈1,𝑁+12, (7) where 𝑣(𝑟) and 𝜑_𝑖 have 𝑁_c+𝑁_e components for each variable of SP_𝑁 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 (𝑁+1)/2SP_𝑁 equations (3) can be written as follows: _elem𝑚_1,𝜑1𝑚_1,𝜑2⋯𝑚_{1,𝜑(𝑁+1)/2}𝑚_2,𝜑1𝑚_2,𝜑2⋯𝑚_{2,𝜑(𝑁+1)/2}𝑚⋮⋮⋮⋮_{(𝑁+1)/2,𝜑1}𝑚_{(𝑁+1)/2,𝜑2}⋯𝑚_{(𝑁+1)/2,𝜑(𝑁+1)/2}×𝜑₁𝜑₂⋮𝜑_(𝑁+1)/2=_elem𝑏₁𝑏₂⋮𝑏_(𝑁+1)/2, (8) where elem denotes all the elements in the region, for the (𝑁+1)/2 variables 𝜑 and 𝑏, each of which has 𝑁_c+𝑁_e components in total according to (7). 𝑚_{𝑖,𝜑𝑗} in the left coefficient matrix is the small element matrix and can be written in details: 𝑚_{𝑖,𝜑𝑗}=∫_Ωe𝐶_{𝑖,∇𝜑𝑖}(𝑟)⋅∇𝑣_𝑝⋅∇𝑣_𝑞𝑑𝑟+∫_Ω𝑒𝐶_{𝑖,𝜑𝑗}(𝑟)𝑣_𝑝⋅𝑣_𝑞𝑑𝑟−∫_𝜕Ω𝐶_{𝑖,∇𝑓𝜑𝑗}(𝑟)𝑓_{𝑛⋅𝜑𝑖}𝑣_𝑝⋅𝑣_𝑞∫𝑑𝑟,𝑖=𝑗._Ωe𝐶_{𝑖,𝜑𝑗}(𝑟)𝑣_𝑝⋅𝑣_𝑞𝑑𝑟−∫_𝜕Ω𝐶_{𝑖,∇𝑓𝜑𝑗}(𝑟)𝑓_{𝑛⋅𝜑𝑖}𝑣_𝑝⋅𝑣_𝑞𝑑𝑟,𝑖≠𝑗, (9) where 𝑣_𝑝 and 𝑣_𝑞 are the corresponding matrix elements of 𝑣(𝑟) in (7) in the finite element Ω_e. 𝑝,𝑞 are the number marks of the four points in element Ω_e. Similarly, 𝑏_𝑖 in the right is 𝑏_𝑖=_Ω𝐶_𝑖,𝑄(𝑟)𝑄_𝑝⋅𝑣_𝑞𝑑𝑟+_𝜕Ω𝐶_{𝑖,𝑓𝑆𝑖}(𝑟)𝑓_𝑆𝑖𝑆_𝑖(𝑟)⋅𝑣_𝑞𝑑𝑟. (10)

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 Ω_e is _Ωe𝐹(𝑟)𝑑𝑟≈⁴_𝑖=1𝐺_𝑖𝐹𝜉_𝑖, (11) where 𝐹(𝑟) is the general integrand, 𝐺_𝑖 is the coefficient and all is 0.25, 𝜉_𝑖 is the Gaussian quadrature point as showed in Figure 2.

Figure 2 

Gaussian quadrature in an enriched element.

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: ∇𝜓_𝑗(𝑟)|_{𝑟=𝜉𝑖}=∇𝑣_𝑗(𝑟)⋅|𝑁(𝑟)|−𝑁_𝑗+𝑣_𝑗(𝑟)⋅∇|𝑁(𝑟)|−𝑁_𝑗|_{𝑟=𝜉𝑖}.𝑁(𝑟)|_{𝑟=𝜉𝑖}=⁴_𝑘=1𝑣_𝑘𝜉_𝑖𝑁_𝑘.𝑣_𝑗𝜉_𝑖𝑣=0.5854,𝑖=𝑗._𝑗𝜉_𝑖=0.1382,𝑖≠𝑗. (12)

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 SP_𝑁 equations is established and can be rewritten in the matrix form: 𝑀_{(𝑁c+𝑁e)∗(𝑁+1)/2,(𝑁c+𝑁e)∗(𝑁+1)/2}Φ_{(𝑁c+𝑁e)∗(𝑁+1)/2,1}=𝐵_{(𝑁c+𝑁e)∗(𝑁+1)/2,1}, (13) where 𝑀 is a matrix including (𝑁_c+𝑁_e)∗(𝑁+1)/2 rows and (𝑁_c+𝑁_e)∗(𝑁+1)/2 columns, Φ is the unknowns including (𝑁_c+𝑁_e)∗(𝑁+1)/2 components, and 𝐵 is the source term including (𝑁_c+𝑁_e)∗(𝑁+1)/2 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 CORR(𝐽⁺_𝑋/𝐹,𝐽⁺_Fine) and the well known mean relative numerical error MRNE_Fine 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: CORR𝐽⁺_𝑋/𝐹,𝐽⁺_Fine=∑^𝑁s_𝑖=1𝐽⁺_{𝑖,𝑋/𝐹}−𝐽⁺_{𝑖,𝑋/𝐹}𝐽⁺_𝑖,Fine−𝐽⁺_𝑖,Fine∑^𝑁s_𝑖=1𝐽⁺_{𝑖,𝑋/𝐹}−𝐽⁺_{𝑖,𝑋/𝐹}²∑_𝑖𝐽⁺_𝑖,Fine−𝐽⁺_𝑖,Fine²,MRNE_Fine=∑^𝑁s_𝑖=1𝐽⁺_{𝑖,𝑋/𝐹}−𝐽⁺_𝑖,Fine/𝐽⁺_𝑖,Fine𝑁_s, (14)(15) where 𝐽⁺_{𝑖,𝑋/𝐹} is the 𝐽⁺ compute by XFEM or FEM, 𝐽⁺_𝑖,Fine is compute by FEM on the fine mesh. 𝑁_s is the number of the sampling point, 𝐽⁺_{𝑖,𝑋/𝐹} and 𝐽⁺_𝑖,Fine is mean value of 𝐽⁺_{𝑖,𝑋/𝐹}𝐽⁺_𝑖,Fine in (14). CORR(𝐽⁺_𝑋/𝐹,𝐽⁺_Fine) = 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 𝐽⁺_𝑖,MC for the 𝐽⁺_𝑖,Fine in (14) and (15), the mean relative numerical error MRNE_MC and correlation coefficient CORR(𝐽⁺_𝑋/𝐹/MC,𝐽⁺_MC) 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³ 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 SP_𝑁 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: 𝜑−𝜑^ℎ→0asℎ→0, (16) 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).

Table 1 

Optical parameters of phantom.

(a)

(b)

(c)

Figure 3 

(a) Heterogeneous phantom, (b) the element mesh of FEM, and (c) element mesh of XFEM without internal boundary generation, red region is the enriched region.

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, SP₃, SP₇ 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 CORR(𝐽⁺_𝑋/𝐹,𝐽⁺_Fine) of XFEM is closer to 1 for DA, SP₃, and SP₇ 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 MRNE_Fine of XFEM is 14%, 15%, and 16%, which is much smaller than FEM, whose MRNE_Fine is 31%, 22%, and 22% for DA, SP₃, and SP₇ equations respectively. The results listed in Table 2 using the two evaluation indexes demonstrate the validity of the proposed method for SP_𝑁 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 SP_𝑁 equations.

Table 2 

Results of XFEM and FEM for SP_𝑁 approximations.

(a)

(b)

(c)

(d)

Figure 4 

Surface detector readings 𝐽+ for SP𝑁 equations around the regular heterogeneous phantom using the two methods on the coarse mesh compared with the FEM on the fine mesh. (a) DA equations. (b) SP3 equations. (c) SP7 equations. (d) Comparison among the results of XFEM for DA, SP3, and SP7 equations.

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³ 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 SP_𝑁 equation and MC to depict the light propagation, the XFEM is compared with the FEM and MC simultaneously.

Table 3 

Optical parameters of the mouse organs.

(a)

(b)

(c)

Figure 5 

(a) The digital mouse model. (b) The external and internal boundary mesh for FEM. (c) The external boundary mesh and enriched region for XFEM.

As shown in Figure 4(d), it is clear that the blue and green curves agree well. Considering the difference of the results between SP₃ and SP₇ equations are very small, for simplicity, only the SP₃ 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 MRNE_MC and correlation coefficient CORR(𝐽⁺_𝑋/𝐹/MC,𝐽⁺_MC) of XFEM and FEM with respect to MC method are obtained and shown in Table 4. The MRNE_MC and CORR(𝐽⁺_𝑋/𝐹/MC,𝐽⁺_MC) 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 time-efficiency and this makes it more practical in imaging process.

Table 4 

Comparison among the FEM, XFEM and MC method.

(a)

(b)

(c)

Figure 6 

The surface detector readings with SP3 approximation (a) using FEM on the fine mesh, (b) using XFEM on the coarse mesh, and (c) using MC method.

Figure 7 

Surface detector readings 𝐽+ using MC method, XFEM, and FEM at the sampling points. Green line denotes the results of MC method, red and blue lines are that of FEM and XFEM with SP3 approximation.

4. Conclusion

We have derived the extended finite element method with SP_𝑁 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 SP_𝑁 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 SP_𝑁 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 SP_𝑁 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.