OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 12, Iss. 22 — Nov. 1, 2004
  • pp: 5402–5417
« Show journal navigation

Adaptive finite element based tomography for fluorescence optical imaging in tissue

Amit Joshi, Wolfgang Bangerth, and Eva M. Sevick-Muraca  »View Author Affiliations


Optics Express, Vol. 12, Issue 22, pp. 5402-5417 (2004)
http://dx.doi.org/10.1364/OPEX.12.005402


View Full Text Article

Acrobat PDF (1727 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A three-dimensional fluorescence-enhanced optical tomography scheme based upon an adaptive finite element formulation is developed and employed to reconstruct fluorescent targets in turbid media from frequency-domain measurements made in reflectance geometry using area excitation illumination. The algorithm is derived within a Lagrangian framework by treating the photon diffusion model as a constraint to the optimization problem. Adaptively refined meshes are used to separately discretize maps of the forward/adjoint variables and the unknown parameter of fluorescent yield. A truncated Gauss-Newton method with simple bounds is used as the optimization method. Fluorescence yield reconstructions from simulated measurement data with added Gaussian noise are demonstrated for one and two fluorescent targets embedded within a 512ml cubical tissue phantom. We determine the achievable resolution for the area-illumination/area-detection reflectance measurement geometry by reconstructing two 0.4cm diameter spherical targets placed at at a series of decreasing lateral spacings. The results show that adaptive techniques enable the computationally efficient and stable solution of the inverse imaging problem while providing the resolution necessary for imaging the signals from molecularly targeting agents.

© 2004 Optical Society of America

1. Introduction

Molecular Imaging is a rapidly advancing research area with the potential of providing early diagnosis and identification of the underlying biochemical causes of human diseases [1

1. M. G. Pomper, “Molecular Imaging: An Overview,” Acad. Radiol. 8, 1141–1153 (2001). [CrossRef] [PubMed]

]. As near infrared (NIR) light can travel several centimeters in tissue, fluorescence enhanced NIR optical imaging promises to open new pathways for the characterization of biological processes in living animals at cellular and molecular levels.

In the past decade, several approaches have been proposed for fluorescence enhanced optical tomography involving the determination of the fluorophore yield and/or fluorescence lifetime distribution in the tissue from a finite number of boundary measurements. Due to the diffusive nature of photon propagation in tissue, fluorescence-enhanced tomography represents a non-trivial inverse problem. Initial efforts did not involve classical optimization but instead focused on approaches such as: perturbative localization [2

2. M. A. O’Leary, D. A. Boas, B. Chance, and A. Yodh, “Reradiation and imaging of diffuse photon density waves using fluorescent inhomogeneities,” J. Luminescence 60, 281–286 (1994). [CrossRef]

, 3

3. J. Wu, Y. Wang, L. Perleman, I. Itzkan, R. R. Desai, and M. S. Feld, “Time resolved multichannel imaging of fluorescent objects embedded in turbid media,” Opt. Lett. 20, 489–491 (1995). [CrossRef] [PubMed]

, 4

4. E. L. Hull, M. G. Nichols, and T. H. Foster, “Localization of Luminescent Inhomogeneities in Turbid Media with Spatially Resolved Measurements of CW Diffuse Luminescence Emittance,” Appl. Opt. 37, 2755–2765 (1998). [CrossRef]

], backprojection [5

5. J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc. Am. A 14(275–279) (1997).

], Born Approximation [6

6. M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 20, 426–428 (1996). [CrossRef]

, 7

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

, 8

8. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26(12), 893–895 (2001). [CrossRef]

], random walk theory [9

9. V. Chernomordik, D. Hattery, I. Gannot, and A. H. Gandjbakhche, “Inverse method 3-D reconstruction of localized in vivo fluorescence-application to Sjøgren syndrome,” IEEE J. Sel. Top. Quantum Electron. 54, 930–935 (1999).

], and more recently, fast fluorescence localization [10

10. H. Quan and Z. Guo, “Fast 3-D Optical Imaging With Transient Fluorescence Signals,” Opt. Express 12(3), 449–457 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-3-449. [CrossRef]

]. These approaches were successful in locating small fluorescent targets in otherwise homogeneous, small sample volumes. Other approaches to fluorescence tomography cast the image reconstruction problem as an optimization problem in which a least squares type minimization is performed in order to determine the fluorescence map which best predicts the measured boundary fluorescence distribution. Optimization approaches are more general in their scope and can handle heterogeneous backgrounds as well as large sample volumes albeit at increased computational cost. These approaches include algorithms based on Newton’s or Newton-type optimization methods [11

11. D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36(10), 2260–2272 (1997). [CrossRef]

, 12

12. R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s Optimization Schemes for Absorption and Fluorescence Optical Tomography: part(1) theory and formulation,” Opt. Express 4(10), 353–371 (1999). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353. [CrossRef]

, 13

13. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,” Appl. Opt. 37(22), 5337–5343 (1998). [CrossRef]

] and Bayesian nonlinear least squares approaches [14

14. M. J. Eppstein, D. E. Dougherty, T. L. Troy, and E. M. Sevick-Muraca, “Biomedical optical tomography using dynamic parametrization and Bayesian conditioning on photon migration measurements,” Appl. Opt. 38(10), 2138–2150 (1998).

, 15

15. M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three dimensional near infrared fluorescence tomography with Bayesian methodologies for image reconstruction from sparse and noisy data sets,” Proc. Nat. Acad. Sci. 99, 9619–9624 (2002). [CrossRef] [PubMed]

, 16

16. A. Godavarty, M. J. Eppstein, C. Zhang, S. Theru, A. B. Thompson, M. Gurfinkel, and E. M. Sevick-Muraca, “Fluorescence-enhanced optical imaging in large tissue volumes using a gain-modulated ICCD camera,” Phys. Med. Biol. 48, 1701–1720 (2003). [CrossRef] [PubMed]

, 17

17. A. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. Boas, and R. P. Milane, “Fluorescence Optical Diffusion Tomography,” Appl. Opt. 42(16), 3061–3094 (2003).

].

The sensitivity of fluorescence enhanced imaging potentially rivals that of the conventional, but “gold-standard” molecular imaging using radiotracers [18

18. J. P. Houston, S. Ke, W. Wang, C. Li, and E. M. Sevick-Muraca, “Optical and nuclear image quality analysis with invivo NIR fluorescence and conventional gamma images acquired using a dual labeled tumor targeting probe,” J. Biomed. Opt. (submitted) (2004).

]. Yet, the achievable resolution for fluorescence tomography is determined first by the signal to noise ratio, and secondly by the level of discretization. To date, the discretization level is selected a priori based on knowledge of the domain and/or computational constraints. Image quality can be improved by uniformly refining the level of discretization throughout the domain. However, this global refinement further increases the ill-posedness of the problem and results in insurmountable computational requirements by increasing the number of unknowns. For example, to achieve a resolution of one millimeter in a volume of one liter would require the use of 10 9 mesh points, a number that is clearly not achievable with today’s computational technologies. In contrast, adaptive mesh refinement provides fine mesh resolution around target locations with coarser resolution in other regions to improve image quality, while maintaining solution stability and computational economy.

The outline of this article is as follows: In Section 2 the fluorescence tomography algorithm is described in a continuous function space setting and its discrete implementation with an adaptive mesh strategy is presented. Sections 3 and 4 detail the computational experiments conducted to demonstrate the adaptive tomography algorithm for reconstructing single and resolving dual fluorescence targets embedded in turbid media with a target to background ratio of 100 : 1. Finally, the results are summarized in Section 5 and we show that Rayleigh resolution (i.e. the minimum distance between which two target can be resolved from one another) can approach the continuum limit using the adaptive strategy. To the best of the authors’ knowledge, this contribution represents the first time that adaptivity has been used in optical tomography to address the issue of reconstructed image quality.

2. Methodology

In this section, the formulation for the nonlinear inverse problem of fluorescence tomography is developed in a continuous function space setting, allowing separate and independent discretization of the parameter map and the finite element mesh used to solve state/adjoint problems for the nonlinear update steps toward the optimal solution. The mesh refinement criteria and implementation described herein are based on the general framework developed in [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 20

20. W. Bangerth, “A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,” Tech. Rep. 04–39, ICES, University of Texas at Austin (2004).

].

2.1. Formulation

Under conditions of multiple scattering, the generation and propagation of diffuse fluorescence photon density waves from modulated, time-periodic sources can be described by the following coupled system of diffusion equations [21

21. E. M. Sevick-Muraca, E. Kuwana, A. Godavarty, J. P. Houston, A. B. Thompson, and R. Roy, Near Infrared Fluorescence Imaging and Spectroscopy in Random Media and Tissues, chap. 33, Biomedical Photonics Handbook (CRC Press, 2003).

]:

.[Dx(r)u(r,ω)]+kxu(r,ω)=0,
(1)
[Dm(r)v(r,ω)]+kmv(r,ω)=βxmu(r,ω),
(2)

whereDx,m=13(μax,mi+μax,mf+μsx,m), kx,m=iωc+μax,mi(r)+μax,mf(r), βxm=ϕμaxf1iωτ(r).

2Dxun+γu+S(r)=0,2Dmvn+γv=0,
(3)

where n denotes the outward normal to the surface and γ is a constant depending on the optical reflective index mismatch at the boundary [22

22. A. Godavarty, D. J. Hawrysz, R. Roy, E. M. Sevick-Muraca, and M. J. Eppstein, “Influence of the refractive index-mismatch at the boundaries measured in fluorescenceenhanced frequency-domain photon migration imaging,” Opt. Express 10(15), 650–653 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-15-653.

]. S(r) is the excitation boundary source. There is no source term for the emission boundary condition. Note that, here the NIR excitation source is modeled as a boundary condition. This is an approximation as the true isotropic source will be one scattering length below the illumination surface. This approximation is justified for a diffuse area illumination scheme [23

23. A. B. Thompson and E. M. Sevick-Muraca, “NIR fluorescence contrast enhanced imaging with ICCD homodyne detection: measurement precision and accuracy,” J. Biomed. Opt. 8, 111–120 (2002). [CrossRef]

] instead of the traditional point illumination by fiber optics and the finite element solutions of equations (1)(3) match the experimentally observed boundary fluorescence [24

24. A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element methods for fluorescence enhanced frequency domain optical tomography: Forward imaging problem,” in International Symposium on Biomedical Imaging, pp. 1103–1106 (IEEE, 2004).

]. The goal of fluorescence tomography is the reconstruction of a spatial map of coefficients µaxf (r) and/or τ(r) from measurements of the fluences u,v on the boundary.

Instead of the (strong) formulation as a PDE above, equations (1)(3) are posed in a weak variational form as is usual in finite element applications. The two equations are multiplied with arbitrary test functions ζ,ξ ∈ H 1, integrated over Ω, and terms with second derivatives are integrated by parts. Here, H 1 is the Sobolev space of functions with integrable (weak) first derivatives (see Adams [25

25. R. A. Adams, Sobolev Spaces (Academic Press, 1975).

]). With this, the resulting variational equation reads:

A(q;[u,v])([ζ,ξ])=0ζ,ξH1,
(4)

where the semilinear form A is nonlinear in its first set of arguments and linear in the test functions. q denotes the set of unknown parameters, i.e. µ axf and/or τ. In this work, we only consider q=µaxf . With (·,·) denoting the L 2 inner product, the definition of A reads

A(q;[u,v])([ζ,ξ])=(Dxu,ζ)Ω+(kxu,ζ)Ω+γ2(u,ζ)Ω+12(S,ζ)Ω
+(Dmv,ξ)Ω+(kmv,ξ)Ω+γ2(v,ξ)Ω(βxmu,ξ)Ω.
(5)

Note that Dx,m =Dx,m (q),kx,m =kx,m (q),βxm =βxm (q).

The goal of the parameter identification problem is to find the set of parameters q for which the predicted boundary fluorescence measurements resulting from equation (4) best match actual measurements. The quantity to be minimized (i.e. the misfit between prediction and measurement) is then

J(q,v)=12vzΣ2+βr(q),
(6)

where the L 2 norm of the difference between the actual measurements, z, and the prediction of the emission fluence v on a part ∑ of the boundary Ω is minimized. In practice, z is interpolated between the pixels of an area detection system such as the gain modulated CCD camera. r(q) is a Tikhonov regularization functional which penalizes certain undesirable aspects of solutions [26

26. A. N. Tikhonov and V. Y. Arsenin, eds., Solution of Ill-Posed Problems (Winston, Washington, DC, 1977).

] and β is the regularization parameter. Information about measurement noise can be incorporated by using a different or weighted norm of the misfit [21

21. E. M. Sevick-Muraca, E. Kuwana, A. Godavarty, J. P. Houston, A. B. Thompson, and R. Roy, Near Infrared Fluorescence Imaging and Spectroscopy in Random Media and Tissues, chap. 33, Biomedical Photonics Handbook (CRC Press, 2003).

].

In the past, optical tomography approaches have followed the output least squares formulation where the state variables u,v are taken to be dependent on the parameters q [27

27. S. R. Arridge and M. Schweiger, “A Gradient Based Optimization Scheme for Optical Tomography,” Opt. Express 2(6), 213–225 (1998). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-2-6-213. [CrossRef]

, 12

12. R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s Optimization Schemes for Absorption and Fluorescence Optical Tomography: part(1) theory and formulation,” Opt. Express 4(10), 353–371 (1999). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353. [CrossRef]

]. In this work, however, we minimize the error functional J by treating {q,u,v} as independent variables where their relationship is enforced by including the state equation as a constraint to the optimization problem [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 20

20. W. Bangerth, “A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,” Tech. Rep. 04–39, ICES, University of Texas at Austin (2004).

, 28

28. L. Beilina, “Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,” Ph.D. thesis, Chalmers University of Technology (2002).

, 29

29. R. Becker, H. Kapp, and R. Rannacher, “Adaptive Finite Element Methods for Optimal Control of Partial Differential Equations: Basic Concept,” SIAM J. Contr. Optim. 39, 113–132 (2000). [CrossRef]

, 30

30. R. Luce and S. Perez, “Parameter identification for an elliptic partial differential equation with distributed noisy data,” Inverse Problems 15, 291–307 (1999). [CrossRef]

]. We also deviate from the larger part of the literature in that we formulate the following steps in function spaces, rather than choosing the route of most optical tomography papers to first discretize state equation and objective function and then stating the optimization problem in the finite-dimensional space of matrices and vectors. While our approach leads to equations that at times look cumbersome, they are exactly equivalent to the matrix formulation if a fixed grid is used. However, they afford for the possibility of using different grids for different nonlinear iterations; in other words, we do not yet want to settle on a fixed discretization and therefore cannot choose the more usual matrix formulation.

With these considerations, the minimization problem reads:

minq,u,vJ(q,v)subjecttoA(q;[u,v])([ζ,ξ])=0.
(7)

This constrained optimization problem can be considered in a Lagrangian framework where the Lagrangian is defined as:

L([u,v],[λex,λem],q)=J(q,v)+A(q;[u,v])([λex,λem]).
(8)

Here, λexem are the Lagrange multipliers corresponding to the excitation and emission diffusion equation constraints, respectively. For simplicity, the abbreviation of x={u,v,λ exem,q} is introduced so that the Lagrangian functional can be written as L(x). The optimum is then characterized by a stationary point of the Lagrangian, i.e. the first order conditions:

Lx(x)(y)=0y={φex,φem,ψex,ψem,χ},
(9)

where Lx (x)(y) is the Frèchet differential [31

31. D. G. Luenberger, Optimization by Vector Space Methods (John Wiley, 1969).

] of L(x) and y denotes possible test functions. Eq. (9) can be expanded for all components of x:

Lu(x)(φex)=Au(q;[u,v])(φex)([λex,λem])=0,
(10)
Lv(x)(φem)=Jv(q,v)(φem)+Av(q;[u,v])(φem)([λex,λem])=0,
(11)
Lλex(x)(ψex)=A(q;[u,v])([ψex,0])=0,
(12)
Lλem(x)(ψem)=A(q;[u,v])([0,ψem])=0,
(13)
Lq(x)(χ)=Jq(q,v)(χ)+Aq(q;[u,v])(χ)([λex,λem])=0.
(14)

Subscripts q,u,v,λ ex , and λ em indicate first order partial Fréchet derivatives of J or A. Equations (12)(13) are the state equations in variational form. Equation (12) can be solved to provide the excitation fluence u which is then used to solve equation (13) to obtain the emission fluence v. Equations (10)(11) are the adjoint equations defining the Lagrange multipliers [λ exem ]. Finally, equation (14) is the control equation.

The above set of coupled nonlinear equations is solved by Newton’s method. The update direction for the kth iteration, δxk ={δuk,δvk, δλkex, δλkem,δqk }, is determined from

Lxx(xk)(δxk,y)=Lx(xk)(y)y,
(15)

where Lxx (xk ) is the Hessian matrix of second derivatives of L at point xk . These equations represent one condition for each variable in δxk and when expanded read as follows [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 31

31. D. G. Luenberger, Optimization by Vector Space Methods (John Wiley, 1969).

]:

Au(qk;[uk,vk])(ϕex)([δλkex,0])+Au(qk;[uk,vk])(ϕex)([0,δλkem])+Auq(qk;[uk,vk])(ϕex,δqk)([λex,λem])=Lu(xk)(ϕex),
Jvv(qk,vk)(δvk,ϕem)+Av(qk;[uk,vk])(ϕem)([δλkex,0])+Av(qk;[uk,vk])(ϕem)([0,δλkem])+Jvq(qk,vk)(δqk,ϕem)+Avq(qk;[uk,vk])(ϕem,δqk)([λex,λem])=Lv(xk)(ϕem),
Au(qk;[uk,vk])(δuk)([ψex,0])+Av(qk;[uk,vk])(δvk)([ψex,0])+Aq(qk;[uk,vk])(δqk)([ψex,0])=Lλex(xk)(ψex),
Au(qk;[uk,vk])(δuk)([0,ψem])+Av(qk;[uk,vk])(δvk)([0,ψem])+Aq(qk;[uk,vk])(δqk)([0,ψex])=Lλem(xk)(ψem),
Aqu(qk;[uk,vk])(δuk,χ)([λex,λem])+Aqv(qk;[uk,vk])(δvk,χ)([λex,λem])+Jqv(qk,vk)(δvk,χ)+Jqq(qk,vk)(δqk,χ)+Aq(qk;[uk,vk])(χ)([δλex,0])+Aq(qk;[uk,vk])(χ)([0,δλem])=Lq(xk)(χ).

From Eq.s (10)(11) one can infer that the Lagrange multipliers [λ exem ] are proportional to Jv (q;v). As a consequence, all terms involving [λ exem ] become negligible near the optimal solution under conditions of low noise and can be dropped from the Hessian of the Lagrange functional, resulting in the simplified Gauss-Newton method. The equations for the Gauss-Newton method are then the same as above with the exception that (i) the last terms on the left hand side of the first and second equations as well as (ii) the first and second term of the last equation are eliminated.

xk+1=xk+αkδxk.
(16)

2.2. Discretization

[M0PT0RCTPC0][δpkδqkδdk]=[F1F2F3],
(17)

where the updates for the primal and dual variables are abbreviated as δpk =[δuk,δdvk ]T, δdk =[δλkex,δλkem ]T. The blocks of the matrix are defined as

M=[000(φi,φj)Σ]ij,R=[βr"(qk,χi,χj)]ij,PT=[AglobalexBglobalexem0Aglobalem].
(18)

Here, P is the representation of the discrete forward diffusion model; Aglobalex and Aglobalem are the global stiffness matrices associated with the excitation and emission diffusion equations; and Bglobalexem is the matrix which couples the excitation and emission linear systems. The subscripts i, j in the preceding equations iterate over all degrees of freedom. The matrix C is defined as CT =[C 1,C 2] with its components obtained by differentiating the semilinear form A in Eq. (5) with respect to the parameter q:

C1=(Dx(qk)qukψi,χj)ij+(kx(qk)qukψi,χj)ij,
C2=(Dm(qk)qvkψi,χj)ij+(km(qk)qvkψi,χj)ij(βxm(qk)qukψi,χj)ij.
(19)

The right hand side terms in Eq. (17) are the discretized components of the negative gradient of the Lagrangian:

F1=[(Dx(qk)λkex,φi)(kx(qk)λkex,φi)γ2(λkex,φi)Ω+(βxm(qk)λkem,φi)(vkz,φi)Σ(Dm(qk)λkem,φi)(km(qk)λkem,φi)γ2(λkem,φi)Ω]i,
F2=[βr(qk,χi)(Dx(qk)qukλkex,χi)(kx(qk)qukλkem,χi)(Dm(qk)qvk.λkem,χi)(km(qk)qvkλkem,χi)+(βxm(qk)qukλkem,χi)]i,
F3=[Dx(qk)ψi,uk)(kx(qk)ψi,uk)12(Sx,ψi)Ωγ2(ψi,uk)ΩDm(qk)ψi,vk)(km(qk)ψi,vk)γ2(ψi,vk)Ω+(βxm(qk)ψi,uk)]i.
(20)

For large domains and fine discretization, the linear system Eq. (17) can be as large as several 100,000 unknowns. Furthermore, the matrix is generally indefinite, restricting the choice of available iterative linear solvers for determining the Gauss-Newton update directions. Thus, instead of solving Eq. (17) directly, we use an efficient solver based on the Schur complement developed in [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 20

20. W. Bangerth, “A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,” Tech. Rep. 04–39, ICES, University of Texas at Austin (2004).

]. By block elimination, Eq. (17) is reduced to the following sequence of three equations:

{R+CTPTMP1C}δqk=F2CTPTF1+CTPTMP1F3,
(21)
Pδpk=F3Cδqk,
(22)
PTδdk=F1Mδpk.
(23)

The matrix R+CTP-TMP-1C is the Schur complement matrix of the Gauss-Newton system. Under practical conditions it is symmetric positive definite [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

] and can thus be inverted efficiently by iterative solvers such as the conjugate gradient method. Furthermore, the Schur complement matrix is comparatively small (at most a few 1,000 to 10,000) with its size equivalent to the number of discretized parameters, rather than the number of parameters, state, and adjoint variables combined. Since the Newton method updates only approximates the direction to the optimal solution, it is not necessary to solve Eq. (21) exactly for each Newton step. The conjugate gradient iteration is thus stopped once the l 2 residual falls below a certain tolerance, for example 10-3 times its initial value. Consequently, this method is a variant of truncated or inexact Gauss-Newton methods [32

32. J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research (Springer, New York, 1999). [CrossRef]

].

2.3. Adaptive Mesh Refinement

The accuracy of finite element solutions of partial differential equations depends on the mesh width. This is reflected in a priori error estimates comparing the exact solution u and the numerically computed finite element solution uh. For the simple example of the Laplace equation, an a priori error estimate would be [33

33. S. C. Brenner and R. L. Scott, The Mathematical Theory of Finite Elements (Springer, Berlin-Heidelberg-New York, 1994).

]:

uuhC(u)h2,
(24)

where h is the maximum mesh size and C(u) is a constant that depends on the exact (and unknown) solution, the domain, and element shapes, but not upon the mesh width. Thus, the accuracy of the finite element solution can be increased by decreasing the maximum mesh size h. Similar estimates may be derived for the system in Eq.s (1)(2). The drawback of a priori estimates is that the exact solution, and thus the numerical value of C(u), is unknown, so that no quantitative bound on the error can be computed.

While Eq. (24) shows that it can be guaranteed that global mesh refinement reduces the error in the solution, this does not represent an efficient strategy in general. A fine mesh is only necessary where the solution varies greatly and where small cells are needed to accurately capture this variation. Since it is in general unknown in advance where these locations are, a posteriori error estimates have been derived in the mathematical community to provide criteria for local mesh refinement. These estimates use the computed solution u h to not only provide a bound on the error ‖u-uh ‖ without requiring knowledge of the exact solution u, but also to indicate on which cells the contribution to this error is largest. Thus, these estimates indicate the cells for which mesh refinement will be most beneficial, and conversely for which cells, mesh refinement will not yield a significant contribution to the reduction of the error. Using this process, the mesh on which uh was computed can be appropriately refined, and the solution is computed again on the refined mesh; this process is iterated until either the error estimate indicates that the requested accuracy is obtained, or computational resources are exhausted. The method in which only selected cells are repeatedly refined based upon an error indicator is commonly referred to as adaptive mesh refinement. For nonlinear problems, an additional advantage is realized by performing the initial Gauss-Newton updates on coarse meshes: finer and thus computationally more expensive meshes are employed only near the solution which reduces both the ill-posedness of the problem and the computational work in the initial steps.

A posteriori error estimates and adaptive mesh refinement have been extensively studied in the last two decades. See [34

34. R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques (Wiley/Teubner, New York, Stuttgart, 1996).

, 35

35. W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations (Birkhäuser Verlag, 2003).

] and the references therein for an overview. However, most of the previouswork is tailored to model equations rather than parameter estimation problems with the exception of Molinari et al. [36

36. M. Molinari, S. J. Cox, B. H. Blott, and G. J. Daniell, “Adaptive Mesh Refinement techniques for Electrical Impedence Tomography,” Physiological Measurement 22, 91–96 (2001). [CrossRef] [PubMed]

, 37

37. M. Molinari, B. H. Blott, S. J. Cox, and G. J. Daniell, “Optimal Imaging with Adaptive Mesh Refinement in Electrical Impedence Tomography,” Physiological Measurement 23, 121–128 (2002). [CrossRef] [PubMed]

] for adaptive mesh refinement in electrical impedance tomography (EIT); see also [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 20

20. W. Bangerth, “A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,” Tech. Rep. 04–39, ICES, University of Texas at Austin (2004).

, 29

29. R. Becker, H. Kapp, and R. Rannacher, “Adaptive Finite Element Methods for Optimal Control of Partial Differential Equations: Basic Concept,” SIAM J. Contr. Optim. 39, 113–132 (2000). [CrossRef]

, 38

38. R. Becker and R. Rannacher, “An optimal control approach to error estimation and mesh adaptation in finite element methods,” Acta Numerica 10, 1–102 (2001). [CrossRef]

, 28

28. L. Beilina, “Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,” Ph.D. thesis, Chalmers University of Technology (2002).

, 39

39. H. Ben Ameur, G. Chavent, and J. Jaffré, “Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities,” Inverse Problems 18, 775–794 (2002). [CrossRef]

, 40

40. A.-A. Grimstad, H. Krüger, T. Mannseth, G. Nævdal, and H. Urkedal, “Adaptive selection of parameterization for reservoir history matching,” in ECMOR VIII: 8th European Conference on the Mathematics of Oil Recovery, Freiberg, Germany, pp. E–46 (European Association of Geoscientists and Engineers (EAGE), 2002).

, 41

41. R. Li, W. Liu, H. Ma, and T. Tang, “Adaptive finite element approximation for distributed elliptic optimal control problems,” SIAM J. Control Optim. 41, 1321–1349 (2002). [CrossRef]

] for approaches in other fields.

In this work, we use two separate meshes that are adaptively refined. The first one is used for the discretization of state and adjoint variables u,v and λexem , while the second one is used for the discretization of the parameter field q. The state/adjoint mesh is finer than the parameter mesh in order to avoid stability problems with the saddle point problem Eq. (17). In addition, and in accordance with the regularity requirements of the equation, we discretize state and adjoint variables using piecewise tri-linear finite elements, and the parameter with piecewise constant functions. Whenever Gauss-Newton iterations on these meshes have reduced the error function by a significant amount, both meshes are refined using a posteriori refinement criteria. In this work, the state and adjoint mesh is refined using a variation of the refinement criterion first derived by Kelly et al. [42

42. D. W. Kelly, J. P. d. S. R. Gago, O. C. Zienkiewicz, and I. Babuška, “A posteriori error analysis and adaptive processes in the finite element method: Part I-Error Analysis,” Int. J. Num. Meth. Engrg. 19, 1593–1619 (1983). [CrossRef]

]:

ηKu=h24nuhK2,ηKv=h24nvhK2,ηK=αηKu+(1α)ηKv,
(25)

where [nuh ] is the jump of the normal derivative of the finite element solution u h across the element boundary K , and α is chosen so that the errors in the two variables are roughly weighted equally. In our implementation, we refine those 35% of elements with highest error indicator ηK and coarsen those 5% of elements with lowest errors in each cycle. Some other elements are also refined to maintain numerical stability of the solution.

The mesh for q is refined by computing, for each cell, a discrete approximation to the gradient of q weighted by the local mesh width [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

]. Note that the choice of two separate meshes means that the first mesh can be fine close to the source where the excitation fluence greatly varies, while the second mesh will only be fine close to the fluorescent target and coarse everywhere else. The detailed derivation of error estimates can be found elsewhere [19

19. W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

, 28

28. L. Beilina, “Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,” Ph.D. thesis, Chalmers University of Technology (2002).

, 35

35. W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations (Birkhäuser Verlag, 2003).

].

2.4. Software Implementation

A simplified block diagram of the tomography algorithm is shown in Fig. 1. At the beginning of the algorithm, all variables are initialized on coarse meshes (mesh level m=1): state/adjoint variables are set to zero, while the parameter is set to the lower bound q 0 (the minimal background concentration). In each subsequent iteration, Gauss-Newton update directions are determined from Eq.s (21)(23) and all variables are updated after computating the step length. The program is terminated if the number of Gauss-Newton iterations exceeds the maximum number of iterations k max or if the model misfit 12vzΣ2 is reduced below a prespecified threshold ε. For the results reported in this paper we used k max=40 and ε=10-15.

Mesh refinement is triggered if either (i) the Gauss-Newton step length α k is less than a prespecified minimum step length α min, or (ii) the nonlinear residual, rk =‖Lx (xk )(·)‖, of the optimality condition (9) has been sufficiently reduced on the current mesh m, i.e. rkεmeshrkm0, where km0 denotes the first iteration on mesh level m. In this work α min=0.15 and the error reduction threshold is ε mesh=10-4. Mesh refinement is followed by re-computation of synthetic measurements z on the new, finer mesh.

This scheme for fluorescence tomography with adaptive mesh refinement was implemented in C++ based on the deal.II finite element library [43

43. W. Bangerth, R. Hartmann, and G. Kanschat, deal.II Differential Equations Analysis Library, Technical Reference (2004). http://www.dealii.org/.

]. Deal.II provides advanced object oriented design techniques and support for the complex data structures needed for adaptive finite element applications. While the current version of the software is used for an area excitation illumination geometry as previously published [23], the implementation is also suitable for point illumination geometries. Reconstructions were carried out on a dual processor (750MHz) Sun Sparc workstation with 3GB of RAM.

3. Computational Experiments

Two computational experiments were performed to test the efficacy of the proposed algorithm: image reconstruction of (i) a single fluorescent target, and (ii) of two closely spaced fluorescence targets of varying separation distance. The synthetic frequency-domain fluorescence data was generated on a simulated 8×8×8 cm3 cube illuminated by a simulated expanded laser beam with a Gaussian profile at the x=0 plane, as shown in Fig. 2. This measurement geometry and phantom corresponds to experimental measurements previously reported [23

23. A. B. Thompson and E. M. Sevick-Muraca, “NIR fluorescence contrast enhanced imaging with ICCD homodyne detection: measurement precision and accuracy,” J. Biomed. Opt. 8, 111–120 (2002). [CrossRef]

]. Fluorescence phase and amplitude measurements were generated at the illumination surface in order to mimic the actual experimental data collection. Random Gaussian noise of zero mean and specified half width is then applied to these synthetic measurements.

Fig. 1. Adaptive tomography algorithm. GN stands for Gauss-Newton; see Section 2.4 for a description of symbols.

3.1. Single Fluorescent Target

Data was generated for a 0.5cm diameter spherical target at a depth of 2.15cm from the illumination surface at an off-center position (x=2.15cm, y=3.15cm, z=3.15cm). The absorption and isotropic scattering properties of 1% Liposyn solution were chosen to mimic the background of the phantom, corresponding to actual experimental measurements using a gain modulated ICCD camera system [23

23. A. B. Thompson and E. M. Sevick-Muraca, “NIR fluorescence contrast enhanced imaging with ICCD homodyne detection: measurement precision and accuracy,” J. Biomed. Opt. 8, 111–120 (2002). [CrossRef]

]. For the simulated background absorption coefficient, we chose µaxi =0.023cm -1 and µami =0.0289cm -1 [44

44. A. Thompson, “Development of a new optical imaging modality for detection of fluorescence enhanced disease,” Ph.D. thesis, Texas A & M University (2003).

]. The absorption coefficient due to fluorophore at the excitation wavelength was set to µaxf =0.5cm -1 in the target and 0.005cm -1 in the background. The emission wavelength absorption coefficient was µ amf =0.0506cm -1 in the target and 0.00506cm -1 in the background. The lifetime of the fluorophore was taken to be τ=0.56ns and the quantum efficiency was ϕ=0.016 to match the corresponding properties of Indocyanine Green (ICG) dye used in experiments. The excitation wavelength for ICG is 785nm and the emission data is collected at 830nm. The reduced scattering coefficient was taken as µ′s =9.84cm -1 in both the target and the background, and for this study was taken to be the same at excitation and emission wavelengths. Two percent random Gaussian noise was added to real and imaginary parts of the synthetic emission fluence solution at the measurement surface at x=0.

Fig. 2. Area illumination and area detection geometry employed by Thompson et al. [23]
Fig. 3. Single target reconstruction: A black wire-frame depicts the actual target and colored blocks represent the reconstruction. Top 10% of the contour levels of µaxf are shown.

3.2. Two Fluorescent Targets

Tomographic reconstructions were performed from synthetic data generated from the solution of Eq.s (1)(3) for two closely spaced fluorescent targets to test the resolving ability of the adaptive inverse algorithm. Two spherical targets with diameter 0.4cm were placed at an edge to edge distance varying from 1cm to 0.1cm, the latter of which is the continuum limit of the diffusion equation for the chosen optical properties. The depth of both targets was simulated to be 1.2cm from the illumination and measurement plane at x=0. The optical properties of targets and background as well as the noise level remained the same as in the single target study.

4. Results

4.1. Single Fluorescent Target

Figure 3 illustrates the tomographic reconstruction for the single target case after 22 Gauss-Newton update steps utilizing approximately 2hrs of CPU time. The location of the target is reconstructed accurately. However, the magnitude 0.0176cm -1 of reconstructed µaxf is lower than the actual value 0.5cm -1, which is an artifact of the L 2 regularization used r(q)=12q2, with a fixed value β=10-3). Figure 4 shows the evolution of state/adjoint and parameter meshes during the reconstruction process. The algorithm started with coarse initial meshes with 64 hexahedral elements. Five automatic mesh refinements were carried out by the algorithm. The final state/adjoint mesh consisted of 116,936 nodes and was refined predominantly on the illumination plane to accurately resolve the Gaussian NIR excitation source. The final parameter mesh had 1016 elements, mostly located around the reconstructed fluorescent target. A uniformly refined mesh with the same parameter mesh resolution as the adaptive mesh surrounding the target would have 323=32,768 unknowns. The advantages of adaptive refinement for the reduction of the total number of unknowns and the improvement of resolution are obvious.

Fig. 4. Adaptive mesh evolution for state/adjoint (left) and parameter discretization (right). Meshes are shown at 1st, 11th and 22nd Gauss-Newton iterations.

Table 1. Summary of results for dual fluorescent target reconstructions. d is the edge to edge target separation in cm; Iter. is the Gauss-Newton iteration for which the other results are reported; ‖q-qtrue‖2 is the error in reconstructed parameter; 12vzΣ2 is the meausurement error; Nq is the number of elements (unknowns) in the parameter mesh.

table-icon
View This Table

4.2. Two Fluorescent Targets

Figure 5 shows the tomographic reconstructions for targets separated by 1.0142, 0.6607, 0.3071, and 0.1657cm. The reconstruction process was started with the same coarse state and parameter meshes as used in the single target case. A heuristic criterion was used to assess the reconstructions of the two targets: the top 10% contour levels of the µaxf map were plotted and if two distinct maximum were visible, then the two targets were considered to be identified separately. The centroid of these two maxima was computed and treated as the reconstructed centroids. The centroids of the two targets was reconstructed accurately for target separations of 0.1657cm and greater, when the reconstructed targets begin to appear as one.

At a target spacing of 0.1cm, the reconstruction is biased towards one of the targets, and no maximum associated with a second target is resolved at this distance (see Fig. 6). This is consistent with the fact that this separation is approximately the mean isotropic scattering length, below which the continuum assumption necessary for the validity of the diffusion approximation is violated. No image resolution below this threshold can be expected.

Table 1 summarizes the computational results for the dual target reconstructions, demonstrating the efficacy of adaptive mesh refinement for determining the maximum achievable resolution for a given experimental configuration in the presence of measurement noise. The reconstructed targets appear slightly closer than the actual separation as the target spacing is reduced. This can be explained by the merging of the fluorescence emission fields with decreasing distance when using only the reflectance data for image reconstruction as we have done here.

5. Conclusions

In this work, we have demonstrated an efficient adaptive finite element algorithm for fluorescence optical tomography. The algorithm is formulated in function spaces independent of any a priori discretization in order to allow for meshes to change as nonlinear (Gauss-Newton) iterations progress. Independent meshes have been used for the state/adjoint variables and the parameter map, and physical information about upper and lower bounds on the unknown parameters have been incorporated to improve the stability of the algorithm. The choice of separately adapted meshes (with the parameter mesh being coarser), different shape functions for the state and parameter variables, as well as the inclusion of bounds on the parameter and a Tikhonov regularization term yield an algorithm that is able to cope with the well known ill-posedness of optical tomography for comparable high resolution imaging with acceptable numerical effort in the presence of noise. The discretization scheme can be contrasted with that recently proposed by Huang et al. [45

45. M. Huang and Q. Zhu, “Dual-mesh optical tomography reconstruction method with a depth correction that uses a priori ultrasound information,” Appl. Opt. 43(8), 1654–1662 (2004). [CrossRef]

], in which prior information about the subdivision (discretization) of the domain is obtained from ultrasound images. While their strategy requires that the target is structurally dissimilar to the surrounding tissue to obtain meaningful information via ultrasound or MRI, our algorithm is able to find high resolution discretizations automatically by using the solutions on coarse grids to generate finer grids where necessary. To reiterate, the final number of unknowns (between 729 and 1359) is 25 to 40 times less than that necessary with the use of conventional finite element based optical tomography schemes for large tissue volumes while still retaining the resolution in the regions of interest. In addition, the reduced number of unknowns in our adaptive finite element scheme greatly reduces the ill-posedness of the inverse problem and thus aids in the quality of the reconstruction. Jiang et al. [46

46. X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. 30(5), 861–869 (2003). [CrossRef]

] have also proposed a dual and adaptive meshing scheme for diffuse optical tomography, However their scheme doesn’t implement automatic mesh adaptation with a posteriori error estimates.

Fig. 5. Dual target reconstructions: A black wire-frame depicts the actual targets and colored blocks represent the reconstruction. Top 10% of the contour levels of µaxf are shown. Edge to edge spacing: (a) 1.0142cm, (b) 0.6607cm, (c) 0.3071cm, and (d) 0.1657cm.
Fig. 6. Dual target reconstruction for 0.1cm target separation.

We have demonstrated image reconstruction from frequency-domain fluorescence reflectance data arising from area-illumination. In comparison to circumferential point illumination measurements, reflectance measurements possess reduced information content for image reconstruction, but make the approach practical for patient imaging. This reflectance geometry and the associated 3-D adaptive finite element reconstruction algorithm could be pertinent for sentinel lymph node mapping as well as intraoperative localization of metastatic lesions (as used to track the progress of breast cancer). Our results provide the first report of a practical approach which may lead to high resolution molecular imaging in the clinic.

Acknowledgments

The authors acknowledge the support of the National Institutes of Health (NIH grant no. R01CA67176) and the Center for Subsurface Modeling at the Institute for Computational Engineering and Sciences, University of Texas at Austin.

References and links

1.

M. G. Pomper, “Molecular Imaging: An Overview,” Acad. Radiol. 8, 1141–1153 (2001). [CrossRef] [PubMed]

2.

M. A. O’Leary, D. A. Boas, B. Chance, and A. Yodh, “Reradiation and imaging of diffuse photon density waves using fluorescent inhomogeneities,” J. Luminescence 60, 281–286 (1994). [CrossRef]

3.

J. Wu, Y. Wang, L. Perleman, I. Itzkan, R. R. Desai, and M. S. Feld, “Time resolved multichannel imaging of fluorescent objects embedded in turbid media,” Opt. Lett. 20, 489–491 (1995). [CrossRef] [PubMed]

4.

E. L. Hull, M. G. Nichols, and T. H. Foster, “Localization of Luminescent Inhomogeneities in Turbid Media with Spatially Resolved Measurements of CW Diffuse Luminescence Emittance,” Appl. Opt. 37, 2755–2765 (1998). [CrossRef]

5.

J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc. Am. A 14(275–279) (1997).

6.

M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 20, 426–428 (1996). [CrossRef]

7.

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

8.

V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26(12), 893–895 (2001). [CrossRef]

9.

V. Chernomordik, D. Hattery, I. Gannot, and A. H. Gandjbakhche, “Inverse method 3-D reconstruction of localized in vivo fluorescence-application to Sjøgren syndrome,” IEEE J. Sel. Top. Quantum Electron. 54, 930–935 (1999).

10.

H. Quan and Z. Guo, “Fast 3-D Optical Imaging With Transient Fluorescence Signals,” Opt. Express 12(3), 449–457 (2004). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-3-449. [CrossRef]

11.

D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36(10), 2260–2272 (1997). [CrossRef]

12.

R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s Optimization Schemes for Absorption and Fluorescence Optical Tomography: part(1) theory and formulation,” Opt. Express 4(10), 353–371 (1999). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353. [CrossRef]

13.

H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,” Appl. Opt. 37(22), 5337–5343 (1998). [CrossRef]

14.

M. J. Eppstein, D. E. Dougherty, T. L. Troy, and E. M. Sevick-Muraca, “Biomedical optical tomography using dynamic parametrization and Bayesian conditioning on photon migration measurements,” Appl. Opt. 38(10), 2138–2150 (1998).

15.

M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three dimensional near infrared fluorescence tomography with Bayesian methodologies for image reconstruction from sparse and noisy data sets,” Proc. Nat. Acad. Sci. 99, 9619–9624 (2002). [CrossRef] [PubMed]

16.

A. Godavarty, M. J. Eppstein, C. Zhang, S. Theru, A. B. Thompson, M. Gurfinkel, and E. M. Sevick-Muraca, “Fluorescence-enhanced optical imaging in large tissue volumes using a gain-modulated ICCD camera,” Phys. Med. Biol. 48, 1701–1720 (2003). [CrossRef] [PubMed]

17.

A. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. Boas, and R. P. Milane, “Fluorescence Optical Diffusion Tomography,” Appl. Opt. 42(16), 3061–3094 (2003).

18.

J. P. Houston, S. Ke, W. Wang, C. Li, and E. M. Sevick-Muraca, “Optical and nuclear image quality analysis with invivo NIR fluorescence and conventional gamma images acquired using a dual labeled tumor targeting probe,” J. Biomed. Opt. (submitted) (2004).

19.

W. Bangerth, “Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,” Ph.D. thesis, University of Heidelberg (2002).

20.

W. Bangerth, “A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,” Tech. Rep. 04–39, ICES, University of Texas at Austin (2004).

21.

E. M. Sevick-Muraca, E. Kuwana, A. Godavarty, J. P. Houston, A. B. Thompson, and R. Roy, Near Infrared Fluorescence Imaging and Spectroscopy in Random Media and Tissues, chap. 33, Biomedical Photonics Handbook (CRC Press, 2003).

22.

A. Godavarty, D. J. Hawrysz, R. Roy, E. M. Sevick-Muraca, and M. J. Eppstein, “Influence of the refractive index-mismatch at the boundaries measured in fluorescenceenhanced frequency-domain photon migration imaging,” Opt. Express 10(15), 650–653 (2002). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-15-653.

23.

A. B. Thompson and E. M. Sevick-Muraca, “NIR fluorescence contrast enhanced imaging with ICCD homodyne detection: measurement precision and accuracy,” J. Biomed. Opt. 8, 111–120 (2002). [CrossRef]

24.

A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element methods for fluorescence enhanced frequency domain optical tomography: Forward imaging problem,” in International Symposium on Biomedical Imaging, pp. 1103–1106 (IEEE, 2004).

25.

R. A. Adams, Sobolev Spaces (Academic Press, 1975).

26.

A. N. Tikhonov and V. Y. Arsenin, eds., Solution of Ill-Posed Problems (Winston, Washington, DC, 1977).

27.

S. R. Arridge and M. Schweiger, “A Gradient Based Optimization Scheme for Optical Tomography,” Opt. Express 2(6), 213–225 (1998). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-2-6-213. [CrossRef]

28.

L. Beilina, “Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,” Ph.D. thesis, Chalmers University of Technology (2002).

29.

R. Becker, H. Kapp, and R. Rannacher, “Adaptive Finite Element Methods for Optimal Control of Partial Differential Equations: Basic Concept,” SIAM J. Contr. Optim. 39, 113–132 (2000). [CrossRef]

30.

R. Luce and S. Perez, “Parameter identification for an elliptic partial differential equation with distributed noisy data,” Inverse Problems 15, 291–307 (1999). [CrossRef]

31.

D. G. Luenberger, Optimization by Vector Space Methods (John Wiley, 1969).

32.

J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research (Springer, New York, 1999). [CrossRef]

33.

S. C. Brenner and R. L. Scott, The Mathematical Theory of Finite Elements (Springer, Berlin-Heidelberg-New York, 1994).

34.

R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques (Wiley/Teubner, New York, Stuttgart, 1996).

35.

W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations (Birkhäuser Verlag, 2003).

36.

M. Molinari, S. J. Cox, B. H. Blott, and G. J. Daniell, “Adaptive Mesh Refinement techniques for Electrical Impedence Tomography,” Physiological Measurement 22, 91–96 (2001). [CrossRef] [PubMed]

37.

M. Molinari, B. H. Blott, S. J. Cox, and G. J. Daniell, “Optimal Imaging with Adaptive Mesh Refinement in Electrical Impedence Tomography,” Physiological Measurement 23, 121–128 (2002). [CrossRef] [PubMed]

38.

R. Becker and R. Rannacher, “An optimal control approach to error estimation and mesh adaptation in finite element methods,” Acta Numerica 10, 1–102 (2001). [CrossRef]

39.

H. Ben Ameur, G. Chavent, and J. Jaffré, “Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities,” Inverse Problems 18, 775–794 (2002). [CrossRef]

40.

A.-A. Grimstad, H. Krüger, T. Mannseth, G. Nævdal, and H. Urkedal, “Adaptive selection of parameterization for reservoir history matching,” in ECMOR VIII: 8th European Conference on the Mathematics of Oil Recovery, Freiberg, Germany, pp. E–46 (European Association of Geoscientists and Engineers (EAGE), 2002).

41.

R. Li, W. Liu, H. Ma, and T. Tang, “Adaptive finite element approximation for distributed elliptic optimal control problems,” SIAM J. Control Optim. 41, 1321–1349 (2002). [CrossRef]

42.

D. W. Kelly, J. P. d. S. R. Gago, O. C. Zienkiewicz, and I. Babuška, “A posteriori error analysis and adaptive processes in the finite element method: Part I-Error Analysis,” Int. J. Num. Meth. Engrg. 19, 1593–1619 (1983). [CrossRef]

43.

W. Bangerth, R. Hartmann, and G. Kanschat, deal.II Differential Equations Analysis Library, Technical Reference (2004). http://www.dealii.org/.

44.

A. Thompson, “Development of a new optical imaging modality for detection of fluorescence enhanced disease,” Ph.D. thesis, Texas A & M University (2003).

45.

M. Huang and Q. Zhu, “Dual-mesh optical tomography reconstruction method with a depth correction that uses a priori ultrasound information,” Appl. Opt. 43(8), 1654–1662 (2004). [CrossRef]

46.

X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. 30(5), 861–869 (2003). [CrossRef]

OCIS Codes
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(170.5280) Medical optics and biotechnology : Photon migration

ToC Category:
Research Papers

History
Original Manuscript: September 7, 2004
Revised Manuscript: October 18, 2004
Published: November 1, 2004

Citation
Amit Joshi, Wolfgang Bangerth, and Eva Sevick-Muraca, "Adaptive finite element based tomography for fluorescence optical imaging in tissue," Opt. Express 12, 5402-5417 (2004)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-12-22-5402


Sort:  Journal  |  Reset  

References

  1. M. G. Pomper, �??Molecular Imaging: An Overview,�?? Acad. Radiol. 8, 1141�??1153 (2001). [CrossRef] [PubMed]
  2. M. A. O�??Leary, D. A. Boas, B. Chance, and A. Yodh, �??Reradiation and imaging of diffuse photon density waves using fluorescent inhomogeneities,�?? J. Luminescence 60, 281�??286 (1994). [CrossRef]
  3. J. Wu, Y. Wang, L. Perleman, I. Itzkan, R. R. Desai, and M. S. Feld, �??Time resolved multichannel imaging of fluorescent objects embedded in turbid media,�?? Opt. Lett. 20, 489�??491 (1995). [CrossRef] [PubMed]
  4. E. L. Hull, M. G. Nichols, and T. H. Foster, �??Localization of Luminescent Inhomogeneities in Turbid Media with Spatially Resolved Measurements of CW Diffuse Luminescence Emittance,�?? Appl. Opt. 37, 2755�??2765 (1998). [CrossRef]
  5. J. C. Schotland, �??Continuous wave diffusion imaging,�?? J. Opt. Soc. Am. A 14(275�??279) (1997).
  6. M. A. O�??Leary, D. A. Boas, B. Chance, and A. G. Yodh, �??Fluorescence lifetime imaging in turbid media,�?? Opt. Lett. 20, 426�??428 (1996). [CrossRef]
  7. E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, �??A submillimeter resolution fluorescence molecular imaging system for small animal imaging,�?? Med. Phys. 30, 901�??911 (2003). [CrossRef] [PubMed]
  8. V. Ntziachristos and R. Weissleder, �??Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,�?? Opt. Lett. 26(12), 893�??895 (2001). [CrossRef]
  9. V. Chernomordik, D. Hattery, I. Gannot, and A. H. Gandjbakhche, �??Inverse method 3-D reconstruction of localized in vivo fluorescence-application to Sjøgren syndrome,�?? IEEE J. Sel. Top. Quantum Electron. 54, 930�??935 (1999).
  10. H. Quan and Z. Guo, �??Fast 3-D Optical Imaging With Transient Fluorescence Signals,�?? Opt. Express 12(3), 449�??457 (2004). <a href=" http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-3-449">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-3-449</a> [CrossRef]
  11. D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, �??Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,�?? Appl. Opt. 36(10), 2260�??2272 (1997). [CrossRef]
  12. R. Roy and E. M. Sevick-Muraca, �??Truncated Newton�??s Optimization Schemes for Absorption and Fluorescence Optical Tomography: part(1) theory and formulation,�?? Opt. Express 4(10), 353�??371 (1999). <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353">http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-353</a> [CrossRef]
  13. H. Jiang, �??Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,�?? Appl. Opt. 37(22), 5337�??5343 (1998). [CrossRef]
  14. M. J. Eppstein, D. E. Dougherty, T. L. Troy, and E. M. Sevick-Muraca, �??Biomedical optical tomography using dynamic parametrization and Bayesian conditioning on photon migration measurements,�?? Appl. Opt. 38(10), 2138�??2150 (1998).
  15. M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, �??Three dimensional near infrared fluorescence tomography with Bayesian methodologies for image reconstruction from sparse and noisy data sets,�?? Proc. Nat. Acad. Sci. 99, 9619�??9624 (2002). [CrossRef] [PubMed]
  16. A. Godavarty, M. J. Eppstein, C. Zhang, S. Theru, A. B. Thompson, M. Gurfinkel, and E. M. Sevick-Muraca, �??Fluorescence-enhanced optical imaging in large tissue volumes using a gain-modulated ICCD camera,�?? Phys. Med. Biol. 48, 1701�??1720 (2003). [CrossRef] [PubMed]
  17. A. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. Boas, and R. P. Milane, �??Fluorescence Optical Diffusion Tomography,�?? Appl. Opt. 42(16), 3061�??3094 (2003).
  18. J. P. Houston, S. Ke,W.Wang, C. Li, and E. M. Sevick-Muraca, �??Optical and nuclear image quality analysis with invivo NIR fluorescence and conventional gamma images acquired using a dual labeled tumor targeting probe,�?? J. Biomed. Opt. (submitted) (2004).
  19. W. Bangerth, �??Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,�?? Ph.D. thesis, University of Heidelberg (2002).
  20. W. Bangerth, �??A framework for the adaptive finite element solution of large inverse problems. I. Basic techniques,�?? Tech. Rep. 04-39, ICES, University of Texas at Austin (2004).
  21. E. M. Sevick-Muraca, E. Kuwana, A. Godavarty, J. P. Houston, A. B. Thompson, and R. Roy, Near Infrared Fluorescence Imaging and Spectroscopy in Random Media and Tissues, chap. 33, Biomedical Photonics Handbook (CRC Press, 2003).
  22. A. Godavarty, D. J. Hawrysz, R. Roy, E. M. Sevick-Muraca, and M. J. Eppstein, �??Influence of the refractive index-mismatch at the boundaries measured in fluorescenceenhanced frequency-domain photon migration imaging,�?? Opt. Express 10(15), 650�??653 (2002).<a href=" http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-15-653"> http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-15-653</a>
  23. A. B. Thompson and E. M. Sevick-Muraca, �??NIR fluorescence contrast enhanced imaging with ICCD homodyne detection: measurement precision and accuracy,�?? J. Biomed. Opt. 8, 111�??120 (2002). [CrossRef]
  24. A. Joshi, W. Bangerth, and E. Sevick-Muraca, �??Adaptive finite element methods for fluorescence enhanced frequency domain optical tomography: Forward imaging problem,�?? in International Symposium on Biomedical Imaging, pp. 1103�??1106 (IEEE, 2004).
  25. R. A. Adams, Sobolev Spaces (Academic Press, 1975).
  26. A. N. Tikhonov and V. Y. Arsenin, eds., Solution of Ill-Posed Problems (Winston, Washington, DC, 1977).
  27. S. R. Arridge and M. Schweiger, �??A Gradient Based Optimization Scheme for Optical Tomography,�?? Opt. Express 2(6), 213�??225 (1998). <a href="http://www.opticsexpress.org/abstract.cfm?URI=OPEX-2-6-213"> http://www.opticsexpress.org/abstract.cfm?URI=OPEX-2-6-213</a> [CrossRef]
  28. L. Beilina, �??Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,�?? Ph.D. thesis, Chalmers University of Technology (2002).
  29. R. Becker, H. Kapp, and R. Rannacher, �??Adaptive Finite Element Methods for Optimal Control of Partial Differential Equations: Basic Concept,�?? SIAM J. Contr. Optim. 39, 113�??132 (2000). [CrossRef]
  30. R. Luce and S. Perez, �??Parameter identification for an elliptic partial differential equation with distributed noisy data,�?? Inverse Problems 15, 291�??307 (1999). [CrossRef]
  31. D. G. Luenberger, Optimization by Vector Space Methods (John Wiley, 1969).
  32. J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research (Springer, New York, 1999). [CrossRef]
  33. S. C. Brenner and R. L. Scott, The Mathematical Theory of Finite Elements (Springer, Berlin-Heidelberg-New York, 1994).
  34. R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques (Wiley/Teubner, New York, Stuttgart, 1996).
  35. W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations (Birkhäuser Verlag, 2003).
  36. M. Molinari, S. J. Cox, B. H. Blott, and G. J. Daniell, �??Adaptive Mesh Refinement techniques for Electrical Impedence Tomography,�?? Physiological Measurement 22, 91�??96 (2001). [CrossRef] [PubMed]
  37. M. Molinari, B. H. Blott, S. J. Cox, and G. J. Daniell, �??Optimal Imaging with Adaptive Mesh Refinement in Electrical Impedence Tomography,�?? Physiological Measurement 23, 121�??128 (2002). [CrossRef] [PubMed]
  38. R. Becker and R. Rannacher, �??An optimal control approach to error estimation and mesh adaptation in finite element methods,�?? Acta Numerica 10, 1�??102 (2001). [CrossRef]
  39. H. Ben Ameur, G. Chavent, and J. Jaffrè, �??Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities,�?? Inverse Problems 18, 775�??794 (2002). [CrossRef]
  40. A.-A. Grimstad, H. Krüger, T. Mannseth, G. Naevdal, and H. Urkedal, �??Adaptive selection of parameterization for reservoir history matching,�?? in ECMOR VIII: 8th European Conference on the Mathematics of Oil Recovery, Freiberg, Germany, pp. E�??46 (European Association of Geoscientists and Engineers (EAGE), 2002).
  41. R. Li, W. Liu, H. Ma, and T. Tang, �??Adaptive finite element approximation for distributed elliptic optimal control problems,�?? SIAM J. Control Optim. 41, 1321�??1349 (2002). [CrossRef]
  42. D. W. Kelly, J. P. d. S. R. Gago, O. C. Zienkiewicz, and I. Babuška, �??A posteriori error analysis and adaptive processes in the finite element method: Part I�??Error Analysis,�?? Int. J. Num. Meth. Engrg. 19, 1593�??1619 (1983). [CrossRef]
  43. W. Bangerth, R. Hartmann, and G. Kanschat, deal.II Differential Equations Analysis Library, Technical Reference (2004). <a href= "http://www.dealii.org/">http://www.dealii.org/</a>
  44. A. Thompson, �??Development of a new optical imaging modality for detection of fluorescence enhanced disease,�?? Ph.D. thesis, Texas A & M University (2003).
  45. M. Huang and Q. Zhu, �??Dual-mesh optical tomography reconstruction method with a depth correction that uses a priori ultrasound information,�?? Appl. Opt. 43(8), 1654�??1662 (2004). [CrossRef]
  46. X. Gu, Y. Xu, and H. Jiang, �??Mesh-based enhancement schemes in diffuse optical tomography,�?? Med. Phys. 30(5), 861�??869 (2003). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited