## Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue

Optics Express, Vol. 15, Issue 11, pp. 6955-6975 (2007)

http://dx.doi.org/10.1364/OE.15.006955

Acrobat PDF (1364 KB)

### Abstract

We have developed fluorescence enhanced optical tomography based upon fully adaptive finite element method (FEM) using tetrahedral dual-meshing wherein one of the two meshes discretizes the forward variables and the other discretizes the unknown parameters to be estimated. We used the 8-subtetrahedron subdivision scheme to create the nested dual-mesh in which each are independently refined. However, two tetrahedrons from the two different meshes pose an intersection problem that needs to be resolved in order to find the common regions that the forward variables (the fluorescent diffuse photon fluence fields) and the parameter estimates (the fluorescent absorption coefficients) can be mutually assigned. Using an efficient intersection algorithm in the nested tetrahedral environments previously developed by the authors, we demonstrate fully adaptive tomography using *a posteriori* error estimates. Performing the iterative reconstructions using the simulated boundary measurement data, we demonstrate that small fluorescent targets embedded in the breast simulating phantom in point illumination/detection geometry can be resolved at reasonable computational cost.

© 2007 Optical Society of America

## 1. Introduction

1. K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, and J. M. Sullivan, Jr., “A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans. Med. Imaging **14**, 504–514 (1995). [CrossRef] [PubMed]

1. K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, and J. M. Sullivan, Jr., “A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans. Med. Imaging **14**, 504–514 (1995). [CrossRef] [PubMed]

*a priori*information on the spatial distribution of the unknown parameters, the discretization of forward and inverse meshes is arbitrarily and nonoptimally chosen. In order to alleviate this difficulty, adaptive mesh refinement or coarsening strategy both for forward and inverse meshes is required.

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

7. A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Non-contact fluorescence optical tomography with scanning patterned illumination,” Opt. Express **14**, 6516–6534 (2006). [CrossRef] [PubMed]

## 2. Algorithms

### 2.1 Nested tetrahedral dual-meshing

9. A. Liu and B. Joe, “Quality local refinement of tetrahedral meshes based on 8-subtetrahedron subdivision,” Math. Comput. **65**, 1183–1200 (1996). [CrossRef]

_{8}, SUB

_{2}, SUB

_{4a}, and SUB

_{4b}, respectively. The regular subdivision SUB

_{8}creates eight subtetrahedrons and the resulting subtetrahedrons are called regular and labeled as

*S*

_{8}. SUB

_{2}, SUB

_{4a}, or SUB

_{4b}are applied only to the boundary of the regularly refined regions to ensure mesh conformity and are performed in the nonregular fashion. The subtetrahedrons produced by SUB

_{2}, SUB

_{4a}, and SUB

_{4b}are labeled as

*S*

_{2},

*S*

_{4a}, and

*S*

_{4b}, respectively and called nonregular. The nonregular tetrahedrons can exist only at leaf levels of the tetrahedron tree built upon parent-progeny relations. All tetrahedrons in the initial mesh become roots of the subdivision level 0 for building their own trees which are also treated as regular and labeled as S

_{8}; for details, see [9

9. A. Liu and B. Joe, “Quality local refinement of tetrahedral meshes based on 8-subtetrahedron subdivision,” Math. Comput. **65**, 1183–1200 (1996). [CrossRef]

_{8}. In the dual-mesh scheme, the forward variables and the unknown parameters are decoupled by discretizing the forward variables on the forward mesh while discretizing the parameters on the inverse mesh. In the fully adaptive scheme, the two meshes are independently refined or coarsened according to 8-subtetrahedron subdivision scheme, therefore posing intersection problems between tetrahedrons belonging to different meshes. The intersection must be resolved between the forward domain and the inverse domain tetrahedral elements in order to obtain the common region where the finite element assembly will be performed.

### 2.2 Intersection algorithm

*partnership*,

*parent-progeny*relations, and

*weights*that are linear finite element basis function values [8]. The unique feature of FINT lies in that it uses the weights for resolving intersection and for obtaining the inter-element weight maps. This feature leads to a robust implementation of intersection computations which includes the ability to appropriately approximate the intersection between tetrahedrons having facets on the curved boundary. As the final output, FINT provides the common region decomposed into disjoint tetrahedron pieces including all necessary information for the finite element assembly.

*partnership*. Hereby, each tetrahedron on the forward mesh has its partner on the inverse mesh for the initial dual-mesh. The partnership can be built only for the

*S*

_{8}-tetrahedrons that are regular because only SUB

_{8}provides the subdivision in the regular pattern. Whenever any tetrahedron on one of the two meshes is subdivided by SUB

_{8}, we can efficiently search and assign partners (if they exist on the other mesh) using the parent-progeny relation together with the partnership. The searching procedure is as follows: (i) for any newly created

*S*

_{8}-subtetrahedron

**T**of the subdivision level

*L*, the

*closest*ancestor that has its partner is searched using the parent-progeny relation. (ii) If such ancestor is not found,

**T**has no partner. Otherwise, we move to the other mesh using the bridge provided by the ancestor’s partner and search the offspring of the ancestor’s partner by tree traversal down to the subdivision level

*L*. (iii) If

**T**’s twin

**T**´ is found, the partnership (that is, the bridge) is built mutually between

**T**and

**T**´. Based upon the above rule, the partnership is updated whenever the nested dual-meshing environment changes.

*nonregular*, that is, of the type

*S*

_{2},

*S*

_{4a}, or

*S*

_{4b}. Because the tetrahedrons created by the subdivision are nested, any intersecting configurations are confined inside the volume occupied by some regular tetrahedron which has been named the

*container*that can be found by using the partnership and the parent-progeny relations. Once the container is found,

*weights*play a role in determining and resolving intersection in the recursive way, providing disjoint tetrahedron pieces completely filling the container. Each disjoint tetrahedron piece contains the complete information required for finite element assembly such as inter-element weight maps, nodal indices, etc. For details, see [8]. Figure 2 illustrates the disjoint tetrahedron pieces obtained by FINT, where the nontrivial intersections have been resolved for the given nested tetrahedrons embedded inside three different container tetrahedrons with nonregular children on the other mesh.

### 2.3 Reconstruction algorithm

#### 2.3.1 Overview

*μ*, is the parameter to reconstruct and the excitation and fluorescent emission fluence fields

_{axf}**Φ**

_{x}and

**Φ**

_{m}are the forward variables. Therefore we distribute the variables

**Φ**

_{x}and

**Φ**

_{m}to the forward mesh and the unknown parameter

*μ*to the inverse mesh. Because optical tomography is highly nonlinear, the spatial distribution of the unknown

_{axf}*μ*is reconstructed iteratively. In this study, we use the optimization based iterative reconstruction scheme and minimize the cost function given as the model-mismatch in the least squares sense:

_{axf}*N*is the number of sources;

_{S}*N*is the number of detectors; and

_{D}*N*is the number of nodes used for the discretization of

_{μ}*μ*. The term

_{axf}*y*denotes the fluorescent emission field measured at

_{ji}*j*th detector for

*i*th source illumination;

**Φ**

*and*

^{i}_{x}**Φ**

*are the nodal solution vectors of the excitation and the emission fields, respectively, predicted by the model for the*

^{i}_{m}*i*th source illumination;

**P**

_{j}is the projection vector whose inner product with

**Φ**

*maps*

^{i}_{m}**Φ**

*onto the measured fluorescence at they*

^{i}_{m}*j*th detector location for the

*i*th source illumination; and

**p**is the array of nodal values of

*p*=

*μ*.

_{axf}*D*=1/[3(

_{x,m}*μ*+

_{axi,ami}*μ*+

_{axf,amf}*μ´*)],

_{sx,sm}*k*=i

_{x,m}*ω*/

*c*+

*μ*(

_{axi,ami}**r**)+

*μ*(

_{axf,amf}**r**), and

*β*=

_{xm}*qμ*/[1-

_{axf}*iωτ*(

**r**)]. An index

*x*denotes the excitation field and

*m*denotes the emission field;

*D*are photon diffusion coefficients;

_{x,m}*μ*is the absorption coefficient due to endogenous chromophores;

_{axi,ami}*μ*is the absorption coefficient due to exogenous fluorophores;

_{axf,amf}*μ´*is the reduced scattering coefficient;

_{sx,sm}*ω*=2

*πf*where

*f*is the modulation frequency;

*q*is the quantum efficiency of the fluorophore; and finally,

*τ*is the lifetime of fluorescence. Eq. (3) is solved under the Robin-type boundary conditions:

**n**is the unit outward surface normal vector,

*γ*is a constant depending only upon the refractive index mismatch on the boundary, and

*Q*(

**r**) is the excitation boundary source distribution.

*μ*in the inverse mesh is passed onto the forward mesh and finite element assembly is performed; (ii)

_{axf}**Φ**

_{x}is solved; (iii)

**Φ**

_{m}is solved; (iv) new estimate of

*μ*is obtained by reducing the cost function; and (v) finally forward/inverse meshes are adapted according to some criterion. For optimization methods, we use the trust-region Gauss-Newton method in the active set method framework [12

_{axf}12. J. Nocedal and S. J. Wright, *Numerical Optimization* (Springer-Verlag, New York, 1999). [CrossRef]

#### 2.3.2 Finite element assembly

**T**and

**T**´ from the forward mesh

*M*and the inverse mesh

*M*´, respectively. Suppose that {

*φ*(

_{n}**r**)∣

*n*=0, 1, 2,⋯} and {

*ψ*(

_{n}**r**)∣

*n*=0, 1, 2,⋯} are linear finite element basis functions in

*M*and

*M*´, respectively, where

**r**is a real space coordinate vector. Each vertex a of a tetrahedron piece

**Δ**in the intersection of

**T**and

**T**´ is associated with two weight arrays provided by FINT:

*i*,

*j*,

*k*,

*l*} and {

*i*´,

*J*´,

*k*´,

*I*´} are nodes of

**T**and

**T**´, respectively, and

**r**

_{α}is the coordinate vector of

*α*. For each vertex

*α*of the piece

*Δ*, a

*virtual*linear basis function

*π*(

_{α}**r**) is introduced.

*π*(

_{α}**r**) subtends

**Δ**and behaves like a usual basis function with

*p*,

*q*,

*r*,

*s*} and

*V*

_{Δ}are the vertices and the volume of

**Δ**, respectively, that is provided by FINT; that is,

*π*is the volume coordinates with respect to

_{α}**Δ**. The linear property of (

*ψ*(

_{n}**r**),

*ψ*(

_{n}**r**), and

*π*(

_{α}**r**) gives rise to the relationships,

*ψ*(

_{n}**r**)=∑

_{α}

*ψ*(

_{n}**r**

_{α})

*π*(

_{α}**r**) and

*ψ*(

_{n}**r**)=∑

_{α}

*ψ*(

_{n}**r**

_{α})

*π*(

_{α}**r**), enabling two 4×4 transformation matrices

**U**and

**V**to be introduced:

**U**and

**V**perform the change of basis from the virtual basis

*π*to the actual basis

*φ*and

*ψ*, respectively, and inside

**Δ**result in:

*φ*,

*ψ*, or both are involved in are performed over

**Δ**for finite element assembly by treating

*π*as the usual linear finite element basis function.

*μ*and that

_{axf}*μ*has a linear relationship with

_{amf}*μ*given by

_{axf}*μ*=

_{amf}*ζ*, only

_{axf}*D*is affected. Therefore we need to discretize both

_{x,m}*D*and

_{x,m}*μ*on the inverse mesh

_{axf}*M*´ giving rise to:

**Φ**

_{x}and

**Φ**

_{m}are discretized on the forward mesh into:

**Z**,

**Θ**,

**E**, and

**T**are given by:

**Γ**and

**Q**can be assembled directly over the boundary triangle elements in the forward mesh

*M*. On the other hand, the finite element assembly of the volume terms

**Z**,

**Θ**, and

**E**are performed over the tetrahedron piece provided by FINT, that is:

**Z**

^{Δ},

**Θ**

^{Δ}, and

**E**

^{Δ}are given by:

**E**can be assembled directly over the tetrahedrons in the forward mesh

*M*since only {

*φ*} is involved. Therefore,

_{i}**K**

_{x,m}and

**B**

_{x,m}are assembled by:

*k*runs through nodal indices of a tetrahedron

**T**´ on the inverse mesh

*M*´ providing

**Δ**by intersection with a tetrahedron

**T**on the forward mesh

*M*.

*D*and

_{x,m}*μ*inside

_{axf}**T**´ by taking the average of the values at the four nodes of

**T**´, that is:

#### 2.3.3 Jacobian matrix computation

*i*th source can be derived from Eq. (12) and the adjoint theorem. With

*p*=

*μ*the derivative of Eq. (12) with respect to

_{axf}*p*at the

_{k}*k*th node in the inverse mesh can be arranged into:

**Ψ**

_{x}and

**Ψ**

_{m}and association after making inner product with Eq. (31) gives rise to:

**K**

_{x,m}and

**B**

_{xm}are symmetric. The sensitivity of the measured fluorescence at the

*j*th detector for the

*i*th source illumination with respect to the variation of

*p*reads:

_{k}**Ψ**

_{x}and

**Ψ**

_{m}satisfy:

**Ψ**

*and*

^{i}_{x}**Ψ**

*are the solutions of Eq. (34). For the point detection geometry, when the*

^{j}_{m}*j*th detector location is

**r**

_{j}, the

*k*th entry

*P*of the projection vector

^{k}_{j}**P**

_{j}becomes:

**J**is assembled using Eq. (29) and Eq. (30). Inserting Eq. (29) and Eq. (30) into Eq. (37) leads to:

**J**

^{Δ}is given by:

*γ*runs through nodal indices of the tetrahedral element

**T**´ on the inverse mesh

*M*´ and {

*α*,

*β*} runs through nodal indices of the tetrahedral element

**T**on the forward mesh

*M*for a pair {

**T**,

**T**´} providing

**Δ**by intersection. In this manner, the assembly is performed over the tetrahedron pieces provided by FINT for

**J**as well.

#### 2.3.4 Adaptation strategy

**Φ**

_{x},

**Φ**

_{m}} and the inverse parameter of

*μ*, are decoupled using the separate meshes, the two meshes can be independently refined. The forward and the inverse meshes are refined using

_{axf}*a posteriori*error estimates upon the spatial variation of {

**Φ**

_{x},

**Φ**

_{m}} and

*μ*, respectively. For

_{axf}*a posteriori*error estimates, we use the Kelly’s flux jump criterion [10

10. D. W. Kelly, J. P. De S. R. Gago, O. C. Zienkiewicz, and I. Babuska, “A *posteriori* error analysis and adaptive processes in the finite element method: Part I-Error analysis,” Int. J. Numer. Meth. Engng. **19**, 1593–1619 (1983). [CrossRef]

*f*for a tetrahedral element

**T**is estimated by:

*∂f/∂n*∣

_{±}is the jump of the normal derivative of

*f*across the triangle facets of the tetrahedral element

**T**and

*h*is the diameter of

**T**. Note that Eq. (40) is computed directly over the tetrahedral element

**T**and the tetrahedron pieces provided by FINT play no role.

*forward element*), requires Eq. (40) to be accumulated for all source/detector configurations with properly weighted combination of error estimates in different fields. We reduce to the error estimate

*e*

^{T}_{Φ}for the forward element

**T**as:

*w*,

_{x,s}*w*,

_{m,s}*w*, and

_{xa,d}*w*are normalizing factors chosen to be squares of the inverse of maximum nodal values in the amplitudes in

_{ma,d}**Φ**

*,*

^{s}_{x}**Φ**

*,*

^{s}_{m}**Φ**

*, and*

^{d}_{x}**Φ**

*, respectively. For a tetrahedral element*

^{d}_{m}**T**in the inverse mesh (

*inverse element*), we simply use Eq. (40) with

*f*=

*μ*for the element error estimate

_{axf}*e*as:

^{T}_{μ}**T**is chosen for refinement in the corresponding meshes if:

*c*,

*d*} are endpoints of an edge of

**T**’s parent which has a mid-point

*m*on it and {

*p*,

_{max}*p*} are the maximum and the minimum values of the parameter

_{min}*p*obtained at the most recent update. The first refinement to the initial forward/inverse meshes is always triggered. For the next refinements, only when there is an inverse element of the nonzero subdivision level satisfying the condition:

**T**is of the nonzero subdivision level, it is chosen for refinement when it fulfills both Eq. (43) and Eq. (45). Note that the condition Eq. (45) is introduced to trigger refinements only when the sufficient amount of spatial variations in the parameters has been accumulated.

**Z**

^{Δ},

**E**

^{Δ}} are updated by running FINT after refinements for both meshes are finished. Although the inverse mesh can be refined independently of the forward mesh, we strictly avoid generations of

*floating*nodes in the inverse mesh: If no node is present on the forward mesh at the location where any new node supposed to be introduced by edge splitting due to the subdivision of an inverse element

**T**, then the subdivision of

**T**is rejected. Hereby, we build the nested dual-mesh satisfying the following rule for any leaf inverse element

**T**: (i) if

**T**is regular (of the type

*S*

_{8}),

**T**’s partner must exist in the forward mesh, and otherwise (ii) if

**T**is nonregular (of the type

*S*

_{2},

*S*

_{4a}, or

*S*

_{4b}), the partner of

**T**’s parent having 8 child tetrahedrons (of the type

*S*

_{8}) must be present in the forward mesh.

#### 2.3.5 Optimization strategy

*precede*the minimization of the object function. The procedure is as follows: (i) the

*reduced*Jacobian matrix is assembled over the tetrahedron pieces provided by FINT that are associated with only

*free*parameters

**x**=

**p**

_{f}. (ii) The

*unconstrained*function minimizing Gauss-Newton direction

**d**is found by the Steihaug-CG inside the

*spherical*trust region with the radius

*D*[12

12. J. Nocedal and S. J. Wright, *Numerical Optimization* (Springer-Verlag, New York, 1999). [CrossRef]

**d**∥

_{2}<

*ε*for the sufficiently small number

_{c}*ε*, then the whole optimization process is considered to be converged and thus terminated. (iii) Next, the step length is controlled and the sets of the

_{c}*free*and the

*bound*parameters are updated; the maximum step length

*λ*is initialized to 1. For a free index

_{max}*i*satisfying

*x*+

_{i}*d*<

_{i}*ε*, if

*x*<

_{i}*ε*and

*d*<0 we set

_{i}*d*=0 and delete

_{i}*i*from the set of

*free*indices and add

*i*to the set of the

*bound*indices and otherwise compute the maximum step length

*λ*satisfying Eq. (2) and update

_{i}*λ*=min{

_{max}*λ*,

_{i}*λ*}. (iv) Backtracking line search is performed for the interval (0,

_{max}*λ*] and the ratio

_{max}*ρ*of the actual to the predicted reduction in the object function is computed. (v) If

*ρ*<1/4, then the trust region radius

*D*is reduced to

*D*=min{

*D*/4,

*D*} and

_{min}**x**from the previous iteration remains unchanged. Otherwise if

*ρ*>3/4,

**x**is updated and the trust region radius

*D*is increased to

*D*=max{2

*D*,

*D*}. Else,

_{max}**x**is updated and

*D*remains unchanged.

*ε*

^{T}_{Φ}satisfying Eq. (43) and the inverse elements satisfying Eq. (43) in

*ε*or both Eq. (43) and Eq. (45) under the restrictions already mentioned are chosen for refinement. For the inverse mesh, the parameter values at the new mid-points introduced by edge splitting are assigned by linear interpolation from the existing values at the endpoints of the parent edge. If the parameters for the new nodes are nonzero, they are added to the set of free parameters and otherwise added to the set of bound parameters. Figure 3 illustrates a flow diagram of the algorithm and decision points.

^{T}_{μ}*free*status and then the iterative reconstructions begin.

## 3. Results

*μ*=0.02483 cm

_{axi}^{-1},

*μ´*=10.8792 cm

_{sx}^{-1},

*μ*=0.0322 cm

_{ami}^{-1},

*μ´*=9.8241 cm

_{sm}^{-1}, and

*ζ*=

*μ*/

_{amf}*μ*=0.1692. The optical parameters specific to fluorescent targets are

_{axf}*μ*=0.1 cm

_{axf}^{-1},

*q*=0.016, and

*τ*=0.56 ns. The refractive index of the medium is

*n*=1.33 and

*f*=100 MHz amplitude modulation is assumed. The locations of the targets are illustrated in Fig. 5. We considered two different configurations of targets with varying the gap lengths as shown in Fig. 5. The single target is located at (2.2 cm, 0.0, 2.2 cm). The gap between the two targets considered is 1 cm, 5 mm, or 2.5 mm. The two target locations are summarized in Table. 1. To obtain the simulated measurement data for each configuration, we employed the dual-mesh FEM illustrated in Figs. 6(a) and 6(b). The input coarse mesh shown in Fig. 6(b) has 1,035 nodes and is used as the initial mesh discretizing

*μ*. The mesh shown in Fig. 6(a) was obtained by refining the boundary of the breast part once and is used as the initial mesh discretizing the fields Φ

_{axf}_{x}and Φ

_{m}. To obtain the accurate fluence fields and the target locations with the specified diameter and gaps for

*μ*, the dual-mesh was refined to the fifth subdivision level using

_{axf}*η*

_{Φ,μ}=0.1 and 7 refinement iterations.

*η*

_{Φ,μ}=0.5 and

*θ*=0.25. The maximum subdivision level for the nested dual-meshes was set as 4, enabling the discretization down to the theoretical diffusion limit. We decided whether to trigger refinements for the dual-mesh by checking the condition in Eq. (45) regularly after five successive iterations. For the trust region, we set

*D*=0.001 and

_{max}*Δ*=0.1. The convergence tolerance for the 2-norm of the Gauss-Newton direction

_{max}**d**is chosen as

*ε*=10

_{c}^{-16}and the tolerance for the bound constraint is chosen as

*ε*=10

^{-4}. We assume no

*a priori*information is available on the sizes of the fluorescent targets. Therefore, totally blind reconstructions were performed in this work, and attempts for localizing and characterizing the unknown targets were made to the length scale of the diffusion limit with progressive refinements relying only on

*a posteriori*basis using Eq. (43). We allowed a maximum 300 iterations to estimate the number of iterations necessary for the sufficient target recovery in different target configurations. All computations were done using the PC with 3.6 GHz Pentium 4 processor and 3.25 GB RAM.

### 3.1 Single target reconstruction

*μ*are 0.0164, 0.0240, 0.0377, 0.0706, and 0.2028 cm

_{axf}^{-1}for IT=30, 60, 100, 150, and 261, respectively. The iterations were terminated at IT=261 since the condition ∥

**d**∥

_{2}<

*ε*was encountered. 30 iterations suffice to locate the target. Quantitatively, 99.7 % of the target

_{c}*μ*value was recovered at IT=172 with 14319/1218 forward/inverse nodes. The accumulated computing times are 0.067, 0.276, 0.906, 1.712, and 2.953 hours for IT=30, 60, 100, 150, and 261, respectively.

_{axf}### 3.2 Two target reconstruction with 1 cm gap

*μ*distribution map are 0.0107, 0.0248, 0.0345, 0.0552, and 0.1790 cm

_{axf}^{-1}for IT=30, 60, 100, 150, and 300, respectively. The two target separation is clear at IT=30. Localization to the correct positions occurs approximately in 150 iterations. Quantitatively, 99.9 % of the target

*μ*value was recovered at IT=201 with 21264/1482 forward/inverse nodes. The accumulated computing times are 0.0857, 0.384, 1.527, 3.474, and 8.268 hours for IT=30, 60, 100, 150, and 300, respectively.

_{axf}### 3.3 Two target reconstruction with 5 mm gap

*μ*distribution map are 0.0129, 0.0204, 0.0320, 0.0553, and 0.1499 cm

_{axf}^{-1}for IT=30, 60, 100, 150, and 300, respectively. The two targets are not separated at IT=30 and seems to be separated at IT=60. Although the separation becomes clear at IT=100, the locations are incorrect. Localization to the correct positions occurs approximately in 150 iterations. Quantitatively, 99.9 % of the target

*μ*value was recovered at IT=230 with 20013/1538 forward/inverse nodes. The accumulated computing times are 0.076, 0.273, 0.853, 2.293, and 6.454 hours for IT=30, 60, 100, 150, and 300, respectively.

_{axf}### 3.4 Two targets with 2.5 mm gap

*μ*distribution map are 0.0241, 0.0356, 0.0431, 0.0689, and 0.2189 cm

_{axf}^{-1}for IT=30, 70, 100, 150, and 300, respectively. The two targets are not separated at IT=30. Although the separation is visible at IT=70 and becomes clear at IT=100, the locations are incorrect. Localization to the correct positions occurs approximately in 150 iterations. Quantitatively, 99.4 % of the target

*μ*value was recovered at IT=197 with 15998/1543 forward/inverse nodes. The accumulated computing times are 0.065, 0.312, 0.685, 1.736, and 4.037 hours for IT=30, 70, 100, 150, and 300, respectively.

_{axf}## 4. Discussions

*μ*distribution map. The number of nodes in the forward mesh reaches around 15,000 to 20,000 depending upon the target configurations. The number of nodes in the inverse mesh reaches around 1,200 to 1,600. Therefore, allocating the Jacobian matrix is not expensive and the overall computational cost is dominated by the forward solution cost.

_{axf}*μ*values with the reduced full-width-half-maximum (FWHM). These artifacts typically trigger superfluous additional refinements. Since termination by our convergence is rare, we require a trade-off between localization and the number of the iterations. We found that the forward meshes and probably the inverse meshes at IT=100 illustrated in Fig. 7 to Fig. 10 are sufficiently dense in order to obtain accurate forward solutions and inverse parameter resolutions, respectively. In all cases considered in this work, all targets are seen to be localized near the correct location in 1 hour of the total computing time which corresponds to 80 to 120 iterations as shown in Fig. 11. Furthermore, since the measured or simulated data contain noise, the increased accuracy of the forward solution for reconstruction purposes is probably not needed at the length scale of the theoretical diffusion limit.

_{axf}*blind*reconstructions in the sense that no

*a priori*information is assumed on the location and the size of the target. Actually, we have assumed that the intrinsic optical properties of the domain at the emission and excitation wavelengths are known, while reconstructing for absorption due to the fluorophore. Therefore, one might ask how this assumption will affect our results. Fluorescence optical tomography for determination of absorption due to fluorophore is robust against small perturbations in values of intrinsic optical properties. The reason can be intuitively understood by interpreting the Jacobian sensitivity matrix. Sensitivity of a given measurement pair to

*μ*depends upon the product of forward and adjoint diffusion equation solutions corresponding to that source detector pair. In the asymptotic limit, the effect of small perturbations in intrinsic optical properties on the diffusion solution gets squared in the Jacobian, and thus does not affect the

_{axf}*μ*reconstruction to a substantial degree. Sahu and coworkers [14

_{axf}14. A. K. Sahu, R. Roy, A. Joshi, and E. M. Sevick-Muraca, “Evaluation of anatomical structure and non-uniform distribution of imaging agent in near-infrared fluorescence-enhanced optical tomography,” Opt. Express **13**, 10182–10199 (2005) [CrossRef] [PubMed]

## 5. Summary and Conclusions

*a priori*information on the length scale for the target sizes and their separation distances in the numerical experiments performed in this work. In these

*blind*situations, the fully adaptive tomography scheme might be the right approach for reconstruction. Using the different meshes to obtain the simulated data and the inverse solutions, we have strictly avoided the inverse crime as well. We have demonstrated that the proposed fully adaptive FEM technique using tetrahedral dual-meshes is able to resolve small fluorescent targets without the large memory overhead.

## Acknowledgments

## References and links

1. | K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, and J. M. Sullivan, Jr., “A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans. Med. Imaging |

2. | H. Jiang, K. Paulsen, U. Osterberg, and M. Patterson, “Frequency-domain near-infrared photo diffusion imaging: initial evaluation in multi-target tissue-like phantoms,” Med. Phys. |

3. | M. Schweiger and S. R. Arridge, “Optical tomographic reconstruction in a complex head model using a priori region boundary information,” Phys. Med. Biol. |

4. | H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, “Three dimensional optical tomography: resolution in small-object imaging,” Appl. Opt. |

5. | M. Huang and Q. Zhu, “Dual-mesh tomography reconstruction method with a depth correction that uses a priori ultrasound information,” Appl. Opt. |

6. | A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence enhanced optical imaging in tissue,” Opt. Express |

7. | A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Non-contact fluorescence optical tomography with scanning patterned illumination,” Opt. Express |

8. | J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, “Fast intersections on nested tetrahedrons (FINT): An algorithm for adaptive finite element based distributed parameter estimation,” ( |

9. | A. Liu and B. Joe, “Quality local refinement of tetrahedral meshes based on 8-subtetrahedron subdivision,” Math. Comput. |

10. | D. W. Kelly, J. P. De S. R. Gago, O. C. Zienkiewicz, and I. Babuska, “A |

11. | 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,” in |

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

13. | M. J. Epstein, F. Fedele, J. Laible, C. Zhang, A. Godavarty, and E. M. Sevick-Muraca, “A comparison of exact and approximate adjoint sensitivities in fluorescence tomography,” IEEE Trans. Med. Imaging , |

14. | A. K. Sahu, R. Roy, A. Joshi, and E. M. Sevick-Muraca, “Evaluation of anatomical structure and non-uniform distribution of imaging agent in near-infrared fluorescence-enhanced optical tomography,” Opt. Express |

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

Medical Optics and Biotechnology

**History**

Original Manuscript: March 8, 2007

Revised Manuscript: May 17, 2007

Manuscript Accepted: May 17, 2007

Published: May 22, 2007

**Virtual Issues**

Vol. 2, Iss. 6 *Virtual Journal for Biomedical Optics*

**Citation**

Jae Hoon Lee, Amit Joshi, and Eva M. Sevick-Muraca, "Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue," Opt. Express **15**, 6955-6975 (2007)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-11-6955

Sort: Year | Journal | Reset

### References

- K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, J. M. Sullivan, Jr., "A dual mesh scheme for finite element based reconstruction algorithms," IEEE Trans. Med. Imaging 14, 504-514 (1995). [CrossRef] [PubMed]
- H. Jiang, K. Paulsen, U. Osterberg, and M. Patterson, "Frequency-domain near-infrared photo diffusion imaging: initial evaluation in multi-target tissue-like phantoms," Med. Phys. 25, 183-193 (1998) [CrossRef] [PubMed]
- M. Schweiger and S. R. Arridge, "Optical tomographic reconstruction in a complex head model using a priori region boundary information," Phys. Med. Biol. 44, 2703-2721 (1999). [CrossRef] [PubMed]
- H. Dehghani, B. W. Pogue, J. Shudong, B. Brooksby, and K. D. Paulsen, "Three dimensional optical tomography: resolution in small-object imaging," Appl. Opt. 42, 3117-3129 (2003). [CrossRef] [PubMed]
- M. Huang and Q. Zhu, "Dual-mesh tomography reconstruction method with a depth correction that uses a priori ultrasound information," Appl. Opt. 43, 1654-1662 (2004). [CrossRef] [PubMed]
- A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, "Adaptive finite element based tomography for fluorescence enhanced optical imaging in tissue," Opt. Express 12, 5402-5417 (2004). [CrossRef] [PubMed]
- A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, "Non-contact fluorescence optical tomography with scanning patterned illumination," Opt. Express 14, 6516-6534 (2006). [CrossRef] [PubMed]
- J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, "Fast intersections on nested tetrahedrons (FINT): An algorithm for adaptive finite element based distributed parameter estimation," (submitted).
- A. Liu and B. Joe, "Quality local refinement of tetrahedral meshes based on 8-subtetrahedron subdivision," Math. Comput. 65, 1183-1200 (1996). [CrossRef]
- D. W. Kelly, J. P. De S. R. Gago, O. C. Zienkiewicz, and I. Babuska, "A posteriori error analysis and adaptive processes in the finite element method: Part I-Error analysis," Int. J. Numer. Methods Eng. 19, 1593-1619 (1983). [CrossRef]
- 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," in Biomedical Photonics Handbook (CRC Press, 2003), Chap. 33.
- J. Nocedal and S. J. Wright, Numerical Optimization (Springer-Verlag, New York, 1999). [CrossRef]
- M. J. Epstein, F. Fedele, J. Laible, C. Zhang, A. Godavarty, and E. M. Sevick-Muraca, "A comparison of exact and approximate adjoint sensitivities in fluorescence tomography," IEEE Trans. Med. Imaging 22, 1215-1223 (2003). [CrossRef]
- A. K. Sahu, R. Roy, A. Joshi, and E. M. Sevick-Muraca, "Evaluation of anatomical structure and non-uniform distribution of imaging agent in near-infrared fluorescence-enhanced optical tomography," Opt. Express 13, 10182-10199 (2005). [CrossRef] [PubMed]

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

### Supplementary Material

» Media 1: AVI (161 KB)

» Media 2: AVI (82 KB)

» Media 3: AVI (68 KB)

» Media 4: AVI (366 KB)

» Media 5: AVI (681 KB)

» Media 6: AVI (182 KB)

» Media 7: AVI (90 KB)

» Media 8: AVI (428 KB)

» Media 9: AVI (477 KB)

» Media 10: AVI (1399 KB)

» Media 11: AVI (179 KB)

» Media 12: AVI (88 KB)

» Media 13: AVI (765 KB)

» Media 14: AVI (570 KB)

» Media 15: AVI (1697 KB)

» Media 16: AVI (176 KB)

» Media 17: AVI (89 KB)

» Media 18: AVI (572 KB)

» Media 19: AVI (422 KB)

» Media 20: AVI (1091 KB)

« Previous Article | Next Article »

OSA is a member of CrossRef.