## Abstract

As an important small animal imaging technique, optical imaging has attracted increasing attention in recent years. However, the photon propagation process is extremely complicated for highly scattering property of the biological tissue. Furthermore, the light transport simulation in tissue has a significant influence on inverse source reconstruction. In this contribution, we present two Galerkin-based meshless methods (GBMM) to determine the light exitance on the surface of the diffusive tissue. The two methods are both based on moving least squares (MLS) approximation which requires only a series of nodes in the region of interest, so complicated meshing task can be avoided compared with the finite element method (FEM). Moreover, MLS shape functions are further modified to satisfy the delta function property in one method, which can simplify the processing of boundary conditions in comparison with the other. Finally, the performance of the proposed methods is demonstrated with numerical and physical phantom experiments.

© 2008 Optical Society of America

## 1. Introduction

Molecular imaging is a very promising and rapidly developing biomedical research field in which the modern tools and methods are being married to represent *in vivo* cellular and molecular processes directly, sensitively, non-invasively and specifically, such as monitoring protein-protein interactions, gene expression, cell trafficking and engraftment [1, 2]. Among molecular imagingmodalities, optical imaging, especially fluorescence and bioluminescence imaging, has become a research focus over the past years for its excellent performance, non-radiativity and high cost-effectiveness compared with traditional imaging techniques like X-ray computed tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET) and single photon emission computed tomography (SPECT) [3, 4, 5]. In fluorescence technology, the fluorophore probe inside the tissue absorbs the incident excitation photons produced by the external light source, and then decays to its ground state, emitting the photons synchronously [6, 7]. Whereas bioluminescence imaging employs luciferase enzymes, which can catalyze the biochemical reactions of substrate luciferin to generate bioluminescent photons in the presence of oxygen, ATP and Mg^{2+} [8, 9]. Although the photon generation schemes are various in different optical imaging modality, the emission spectrum in the optical engineering field is generally in the so-called near-infrared light window of the biological tissue ([650–900]*nm*), which can travel several centimeters in tissue due to the low photon absorption in the above spectral window [7, 10, 11].

The propagation of the emission photons in tissue can be accurately represented by the radiative transfer equation (RTE) approximated from Maxwell’s equations, but it is extremely computationally expensive for its integro-differential nature [12, 13]. Therefore, the commonly used mathematical model in optical imaging field is the diffusion equation derived from RTE in view of highly scattering property of the biological tissue [12, 14, 15]. Furthermore, many algorithms have been presented to simulate light transport in the turbid tissue and predict the diffuse light flux on the surface of the small animal based on diffusion model [16, 17]. The corresponding simulation results can be employed not only to verify the correctness of the physical model, but also to generate a sensitivity matrix which relates the surface measurements to the internal optical properties and will be employed in the inverse light source reconstruction [12].

It is well known that analytical, statistical and numerical techniques are three kinds of methods to solve the aforementioned diffusion approximation model [12, 18, 19, 20]. Among these methods, numerical techniques are studied widely because of its high efficiency and good applicability, and finite element method (FEM) is one of the most typical and successful algorithms. For example, Cong *et al.* [17, 21] employed FEM to obtain the photon flux density on the boundary of the homogeneous and heterogeneous phantoms. Lv *et al.* [22] used the adaptive FEM to compute the photon energy distribution on the phantom surface. In addition, an improved FEM was presented to handle light propagation in nonscattering regions within diffusing domains by Arridge *et al.* [23]. However, finite element mesh generation and data pre-processing are difficult and time-consuming, especially for three-dimensional irregular objects with complex internal structure like the heterogeneous biological tissue. For instance, the human head model was discretized into 155915 elements in the literature [24]. What is more, Shephard *et al.* [25] reported a mesh with two million tetrahedral elements! Although remarkable progress has been made in generating the three-dimensional structured meshes for FEM analysis of solids and structures, it is generally recognized that the development of fast and robust meshing techniques of three-dimensional objects with complex geometrical shape and internal structure is still a challenge in practice. Furthermore, for the linear FEManalysis, mesh generation and data preparation before calculation often need much more time than the assembly and solution of the FEM equations. Therefore, meshless methods are explored as a novel numerical analysis approach which can avoid or greatly simplify meshing task, and they have been successfully applied to solve problems of solid mechanics, heat transfer, fluid flow, *etc.* [26, 27, 28].

In this contribution, two Galerkin-based meshless methods (GBMM) are proposed and developed for simulating the photon propagation process in the biological tissue. Compared with FEM, GBMM uses only a set of discretized points and does not require any node connectivity or element information, which helps not only to avoid the burdensome meshing but also to describe complex inhomogeneous domains more accurately. In GBMM algorithms, moving least squares (MLS) approximation plays an important role, and two kinds of GBMM methods are provided in this paper according to whether the MLS shape functions are modified. Then, a linear matrix form linking the source distribution and the photon flux density can be established based on Galerkin approach and Gauss theory with the diffusion equation and Robin boundary condition. Moreover, our two proposed algorithms should incorporate *a priori* knowledge of the tissue optical parameters, which can be assigned on the basis of available data in literature [29] or *in vivo* diffuse optical tomography (DOT) measurements. The paper is organized as follows. The next section presents the details of the proposed GBMM, diffusion model based methods for photon propagation in the biological tissue. In the third section, the performance of the two methods is validated using numerical and physical phantoms and compared with the simulation data based on FEM or MC and the measured photon flux density by a cooled CCD camera. In the last section, relevant issues are discussed and conclusions are provided.

## 2. Methodology

#### 2.1. Diffusion approximation and boundary condition

Under the assumption that light scattering dominates over absorption, the propagation of the emitting photons in the biological tissue can be represented by steady-state diffusion equation when a continuous wave external excitation light source is used in fluorescence imaging experiment or a bioluminescent light source is employed in bioluminescence technology [7, 14]:

where Ω is the region of interest; Φ(x) represents the photon flux density at location x [*Watts*/*mm*
^{2}]; *S*(**x**) denotes the internal source density [*Watts*/*mm*
^{3}]; *µ*
* _{a}*(

**x**) is the absorption coefficient [

*mm*

^{-1}];

*D*(

**x**)=(3(

*µ*

*(*

_{a}**x**)+(1-

*g*)

*µ*

*(*

_{s}**x**)))

^{-1}is the optical diffusion coefficient,

*µ*

*(*

_{s}**x**) the scattering coefficient [

*mm*

^{-1}] and

*g*the anisotropy parameter.

In order to eliminate the diffusion approximation error near the surface where light does not propagate diffusively, an appropriate boundary condition should be specified [12, 30]. When the optical imaging experiment is carried out in a totally dark environment, Robin boundary condition can be employed [21, 30]:

where *∂*Ω is the corresponding boundary of Ω; *v*(**x**) refers to the unit normal vector outward to the boundary *∂*Ω; *A*(*x*;*n*,*n*
^{′}) is a function to incorporate the mismatch between the refractive indices *n* within Ω and *n*
^{′} in the surrounding medium. Furthermore, *A*(**x**;*n*,*n*
^{′}) can be approximately expressed as the following formula when the imaging experiment is performed in air, for which *n*
^{′}≈1.0:

where *R*(**x**) is a parameter relative to the internal reflection at *∂*Ω. According to the reference [31], *R*(**x**) can be approximated by *R*(**x**)≈-1.4399*n*
^{-2}+0.7099*n*
^{-1}+0.6681+0.0636*n*. In our study, the measured outgoing flux density *Q*(**x**) on *∂*Ω is:

## 2.2. Galerkin-based meshless methods

#### 2.2.1. Moving least squares approximation

The MLS approximation is the basis of the proposed GBMM algorithms. According to the literatures [32] and [33], the photon flux density Φ(**x**) at node **x** can be approximated by Φ* ^{h}*(

**x**) in the region of interest Ω:

where *p*
* _{j}*(

**x**) is the monomial basis function of the spatial coordinates, and

*m*is the number of the basis functions;

*a*

*(*

_{j}**x**) is the non-constant coefficient which can be determined by minimizing the following weighted discrete

**L**

_{2}norm

*J*(

**x**):

where *N*
* _{n}* is the number of the nodes in the region of interest Ω;

*w*(

**x**-

**x**

*) denotes the weight function related to the node*

_{i}**x**, and

**x**

*is a point in the support domain of*

_{i}**x**for which

*w*(

**x**-

**x**

*) 6=0; Φ*

_{i}*is the approximation to the value Φ(*

_{i}**x**) at the node

*x*

*, which is called as the generalized photon flux density.*

_{i}After the minimization of *J*(**x**), the following linear equation can be obtained:

where

Solving **a**(**x**) from Eq. (7) and inserting it into Eq. (5), we can obtain the following MLS approximation form:

The partial derivatives of *N*
* _{i}*(

**x**) are given by

where **B**
* _{i}*(

**x**) stands for the ith column of the matrix

**B**(

**x**).

The partial derivatives of *N*
* _{i}*(

**x**) are given by

where *s* represents the space variable *x*, *y* or *z*, and the comma indicates the partial derivative with regard to the spatial coordinate that follows.

From Eq. (12), we can see that the performance of the MLS approximation is governed by the basis function and the weight function. In this paper, the quadratic basis function

and the following quartic spline weight function are used in the three-dimensional case according to the references [33] and [34].

where *r*=‖**x**-**x**
* _{i}*‖/

*d*is the ratio of the Euclidean distance between the evaluation point

**x**and the sampling node

**x**

*to the radius of the support domain*

_{i}*d*. Removing all the variables correlated with

*z*-component in three-dimensionalGBMM procedure, we can obtain two-dimensional GBMM program easily.

## 2.2.2. Modified MLS shape function

Substituting **x**=**x**
* _{k}* back into Eq. (11), we have

where Φ designates the generalized photon flux density at all nodes in Ω, and *N*
* _{k}*=[

*N*

_{1}(

**x**

*),*

_{k}*N*

_{2}(

**x**

*), …,*

_{k}*N*

_{N}*(*

_{n}**x**

*)]*

_{k}**. And then, Eq. (16) can be furtherwritten as the followingmatrix equation:**

^{T}where $\widehat{\Phi}$ denotes the nodal photon flux density, and Λ is referred to as the transformationmatrix [35, 36]. They can be expressed as:

Therefore, the generalized photon flux density can be obtained from Eq. (17):

and

where Λ^{-1} is the inverse matrix of L. Incorporating Eq. (21) with Eq. (11), we have

where ${M}_{l}\left(\mathbf{x}\right)=\sum _{i=1}^{{N}_{n}}{N}_{i}\left({\mathbf{x}}_{k}\right){N}_{l}{\left({\mathbf{x}}_{i}\right)}^{-1}$ is called modified MLS shape function, and it satisfies the Kronecker delta function property:

## 2.2.3. Numerical implementation

According to whether modifying the MLS shape function, two meshless methods based on Galerkin approach are developed in this article.

Using Galerkin method and Gauss theory, Eqs. (1) and (2) can be transformed to the following weak form [30, 37]:

$$+{\int}_{\partial \Omega}\frac{1}{2A(\mathbf{x};n,n\text{'})}\Phi \left(\mathbf{x}\right)\Psi \left(\mathbf{x}\right)d\mathbf{x}={\int}_{\Omega}S\left(\mathbf{x}\right)\Psi \left(\mathbf{x}\right)d\mathbf{x}$$

where Ψ(**x**) is a test function from Sobolev space. The aforementioned Eqs. (11) and (22) can be rewritten in the following unified form:

where γ* _{l}*(

**x**) represents

*N*

*(*

_{l}**x**) in Eq. (11) or

*M*

*(*

_{l}**x**) in Eq. (22); Γ denotes the generalized or nodal photon flux density.

Substituting Eq. (25) into Eq. (24) and using γ* _{l}*(

**x**) as the test function, and then we can obtain the matrix equation as follows:

where the components of the matrices **K**, **C**, **F** and the vector **S** are given by

Since the matrix **G** is symmetric and positive definite [34], Γ can be uniquely determined from

When the MLS shape functions are employed, Γ that solved from Eq. (28) are only the generalized photon flux density Φ because the MLS approximation does not pass through the nodal function values according to the preceding subsection 2.2.1. In order to get the actual photon flux density, that is the nodal photon flux density $\widehat{\Phi}$ , at any point in the given domain W, we need to use the MLS approximation again with the formula:

where $\widehat{\Phi}$
* _{k}* is the nodal photon flux density at point

**x**

*;*

_{k}*N*

_{1}(

**x**

*),*

_{k}*N*

_{2}(

**x**

*), ‖,*

_{k}*N*

_{N}*(*

_{n}**x**

*) are the MLS shape functions without modification. However, the nodal photon flux density can be obtained directly through solving Eq. (28) when we use the modified MLS approximation.*

_{k}To sum up, the flowchart of the above two GBMM algorithms is shown in Fig. 1.

## 3. Experiments and results

To evaluate the developed GBMM algorithms in this paper, numerical and physical phantom experiments were performed respectively. Furthermore, the computational results based on GBMM were compared with the simulation data by FEM or MC and the measured photon flux density using a CCD camera. For the sake of convenience, GBMM1 represents the GBMM algorithm without modification, and GBMM2 indicates the other method in this section.

## 3.1. Numerical phantom experiments

#### 3.1.1. Homogeneous numerical phantom experiments

Firstly, a homogeneous tissue-like phantom was used in this experiment, and a light source with a total power of 0.125*nano*–*Watts* was placed in the phantom with its center at (6.5358,7.1537,15.0729). In GBMM study, 2178 random points were distributed in the above phantom, as showed in Fig. 2(a). In addition, *a priori* optical parameters were specified as *µ*
* _{a}*=0.035,

*µ*

*=6.0,*

_{s}*g*=0.9 and

*n*=1.37 to represent the diffusive biological tissue, which could be obtained from the literature [29]. Finally, the light exitance map on the phantom surface was solved using GBMM1 and GBMM2 procedures respectively, as presented in Fig. 2(b) and 2(c). Furthermore, the computational result of GBMM1 algorithm is completely identical to that of GBMM2.

In order to demonstrate the accuracy of GBMM, FEM was also employed to calculate the surface photon flux density. In FEM framework, the photon flux densities at 24654 discretized nodes were simulated, and then we used the interpolation method to determine the photon fluence at the arranged points in Fig. 2(a). Figure 3(a) and 3(b) give the volumetric mesh used in FEM and the corresponding simulation result respectively. The computational photon flux density based on GBMM was in good agreement with the numerical result by FEM, with the average relative error being about 0.9%, as showed in Fig. 3(c).

As we all know, MC method, regarded as the “*gold standard*”, is rigorous, flexible and powerful to study photon transport phenomena in the biological tissue, with which other numerical techniques are often compared [12, 18, 19, 20, 30]. Therefore, MC approach was also used to verify the performance of the above two GBMM algorithms in this paper. In the simulation experiment, a cubic light source of 2.0*mm* side length and 1.0*nano* – *Watts*/mm^{3} power density was placed in (11.0,11.0,11.0), and four cube phantoms centered at (10.0,10.0,10.0) with different side length from 5*mm* to 15*mm* were employed to obtain their respective surface photon fluxes. Table 1 lists the corresponding comparative data between GBMM1, GBMM2, FEM and MC method. RE1, RE2 and RE3 in Table 1 are the relative errors (RE) between the computational results by GBMM1, GBMM2, FEM and the simulation data using MC method respectively. Comparing the GBMM, FEMsolution with theMC simulation results, it has found that they have the same tendency, while the GBMM results are in better agreement with MC simulation data than the FEM solution. Furthermore, two GBMM methods in this paper have the same calculation precision.

## 3.1.2. Heterogeneous numerical phantom experiments

A cubical heterogeneous phantom with 8000*mm*
^{3} volume was utilized to further test the performance of the GBMM algorithms, and two cube light sources of 2.5*mm* side length and 1.0*nano* – *Watts*/*mm*
^{3} power density were located at (6.25,6.25,11.25) and (13.75,13.75,11.25) in the above phantom respectively. The optical parameters were assigned to each of the two components, as listed in Table 2. In Fig. 4 and 5, we illustrate the ability of our developed GBMM algorithms to simulate photon transport in tissue in the multiple light sources case. We can see that the simulation data by GBMM1 and GBMM2 is identical. Furthermore, the average relative error of the results obtained using GBMM and FEM is about 2.40%. In order to better demonstrate the performance of GBMM1 and GBMM2, we compared GBMM photon flux density with FEM calculational result along the detection square on the phantomsurface at heights 5*mm*, 10*mm* and 15*mm* from the bottom of the geometricmodel, as showed in Fig. 5(c)–5(e) respectively.

Similarly, four heterogeneous cubic phantomswith different side length from 11*mm* to 20*mm* were set up for numerical simulation usingMC method. A cube light source with a total volume of 8.0*mm*
^{3} and a power density of 1.0*nano*–*Watts*/*mm*
^{3} was embedded into the phantomswith its center at (11.0,11.0,11.0). Finally, the photon fluxes on the boundary of the phantoms were solved by GBMM1, GBMM2, FEM and MC respectively, see Table 3. There are no any difference between the computational results of GBMM1 and those of GBMM2, and the GBMM simulation data have better agreement with MC results than FEM solution. Furthermore, from Table 1 and Table 3, it can be seen that the relative error of photon flux between numerical techniques (GBMM, FEM) and statistical method (MC) on the phantomsurface decreases gradually with increase of the side length of the phantom, which is caused by that light propagation near the source is anisotropic [30].

In addition, finite element mesh generation before the assembly and solution of the FEM equations is an extremely burdensome task required by FEM in the experiment. For example, the time cost of discretizing the cubic phantom into 74861 nodes and 408561 elements using free meshing software Netgen [38] is 439.141*s*. Comparatively, arranging 19902511 nodes in the phantom only needs 0.4375*s* in our procedure, which shows a very faithful performance of the proposed algorithms. In FEM analysis, the corresponding time cost of FEM computation was about 57.0*s* for 61067 nodes and 331376 tetrahedral elements. However, for 68921 discretized points in the same phantom, GBMM1 spent approximately 25.0*s* for calculating the photon flux density based on the MLS interpolation because the computational result has been already smooth. Unfortunately, for large three-dimensional problems, the time cost of GBMM2 will increase largely because a matrix inversion step and two times matrix multiplication are added to modify the MLS shape function to satisfy the Kronecker delta function property. Using our GBMM2 procedure, 382.0*s* were required to determine the light power density in the phantom for the same meshless node distribution. Finally, it is worth mentioning that all the computational results in the section 3.1 have not been processed by any algorithm like normalization.

## 3.2. Physical phantom experiments

#### 3.2.1. Experimental preparation

Besides the above simulation studies using numerical phantoms, a series of physical phantom experiments were carried out for demonstrating the feasibility of the proposed methods in this contribution. In the optical imaging experiments, two cube resinous phantoms with 20*mm* side length were designed and manufactured. The two phantoms are both made from polyoxymethylene, and one or two small holes of 1.25*mm* radius were drilled in the phantoms to emplace the light source, as shown in Fig. 6. According to luminescence principle of luminescent light stick, peroxide, ester compound solutions and fluorescent dye were injected into the holes in the phantoms after thorough mix, and then the red light whose central wavelength located about 650*nm* was emitted due to the chemical reaction of the mixed resolutions. The aforementioned red light served as the internal light source in the physical phantomexperiments. However, the detectable photon flux density of the above light source on the phantom surface was very weak, reaching *pico*–*Watts*/*mm*
^{2} level. Therefore, a liquid-nitrogen-cooled back-illuminated CCD camera (Princeton Instruments VersArray 1300B, Roper Scientific, Trenton, NJ) was adopted to collect the transmitted and scattered near-infrared photons on surface of the phantom, and the CCD array can generate 1340×1300 pixels and 16 bits dynamic range images with 20*µm*×*20µm* sized pixel. The dark current of the camera is reduced largely through cooling the CCD chip down to -110°C using liquid nitrogen for long exposures. Furthermore, quantum efficiency (QE) of the CCD camera is greater than 80% for the wavelength range between 500*nm* and 700*nm*. In addition, the optical parameters of the phantomwere determined by a time-correlated single photon counting (TCSPC) system specifically constructed for the optical properties of the turbid medium [39]. The final measured optical parameters of the phantom at the wavelength around 660*nm* were listed as follows: the absorption coefficient *µ*
* _{a}*≈0.0002

*mm*

^{-1}and the reduced scattering coefficient

*µ*

^{′}

*=(1-*

_{s}*g*)

*µ*

*≈1.0693*

_{s}*mm*

^{-1}.

## 3.2.2. Data acquisition and analysis

All the physical phantomexperimentswere performed in a light-tight imaging chamber to avoid external disturbance and light leaking. Under the computer control, a motorized rotation stage was used to rotate the phantom for acquisition of the photon flux density on the four sides of the phantom, as illustrated in Fig. 6(b) and 6(d). Firstly, the phantom with a single light source was mounted on the sample stage, and the exposure time of each data acquisition was set to 30 seconds. Figure 7(a)–7(d) depict the collected luminescent light distribution from each side of the cubic phantom when the phantom was rotated three 90°, and the corresponding computational photon fluence by GBMM is showed in Fig. 7(e)–7(h). Because the GBMM1 and GBMM2 results were completely identical, which is demonstrated in Fig. 8, the simulated luminescent imageswere portrayed using alternative set of GBMMdata. As a result, the average relative error between the GBMM solution and the measured photon flux density was about 3.839%, which have similar tendency and good agreement.

Subsequently, the cube phantom with two light sources showed in Fig. 6(c) was utilized for luminescence imaging experiment. The photon energy distribution on the phantom surface was recorded using the CCD camera with an exposure time of 20 seconds, as orderly shown in Fig. 9(a)–9(d) along four directions separated by 90° presented in Fig. 6(d). As seen in Fig. 9(e)–9(h), the computed photon flux density based on GBMM was in good agreement with the experimental counterpart. Moreover, the computational results by GBMM1 and GBMM2 were the same. As it is expected, GBMM worked well with the average relative error about 5.092%, as shown in Fig. 10.

## 4. Discussion and conclusion

Optical molecular imaging has been rapidly developed and used to study physiological and pathological processes *in vivo* at the cellular and molecular levels for its unique advantages in sensitivity, specificity and directness. The research of photon propagation in the biological tissue is necessary because it helps not only to improve the spatial resolution of inverse source reconstruction but also to test the accuracy of the physical model. However, the application scope of the existing analytical and statistical techniques is restricted for their inherent characteristic. Therefore, numerical techniques should be studied deeply for better applicability and higher efficiency, especially for the more complex geometries. In this paper, we have presented two kinds of Galerkin-based meshless methods to solve photon transport problem in the heterogeneous biological tissue. In comparison with conventional FEM analysis or MC approach, GBMM algorithms have the following novel features.

First, although many excellent and powerful automatic mesh generators are available, meshing is computationally burdensome and expensive for three dimensional arbitrarily shaped domains with complex internal structure, especially when the geometric model is discretized into more than 10^{6} elements. Moreover, the transformation fromthe geometricmodel to the finite element data is also computationally arduous. However, two kinds of GBMMmethods developed in this paper do not require any mesh to discretize the problem, in which the photon flux density is approximated entirely in terms of a group of scatter nodes, and the discrete equations are constructed without consideration of element connectivity information and interrelationship among the nodes. Therefore, GBMM methods can avoid the time-consuming three-dimensional finite element mesh generation compared with FEM. Furthermore, irregular complex objects with heterogeneous internal structure can be represented more accurately, especially for the inhomogeneous distribution of optical properties in the complicated biological tissue. In addition, the adaptive GBMM method can be realized easily through inserting some nodes in the high gradient region of the field function or deleting several points in the low gradient area, but the adaptive FEM needs to remesh the domain of the problem based on error estimation after every iteration.

Secondly, when *in vivo* small animal experiment is performed, the living small animal is a four-dimensional object because it is impossible to maintain absolutely stationary state, that is, there may be some surface shape variations and involuntary internal structure motions, such as the cardiac and respiratory motions. For this moving boundary conditions case, FEM analysis can not monitor the time-course of the above motions for the expensive computational cost of remeshing the problem domain, which will lead to discontinuity as well as a degradation of accuracy. On the contrary, GBMM can simplify the problem with moving boundary conditions because no element connectivity is needed.

Lastly, two kinds of GBMM algorithms are developed for photon propagation in the biological tissue in this paper. Both of them are based on MLS approximation. However, the MLS approximated function does not pass through the nodal values, so we need to use MLS approximation again to obtain the actual nodal photon flux density. In other words, the MLS interpolation lacks the Kronecker delta function property of the conventional FEM shape functions. To ameliorate this difficulty, the MLS shape functions are modified in one of the aforementioned GBMM algorithms (GBMM2), with which we can not only calculate the nodal photon flux density directly but also impose Dirichlet and Robin boundary condition reported in the literature [30] accurately. Unfortunately, for large three-dimensional problems, the time and memory cost of the modified GBMM algorithm will increase slightly because a matrix inversion step and two times matrix multiplication are superinduced. Moreover, the GBMM results are smooth due to the high order smoothness of the MLS shape functions, which can simplify post-processing largely. However, some post-processing techniques are necessary to obtain the smooth results in the framework of FEM analysis.

In conclusion, we have presented two kinds of GBMM algorithms to simulate photon transport in the biological tissue. Furthermore, the feasibility of the proposed methods has been demonstrated with a series of numerical and physical phantom experiments. The preliminary results show that meshless methods are of huge potential in the optical engineering field. Our future work will focus on employing the GBMM algorithms to reconstruct the internal light source. The corresponding results will be reported later.

## Acknowledgments

This paper is supported by NBRPC (2006CB705700), PCSIRT (IRT0645), CAS HTP, CAS SREDP (YZ0642, YZ200766), 863 Program (2006AA04Z216), JRFOCYS (30528027),NSFC (30672690, 30600151, 30500131, 60532050, 60621001) and BNSF (4071003) in China.

## References and links

**1. **V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole body photonic imaging,” Nat. Biotechnol. **23**, 313–320 (2005). [CrossRef] [PubMed]

**2. **G. Wang, W. Cong, H. Shen, X. Qian, M. Henry, and Y. Wang, “Overview of bioluminescence tomography–a new molecular imaging modality,” Front. Biosci. **13**, 1281–1293 (2008). [CrossRef]

**3. **S. Bhaumik and S. S. Gambhir, “Optical imaging of Renilla luciferase reporter gene expression in living mice,” Proc. Natl. Acad. Sci. USA **99**, 377–382 (2002). [CrossRef]

**4. **T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. **17**, 545–580 (2003). [CrossRef] [PubMed]

**5. **W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt. **6**, 432–440 (2001). [CrossRef] [PubMed]

**6. **E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. **30**, 901–911 (2003). [CrossRef] [PubMed]

**7. **A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express **13**, 9847–9857 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-24-9847. [CrossRef] [PubMed]

**8. **C. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. **4**, 235–260 (2002). [CrossRef] [PubMed]

**9. **G. Wang, H. Shen, W. Cong, S. Zhao, and G. W. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express **14**, 7852–7871 (2006), http://www.opticsinfobase.org/abstract.cfm? URI=oe-14-17-7852. [CrossRef] [PubMed]

**10. **V. Y. Soloviev, “Tomographic bioluminescence imaging with varying boundary conditions,” Appl. Opt. **46**, 2778–2784 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=ao-46-14-2778. [CrossRef]

**11. **Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finite element analysis: methodology and simulation,” Phys. Med. Biol. **52**, 4497–4512 (2007). [CrossRef] [PubMed]

**12. **A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. **50**, R1–R43 (2005). [CrossRef] [PubMed]

**13. **W. Cong, A. Cong, H. Shen, Y. Liu, and G. Wang, “Flux vector formulation for photon propagation in the biological tissue,” Opt. Lett. **32**, 2837–2839 (2007), http://www.opticsinfobase.org/abstract. cfm?URI=ol-32-19-2837. [CrossRef] [PubMed]

**14. **Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express **14**, 8211–8223 (2006), http://www.opticsinfobase. org/abstract.cfm?URI=oe-14-18-8211. [CrossRef] [PubMed]

**15. **A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express **12**, 5402–5417 (2004), http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-22-5402. [CrossRef] [PubMed]

**16. **A. Joshi, W. Bangerth, A. B. Thompson, and E. M. Sevick-Muraca, “Adaptive finite element methods for fluorescence enhanced frequency domain optical tomography: forward imaging problem,” IEEE International Symposium on Biomedical Imaging (ISBI 2004) **2**, 1103–1106 (2004).

**17. **W. Cong, D. Kumar, Y. Liu, A. Cong, and G. Wang, “A practical method to determine the light source distribution in bioluminescent imaging,” Proc. SPIE **5535**, 679–686 (2004). [CrossRef]

**18. **L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML-Monte Carlo modeling of photon transport in multi-layered tissues,” Comput. Meth. Prog. Biomed. **47**, 131–146 (1995). [CrossRef]

**19. **D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express **10**, 159–169 (2002), http://www.opticsinfobase.org/abstract.cfm?URI=oe-10-3-159. [PubMed]

**20. **H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with Monte Carlo method,” Acad. Radiol. **11**, 1029–1038 (2004). [CrossRef] [PubMed]

**21. **W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express **13**, 6756–6771 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-18-6756. [CrossRef] [PubMed]

**22. **Y. Lv, J. Tian, H. Li, J. Luo, W. Cong, G. Wang, and D. Kumar, “Modeling the forward problem based on the adaptive FEMs framework in bioluminescence tomography,” Proc. SPIE **6318**, 63180I (2006). [CrossRef]

**23. **S. R. Arridge, H. Dehghani, M. Schweiger, and E. Okada, “The finite element model for the propagation of light in scattering media: A direct method for domains with nonscattering regions,” Med. Phys. **27**, 252–264 (2000). [CrossRef] [PubMed]

**24. **R. H. Bayford, A. Gibson, A. Tizzard A, T. Tidswell, and D. S. Holder, “Solving the forward problem in electrical impedance tomography for the human head using IDEAS (integrated design engineering analysis software), a finite element modelling tool,” Physiol. Meas. **22**, 55–64 (2001). [CrossRef] [PubMed]

**25. **S. J. Koopman, A. C. Harvey, J. A. Doornick, and N. Shephard, *Stamp 5.0: structural time series analyser, modeller and predictor*, (The Manual. Chapman & Hall, London, 1995).

**26. **I. V. Singh, K. Sandeep, and R. Prakash, “The element free Galerkin method in three dimensional steady state heat conduction,” Int. J. Comput. Eng. Sci. **3**, 291–303 (2002). [CrossRef]

**27. **I. V. Singh, “Parallel implementation of the EFG method for heat transfer and fluid flow problems,” Adv. Eng. Software **34**, 453–463 (2004).

**28. **T. Belytschko, L. Gu, and Y. Y. Lu, “Fracture and crack growth by element-free Galerkin methods,” Modelling Simul. Mater. Sci. Eng. **2**, 519–534 (1994). [CrossRef]

**29. **G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys.Med. Biol. **50**, 4225–4241 (2005). [CrossRef] [PubMed]

**30. **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,” Med. Phys. **22**, 1779–1792 (1995). [CrossRef] [PubMed]

**31. **W. G. Egan and T. W. Hilgeman, *optical properties of inhomogeneous materials*, (Academic, New York, 1979).

**32. **T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin method,” Int. J. Numer. Methods Eng. **37**, 229–256 (1994). [CrossRef]

**33. **J. Dolbow and T. Belytschko, “An introduction to programming the meshless element free Galerkin method,” Arch. Comput. Methods Eng. **5**, 207–241 (1998). [CrossRef]

**34. **X. Zhang and Y. Liu, *Meshless methods*, (Tsinghua University Press, Beijing, 2004).

**35. **J. S. Chen and H. P. Wang, “New boundary condition treatments in meshfree computation of contact problems,” Comput. Methods Appl. Mech. Eng. **187**, 441–468 (2000). [CrossRef]

**36. **S. Li, W. Hao, and W. K. Liu, “Numerical simulations of large deformation of thin shell structures using meshfree methods,” Comput. Mech. **25**, 102–116 (2000). [CrossRef]

**37. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. **20**, 299–309 (1993). [CrossRef] [PubMed]

**38. **J. Schöberl, “Netgen an advancing front 2D/3D-mesh generator based on abstract rules,” Comput. Visual. Sci. **1**, 41–52 (1997). [CrossRef]

**39. **D. Qin, H. Zhao, Y. Tanikawa, and F. Gao, “Experimental determination of optical properties in turbid medium by TCSPC technique,” Proc. SPIE **6434**, 64342E (2007). [CrossRef]