Discretization Error Analysis and Adaptive Meshing ... - IEEE Xplore

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

217

Discretization Error Analysis and Adaptive Meshing Algorithms for Fluorescence Diffuse Optical Tomography: Part I Murat Guven, Laurel Reilly-Raska, Lu Zhou, and Birsen Yazıcı*, Senior Member, IEEE

Abstract—For imaging problems in which numerical solutions need to be computed for both the inverse and the underlying forward problems, discretization can be a major factor that determines the accuracy of imaging. In this work, we analyze the effect of discretization on the accuracy of fluorescence diffuse optical tomography. We model the forward problem by a pair of diffusion equations at the excitation and emission wavelengths and consider a finite element discretization method for the numerical solution of the forward problem. For the inverse problem, we use an optimization framework which allows incorporation of a priori information in the form of zeroth- and first-order Tikhonov regularization terms. Next, we convert the inverse problem into a variational problem and use Galerkin projection to discretize the inverse problem. Following the discretization, we analyze the error in reconstructed images due to the discretization of the forward and inverse problems and present two theorems which point out the factors that may lead to high error such as the mutual dependence of the forward and inverse problems, the number of sources and detectors, their configuration and their positions with respect to fluorophore concentration, and the formulation of the inverse problem. Finally, we demonstrate the results and implications of our error analysis by numerical experiments. In the second part of the paper, we apply our results to design novel adaptive discretization algorithms. Index Terms—Adaptive meshing algorithms, error analysis, fluorescence diffuse optical tomography.

I. INTRODUCTION

T

HERE HAS been a growing interest in fluorescence imaging due to its potential for the characterization of biological processes at cellular and molecular levels [1]–[4].

Manuscript received June 07, 2009; revised July 31, 2009; accepted August 25, 2009. Current version published February 03, 2010. This work was supported in part by U.S. Army Medical Research under Grant W81XWH-04-10559 and in part by the Center for Subsurface Sensing and Imaging Systems under the Engineering Research Centers Program of the National Science Foundation under Award EEC-9986821). Asterisk indicates corresponding author. M. Guven is with Intel Corporation, Santa Clara, CA 95054 USA (e-mail: [email protected]). L. Reilly-Raska is with Lincoln Laboratory, Massachusetts Institute of Technology, Lexington, MA 02421 USA (e-mail: [email protected]). L. Zhou is with the Departments of Electrical, Computer and Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180 USA (e-mail: [email protected]). *B. Yazıcı is with the Departments of Electrical, Computer and Systems Engineering and Biomedical Engineering, Rensselaer Polytechnic Institute, Troy, NY 02180 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMI.2009.2031492

In particular, fluorescence diffuse optical tomography (FDOT) offers the quantification, 3-D imaging, and depth retrieval of fluorescence activity with high sensitivity, which can be used for functional and molecular characterization of normal and diseased tissues [5]–[7]. Like its analogue diffuse optical tomography, (see [8] and [9], and the references therein), FDOT poses a computationally intense imaging problem. This stems from the necessity of solving a forward problem comprised of a pair of coupled partial differential equations (PDEs) and an ill-posed integral equation resulting from the linearization of the inverse problem whose formulation relies on the solutions of the forward problem. In general, no closed form solutions can be computed except for specific domain geometries and optical medium properties [10], [11]. Hence, one has to consider numerical solutions of the forward and inverse problems. However, numerical solutions are merely approximations to the actual solutions and they possess error as a result of the discretization involved in the process of solving them. Due to the interdependence of the forward and inverse problems, discretization of each problem typically leads to error in the reconstructed images. While the effect of inverse problem discretization on the imaging accuracy can be deduced rather intuitively [12], errors due to the discretization error in the forward problem solutions can also result in severe imaging artifacts in optical tomography [12], [13]. There is a vast degree of work on the estimation and analysis of discretization error in the solutions of partial differential equations (PDEs) [14]–[19]. A different approach is followed in [20], [21] in which error in quantities of interest is related to the discretization of the second-order elliptic PDEs. In the area of parameter estimation problems governed by PDEs, relatively little has been published. See for example [22], [23] for an a posteriori error estimate for the Lagrangian in the inverse scattering problem for the time-dependent acoustic wave equation and [24] for a similar approach, and [25] for a posteriori error estimates for distributed elliptic optimal control problems. In the area of DOT, it was numerically shown that the approximation errors resulting from the discretization of the forward problem can lead to significant errors in the reconstructed optical images [13]. In [26], we presented an approach for the analysis of the error in reconstructed optical absorption images due to discretization, which led to the development of new adaptive mesh generation techniques [12]. In these studies [12], [26], the ill-posed nature of the inverse problem was addressed by zeroth-order Tikhonov regularization [27] and a collocation method was used for the discretization of the inverse problem.

0278-0062/$26.00 © 2010 IEEE

218

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

The main premise in the first part of this work is to analyze the error in fluorescence optical imaging due to discretization. In this respect, we identify the key factors specific to the imaging problem that show how discretization of the forward and inverse problems impacts the accuracy of the reconstructed optical image. In particular, we focus on the estimation of the fluorophore absorption coefficient in a bounded optical domain where the light propagation at the excitation and emission wavelengths is modeled by a pair of coupled frequency-domain diffusion equations with appropriate boundary conditions. We consider an iterative linearization algorithm based on first-order Fréchet derivatives to address the nonlinearity of the inverse problem, and use the gradient descent algorithm to solve the resulting optimization problem. Following each linearization step, we formulate the inverse problem in the optimization framework which enables us to incorporate a priori information in the form of zeroth- and first-order Tikhonov regularization terms. For the discretization of the forward problem, we use finite elements with first-order Lagrange basis functions. Before we discretize the inverse problem, we first convert the optimization problem into a variational formulation which can be regarded as a boundary value problem with the assignment of suitable boundary conditions. Then we use projection by Galerkin method with first-order Lagrange basis functions to discretize the resulting inverse problem. Following the discretization of the forward and inverse problems, we derive two new error estimates which show respectively the effect of forward and inverse problem discretizations on the accuracy of reconstructed optical absorption coefficient of the fluorophore. Next we discuss the major implications established by these estimates. The error analysis presented in this work motivates the development of novel adaptive discretization schemes. In the second part of this work, we propose two novel adaptive mesh generation algorithms for the discretization of the forward and inverse problems [28], and discuss and validate the potential computational savings and improvements in the accuracy of optical fluorescence imaging. The outline of this paper is as follows. We describe the forward model in the next section. In Section III, we state the inverse problem formulation under iterative linearization and discuss regularization to address ill-posedness. We next pose the inverse problem as an optimization problem incorporating a priori information and derive the variational formulation from this optimization problem. Section IV details the discretization of forward and inverse problems, respectively. In Section V, we present the theorems which describe the error in fluorescence imaging arising from these discretizations. In Section VI, we demonstrate the dependency of the error to a number of specific factors in numerical experiments. Finally, in Section VII, we conclude our discussion.

TABLE I DEFINITION OF FUNCTION SPACES AND NORMS

to denote vectorized quantities such as , . Table I provides a summary of key variables and function spaces and norms used throughout the paper. B. Forward Problem Derivation We start with the coupled diffusion equations which describe the light transport in a fluorescent medium of a bounded domain with Lipschitz boundary

(1)

(2) , is the source operating frequency, where denote the excitation and emission wavelengths, subscripts represents the optical fields, represents the isotropic is the th excitation source. We diffusion coefficients, and assume that the diffusion coefficients are known and they are identical at both the excitation and emission wavelengths in , the closed domain; this implies . The quantum efficiency is denoted by ; is the absorption coefficient of the fluorophore and is the lifetime of the fluorophore. For the sake of exposition, we make . the following simplifying assumption that the frequency Subsequent developments can be extended to include multiple and repfrequencies where is known. The quantities resent the absorption coefficient of the medium at the excitation and emission wavelengths, respectively. Typically these are represented as (3) (4)

II. FORWARD PROBLEM A. Notational Conventions In this paper, we denote operators by capital cursive Latin and matrices by bold capital Latin letters . Funcletters tions are denoted by lowercase Latin or Greek letters. For a function , denotes its finite element approximation. We use bold

where the subscript denotes endogenous properties and denotes exogenous properties. Without loss of generality, we asand are nonnegative and bounded sume that both on . be the number of point sources at position for Let . Based on the assumptions stated above, we use

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

the following boundary value problem to model NIR light propagation at the excitation wavelength due to the th source at

where the

219

th measurement satisfies the following model: (14)

(5) (6) where

. The Robin-type boundary conditions are (7) (8)

, is a parameter governing the internal reflection where , and denotes the directional derivative at the boundary along the unit normal vector on the boundary. In this work, represents the th point source which is modeled by a Gaussian function centered at source position [26]. In order to simplify the analysis of later sections, we make use be of the adjoint problem associated with (6) and (8). Let the number of detectors. Then, for a detector located at ,

Our objective is to recover the quantity using the measurement vector based on the nonlinear integral equation (14) for th pair. In this model, we assume that the measureeach ments are noise-free. This assumption allows us to eliminate the effect of noise in our error analysis, and to focus primarily on the effect of the discretization error. In the next subsection, to address the problem of nonlinearity in (14), we select an iterative linearization scheme based on first-order Fréchet derivatives. Next, to address the ill-posedness, we discuss regularization in an optimization framework and incorporation of a priori information about the unknown image . Then by taking the derivative of the resulting optimization problem and defining appropriate boundary conditions for it, we convert this optimization problem into a boundary value problem. In the final subsection, we show the variational formulation of the boundary value problem and comment on the existence and uniqueness of the solution. A. Iterative Linearization

(9)

Consider an infinitesimal perturbation on

, [3], [29]

(10) (15) where is the adjoint source. For a point adjoint source located at the detector position , the following holds [8]: (11)

Then at each linearization step, the corresponding perturbation in the emission field at detector position due to the th source at is given by the following linear integral equation:

is the Green’s function of (6) and (8). Note that in this where paper, we model the point adjoint source by a Gaussian function with sufficiently low variance, centered at [26]. The emission field at due to the source at is given by the following nonlinear integral equation: (12)

(16)

and defined in (12) is nonThe relationship between is dependent on which in turn is related linear because to and the dependence of on is clear. In the next section, we formally state the inverse problem and address the nonlinearity by using an iterative linearization scheme based on first-order Fréchet derivatives.

where is the solution to the boundary value problem (9), (10) is the solution to the boundary value problem (5) and and in (5) is replaced by . In (16), (7) where the first integral results from the right-hand side of (6), while the second and third integrals originate from the dependence of and on the unknown fluorophore absorprespectively tion coefficient. We note that the kernels of the second and third integrals are much smaller than the kernel of the first integral. Therefore, the first integral in (16) dominates and the last two terms can be neglected. As a result, inverse problem at each linearization step reads

III. INVERSE PROBLEM Given sources and detectors, we define to be the , due measurement for a detector at position , to a source at , . The individual measurements can be grouped into the vector form

(17) (13)

220

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

To simplify the notation, we introduce which represents the unknown perturbed fluorophore absorption coefficient scaled by the quantum efficiency and assume is known and . Furthermore, noting that that the emission and excitation subscripts are fixed for the duand ration of this analysis, we represent , suppressing the dependence of these functions. to be the differential measurement at the th We define detector due to th source normalized with respect to a known . background fluorophore absorption. Let Using (17) we model as follows:

In the next subsection, we address the ill-posedness in an optimization framework by incorporating regularization terms. B. Inverse Problem as an Optimization Problem and Regularization In this subsection, we address the ill-posedness of (24) through regularization in the optimization framework which provides a suitable means for the incorporation of a priori information about the solution. In this respect, we consider the following minimization problem where we seek a solution (25)

(18) We represent the differential measurements corresponding to individual source-detector pairs as elements of a vector

where the smoothness on the solution is imposed through the use of appropriate regularization terms. The funcand tional in (25) can be decomposed into two parts, as follows: (26) where measures the difference between the predicted and actual measurements

Then (19)

(27)

is a vector of operators whose where th entry acting on corresponds to (18). Although all norms on a finite-dimensional space are equivalent, we select the norm on the range of to be the norm as this proves useful in later analysis. Then an upper bound for the linear operator can be given by

and the regularization term contains the a priori information. In this work, we assume that a priori information on the image denote the a priori inand image gradient is available. Let formation on the image and denote the a priori information on the image gradient. We incorporate the a priori information via zeroth- and first-order Tikhonov regularization terms as follows [31]:

(20) The boundedness and the finite-dimensional range of operator means it is compact [30]. We define the adjoint operator acting on the vector as (21) for

where

and

(28) where is the image gradient and are regularization parameters. Using (27) and (28), the minimization problem (25) can be rewritten as follows:

. Let

, we have (22)

where (23) Then an alternate form of (19) can be expressed as follows: (24) . Note that where Therefore, (24) is ill-posed.

is compact since

is compact.

(29) There are a number of methods in choosing appropriate regularization parameters, see [32]–[36]. In this work, we assume and are properly chosen and focus on deriving disthat cretization error estimates. In the next subsection, after defining appropriate boundary conditions, we consider the equivalent variational formulation of the minimization problem in (29).

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

C. Inverse Problem as Boundary Value Problem and Variational Formulation

A more convenient way to express (34) is through a bilinear form. Thus, we define

In this work, we follow a finite element method for the discretization of the inverse problem. Due to incorporation of the regularization term on the gradient of the solution, a natural step is to formulate the minimization problem as a variational one. In this subsection, we describe the derivation of the variational problem formulation of the inverse problem by first considering the first-order optimality condition for the minimization problem (29). Next, with the aid of properly chosen boundary conditions, we transform the optimization problem into a boundary value problem (BVP), which is followed by the variational formulation of the BVP. Finally, we show that a unique solution exists to the variational formulation of the regularized inverse problem. , where The solution of (29) satisfies is the gradient with respect to the th direction for , 2, 3. In particular, if , the Gâteaux derivative [37] is defined by (30) Taking the Gâteaux derivative of (29) with respect to setting it equal to zero yields

221

and (31)

where (32) Note that is composed of known terms from a priori information and measurements. We consider (31) with the following Neumann boundary condition: (33) where is the directional derivative of along the unit . The boundary condition in normal vector at the boundary (33) implies that no changes in the perturbed fluorophore concentration occur across the boundary. At this point, one can consider a finite difference scheme for the solution of the inverse problem which is posed as a boundary value problem (31)(33). However, as our goal in this paper is to apply a finite element scheme for the discretization of the BVP, we obtain the corresponding variational (weak) problem. Hence, , we multiply both sides of (31) by a test function and integrate over . Applying Green’s first theorem to the last term on the left and using the boundary condition in (33), we obtain

(34)

(35) (36) where the inner product is defined by

Hence, (34) can be expressed as (37) It can be shown that the bilinear form (35) is bounded and co[30]. Thus, by ercive for regularization parameters the Lax–Milgram lemma, a unique solution exists for the reguused in larized inverse problem (37) for each pair of the formulation of the Galerkin problem (34) [30], [38]. For an explicit statement of the Lax–Milgram lemma, see Appendix A. In the following section, we describe the discretization methods selected in this paper for each of the separate forward and inverse discretizations as well as the combined forward and inverse discretization. IV. DISCRETIZATION BY FINITE ELEMENT METHOD OF THE FORWARD AND INVERSE PROBLEMS In the following subsections, we first discuss the variational formulation and finite element discretization of the forward problem. In practice, for arbitrary domains and background optical properties, no analytical solutions exist for the forward problem when it is defined in a variational form. Thus, we discretize the forward problem and obtain finite-dimenand , for , sional approximations of . Next, we use the finite element solutions of the forward problem in the inverse problem formulation, which implies an approximation to the inverse problem. The resulting inverse problem in general does not possess a closed form solution. Therefore, finding the solution calls for numerical techniques. We discuss the discretization of the resulting approximate inverse problem using projection by Galerkin method. A. Discretization of the Forward Problem In this subsection, we discuss the forward problem discretization. We express the coupled PDEs in their variational form in order to apply a finite element method. , To do so, we multiply (5) by a test function and apply Green’s theorem to the second derivative term. Then using the boundary condition in (7) we have

(38) It can be shown that a unique solution for (38) exists and is , bounded [38], [39]. Similarly, for a test function

222

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

the variational form for the adjoint forward problem (9), (10) becomes

where

(39) for which it is possible to show that a unique, bounded solution exists as well. denote the piecewise linear Lagrange basis funcLet , , as the finite tions. We define , . dimensional subspace spanned by Note that are associated with the set of points , , on . Similarly, we define , , as the finite-dimensional subspace spanned by , , which are associated with the set of points , . in (38) and in (39) are reNext, the functions placed by their finite-dimensional counterparts

is a positive constant, and denote the and norms, respectively on ; is the diameter of the smallest ball containing the and in the solution . finite element In the next subsection, these approximate solutions to the forward problem are substituted into the inverse problem operator. The error is estimated based on the resulting operators with approximations. B. Simultaneous Discretization of the Forward and Inverse Problems

We substitute the finite element solutions , and , of the forward problem into in the operators defined by (18) and (21). The resulting and , indiapproximate operators are denoted by tildes cating that the finite element solutions of the forward problem are used. By so doing, we arrive at the approximate variational problem formulation (46)

(40) In (46), (41)

and

are given, respectively, by (47) (48)

The representation is an approximation to the function for each source (detector). This means that for each source and detector, the dimension of the solution can be difcan vary for each and , respecferent; the parameters tively. The finite-dimensional expansions are therefore depenas represented by the superscript. dent on the parameters However, we suppress this cumbersome notation as the dependence is clearly understood. Substitution of (40), (41) into the variational forward problem (38), (39) yields the matrix equations

where

and

(42) (43) and . Here and are the finite element matrices and and are the load vectors resulting from the finite element discretization of the forward problem. denote the set of linear elements used to discretize Let , where is the number of elements (38) for for all . for the th source such that denote the set of linear elements used to Similarly, let where is the number of discretize (39) for for all elements for the th detector such that . A bound for error and between and respectively on each finite element can be given by [38] for

(44) (45)

Next, we discretize the functions and . Let denote a sequence of finite-dimensional subspaces of dimension , spanned by the first-order Lagrange basis functions which are associated with the set of points , , on . We replace and in (37) by and their respective, finite-dimensional counterparts

(49)

(50)

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

where and are unknown coefficients. As it is clear that the finite-dimensional expansions are dependent on the parameter , this dependence is hereafter suppressed. Substituting (49), (50) into (46) we arrive at (51) This can be transformed to a matrix equation (52) where represents the unknown coeffiand are recients in the finite expansion of (50) and spectively the finite element matrix and the load vector resulting from the projection of (46) by Galerkin method. denote the set of linear elements used to discretize Let , where is the number of elements (46) for such that . Note that the inverse problem mesh is independent to the meshes and , which are used to discretize the forward problem. Similar to the forward problem, a traditional error estimate for error between and on each finite element can be given by

223

In the following, we analyze each of the error contributors and norm derive estimates in the form of upper bounds of the of these errors. A. Error in Fluorescence Imaging Due to Forward Problem Discretization of the The following theorem presents a bound for the error , where satisfies the exact inverse problem (37) and satisfies the approximate inverse problem (46). denote the set of linear elements used Theorem 1: Let to discretize (38) for ; such that and is the diameter of the smallest ball that contains the th . Similarly, let element in the solution , for all denote the set of linear elements used to discretize (39) ; such that and is the for diameter of the smallest ball that contains the th element in the solution , for all . Then a bound for the error between the solution of (37) and the solution of (46) due to the approximations and is given by

(53) where is a positive constant, and denote the and norms, respectively on ; and is the diameter in the of the smallest ball containing the finite element . solution V. ANALYSIS OF THE ERROR IN FLUORESCENCE IMAGING DUE TO DISCRETIZATION In this work, we consider the solution of the problem stated in (37) to be exact solution since neither the forward problem nor the inverse problem is discretized. Our objective is to examine the error in fluorescence absorption imaging due to the finite element discretization of the forward and inverse problems. Then the error analysis can be used in the design of adaptive meshes that could reduce the total error in the reconstructed images due to discretization. We have divided this section into two subsections: In the first subsection, we derive an estimate for the error in fluorescence absorption imaging due to the forward problem discretization as described in the previous section. Thus, the first error we find is , where satisfies the exact inverse the difference problem (37) and satisfies the approximate inverse problem (46). Note that the inverse problem is not discretized for this case. In the second subsection, we analyze the error in the reconstructed fluorophore absorption coefficient resulting from the finite element discretization of the inverse problem. Therefore, of (46) and the we examine the error between the solution solution of (51), i.e., . In this case, the error is due entirely to the discretization of the approximate inverse problem (46). Finally, we define the total error as the difference between of (37) and of (51) in terms of the two the exact solution contributors (54)

(55) are regularization parameters. where Proof: The proof of this theorem is given in Appendix D. The error estimate (55) shows the specific effect that the forward problem discretization has on the accuracy of the inverse problem solution. In this respect, for the forward problem, Theorem 1 suggests a discretization criterion for the forward problem which includes the accuracy of the inverse problem solution, in addition to the accuracy of the forward problem solution. First, the forward problem discretization should inand . Second, clude the discretization of each solution of the th element to keep the total error bound low, the has to be small when in solution is large on this element; and the of the th element in solution has to be small when is large on or is large on the elethis element. Note that ment close to the th source or the th detector, respectively, and the value of the terms and depend on the absorption coefficient of the fluorophore and its position with respect to the sources and detectors. Therefore, finer elements near the designated source or detector as well as near the

224

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

heterogeneity may result in a lower error bound. The traditional error estimates (44) and (45) only depend on the smoothness and and the finite-dimensional space of and support of approximating functions. Traditional error estimates require to or small where or is large, keep respectively. However, lower error bound in (44) and (45) only guarantees to reduce the error in the solution of the forward problem, and may not necessarily reduce the error in the solution of inverse problem. Furthermore, it is clear that other factors can lead to a higher error bound. The regularization parameters scale the sum of results in the terms. Thus, choosing smaller values for a higher error estimate. Note, too, that since the error bound is a sum over all sources and detectors, increasing the number of either can have an impact on the error estimate.

and , as well as the spathe forward problem solutions tial relationship among these solutions. This also implies that simply keeping the mesh parameter small over regions where is large, as indicated in (53), may not ensure a lower error bound, thereby a reduction of the error in the reconstructed image because of the dependence of the bound on the location of the heterogeneity with respect to the sources and detectors. Similar to Theorem 1, the regularization parameters and the number of sources and detectors also affect the error bound. Combining results of Theorems 1 and 2 and rearranging the terms, both error estimates can be viewed in a single bound. Let for all and . Then

B. Error in Fluorescence Imaging Due to Inverse Problem Discretization In the following theorem, we present a bound for the norm of the error between the solution of (46) and the solution of (51). Theorem 2: Consider the Galerkin projection of the variational problem (46) on a finite-dimensional subspace using a set of linear finite elements , , whose vertices are at , , such that , and let be the diameter of the smallest ball that contains the th element. Assume that the solution of (46) also satisfies . Then a bound for the error in the solution of (51) with respect to the solution of (46) can be given by (57) We note that the error estimates in Theorem 1 and 2 can be extended to show the effect of noise as an additional error source [26]. C. Error in Iterative Linearization (56)

are regularization parameters. where Proof: The proof of this theorem is given in Appendix E. The error estimate (56) shows that the error bound for the discretization of inverse problem not only depends on the inverse problem solution itself, but also on the solutions of the forward problem. The first term in the bracket of (56) shows that the term is scaled by the finite element solutions of the forward problem . This implies that the error bound is dependent on the location of the heterogeneity with respect to the sources and detectors. The second and third term in the bracket suggest that keeping the mesh size small where is large, can help to lower the error bound, but it also depends on the regularization parameters and . Comparing (56) with the traditional error estimate (53), (56) suggests a discretization criterion based on the inverse problem solution ,

The error bounds described in the previous subsections only address one linearized step of the nonlinear inverse problem. To consider all iterative linearization steps, we use the following method to describe the propagation of errors at each linearization step: and represent the solution to (37) and (51) Let at the th linearization step, respectively. The absorption coefth iteration step ficient as given in (15) at the end of the is , where is an initial guess for the value of the background absorption coefficient, and contains discretization error with respect to the exact so. Three errors are introduced at the th iteration step. lution The first is due to the finite element discretization of the inverse problem. The next is due to finite element discretization in the in operators forward problem by the error and the corresponding measurement vector error in . Finally, and , which are related to each other, appear as the coefficient in the forward problem (5)–(7) and the

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

225

Fig. 1. Setups used for the simulation studies 1 and 2. The squares and triangles denote the detectors and sources, respectively. (a) Optical domain and sourcedetector configuration for Simulation study 1. (b) Optical domain and source-detector configuration for Simulation study 2. The radius of the circles is 3 mm.

adjoint problem (9), (10), respectively. Then, the error in the soat each step will propagate and lead to an additional lution and at the th step in addition to the discretizaerror in tion error analyzed before. If the exact solution (without any discretization error) th iteration step is given by after , then the error in at the th iteration is bounded by

(58)

VI. NUMERICAL EXPERIMENTS To demonstrate the implications of Theorem 1 and 2, we performed a series of numerical experiments. In the first simulation, we considered a series of image reconstructions to show how the error, in reconstructed images due to discretization, changes as the absorption coefficient of the fluorophore increases. In the next set of experiments, we showed the effect of the relative position of the fluorophore concentration with respect to the sources and detectors on the accuracy of reconstructed images. We first describe the numerical phantom, data generation and image reconstruction procedures, and next, describe the results of numerical simulations. A. Description of the Phantom, Data Generation and Reconstruction For both simulation studies, we considered a cubic domain shown in Fig. 1(a). We placed 25 sources and 25 detectors evenly on two 5 5 grids at the bottom and top surface of the domain, respectively. The circular heterogeneity represents the fluorophore concentration with constant absorpand radius 3 mm, embedded in an optically tion coefficient at both exhomogeneous background with citation and emission wavelengths. We set the diffusion coeffifor and the refractive index cient for the boundary. Using the parammismatch parameter eters above, we simulated the fluorescence data by solving the coupled diffusion equations (5) and (6) with their corresponding

boundary conditions (7) and (8) on a fine uniform grid with nodes. To demonstrate the effect of discretization, we considered three image reconstruction scenarios corresponding to three different meshing schemes. 1) To obtain an accurate solution, we solved the forward problem on a very fine uniform mesh with nodes and obtained the corresponding finite element solutions, which we assume has negligible error due to discretization. Next, we formulated the inverse problem with these finite element solutions and discretized the resulting variational inverse problem on the same uniform mesh. As a result of using very fine meshes, we assume that the solution of the inverse problem on this mesh, denoted by , has negligible error that can be attributed as a baseline for to discretization. Therefore, we used comparison. 2) To demonstrate the effect of discretization of the forward problem, we solved the forward problem on relatively coarse mesh shown in Fig. 2(a). Then, we used the fine nodes to compute the inuniform mesh with and verse problem solution. We denote this solution by assume that it possesses error with respect to the baseline due to only the forward problem discretization. image 3) To demonstrate the effect of discretization of the inverse problem, we discretized the forward and inverse problem on two coarse meshes shown in Fig. 2(a) and (b), respectively, and computed the solution of the resulting inverse and assume that it problem. We denote the solution by possesses error with respect to the solution due to only the inverse problem discretization. We note that for each case the forward problem was solved with the same parameters as the data generation with the exception of the mesh size. For the inverse problem, we chose the regularization parameters as small as possible, yet large enough to enable a robust image reconstruction. In this respect, appropriate values for the regularization parameters were empirically selected as and . We assumed that no a priori information were available about the image or its gradient, hence and . we set Note that our simulation study is performed using a C++ finite element library—deal.II [40]. We used hexahedral finite elements associated with trilinear Lagrange basis functions to dis-

226

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

Fig. 2. Coarse uniform meshes used to discretize the forward and inverse problems in the Simulation studies 1 and 2. The mesh is cut through to show the mesh structure inside. (a) The coarse uniform mesh with 25 25 13 nodes used to discretize the forward problem. (b) Coarse uniform mesh with 17 17 9 nodes used to discretize the inverse problem.

2 2

H

( ) NORM

2 2

TABLE II

OF THE ERROR RESULTING FROM FORWARD AND INVERSE PROBLEM DISCRETIZATIONS IN THE RECONSTRUCTED OPTICAL IMAGES IN SIMULATION STUDY 1. ABSORPTION COEFFICIENT OF THE FLUOROPHORE IS GIVEN IN cm



cretize both the forward and inverse problems. Then we used the Gaussian quadrature method to evaluate the integrals in the variational problems given in (38), (39) and (51) [41]. While solving the forward (or inverse) problem, we evaluated the value of the inverse (forward) problem solution at the Gaussian quadrature points associated with the forward (inverse) problem mesh. B. Simulation Results We performed two sets of numerical experiments to show the effect of fluorophore concentration and fluorophore location on the error in reconstructed images due to discretization. In the first simulation study, we consider the geometry shown in Fig. 1(a). To show the effect of fluorophore concentration on the error in reconstructed images resulting form discretization, we chose 5 data sets corresponding to 5 different values . In the for second simulation study, we considered the geometry shown in Fig. 1(b). To show the effect of the fluorophore location on the error due to discretization in reconstructed images, we simulated 5 data sets corresponding to 5 different positions of the circle centered at: (0,0, 0.8), (0,0, 0.4), (0,0,0), (0,0,0.4) and (0,0,0.8), respectively. Note that the center of the square domain is positioned at (0,0,0). The absorption coefficient of the fluo. Tables II rophore concentration was set to and III show the error in the solutions and due to forward and inverse problem discretizations, respectively, in two sets of simulations.

H

TABLE III

( ) NORM OF THE ERROR RESULTING FROM FORWARD AND INVERSE

PROBLEM DISCRETIZATIONS IN THE RECONSTRUCTED OPTICAL IMAGES IN SIMULATION STUDY 2

The results in Tables II and III show the following. i) The error due to both the forward and inverse problem discretization is comparable to the norm of the accurate . This indicates that the error due to dissolution cretization can be a significant factor that determines the accuracy of fluorescence imaging. ii) The error in reconstructed images due to the forward problem discretization is larger as compared to the error due to the inverse problem discretization. This shows the effect of the accuracy of the forward problem solutions in reconstructed images. This result is consistent with Theorem 1 which states the interrelatedness of the forward and inverse problem discretization in terms of the error in reconstructed images. iii) The error due to discretization increases nonlinearly of the with the increasing absorption coefficient fluorophore, and also increases when the fluorophore concentration moves closer to the detectors. The results show the dependency of the discretization error on image and/or geometry specific factors such as the value of fluorophore absorption coefficient and its position. The effect of image and geometry specific factors is not indicated by the traditional error estimate, but clearly spelled out by the error estimates in Theorems 1 and 2. In the second part of this work [28], we present two novel adaptive mesh generation algorithms developed based on the theorems presented in this paper.

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

VII. CONCLUSION In this work, we analyzed the error in fluorescence diffuse optical tomography resulting from the discretization. We presented the results of the error analysis in two theorems which provide two estimates for the error in the reconstructed absorption coefficient of the fluorophore resulting from the discretization of the forward and inverse problems, respectively. These theorems show that the error in the reconstructed optical image due to the discretization of each problem is affected by the absorption coefficient value of the fluorophore, the number of source-detector pairs, their locations with respect to the fluorophore concentration and the regularization parameters that are used to address the ill-posedness of the inverse problem. To demonstrate the dependency of the error due to discretization on various medium parameters, we conducted a series of simulation experiments. The first set of experiments shows the dependence of the error due to discretization on the absorption coefficient of the fluorophore. The second set of experiments show that the location of the fluorophore concentration with respect to the sources and detectors can be an important factor that determines the extent of the error resulting from discretization in the reconstructed optical image. APPENDIX A LAX–MILGRAM LEMMA Given a Hilbert space , a continuous, coercive, bilinear and a continuous linear functional , there exists form a unique such that

APPENDIX D PROOF OF THEOREM 1: ERROR ESTIMATE DUE TO FORWARD PROBLEM DISCRETIZATION from both sides of (46) yields

Subtracting

(63) Adding and subtracting to

on the right hand side of (63) leads

(64) Then following the Lax–Milgram Lemma in Appendix A, the is bounded by error

(65) where Clearly,

is defined in Appendix B.

(66)

(59) The proof may be found in [38]. As a direct consequence of this lemma, an upper bound on may be established as

227

Note that the dual norm of the functional defined by

is

(60) where .

is the coercivity constant and

is the dual space of

APPENDIX B DEFINITION OF DUAL NORM The dual norm of

is defined by [38]

(67) (61)

where space of

Recall parameters

denotes the norm of . APPENDIX C COERCIVITY OF . Thus , we see

Hence, (68)

which is the dual Following [26], we express and from the positive

(62) (69) where

is the coercivity constant.

228

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 29, NO. 2, FEBRUARY 2010

where and errors. Then if we expand the integral on integral on the finite element , , an upper bound for

are the discretization as a summation of the and , becomes

APPENDIX E PROOF OF THEOREM 2: ERROR ESTIMATE DUE TO INVERSE PROBLEM DISCRETIZATION Taking (46) as the starting point, by coercivity we can write

Let

. Then the above inequality is equivalent to

since and the error is orthogonal to the finite-dimensional subspace with respect to the norm induced by the bilinear form. Noting

(75)

(70) it is clear Next, following a similar argument as in (67) and (68), we obtain (71) Using the definition of is given by [26]

in (32), an upper bound for

(76) Cancel

terms

(72) A bound for

can be obtained by using (18)

(77) (73)

, we expand Finally, to compute an upper bound for norm computed on as a summation on the finite the elements , and ,

(74) In the end, using the discretization error estimates (44) and (45) leads to the theorem.

be the interpolant of and be the Let interpolation error. Then the first term in the bound (77) can be expanded as follows:

GUVEN et al.: DISCRETIZATION ERROR ANALYSIS AND ADAPTIVE MESHING ALGORITHMS FOR FLUORESCENCE DIFFUSE OPTICAL TOMOGRAPHY

Then, if we expand the integral on integral on the finite element ,

as a summation of the

The remaining two terms in (77) can be expressed in a straightforward way

Assume that our solution also satisfies . Then a bound for the interpolation error and its gradient on each element can be given by [38] (78) (79) and are respectively where is a positive constant, norms on and is the diameter of the smallest the and ball containing the finite element . Finally, substituting (78) and (79) into (77) proves the theorem. REFERENCES [1] A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt., vol. 42, pp. 3081–3094, 2003. [2] A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Expr., vol. 12, no. 22, pp. 5402–5417, 2004. [3] M. J. Eppstein, F. Fedele, J. P. Laible, C. Zhang, A. Godavarty, and E. M. Sevick-Muraca, “A comparison of exact and approximate adjoint sensitivies in fluorescence tomography,” IEEE Trans. Med. Imag., vol. 22, no. 10, pp. 1215–1223, Oct. 2003. [4] R. Roy, A. Godavarty, and E. M. Sevick-Muraca, “Fluorescence-enhanced optical tomography using referenced measurements of heterogeneous media,” IEEE Trans. Med. Imag., vol. 22, no. 7, pp. 824–836, Jul. 2003. [5] A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: Study of the normalized born ratio,” IEEE Trans. Med. Imag., vol. 24, no. 10, pp. 1377–1386, Oct. 2005. [6] V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized born approximation,” Opt. Lett., vol. 26, no. 12, pp. 893–895, 2001. [7] V. Ntziachristos, G. Turner, J. Dunham, S. Windsor, A. Soubret, J. Ripoll, and H. A. Shih, “Planar fluorescence imaging using normalized data,” J. Biomed. Opt., vol. 10, no. 6, pp. 1–8, 2005. [8] S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl., vol. 15, pp. R41–R93, 1999. [9] A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., vol. 50, pp. R1–R43, 2005. [10] V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A, vol. 19, no. 3, pp. 1–9, 2002. [11] V. A. Markel, V. Mital, and J. C. Schotland, “Inverse problem in optical diffusion tomography. III. Inversion formulas and singular-value decomposition,” J. Opt. Soc. Am. A, vol. 20, no. 5, pp. 890–902, 2003. [12] M. Guven, B. Yazici, K. Kwon, E. Giladi, and X. Intes, “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: II,” Inv. Probl., vol. 23, pp. 1135–1160, 2007.

229

[13] S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl., vol. 22, pp. 175–195, 2006. [14] M. Ainsworth and J. T. Oden, “A unified approach to a posteriori error estimation using elemental residual methods,” Numerische Mathematik, vol. 65, pp. 23–50, 1993. [15] I. Babuˇska and W. C. Rheinboldt, “Error estimates for adaptive finite element computations,” SIAM J. Numer. Anal., vol. 15, pp. 736–754, 1978. [16] I. Babuˇska, O. C. Zienkiewicz, J. Gago, and E. R. de A. Oliveira, Accuracy Estimates and Adaptive Refinements in Finite Element Computations. New York: Wiley, 1986. [17] R. E. Bank and A. Weiser, “Some a posterior error estimators for elliptic partial differential equations,” Math. Comput., vol. 44, pp. 283–301, 1985. [18] T. Strouboulis and K. A. Hague, “Recent experiences with error estimation and adaptivity, Part I: Review of error estimators for scalar elliptic problems,” Comput. Meth. Appl. Mech. Eng., vol. 97, pp. 399–436, 1992. [19] R. Verfurth, A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques. New York: Teubner-Wiley, 1996. [20] J. T. Oden and S. Prudhomme, “Goal-oriented error estimation and adaptivity for the finite element method,” Comput. Math. Applicat., vol. 41, pp. 735–756, 2000. [21] V. Heuveline and R. Rannacher, “Duality-based adaptivity in the hp-finite element method,” J. Numer. Math., vol. 11, pp. 95–113, 2003. [22] L. Beilina and C. Johnson, “A posteriori error estimation in computational inverse scattering,” Math. Models Meth. Appl. Sci., vol. 15, pp. 23–37, 2005. [23] L. Beilina and C. Johnson, “Adaptive finite element/difference method for inverse elastic scattering waves,” Appl. Comput. Math., vol. 2, pp. 158–174, 2003. [24] W. Bangerth, “Adaptive finite element methods for the identification of distributed parameters in partial differential equations,” Ph.D. dissertation, Univ. Heidelberg, Heidelberg, Germany, 2002. [25] R. Li, W. Liu, H. Ma, and T. Tang, “Adaptive finite element approximation for distributed elliptic optimal control problems,” SIAM J. Contr. Optim., vol. 41, pp. 1321–1349, 2002. [26] M. Guven, B. Yazici, K. Kwon, E. Giladi, and X. Intes, “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: I,” Inv. Probl., vol. 23, pp. 1115–1133, 2007. [27] J. Kaipio and E. Somersalo, Computational and Statistical Inverse Problems, ser. Applied Mathematical Sciences. New York: Springer-Verlag, 2005, vol. 160. [28] M. Guven, L. Zhou, L. Reilly-Raska, and B. Yazici, “Discretization error analysis and adaptive meshing algorithms for fluorescence diffuse optical tomography: Part II,” IEEE Trans. Med. Imag., to be published. [29] 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., vol. 187, pp. 597–619, 2003. [30] R. Kress, Linear Integral Equations, ser. Applied Mathematical Sciences, 2nd ed. New York: Springer-Verlag, 1999, vol. 82. [31] M. Guven, B. Yazici, and V. Ntziachristos, “Optical fluorescence tomography with a priori information,” Proc. SPIE, vol. 6431, 2007. [32] N. P. Galatsanos and A. K. Katsaggelos, “Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation,” IEEE Trans. Image Process., vol. 1, pp. 322–336, Jul. 1992. [33] M. Hanke and T. Raus, “A general heuristic for choosing the regularization parameter in ill-posed problems,” SIAM J. Sci. Comput., vol. 17, pp. 956–972, 1996. [34] R. Molina, A. K. Katsaggelos, and J. Mateos, “Bayesian and regularization methods for hyperparameter estimation in image restoration,” IEEE Trans. Image Process., vol. 8, no. 2, pp. 231–246, Feb. 1999. [35] M. Heath, G. Golub, and G. Wahba, “Generalized cross-validation as a method for choosing a good ridge parameter,” Technometrics, vol. 21, p. 215C223, 1979. [36] P. Hansen and D. O’Learly, “The use of the l-curve in the regularization of discrete ill-posed problems,” SIAM J. Comput., vol. 14, p. 1487C1503, 1993. [37] G. Aubert and P. Kornprobst, Mathematical Problems in Image Processing. New York: Springer Verlag, 2002. [38] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods. New York: Springer-Verlag, 2002. [39] D. Daners, “Robin boundary value problems on arbitrary domains,” Trans. Am. Math. Soc., vol. 352, no. 9, pp. 4207–4236, 2000. [40] W. Bangerth, R. Hartmann, and G. Kanschat, “Deal.II—A general-purpose object-oriented finite element library,” ACM Trans. Math. Softw., vol. 33, no. 4, 2007. [41] I. M. Smith and D. V. Griffiths, Programming the Finite Element Method, 4th ed. New York: Wiley, 2004.