## Adaptive finite element based tomography for fluorescence optical imaging in tissue

Optics Express, Vol. 12, Issue 22, pp. 5402-5417 (2004)

http://dx.doi.org/10.1364/OPEX.12.005402

Acrobat PDF (1727 KB)

### 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 512*ml* cubical tissue phantom. We determine the achievable resolution for the area-illumination/area-detection reflectance measurement geometry by reconstructing two 0.4*cm* 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

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]

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]

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]

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]

*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.

## 2. Methodology

### 2.1. Formulation

*x*denotes the excitation light field and

*m*denotes the emission field;

*u,v*are the complex-valued photon fluence fields at excitation and emission wavelengths, respectively;(Note that part of the literature uses the symbols Φ

_{x},Φ

_{m}for these variables;

*u,v*are used to avoid overly complicated expressions with many indices in the next sections.)

*D*

_{x,m}are the photon diffusion coefficients;

*µ*

_{ax,mi}is the absorption coefficient due to endogenous chromophores;

*µ*

_{ax,mf}is the absorption coefficient due to exogenous fluorophores;

*µ′*

_{sx,m}is the reduced scattering coefficient;

*ω*is the modulation frequency;

*ϕ*is the quantum efficiency of the fluorophore; finally,

*τ*is the fluorophore lifetime associated with first order fluorescence decay kinetics. The fluorescence generation mechanism is detailed in [21]. These equations are complemented by Robin-type boundary conditions on the boundary

*∂*Ω of the domain Ω:

*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]

*µ*

_{axf}(

**r**) and/or τ(

**r**) from measurements of the fluences

*u,v*on the boundary.

*ζ*,ξ ∈

*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]). With this, the resulting variational equation reads:

*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

*D*

_{x,m}=

*D*

_{x,m}(

*q*),

*k*

_{x,m}=

*k*

_{x,m}(

*q*),

*β*

_{xm}=

*β*

_{xm}(

*q*).

*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] and

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

*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. 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]

*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, 20, 28, 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. R. Luce and S. Perez, “Parameter identification for an elliptic partial differential equation with distributed noisy data,” Inverse Problems **15**, 291–307 (1999). [CrossRef]

*λ*

^{ex}

*,λ*

^{em}are the Lagrange multipliers corresponding to the excitation and emission diffusion equation constraints, respectively. For simplicity, the abbreviation of

*x*={

*u,v,λ*

^{ex}

*,λ*

^{em}

*,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:

*L*

_{x}(

*x*)(

*y*) is the Frèchet differential [31] of

*L*(

*x*) and

*y*denotes possible test functions. Eq. (9) can be expanded for all components of

*x*:

*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 [

*λ*

^{ex}

*,λ*

^{em}]. Finally, equation (14) is the control equation.

*k*th iteration,

*δx*

_{k}={

*δu*

_{k}

*,δv*

_{k},

*,δq*

_{k}}, is determined from

*L*

_{xx}(

*x*

_{k}) is the Hessian matrix of second derivatives of

*L*at point

*x*

_{k}. These equations represent one condition for each variable in

*δx*

_{k}and when expanded read as follows [19, 31]:

*λ*

^{ex}

*,λ*

^{em}] are proportional to

*J*

_{v}(

*q;v*). As a consequence, all terms involving [

*λ*

^{ex}

*,λ*

^{em}] 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.

*α*

_{k}:

*α*

_{k}can be computed from one of several methods, such as the Goldstein-Armijo backtracking line search [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]

*q*are available. For example, background or maximal uptake concentrations may be known. This information should be used in the inverse problem to stabilize its solution, as well as to enforce physically reasonable solutions. Thus, we incorporate bounds

*q*

_{0}≤

*q*(

**r**)≤

*q*

_{1}using the scheme presented in [19, 20]. The method is a variation of the active set strategy, see [32

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

### 2.2. Discretization

*φ*

_{i}} as the basis functions for the state and adjoint variables

*u,v,λ*

^{ex}, and

*λ*

^{em}, and {

*χ*

_{i}} as the basis for the parameter

*q*. Piecewise linear, continuous shape functions on hexahedral meshes for {

*φ*

_{i}} and piecewise constant, discontinuous functions for {

*φ*

_{i}} are used. The discrete equations resulting from the Gauss-Newton modification of (15) are then represented by the following block system:

*δp*

_{k}=[

*δu*

_{k}

*,δdv*

_{k}]

^{T},

*δd*

_{k}=[

^{T}. The blocks of the matrix are defined as

*P*is the representation of the discrete forward diffusion model;

*i, j*in the preceding equations iterate over all degrees of freedom. The matrix

*C*is defined as

*C*

^{T}=[

*C*

_{1},

*C*

_{2}] with its components obtained by differentiating the semilinear form

*A*in Eq. (5) with respect to the parameter

*q*:

*R*+

*C*

^{T}

*P*

^{-T}

*MP*

^{-1}

*C*is the Schur complement matrix of the Gauss-Newton system. Under practical conditions it is symmetric positive definite [19] 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

*a priori*error estimates comparing the exact solution

*u*and the numerically computed finite element solution

*u*. For the simple example of the Laplace equation, an

_{h}*a priori*error estimate would be [33]:

*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.

*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-u*

_{h}‖ 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

*u*

_{h}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, 35] 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. 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]

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. 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]

*u,v*and

*λ*

^{ex}

*,λ*

^{em}, 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]

*∂*

_{n}

*u*

_{h}] 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.

*q*is refined by computing, for each cell, a discrete approximation to the gradient of

*q*weighted by the local mesh width [19]. 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, 28, 35].

### 2.4. Software Implementation

*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

*ε*. For the results reported in this paper we used

*k*

_{max}=40 and

*ε*=10

^{-15}.

*α*

_{k}is less than a prespecified minimum step length

*α*

_{min}, or (ii) the nonlinear residual,

*r*

_{k}=‖

*L*

_{x}(

*x*

_{k})(·)‖, of the optimality condition (9) has been sufficiently reduced on the current mesh

*m*, i.e.

*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.

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

## 3. Computational Experiments

^{3}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]

### 3.1. Single Fluorescent Target

*cm*diameter spherical target at a depth of 2.15

*cm*from the illumination surface at an off-center position (

*x*=2.15

*cm*,

*y*=3.15

*cm, z*=3.15

*cm*). 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]

*µ*

_{axi}=0.023

*cm*

^{-1}and

*µ*

_{ami}=0.0289

*cm*

^{-1}[44]. The absorption coefficient due to fluorophore at the excitation wavelength was set to

*µ*

_{axf}=0.5

*cm*

^{-1}in the target and 0.005

*cm*

^{-1}in the background. The emission wavelength absorption coefficient was

*µ*

_{amf}=0.0506

*cm*

^{-1}in the target and 0.00506

*cm*

^{-1}in the background. The lifetime of the fluorophore was taken to be

*τ*=0.56

*ns*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 785

*nm*and the emission data is collected at 830

*nm*. The reduced scattering coefficient was taken as

*µ′*

_{s}=9.84

*cm*

^{-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.

### 3.2. Two Fluorescent Targets

*cm*were placed at an edge to edge distance varying from 1

*cm*to 0.1

*cm*, 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.2

*cm*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

*hrs*of CPU time. The location of the target is reconstructed accurately. However, the magnitude 0.0176

*cm*

^{-1}of reconstructed

*µ*

_{axf}is lower than the actual value 0.5

*cm*

^{-1}, which is an artifact of the

*L*

_{2}regularization used

*β*=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 32

^{3}=32,768 unknowns. The advantages of adaptive refinement for the reduction of the total number of unknowns and the improvement of resolution are obvious.

### 4.2. Two Fluorescent Targets

*cm*. 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.1657

*cm*and greater, when the reconstructed targets begin to appear as one.

*cm*, 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.

## 5. Conclusions

*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]

*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]

*a posteriori*error estimates.

## Acknowledgments

## References and links

1. | M. G. Pomper, “Molecular Imaging: An Overview,” Acad. Radiol. |

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 |

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. |

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. |

5. | J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc. Am. A |

6. | M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. |

7. | E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. |

8. | V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. |

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. |

10. | H. Quan and Z. Guo, “Fast 3-D Optical Imaging With Transient Fluorescence Signals,” Opt. Express |

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. |

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 |

13. | H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,” Appl. Opt. |

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. |

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. |

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. |

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. |

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, |

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 |

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. |

24. | A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element methods for fluorescence enhanced frequency domain optical tomography: Forward imaging problem,” in |

25. | R. A. Adams, |

26. | A. N. Tikhonov and V. Y. Arsenin, eds., |

27. | S. R. Arridge and M. Schweiger, “A Gradient Based Optimization Scheme for Optical Tomography,” Opt. Express |

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. |

30. | R. Luce and S. Perez, “Parameter identification for an elliptic partial differential equation with distributed noisy data,” Inverse Problems |

31. | D. G. Luenberger, |

32. | J. Nocedal and S. J. Wright, |

33. | S. C. Brenner and R. L. Scott, |

34. | R. Verfürth, |

35. | W. Bangerth and R. Rannacher, |

36. | M. Molinari, S. J. Cox, B. H. Blott, and G. J. Daniell, “Adaptive Mesh Refinement techniques for Electrical Impedence Tomography,” Physiological Measurement |

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 |

38. | R. Becker and R. Rannacher, “An optimal control approach to error estimation and mesh adaptation in finite element methods,” Acta Numerica |

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 |

40. | A.-A. Grimstad, H. Krüger, T. Mannseth, G. Nævdal, and H. Urkedal, “Adaptive selection of parameterization for reservoir history matching,” in |

41. | R. Li, W. Liu, H. Ma, and T. Tang, “Adaptive finite element approximation for distributed elliptic optimal control problems,” SIAM J. Control Optim. |

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. |

43. | W. Bangerth, R. Hartmann, and G. Kanschat, deal.II |

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. |

46. | X. Gu, Y. Xu, and H. Jiang, “Mesh-based enhancement schemes in diffuse optical tomography,” Med. Phys. |

**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

- M. G. Pomper, �??Molecular Imaging: An Overview,�?? Acad. Radiol. 8, 1141�??1153 (2001). [CrossRef] [PubMed]
- 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]
- 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]
- 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]
- J. C. Schotland, �??Continuous wave diffusion imaging,�?? J. Opt. Soc. Am. A 14(275�??279) (1997).
- 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]
- 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]
- 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]
- 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).
- 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]
- 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]
- 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]
- H. Jiang, �??Frequency-domain fluorescent diffusion tomography: a finite element based algorithm and simulations,�?? Appl. Opt. 37(22), 5337�??5343 (1998). [CrossRef]
- 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).
- 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]
- 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]
- 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).
- 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).
- W. Bangerth, �??Adaptive Finite Element Methods for the Identification of Distributed Coefficients in Partial Differential Equations,�?? Ph.D. thesis, University of Heidelberg (2002).
- 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).
- 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).
- 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>
- 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]
- 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).
- R. A. Adams, Sobolev Spaces (Academic Press, 1975).
- A. N. Tikhonov and V. Y. Arsenin, eds., Solution of Ill-Posed Problems (Winston, Washington, DC, 1977).
- 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]
- L. Beilina, �??Adaptive Hybrid FEM/FDM Methods for Inverse Scattering Problems,�?? Ph.D. thesis, Chalmers University of Technology (2002).
- 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]
- R. Luce and S. Perez, �??Parameter identification for an elliptic partial differential equation with distributed noisy data,�?? Inverse Problems 15, 291�??307 (1999). [CrossRef]
- D. G. Luenberger, Optimization by Vector Space Methods (John Wiley, 1969).
- J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research (Springer, New York, 1999). [CrossRef]
- S. C. Brenner and R. L. Scott, The Mathematical Theory of Finite Elements (Springer, Berlin-Heidelberg-New York, 1994).
- R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh Refinement Techniques (Wiley/Teubner, New York, Stuttgart, 1996).
- W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations (Birkhäuser Verlag, 2003).
- 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]
- 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]
- 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]
- 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]
- 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).
- 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]
- 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]
- 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>
- A. Thompson, �??Development of a new optical imaging modality for detection of fluorescence enhanced disease,�?? Ph.D. thesis, Texas A & M University (2003).
- 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]
- 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.