1. Introduction

This paper focuses on near-infrared diffuse optical tomography and its main application in cancer research and small animal imaging (mice in our work) [1]. More specifically, fluorescence-enhanced Diffuse Optical Tomography (fDOT) enables to estimate three-dimensional (3D) locations and geometries of areas targeted by fluorescent markers, such as tumors in a small animal. The physical parameters estimation, that is called the inverse problem or reconstruction, is usually based on the minimization of the error (in the least-squares sense) between measured outputs and outputs computed through a parameterized set of Partial Differential Equations (PDEs). This latter modelling step, e.g. the resolution of the PDEs, is the forward problem.

In the case of simplified geometries (semi-infinite media for example), the forward problem can be solved with analytical solutions [2]; however these solutions require a complex acquisition protocol, for example by immersing the small animal in an optical index matching liquid. A numerical resolution of the diffusion equation using for example the polynomial Finite Element Method (FEM) [3], is generally used [4], because it has no restriction in terms of neither the geometry nor the homogeneity of the system. However it is well-known to be time-consuming and computation heavy. Then solving the adjoint problem generally speeds up treatment and allows a full 3D reconstruction [5]. Nevertheless, since the accuracy of the results is mesh dependent, the matrices dimensions involved in the discretized version of the PDEs are large and the resolution of the problem remains long.

Our ultimate goal being to process experimental data in real-time, the principal objective of the present work is to study an alternative numerical method in order to speed up the computation of the forward model, that takes much more time than the reconstruction stage (20 or 30 times longer if we choose an accurate method such as the FEM: the inverse problem takes less than 3 minutes if we use a linearization method, the forward problem 40 or 50 minutes in the case of an object whose size is similar to the a mouse). Knowing that the mesh generation in the standard FEM is very costly in time, especially in 3D, a wavelet-Galerkin method is chosen, in particular because it is a meshless method. A regular grid is thus established, and wavelets are chosen as interpolation functions, instead of the polynomial functions of the standard FEM [3]. This method is interesting because it allows an accurate solution comparable to the one obtained with the classical FEM, while saving computation time.

It is well-known that wavelets are an efficient mathematical tool to solve PDEs [6] and the physical applications are numerous. The wavelets can be used in the collocation method [7] or in the Galerkin method [8]: these meshless methods are alternative methods to the standard FEM in the resolution of all kinds of PDEs. A large mathematical theory shows significant advantages of the wavelets compared to the classical polynomial functions of the FEM: G. Oberschmidt et al. [9] proved that wavelets allow an important reduction of the number of non-zero coefficients in the system matrix. Due to the properties of the wavelet basis, wavelet-based numerical methods reach both good spatial and spectral resolution [10]. Another advantage is the multi-scale capability of wavelets [11], that allows a finer resolution in the Regions Of Interest (ROIs). Finally the absence of mesh allows an efficient time reduction compared to the standard FEM.

In the field of tomography, the wavelet-based methods are principally used in the inverse problem. For example, wavelets are frequently used in X-ray tomography. Their application to the Radon transform [12,13] allows for a multi-scale reconstruction to reduce noise and to obtain an approximate reconstructed image without loosing quality. In optical tomography, wavelets are sometimes used for the purpose of denoising, compressing data [14], and regularizing reconstruction: Zhu et al. [15] use the wavelet-multiresolution approach to reconstruct (at coarse approximation) the unimportant regions of the studied domain, and to reconstruct in detail the ROIs. Computation time is thus reduced without quality loss in ROIs. B. Kanmani et al. [16] describe the wavelets’ insertion into the resolution of the inverse problem in optical tomography: the reconstruction equation is projected on to a wavelet basis, which reduces noise and computation time.

Galerkin method has been sometimes used in the domain of fDOT, to solve the forward problem, with Daubechies wavelets [17] (where the wavelet-Galerkin method is used to solve the general diffuse optical tomography equation, for finding the value of physical parameters), splines [18] or exponential functions [19].

In the present work, Daubechies wavelets are chosen [20]. The solutions of the PDEs are approximated by Daubechies scaling functions of order 3, instead of the classical polynomial functions used in the standard FEM. The weak form of the PDEs is then projected in the same basis of wavelets (Galerkin method). The stiffness matrix obtained in the matrix form of the PDEs is then easily computable thanks to the properties of the wavelets, and the PDEs are quickly solved.

The use of wavelets are however possible only on simple geometries, such as a rectangle in 2D and a parallelepiped in 3D. In order to extend the wavelet-Galerkin method to any geometry, the fictitious domains method [21] is chosen. This allows to satisfy the boundary conditions at the same time.

In this paper, we propose to solve the forward problem in fluorescence diffuse optical tomography by a wavelet-Galerkin method. In order to apply this method to any complex geometry, the fictitious domain method is chosen. Some simulations in 2D and 3D illustrate the feasibility and the effectiveness of the method. The results of these simulations indicate a time reduction in the solving of the PDEs, that can reach a factor 70 compared to the FEM.

The paper is organized as follows: section 2 is devoted to the description of the forward problem in fluorescence diffuse optical tomography and its resolution by a general numerical method, and a brief description of the inverse problem is given. In section 3 we present and justify the choice of Daubechies wavelets, then we describe the chosen meshless wavelet-Galerkin method, applied to the diffusion equation implied in the forward problem and to the associated boundary conditions. In section 4 we present the results obtained with simulations in 2D- and 3D-domains, and the comparison with the results obtained with a classical polynomial FEM. A discussion about these results is proposed in section 5. Finally the conclusion is drawn in section 6.

2. Equations in optical tomography and resolution methods

In this section, we begin with a reminder of the basic equations governing light propagation through diffusive media, within the diffusion approximation and more specifically in the case of the fDOT. In the presented theory, discussion is restricted to a continuous wave illumination system, although this work can be extended to other optical imaging approaches. Some details of a general numerical method chosen to solve the forward problem are provided in the second subsection, and the inverse problem is presented on the third subsection.

2.1 Analytical expression for the fluorescence signal in the framework of the fluorescence Diffuse Optical Tomography

The objective of fDOT consists in reconstructing in 3D internal optical property maps of the tissues and the local concentration of fluorescent markers (fluorochromes). This reconstruction is based on an accurate forward model of coupled excitation and emission light propagation through the studied medium, which is highly scattering. Consider a diffusing medium with a total volume Ω, and bounded by the surface ∂Ω. Suppose a point source located at the position r _s inside the diffusing medium with intensity q(r _s). The fluorochromes present within the medium at position r are excited by the incident Diffuse Photon Density Wave (DPDW) u_x at the excitation wavelength λ_x (typically λ_x~700nm). Each excited fluorescent particle then acts as a secondary point source at position r and generates a fluorescent DPDW u_m, propagating with a longer wavelength λ_m. For simplicity we suppose λ_m=λ_x; the light propagation through highly scattering media is commonly modeled, within the Diffusion Approximation (the radiative transfer equation is more adapted for small animals imaging, but this approximation is often used), by the following set of weakly coupled PDEs [4]:

The first equation represents the propagation of the excitation light (subscript x), from the point-like source, located at position r _s, to the fluorochrome at a position r inside the diffusing volume Ω. The second models the generation at position r and the propagation of the emitted light (subscript m), observed in the position r’. The coefficient µ_a is the total absorption in the volume and D is the diffusion coefficient; for simplicity, these two coefficients are supposed to be constant in this work. The parameter β(r) involved in the fluorescence source term is the conversion factor of the excitation light, at a point r, into fluorescence light, and depends on the intrinsic properties of the fluorochromes, such as the quantum yield and also their local concentration. Hence, this parameter is precisely the parameter to be recovered. The other physical parameters µ_a and D are considered to be known. In the time resolved domain, with a more experimental method, one can also recover these parameters by generalizing the technique [22]. Since the parameter β(r) is directly proportional to the local concentration, it allows the recollection of the position of the fluorochromes as well as their quantification.

Eq. (1) are subject to the Robin boundary conditions, which have to be verified at all the boundaries of the system ∂Ω:

where n is the outward unit vector normal to the surface, and ζ_m (supposed to be equal to ζ_m) is a multiplicative factor, self-consistently accounting for index mismatching at the boundaries and depending on the Fresnel reflection coefficients [23]. We note ζ_x=ζ_m=ζ.

By using the adjoint method [5], it can be shown that the observed fluorescent DPDW u^obs_m(r _d,r _s) m u at detector position r _d for a source at position r _s, that is the contribution of all the values u_m(r′=r _d,r,r _s) is given by [22]:

where u_x is the solution of the first equation of the system (1), and g_m is the adjoint solution, which verifies the following equation:

rd being the position of a fictitious source placed at a detector position. The forward problem consists in solving the first equation of (1) to determine ux, then Eq. (4) to determine g_m. Hence these equations will be written with the generic form:

The inverse problem consists in solving Eq. (3) to determine the parameter β(r), knowing the measures u^obs_m(r _d,r _s) and the solutions of the forward problem u_x and g_m.

2.2 Numerical method and matrix form of the PDEs in the forward problem

Usually the medium has a complex inhomogeneous geometry, so a numerical resolution is definitely required for solving (5): the Finite Element Method (FEM) is the method of choice. Within the FEM, the volume object Ω is divided into a limited number N of elements which are joint at vertices. In 2D, these elements can be triangles or rectangles. In 3D the elements can be tetrahedrons or parallelepipeds.

Two steps are considered in the FEM [3]:

- First is the interpolation step, where an approximative solution ũ(ũ∊{ũ_x,g̃_m} in the fDOT case) of the PDE (5) is defined by linear combination of interpolation functions {f^interp_k}_k=1…N :

The numerical solution is the vector $\mathbf{u}=\left(\begin{array}{c}\multicolumn{1}{c}{u_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{u_N}\end{array}\right)$ u corresponding to the weight values for each of the N interpolation functions [24]. q is approximated in the same way:

- Then the projection step consists in projecting the residual of the equation, defined by R=Δ∇.(-D∇ũ)µ _aũ-q̃, which must be minimized, in a projection functions basis {f^proj_l}l=_1…N:

We summarize in the Table 1 the different existing classical numerical methods classified in terms of interpolation/projection functions. The Galerkin projection method [3] is used here to express the equations in a weak form: the basis of projection functions is equal to the basis of interpolation functions:, f^interp_k=^fproj_k, ∀k=1,…,N.

Table 1. Numerical methods by the Finite Elements Methods.

View Table | View all tables in this article

By inserting Eq. (6) and (7) in Eq. (8), the PDE is then reformulated into a linear compact matrix form:

where $\mathbf{q}=\left(\begin{array}{c}\multicolumn{1}{c}{q_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{q_N}\end{array}\right)$ is the vector of weight values given by (7), P _q is the projection matrix whose coefficients are defined by the integrals ∫f ^interp _k f ^proj _l, and A is called the stiffness matrix and accounts for the different optical parameters and contains integrals of the form ${\int }_{\Omega }{f}_k^{{\mathrm{interp}}^{\left(m\right)}}{f}_l^{{\mathrm{proj}}^{\left(m\right)}}$where f ^(m) is the m ^th derivative of f, what we will describe in detail in section 3. A is generally a sparse and quasi-diagonal matrix whose size depends on the number of degrees of freedom in the FEM [3]. The different quantities are generally approximated by piecewise quadratic functions, used in particular by the commercial software used in the present work [25]. In this case, a triangular (or tetraedric in 3D) mesh is often required to obtain a good accuracy.

2.3 The inverse problem

The inverse problem is the resolution of Eq. (3) in order to estimate β (r),∀r. When discretizing the medium into elementary volumes, (3) readily appears as a linear function of the parameters β. In order to introduce expression (3) into a reconstruction scheme, it is usual to reformulate the problem in terms of a weight matrix (discretization numerical implementation):

where the weight matrix W depends on the values of u_x and g_m, which are estimated in the forward problem step. Each value of the vector u ^obs _m represents the measurement on one of the detectors for a given point source. β is the vector of the unknown parameters to be determined. The inversion method used in this work is the Algebraic Reconstruction Technique (ART), which is an iterative method [26].

3. A wavelet-Galerkin method to solve the forward problem

There are a few motivations to replace polynomial functions by wavelet or scaling functions: first the solution will have both good spatial and spectral resolution, not only spatial in the case of the standard FEM, or only spectral in the case of spectral methods; second the wavelet and scaling functions have a multi-scale capability, that allows a faster and finer resolution in the ROIs; third the functions are continuous on the whole domain, so there is no equation inter-element to solve, contrary to the standard polynomial FEM [3]; fourthly they are separable, that implies an easy implementation in 2D and 3D; and finally, if the chosen interpolation/projection functions of a Finite Element Method are wavelets, the nodes of the mesh must be equidistant in order to keep the orthonormality of the wavelet basis, so the solution is approximated in a regular grid (rectangular or parallelepipedic mesh, where all elements have the same size): with the wavelets, that is sufficient to obtain a good accuracy in the solution, and this regular grid has the advantage to gain a lot of time compared to the definition of triangular (or tetraedric in 3D) mesh.

3.1 Choice of appropriate wavelets

3.1.1 General definition and construction of the wavelets

The concept of wavelets, introduced for the first time in the 80’s by the geophysicist J. Morlet [27], has been widely applied for signal decomposition and approximation [11]. Let us briefly give some properties of the wavelet basis:

If L²(R) is the space of finite energy signals, a multiresolution analysis is a succession of subspaces V_j in L²(R) which verify the following conditions:

(i) ∀j∊ℤ V_j+1⊁V_j;

(ii)lim _j→∞ V_j={0} and lim _j→-∞ V_j is dense in L²(R);

(iii) s ∊ V_j ⇔ s(./2) ∊ V _j+1;

(iv) ∃Φ, “scaling function”, such that {Φ(.-k),k∊ℤ} forms an orthonormal basis of V₀. This implies that ${\left\{{\Phi }_{j,k}\right\}}_{k\in \mathbb{Z}}$ is an orthonormal basis of V_j, where

and it exists a sequence {h(n)}_n∊ℤ, the scaling filter, such that:

We have furthermore the normalisation condition: ∫Φ=1.

3.1.2 Choice of Daubechies scaling function

Daubechies scaling function Φ associated with the wavelet function Ψ of order M [20] is often chosen to solve the PDEs [17, 28], because they have the following properties:

- Wavelets are directly issued from the multiresolution analysis [11] and they form an orthonormal basis of L²(R);

- The wavelet functions Ψ have M vanishing moments, that means that it reproduces polynomials up to degree M; the more M will be large, the more Φ will be regular and the solution of the PDE better approximated;

- Φ is compactly supported, due to the construction of the filter h: if h(n)=0 for n<0 and n>2M-1, support(Φ)=support(Ψ)=[0,2M-1]=[0,p].

Daubechies wavelets have the important property to offer the smallest support for a given number M of vanishing moments. In fact, the smaller the support is, the faster the calculations, in this basis, are. However, the greater the number of vanishing moments is, the more regular the wavelet and the scaling function are, but the larger the wavelet support will become. So a trade-off has to be found: Daubechies constructed its basis by optimizing these two properties. Daubechies scaling function with M=3, the chosen order in this work, is presented in Fig. 1: the support is not too large (p=2M-1=5) and three vanishing moments lead to a wavelet function not too irregular.

 

Fig. 1. Daubechies scaling function of order M=3.

Download Full Size | PDF

3.2 Resolution of the diffusion equation by the wavelet-Galerkin method

3.2.1 General case of Galerkin method

Resume the diffusion Eq. (5) and write its integral form (in 1D) obtained from Eq. (6), (7) and (8):

After an integration by parts of the first term and by applying the Green-Ostrogradski theorem, then including the Robin boundary conditions (2) (which have the advantage to insert naturally in the integral form of the PDE [29], because they are necessary and sufficient), we obtain, by placing in a Galerkin scheme (f^nterp_k=f^proj_k=f_k=, ∀k=1,…N), the weak form of the PDE:

This set of N equations can then be reformulated into a linear compact matrix form (9).

3.2.2 Weak form expression with Daubechies scaling function

Now we introduce the Daubechies scaling function Φ and the corresponding basis f_k=Φ_j,k=2-^j/2Φ(./2^j-k) as interpolation/projection functions. Let ũ_j be an approximation at scale 2^j (approximation level j) of the exact solution of (5). Interpolation step is written under the following general form (in 1D):

where

N_x is the meshes number in a 1D regular grid (dimension x), that has to be determined according to the geometry of the studied domain and the chosen scale 2^j: the more j will be large, the more the grid and then the solution will be coarse; p=2M-1 is the wavelet support size, p=5 in the case of Daubechies order 3. Figure 2(c) gives a simple schema of used scaling functions that are included in the studied domain Ω=[0 8] (N_x=8): the solid curves represent the functions which are totally included in Ω, that is Φ_j,k for j=0 and k=0 to N_x-p=3. One talks here about meshless numerical method, because the mesh is a simple regular grid, in order to preserve the basis orthonormality.

Eq. (14) becomes:

At this stage, we see that three types of integrals appear: ∫ΩΦ_j,kΦ_j,l, ∫Ω∇Φ_j,k∇Φ_j,l and ∫∂ΩΦ_j,kΦ_j,l. Such derivative wavelet product integrals, called the connection coefficients, are defined in the unbounded case by the Eq. (18):

where ${\Phi }^{\left(m\right)}=\genfrac{}{}{0.1ex}{}{{\partial }^m\Phi \left(x\right)}{{\partial x}^m}$. They can be accurately calculated from the scaling filter h(n) (Beylkin’s theory) [30]. Coefficients tables are sometimes given for different wavelet orders M and for different derivation orders m ₁, m ₂ [31], and can be applied generally to each approximation level j; in the following, we choose j=0 in order to simplify the notations. The tables that are useful in the diffusion problem are reported in Appendix, for M=3. Note that because of the orthonormal specificity of Daubechies scaling functions, we have:

In the case of optical tomography, studied domains are not infinite, and we need to compute connection coefficients for intervals smaller than the wavelet support:

where n is an integer. Γ⁰⁰ _k,l[0;1]=∫¹ ₀Φ(x-k)Φ(x-l)dx can be calculated in different ways, for example by [31, 28]. Then we can show that (more details in Appendix):

Finally the Eq. (17) can be written on the matrix form Au=P_qq. Here A is the stiffness matrix whose terms depend on the connection coefficients and on the physical parameters, u is the vector of the N unknowns u_j,k, P_q is the projection matrix whose coefficients are defined by the integrals ∫ΩΦ_j,kΦ_j,l, q is the vector of the N coefficients q_j,k of the approximated source function q̃ defined by Eq. (15).

In 2D and 3D, the previous equations can be generalized by tensor products [32]. Thanks to the separability property of the wavelets, the interpolation functions become for example in 2D:

Eq. (17) can be generalized in 2D and 3D. The meshes number becomes N_x×N_y in 2D, N_x×N_x×N_z in 3D, that depends on the geometry of the studied domain and on the meshsize, i.e. the approximation level j.

 

Fig. 2. (a) Schema of a fictitious domain (rectangle Ω’) larger than the studied domain Ω in 2D. (b) Schema of the discretized approximation of edges ∂Ω. (c) Schema of a fictitious domain in 1D and corresponding scaling functions Φ_j,k=2^-j/2Φ)./2^j-k) with j=0 and k=-4 to 7.

Download Full Size | PDF

3.3 Fictitious domain method

Wavelet basis used in PDEs solving have however the inconvenient to be applicable only on simple geometries, that is a rectangle in 2D and a parallelepiped in 3D. A fictitious domain Ω’ larger than the studied domain Ω must be defined (Fig. 2(a)): it is a rectangle in 2D and a parallelepiped in 3D, allowing thus the use of the wavelets. The PDEs are solved into this new domain [21]. The parameters k and l of Eq. (17) are consequently not defined from 0 to N_x-p (in dimension x), but from (-p+1) to N _x-1, that means p-1 additional scaling functions on either side of the studied domain, for each dimension. The Fig. 2(c) shows a schema in 1D of fictitious domain and of the p-1 added functions on the two sides: in this example, the 8 added functions are schematized by the dotted curves, for k=(-p+1)=-4 to (-1) and k=(N_x-p+1)=4 to (N_x-1)=7.

Only the results in the region of interest Ω are considered.

We must also pay attention to the boundary conditions, that have to be verified on ∂Ω and not on ∂Ω′. The boundary conditions are furthermore inserted during solving the PDEs (Eq. (14)), that allows to include the actual physical boundaries into the fictitious domain. The edges ∂Ω of the studied domain have to be approximated by small regular edges (Fig. 2(b)). The first term of Eq. (17) which depends on ∫∂ΩΦ_j,kΦ_j,l is consequently decomposed on each edge. For example, if the domain is a 2D-rectangle delimited by x=0, x=L_x, y=0 and y=_Ly, e.g. 4 edges, this integral becomes:

Eq. (23) is applied on any complex edges system ∂Ω. But we observe that the boundary integrals ∫∂ΩΦ_j,kΦ_j,l on these small edges require a specific calculation, thanks to the connection coefficients defined by (20) and (21), and that allows to reach a good accuracy in the PDEs resolution.

4. Numerical simulations: comparison with the polynomial FEM

Two 2D-simulations were studied: first a simple rectangle of size 4×3 cm², then an approximated irregular geometry, whose edges are small segments. These geometries are inserted in a fictitious rectangular domain (Fig. 3(a) and 3(b)) [33]. The chosen physical parameters are D=0.0327 cm (diffusion coefficient), µ_a=0.2 cm^-1 (absorption coefficient), which are typical diffusion and absorption coefficients in biological tissues with the used wavelength λ_x [34], and the source function is Dirac function: q=1 on the source node, 0 elsewhere. ζ=½ is the multiplicative factor in the boundary conditions (Eq. (2)). These 2D-objects are simulated with 2 fluorescent elliptic zones, simulating the fluorophores positions (Fig. 3(a) and 3(b)).

A third simulation with a simple cube, was chosen (Fig. 3(c)) to illustrate the feasibility in 3D. The corresponding fictitious domain is a larger cube. The fluorescent zone is an ellipsoid zone inside the cube.

 

Fig. 3. (a)(b) Schemas of the 2D geometries: black points: positions of the source; gray points: detectors; ellipses: 2 simulated fluorophores; gray grid: regular grid for a scale 2–³; large rectangle: fictitious larger domain; (c) schemas of the 3D geometry: ellipsoid: simulated fluorophore. (d) Example of result for the forward problem (ux for the source #8 (surrounded in (a)), in the first simulation), by wavelet-Galerkin method, analytical modelling and FEM. (e) (h) Results of the rectangular simulation. (f)(i) Results of the irregular geometry simulation. g) (j) Results of the cube simulation. (e) (f) (g) Results after reconstruction, with a forward modelling by wavelet-Galerkin method and a reconstruction by ART; (h) (i) (j) Results with a forward modelling by standard FEM method and a reconstruction by ART.

Download Full Size | PDF

To obtain the simulated emission signal u^obs_m, the finite element software COMSOL Multiphysics^™ [25] is used. The simulated detectors are located on the physical surface of the object (gray points on Fig. 3(a) and 3(b) for the 2D cases and 3g for the 3D case, placed on the top of the objects). For the source positions simulation, a classical model to describe the interaction between the medium and a laser beam is to consider an isotropic point source at distance $\genfrac{}{}{0.1ex}{}{1}{{\mu }_s^{\prime }}$ (where µ′_s is the reduced scattering parameter [35]) from the boundary ∂Ω, inside the domain.

For simplicity, the virtual sources are placed in these simulations on grid nodes (black points on Fig. 3(a) and 3(b) for the 2D cases and black crosses on Fig. 3(g), placed at the bottom of the objects); but the generalization to another distance is possible by interpolation of the source function q on neighbouring nodes.

For the forward model, the studied domain is divided according to a regular grid. In the case of the irregular 2D geometry, the fictitious domain is the same as the one of the rectangle defined in the first simulation (Fig. 3(a) and 3(b)). In this work, we chose j=-3 in the whole studied domain, so we have 2³=8 nodes for 1 cm.

Then the solutions u_x and g_m of the forward model (5) are computed; the Fig. 3(d) gives an example of u_x values along a cross section represented by the line arrow on Fig. 3(a): 3 curves are presented: first the result obtained by the FEM, second the result obtained by an analytical calculation (with a semi-infinite medium hypothesis), third the result obtained by the wavelet-Galerkin method.

Finally the inverse problem is solved by using the ART method [26]. The reconstructions of the spatial distribution of parameter β, are represented with a gray level scale on Fig. 3(e), 3(f), 3(g), 3(h), 3(i) and 3(j). The results are normalized: the white color represents the value 0, the black color represents the normalized maximal value 1.

The parameters for each of these simulations are summarized in the Table 2. The results are presented with a resolution of the forward model first by a wavelet-Galerkin method (noted W-G), second by a classical polynomial FEM.

Table 2. Description of the 3 simulations

View Table | View all tables in this article

4.1 Simulation with a rectangle

The simple rectangle is of size 4×3 cm², with two elliptic zones simulating the positions of the fluorophores (Fig. 3(a)). The sources are positioned on the bottom horizontal line, the detectors on the top horizontal line (Fig. 3(a)).

The meshsize is determined by the scale 2^j of the scaling function. x takes the values 0, 2^j, 2×2^j, …, 4. So we have in the fictitious domain N_x=(4×2^-j+1)+4. The presented results are obtained with j=-3 in the case of wavelet-Galerkin method (N_x=37 and N_y=29), and with an equivalent mesh size in the case of classical FEM.

Figure 3(e) gives reconstructions obtained with a wavelet-Galerkin model, Fig. 3(h) with a classical polynomial FEM model. The forward model computation takes 2 seconds with a wavelet-Galerkin method, and 140 seconds with FEM.

4.2 Simulation with an irregular geometry

The studied domain is included in the same rectangle as previously. 15 source positions are considered and the 7 detectors are positioned on the top boundaries (Fig. 3(b)). Two elliptic zones simulate the positions of the fluorophores. The resolution level is the same than previously (j=-3).

Figure 3(f) gives reconstructions with a wavelet-Galerkin model, Fig. 3(i) with a classical polynomial FEM model. The forward model computation takes 6 seconds with a wavelet-Galerkin method, and 70 seconds with FEM.

4.3 Simulation with a cube

The cube size is 3×3×3 cm³. 121 source positions are considered at a distance $\genfrac{}{}{0.1ex}{}{1}{{\mu }_s^{\prime }}$ from the bottom surface and 121 detectors are positioned on the top surface.

One ellipsoidal zone simulates the fluorophores. The level is in this case j=-2, so we have a coarser mesh than in the 2D cases (4913 nodes).

Figure 3(g) gives reconstruction by ART with a wavelet-Galerkin model, Fig. 3(j) with a classical FEM model. The forward model computation takes 25 seconds with a wavelet-Galerkin method, and 15 minutes with FEM.

5. Discussion

The simulations presented in the previous section have shown the effectiveness of the wavelet-Galerkin method and its relevance in terms of computation time (a factor of 10 to 70 is gained depending on the studied geometry). The final results obained by a wavelet-Galerkin modelling are comparable in terms of quality with the one obtained by a standard FEM modelling. Consequently we are able to solve the forward problem with the wavelet-Galerkin method, from ~10 to 70 times faster.

We observe in the first simulation that the right hand side fluorophore is reconstructed with a smaller intensity than the left one (in the 2 modelling 2D cases). That is due to the fluorophore horizontal orientation, since there is a lower resolution in the direction perpendicular to the source-detector axis.

An important remark is that if we want more reproducibility in the reconstruction by increasing the number of detectors or the number of source positions, the standard FEM will take much more time on account of the meshing step, for each source position; but in the case of the meshless modelling, this step does not exist, and the calculation time will be the same whatever the number of detectors and the number of source positions.

We see that the fictitious domain and the calculation of connection coefficients do not add significant computational time, and they allow the boundary conditions to be well verified, as well as with the FEM. In wavelet-Galerkin case, the accuracy of the boundary conditions is furthermore improved thanks to the calculation of the connections coefficients on intervals. For example, the Fig. 4 illustrates the differences in the reconstruction of the second simulation, between the exact calculation of connection coefficients on ∂Ω, and an approximate calculation which uses only the global connection coefficients defined on unbounded domain: artifacts on the edges are slightly reduced.

 

Fig. 4. Reconstruction of the second simulation (a) with connection coefficients accurately calculated on ∂Ω; (b) without accurate connection coefficients on ∂Ω.

Download Full Size | PDF

Finally, to illustrate the reconstruction error, the quadratic normalized error is mesured:

where u^obs_m is the observed fluorescence DPDW (given by Eq. (3)) simulated by COMSOL Multiphysics, and u^rec_m the reconstructed fluorescence DPDW, given by u^rec_m(r _d,r _s)=∫Ωg_m(r,r _d)βr^ec (r)u_x (r,r _s)d r, where β^rec is the reconstructed fluorescence parameter. In the case of wavelet-Galerkin model 3D, ε=1.68 %; in the case of FEM model, ε=5.98 %. So we can see that the wavelet-Galerkin method reduces time computation in the forward model, and the quality of the reconstruction is improved compared to the classical FEM-based forward problem solving.

6. Conclusion

In this work we have, as extensively as possible, described a method for solving the diffusion equation appearing in the forward problem in fluorescence diffuse optical tomography: the wavelet-Galerkin method. We have chosen the Daubechies scaling function of order 3, and we have calculated the connection coefficients to verify accurately the boundary conditions. Some simulations in 2D and 3D have validated the accuracy of the method, and have given faster forward problem solving than the classical FEM.

Some improvements and works need to be done: first experiments in 3D with a complex geometry would corroborate the advantages of the wavelet-Galerkin method. Then the multiresolution properties of the wavelets would be applicable to a multi-scale solving system: that would allow a better accuracy in the ROIs of the studied domain, from the forward problem and not uniquely in the reconstruction.

Note that the wavelet-Galerkin method can be applied to any inhomogeneous medium: suppose that the physical parameters are not constant in the medium. For example, $D={\left(D_1\cdots D_{N_x}\right)}^T$ and ${\mu }_a={\left({\mu }_{a1}\cdots {\mu }_{{\mathrm{aN}}_x}\right)}^T$, where N_x is the grid nodes number. Eq. (17) becomes

$\zeta \sum _{k=0}^{N_x-p}{u}_{j,k}{\int }_{\partial \Omega }{\Phi }_{j,k}{\Phi }_{j,l}+\sum _{k=0}^{N_x-p}{u}_{j,k}{\int }_{\Omega }{D}_k\nabla {\Phi }_{j,k}\nabla {\Phi }_{j,l}+\sum _{k=0}^{N_x-p}{u}_{j,k}{\int }_{\Omega }{\mu }_{\mathrm{ak}}{\Phi }_{j,k}{\Phi }_{j,l}$ $=\sum _{k=0}^{N_x-p}{q}_{j,k}{\int }_{\Omega }{\Phi }_{j,k}{\Phi }_{j,l},l=1,\dots ,N_x.$

We suppose for example that µ_ak is piecewise constant; it is then decomposed on a wavelet basis, for example the appropriate Haar’s wavelet: µ_ak=∑_m µ_a(k,m)Θ_j,m, where Θ_j,m is the Haar scaling function. Then we can write ∫µ_akΦ_j,kΦ_j,l=∑m µ_a(k,m)∫ΩΘ_j,mΦ_j,kΦ_j,l. Three wavelets product integrals have to be calculated. S. Dumont [28] for example describes the method to calculate these 3-term connection coefficients. As for the connection coefficients described in Eq. (18), the values are defined in tables.

The method is general and can be applied to a time-resolved diffusion equation or another time-resolved equation. It is very interesting in terms of computation time, because in the case for example of deformable mesh in time, the classical FEM would take a very long time: a wavelet-Galerkin method would be very fast thanks to its meshless property, while keeping the accuracy of the results.

Appendix : Connection coefficients

Table A-1 provides the values of the scaling filter h(n), for n=0 to 2M-1, in the case of Daubechies wavelet of order M=3. Tables A-2 and A-3 present the values of the connection coefficients useful in this work, that means Γ^1,1 _k,l[-∞;+∞], for different k and l, where ${\Gamma }_{k,l}^{m_1{m}_2}[-\infty ;+\infty ]$ is defined by Eq. (18), and Γ^0,0 _k,l[0;1] defined for a bounded interval [0,1], where Γ^0,0 _k,l[0;n] is defined by Eq. (20).

Table A-1. Daubechies scaling filter coefficients, order M=3.

View Table | View all tables in this article

Table A-2. Values of Γ^1,1_k,l[-∞;+∞]=∫+∞-∞Φ⁽¹⁾(x-k) Φ⁽¹⁾(x-l)dx, with Φ the Daubechies scaling function of order M=3. This filter is symmetrical, only the value |k-l| is considered.

View Table | View all tables in this article

Table A-3. Values of Γ^0,0_k,l[0;1]=∫¹₀Φ(x-k)Φ(x-l)dx, with Φ the Daubechies scaling function of order M=3.

View Table | View all tables in this article

The details of Eq. (21) are given by the following equation: ${\Gamma }_{k,l}^{0,0}[0;n]={\int }_0^n\Phi (x-k)\Phi (x-l)\mathrm{dx}$ $={\int }_0^1\Phi (x-k)\Phi (x-l)\mathrm{dx}+{\int }_1^2\Phi (x-k)\Phi (x-l)\mathrm{dx}+\dots +{\int }_{n-1}^n\Phi (x-k)\Phi (x-l)\mathrm{dx}$ $={\int }_0^1\Phi (x-k)\Phi (x-l)\mathrm{dx}+{\int }_0^1\Phi (x+1-k)\Phi (x+1-l)\mathrm{dx}+\dots +{\int }_0^1\Phi (x+n-1-k)\Phi (x+n-1-l)\mathrm{dx}$ $={\Gamma }_{k,l}^{0,0}[0;1]+{\Gamma }_{k-1,l-1}^{0,0}[0;1]+\dots +{\Gamma }_{k-n+1,l-n+1}^{0,0}[0;1]$ $=\sum _{p=0}^{n-1}{\Gamma }_{k-p,l-p}^{0,0}[0;1].$

References and links

1. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nature Med 9(1), 123–128 (2003).

2. S. R. Arridge, M. Cope, and D. T. Delpy, “The Theoretical Basis For The Determination Of Optical Pathlengths In Tissue - Temporal And Frequency-Analysis,” Phys. Med. Bio. 37(7), 1531–1560 (1992). [CrossRef]  

3. D. S. Burnett, Finite Element Concept (Addison-Wesley Publishing Company, 1987).

4. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15(2), R41–R93 (1999). [CrossRef]  

5. G. I. Marchuk, V. I. Agoshkov, and V. P. Shutyaev, “Adjoint equations and perturbation algorithms in nonlinear problems,” CRC Press, Boca Raton, FL (1996).

6. S. Goedecker, Wavelets and their application for the Partial Differential Equations in Physics, PPUR (1998).

7. O. V. Vasilyev, S. Paolucci, and M. Sen, “A Multilevel Wavelet Collocation Method For Solving Partial-Differential Equations In A Finite Domain,” J.Comp. Phys. 120(1), 33–47 (1995). [CrossRef]  

8. S. Jaffard, “Wavelet Methods For Fast Resolution Of Elliptic Problems,” SIAM J. Num. Anal. 29(4), 965–986 (1992). [CrossRef]  

9. G. Oberschmidt, G. Schneider, and A. F. Jacob, “A priori size estimation of wavelet-based Galerkin matrices,” Proceeding of 25th International Conference of Infrared and Millimeter Waves, 241–242, 2000.

10. S. Mallat, A wavelet tour of signal processing (Academic Press, 1999).

11. S. Mallat, “A Theory of Multiresolution Signal Decomposition: the Wavelet Representation,” IEEE Trans. Pattern Anal. Machine Intell. 11(7), 674–693 (1989). [CrossRef]  

12. B. Sahiner and A. E. Yagle, “Image-Reconstruction From Projections Under Wavelet Constraints,” IEEE Trans. SignalProc. 41(12), 3579–3584 (1993). [CrossRef]  

13. S. Bonnet, F. Peyrin, F. Turjman, and R. Prost, “Multiresolution Reconstruction in Fan-Beam Tomography,” IEEE Trans. Image Proc. 11(3), 169–176 (2002). [CrossRef]  

14. S. K. Nath, R. M. Vasu, and M. Pandit, “Wavelet based compression and denoising of optical tomography data,” Opt. Commun. 167(1–6), 37–46 (1999). [CrossRef]  

15. W. Zhu, Y. Wang, Y. N. Deng, Y. Q. Yao, and R. L. Barbour, “A wavelet-based multiresolution regularized least squares reconstruction approach for optical tomography,” IEEE Trans. Med. Imaging 16(2), 210–217 (1997). [CrossRef]  

16. B. Kanmani, P. Bansal, and R. M. Vasu, “Computationally efficient optical tomographic reconstruction through waveletizing the normalized quadratic perturbation equation,” Proc. SPIE 6164, R1–R10 (2006).

17. B. Kanmani, “Wavelet-Galerkin solution to the diffusion equation of diffuse optical tomography,” Internet Report, Saratov Fall Meeting: Optical Technologies in Biophysics ans Medicine IX, Russia (2007).

18. J. C. Baritaux, S. C. Sekhar, and M. Unser, “A Spline-based forward model for Optical Diffuse Tomography,” Biomedical Imaging, ISBI, 384–387 (2008).

19. D. Georges, “A fast Method for the solution of some Tomography Problems,” Decision and Control, CDC IEEE Conf. (2008).

20. I. Daubechies, “Ten lectures on Wavelets,” SIAM, 167–213 (1992).

21. J. Baccou and J. Liandrat, “On coupling wavelets with fictitious domain approaches,” Appl. Math. Lett. 18(12), 1325–1331 (2005). [CrossRef]  

22. F. Fedele, J. P. Laible, and M. J. Eppstein, “Coupled complex adjoint sensitivities for frequency-domain fluorescence tomography: theory and vectorized implementation,” J. Comp. Phys. 187(2), 597–619 (2003). [CrossRef]  

23. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, and M. S. McAdams, “Boundary-Conditions For The Diffusion Equation In Radiative-Transfer,” J. Opt. Soc. Am. A 11(10), 2727–2741 (1994). [CrossRef]  

24. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004). [CrossRef]  

25. http://www.comsol.com/

26. A. C. Kak and M. Slaney, Principles of computerized tomographic imaging (NY, IEEE Press, 1987).

27. A. Grossmann and J. Morlet, “Decomposition Of Hardy Functions Into Square Integrable Wavelets Of Constant Shape,” SIAM J. Math. Anal. 15(4), 723–736 (1984). [CrossRef]  

28. S. Dumont and F. Lebon, “Representation of plane elastostatics operators in Daubechies wavelets,” Comp. Structures 60(4), 561–569 (1996). [CrossRef]  

29. A. Ern and J. Guermond, Eléments Finis: théorie, applications, mise en oeuvre, Chap. 4 (Springer, 2002).

30. G. Beylkin, “On The Representation Of Operators In Bases Of Compactly Supported Wavelets,” SIAM J. Numerical Anal. 29(6), 1716–1740 (1992). [CrossRef]  

31. M. Q. Chen, C. I. Hwang, and Y. P. Shih, “The computation of wavelet-Galerkin approximation on a bounded interval,” Int. J. Num. Meth. Engin. 39(17), 2921–2944 (1996). [CrossRef]  

32. P. Charton and V. Perrier, “Produits rapides matrice-vecteur en bases d’ondelettes: application à la résolution numérique d’équations aux dérivées partielles,” Math. Modelling Numerical Anal. 29(6), 701–747 (1995).

33. R. Glowinski, T. W. Pan, R. O. Wells, and X. D. Zhou, “Wavelet and finite element solutions for the Neumann problem using fictitious domains,” J. Comput. Phys. 126(1), 40–51 (1996). [CrossRef]  

34. T. Vo-Dinh, Biomedical Photonics Handbook (CRC Press, 2003). [CrossRef]  

35. R. Elaloufi, R. Carminati, and J. J. Greffet, “Definition of the diffusion coefficient in scattering and absorbing media,” J. Opt. Soc. Am. A 20(4), 678–685, (2003). [CrossRef]