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

Over the past few decades, diffuse optical and fluorescence enhanced optical tomography using the near-infrared light has been sought as a new biomedical imaging modality. As a light propagation model in tissue, the diffusion approximation has extensively been used, where the underlying tissue properties are described by the scattering and absorption coefficients. Optical tomography reconstructs the spatial map of the unknown scattering or absorption coefficients as the unknown parameters using the diffusion equations as the forward model. Optical tomography is highly nonlinear and iteratively reconstructs the unknown parameters by reducing the mismatch between the measured and the predicted fluence field data on the tissue boundary.

To obtain numerical solutions of the fluence fields in arbitrary geometries, the finite element method (FEM) has widely been used. Accurate numerical solutions of the fluence fields require fine discretizations of the tissue volume. However, when using the same mesh to solve the inverse parameter estimation problem as used to solve the forward problem, large computational resources are required. Consequently, the dual-mesh scheme using two separate meshes was introduced, wherein the forward mesh defines the forward variables and the inverse mesh defines the unknown parameters to reconstruct [1]. In other words, the dual-mesh scheme decouples the forward variables and inverse parameters. Since the forward variables and the inverse parameters are coupled in the governing equations, the dual-mesh scheme exploits the inter-element maps between the two meshes. The unknown parameter recovery using fixed dual-meshes have been proven to be successful using the fine forward and the coarse inverse meshes while maintaining accurate forward solutions and reasonable reconstruction costs [1–5]. However, when there is no *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.

Joshi et al. [6, 7] first developed fully adaptive, three-dimensional fluorescence enhanced optical tomography technique using hexahedral dual-meshing wherein the forward and the inverse meshes were independently refined or coarsened. The approach employed a forward mesh that was refined/coarsened based upon the spatial gradient of excitation and emission fluence fields and of tissue properties while the inverse mesh was independently refined/coarsened based upon the spatial gradient of tissue fluorophore concentration. For a cubic tissue phantom, they demonstrated that fully adaptive tomography using dual-meshing gives rise to computational efficiency while improving the reconstruction resolution of small fluorescent targets to within theoretical limits. The use of hexahedral elements, however, may limit application to simple rectilinear geometries while medical imaging applications may be expected to require curvilinear geometries.

For optimal, three-dimensional adaptive tomography of arbitrarily shaped volumes or curved boundaries, a tetrahedron based dual-mesh scheme is required. However, independent refinement of tetrahedral elements in the two meshes requires solution to a significant intersection problem between two tetrahedrons on the different meshes. Intersection of any two tetrahedrons from the forward and inverse meshes provides the common region through which the forward variables and the inverse parameters can be mutually referenced to each other in order to enable the finite element assembly. If not optimally solved, this problem prohibits deployment of tetrahedral-based adaptivity in iterative parameter estimation algorithms.Unfortunately, no efficient algorithm resolving intersection in the nested dual-meshing environment has been available and consequently there is no demonstrated method enabling fully adaptive tetrahedron based tomography.

Recently, we developed an efficient algorithm named FINT (*F*ast *I*ntersection on *N*ested *T*etrahedron) resolving intersection in nested tetrahedral dual-meshing by 8-subtetrahedron subdivision scheme [8, 9] and enabling fully adaptive tomography using tetrahedral dual-meshing. In this contribution, as the first application of FINT to model-based parameter reconstruction problems, we present fully adaptive fluorescence enhanced optical tomography technique using tetrahedron based dual-meshing and demonstrate its ability to resolve small fluorescent targets immersed in the breast-simulating phantom with the point illumination/detection configuration and the simulated boundary measurement data. To the authors’ best knowledge, this work is the first demonstration of the fully adaptive tomography using tetrahedral dual-meshing.

This paper is organized as follows: In the Algorithms section, the tetrahedral dual-meshing together with the intersection algorithm FINT is first briefly introduced and then the finite element assembly and the optimization based reconstruction algorithm in the adaptive tetrahedral dual-meshing environment are presented. In the Results and Discussions sections, we show the results for reconstruction of fluorescent targets embedded in a breast phantom. Finally, we conclude by commenting upon potential clinical applications.

## 2. Algorithms

#### 2.1 Nested tetrahedral dual-meshing

The 8-similar subtetrahedron subdivision scheme [9] is employed for creating nested tetrahedral environments from a coarse initial mesh. The 8-subtetrahedron subdivision scheme enforces only 0, 1, or 3 split points on each triangle facets of a tetrahedron, resulting in four different subdivision patterns as illustrated in Fig. 1. The four subdivision operation in Fig. 1 (a) through (d) are termed SUB_{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].

For the initial dual-meshing, the input mesh is duplicated to create a twin mesh. The two twin meshes thus obtained constitute the initial dual-mesh. All elements in the two meshes (labeled as the forward and inverse meshes) are initialized by labeling them as S_{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

Previously, we developed a fast and robust intersection handling algorithm named FINT (Fast Intersection on Nested Tetrahedrons) enabled by a data field named *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.

The partnership is introduced to build a bridge between the twin tetrahedrons wherein one lies in the forward and the other lies in the inverse mesh. Through the partnership, one can iteratively make common assignments between the two independently adapted meshes. For example, we consider that the two meshes in the initial dual-mesh are obtained by duplicating the input mesh and that each tetrahedron duplicated from the original is the twin of the original tetrahedron. The twin relationship is termed as *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.

Intersection between the regular tetrahedrons is trivial. Any nontrivial intersection occurs if one of the tetrahedrons is *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

For the fluorescence enhanced optical imaging in tissue, the coefficient of absorption of the excitation near-infrared light by the fluorophore, *μ _{axf}*, is the parameter to reconstruct and the excitation and fluorescent emission fluence fields

**Φ**

_{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}subject to the simple lower bound constraint:

where *N _{S}* is the number of sources;

*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}The forward model for the frequency domain photon migration is given by the well-known coupled diffusion equations [11]:

where *D _{x,m}*=1/[3(

*μ*+

_{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:

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

The minimization of cost function was performed as follows; (i) the current estimate of *μ _{axf}* in the inverse mesh is passed onto the forward mesh and finite element assembly is performed; (ii)

**Φ**

_{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}### 2.3.2 Finite element assembly

FINT provides disjoint tetrahedron pieces that completely cover up the intersection of two leaf tetrahedrons **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:

$${\mathbf{v}}_{\alpha}=\{{\psi}_{i\text{'}}\left({\mathbf{r}}_{\alpha}\right),{\psi}_{j\text{'}}\left({\mathbf{r}}_{\alpha}\right),{\psi}_{k\text{'}}\left({\mathbf{r}}_{\alpha}\right),{\psi}_{l\text{'}}\left({\mathbf{r}}_{\alpha}\right)\},$$

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

where vertices {*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:

The matrices **U** and **V** perform the change of basis from the virtual basis *π* to the actual basis *φ* and *ψ*, respectively, and inside **Δ** result in:

Integrations of products that *φ*, *ψ*, or both are involved in are performed over **Δ** for finite element assembly by treating *π* as the usual linear finite element basis function.

Since we assume that the heterogeneity is present only in *μ _{axf}* and that

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

while the fields **Φ**
_{x} and **Φ**
_{m} are discretized on the forward mesh into:

The standard Galerkin method using {*φ _{i}*,(

**r**)} discretizes and arranges Eq. (3) into:

The entries of the matrices and the column vector become:

where **Z**, **Θ**, **E**, and **T** are given by:

with

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

where **Z**
^{Δ}, **Θ**
^{Δ}, and **E**
^{Δ} are given by:

that can be calculated exactly using Eq. (7) and Eq. (8). Alternatively, **E** can be assembled directly over the tetrahedrons in the forward mesh *M* since only {*φ _{i}*} is involved. Therefore,

**K**

_{x,m}and

**B**

_{x,m}are assembled by:

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

For practical implementation, we approximate *D _{x,m}* and

*μ*inside

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

**T**´, that is:

The above approximation is not inconsistent with the spirit of the finite element method. Hereby, Eq. (22), Eq. (25), and Eq. (26) reduce to:

and Eq. (23) reduces to **E**
^{Δ} by the approximation. The two 4×4 matrices **Z**
^{Δ} and **E**
^{Δ} computed from Eq. (24) and Eq. (28) are stored in the computer memory for each tetrahedron piece **Δ** and they are updated by running FINT whenever the forward or inverse mesh is changed by adaptation.

### 2.3.3 Jacobian matrix computation

The Jacobian matrix in the boundary fluorescence field measurement for the *i*th source can be derived from Eq. (12) and the adjoint theorem. With *p*=*μ _{axf}* the derivative of Eq. (12) with respect to

*p*at the

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

Introducing adjoint solution vectors **Ψ**
_{x} and **Ψ**
_{m} and association after making inner product with Eq. (31) gives rise to:

since **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 _{k}* reads:

If **Ψ**
_{x} and **Ψ**
_{m} satisfy:

then the L.H.S. of Eq. (32) reduces to Eq. (33). Thereby, the combination of Eq. (32) and Eq. (34) results in the Jacobian matrix entry:

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

**Ψ**

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

The third term is negligible [13] and we reduce to:

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

where **J**
^{Δ} is given by:

In Eq. (39), *γ*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

Since the forward variables of {**Φ**
_{x}, **Φ**
_{m}} and the inverse parameter of *μ _{axf}*, are decoupled using the separate meshes, the two meshes can be independently refined. The forward and the inverse meshes are refined using

*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] wherein the error of a physical quantity

*f*for a tetrahedral element

**T**is estimated by:

where ∣*∂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.

When multiple source and detector data is employed, each tetrahedral element in the forward mesh (*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:

where *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}_{μ}Any tetrahedral element **T** is chosen for refinement in the corresponding meshes if:

The most intriguing question in the fully adaptive strategy using dual-meshing is when to *trigger* refinement in the sequence of the iterative reconstructions. This question is indirectly related to the issue of optimal parameter discretization. From the viewpoint of the discretization level, the solution to the optimal parameter discretization problem can be translated into the answer to the question on the degree of approximation on the smoothness of the discretized continuous parameters. Hereby, we define a *smoothness* measure *κ* for a leaf inverse element **T** as:

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

the refinements to the forward/inverse meshes are triggered. Once the refinement is triggered, the forward elements are refined first and then the inverse elements are refined. In any case, any tetrahedron of the subdivision level 0 is refined using the criterion Eq. (43). A restriction is made to the inverse elements of the nonzero subdivision levels: if an inverse element **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.

The partnership is updated whenever either the forward or the inverse mesh is refined, while {**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

In minimizing the object function given by Eq. (1) subject to the simple bound constraint Eq. (2), we used the variant of the active set method where the standard procedure is followed with the exception that the updates of the active set *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]. If ∥**d**∥_{2} < *ε _{c}* for the sufficiently small number

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

The above optimization scheme must be combined with the adaptive scheme. During the iterative reconstructions, we periodically check whether or not to trigger adaptation with a predefined period of the iterations. The first refinement to the initial dual-mesh is always triggered and performed for elements satisfying Eq. (43). The next upcoming adaptations are triggered by Eq. (45). Once triggered, the forward elements with *ε ^{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}_{μ}An intriguing situation that can arise is the successive failing in object function minimization while no mesh adaptation is performed even though *D* is *D _{min}*. Since our optimization scheme is not based upon the Lagrangian formalism, we can rescue the situation by switching the status of some

*bound*parameters back to

*free*on the

*proximity*basis. The proximity of the parameter distribution provides the good candidate for the switching. We search the

*leaf*-level edges having bound and free endpoints together and change the status of the bound endpoints back to free. That is, only the bound parameters

*graph-theoretically connected*to the free parameters are switched back to free. However, the authors’ experience suggests that such situations are not encountered for reasonably chosen

*D*and

_{min}*D*.

_{max}In the meantime, we also use the above proximity based status switching for the first few refinements, enabling the regions of interest to grow or move from place to place until they remain relatively constant. All parameters are initially set to 0 and assigned the *free* status and then the iterative reconstructions begin.

## 3. Results

To assess fully adaptive fluorescence optical tomography using tetrahedral dual-meshing, we chose the breast simulating phantom geometry with the point illumination/detection configurations as illustrated in Fig. 4. The breast phantom has two parts where the chest-simulating part is the bottom cylinder of 20 cm diameter and 3 cm height and the breast-simulating part consists of a hemisphere of 10 cm diameter and a cylinder of 10 cm diameter and 0.5 cm height. We used 27 point sources and 128 point detectors attached on the boundary of the breast-part as shown in Fig. 4.

In the simulations presented in this contribution, one or two fluorescent spherical targets are embedded in the otherwise homogeneous background of the breast phantom where the size of the fluorescent spheres considered is 5 mm in diameter. The optical properties assumed are *μ _{axi}*=0.02483 cm

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

For the iterative reconstructions, random noise of maximum 5 % amplitude and ±2% phase were added to the amplitude and the phase of the simulated boundary measurement data. The initial forward and inverse meshes for the reconstruction are illustrated in Fig. 6(c) and (d), respectively. The initial forward mesh with 1129 nodes shown in Fig. 6(c) is obtained by the boundary refinement of the breast part of the initial inverse mesh with 314 nodes shown in Fig. 6(d). For the reconstruction problem, the refinement tolerances in Eq. (43) and Eq. (45) were chosen as *η*
_{Φ,μ}=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 _{max}*=0.001 and

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

The movies of adapted dual-meshes and reconstructions are illustrated in Fig. 7. The number of nodes in the forward/inverse meshes are 1955/554, 6279/794, 12660/945, 14139/1218, and 14139/1218 for IT=30, 60, 100, 150, and 261, respectively. The maximum values or peaks of the reconstructed *μ _{axf}* are 0.0164, 0.0240, 0.0377, 0.0706, and 0.2028 cm

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

The movies of adapted dual-meshes and reconstructions are illustrated in Fig. 8. The number of nodes in the forward/inverse meshes are 1966/486, 7821/726, 16819/1015, 21246/1075, and 21264/1515 for IT=30, 60, 100, 150, and 300, respectively. The maximum values or peaks of the reconstructed *μ _{axf}* distribution map are 0.0107, 0.0248, 0.0345, 0.0552, and 0.1790 cm

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

The movies of adapted dual-meshes and reconstructions are illustrated in Fig. 9. The number of nodes in the forward/inverse meshes are 1769/489, 4862/737, 12887/1099, 19987/1306, and 20013/1538 for IT=30, 60, 100, 150, and 300, respectively. The maximum values or peaks of the reconstructed *μ _{axf}* distribution map are 0.0129, 0.0204, 0.0320, 0.0553, and 0.1499 cm

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

The movies of adapted dual-meshes and reconstructions are illustrated in Fig. 10. The number of nodes in the forward/inverse meshes are 1686/452, 5879/787, 10254/1013, 15998/1543, and 15998/1543 for IT=30, 70, 100, 150, and 300, respectively. The maximum values or peaks of the reconstructed *μ _{axf}* distribution map are 0.0241, 0.0356, 0.0431, 0.0689, and 0.2189 cm

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

For all cases considered so far, we have successively demonstrated that the targets can be reconstructed at the correct locations with a small number of discretized parameters by the fully adaptive FEM using the nested tetrahedral dual-meshes. Even for the small targets whose separation reaches the diffusion limit, adaptive FEM is able to resolve them without the large memory overhead. We were able to minimize number of nodes in the inverse mesh using the smoothness measure in Eq. (44). Our results indicate that the criterion in Eq. (45) is effective in the adaptation control. Fig. 11 illustrates the statistics in the number of nodes in the forward/inverse meshes, total accumulated computing time, and the changes in maximum value of the reconstructed *μ _{axf}* 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.

However, the 300 iterations performed in this work is excessive and not practical. Once the reconstructed targets have been correctly located, further iterations generally result in more localized peaks eventually increasing above the correct *μ _{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.

As the iterative reconstruction proceeds, the average slope in the maximum *μ _{axf}* value changes around IT=130 as indicated by the dotted circled region in Fig. 11. We observe that the dotted circled region in the plot of the maximum

*μ*value is related with the transition to the narrowing stage from the localization stage and that the FWHM (i.e., the average diameter of 50% isosurfaces) of the peaks in

_{axf}*μ*distribution is approximately the same as the ideal target size. For the iterations above the circled region, peaks are clearly seen to increase with the narrowing FWHM. Detection of such slope-changing feature indicates that overly iterated situations arise and suggest a potential criterion for stoping iterations. The average slope-changing features are also observed in the plot of the number of nodes in the inverse mesh as indicated by the dotted circled regions in Fig. 11 around IT=50, 80, 120, and 140 for the single target, and for the 1 cm, 5 mm, and 2.5 mm separated two targets. For the iterations above the circled region, the increase in the number of nodes in the inverse mesh is retarded. Such retardation in the increase in the number of nodes feature probably indicates that the refinements are achieved to sufficient levels and overly refined situations are arising. However, there are cases in which more than one average slope-changing feature exists. In addition, one could conceive of cases in which even such features are not present. The algorithms on how to implement these ideas of truncating refinements or iterations have not yet been developed and further work is necessary to optimize our tetrahedral dual-mesh based fully adaptive tomography algorithm. Future work will include the development of the algorithms handling the truncations using features mentioned above.

_{axf}One might argue against our claim that we have performed totally *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 *μ _{axf}* 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

*μ*reconstruction to a substantial degree. Sahu and coworkers [14] validated this fact by showing that fluorescence reconstructions are robust for up to 100% random perturbation in the values of background optical properties. Therefore, heterogeneity of background endogenous properties has minimal influence on reconstruction results. In practice, a prior spectroscopy measurement of average background intrinsic optical properties should suffice for accurate fluorescence reconstructions. In future work, the influence of the deviations from average background properties on the reconstructed image quality will be assessed.

_{axf}## 5. Summary and Conclusions

For the first time, we have developed the fully adaptive FEM using tetrahedral dual-meshing and applied to the fluorescence enhanced optical tomography. We have assumed no *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.

A few issues need to be resolved on how and when to truncate the iterative reconstructions in the adaptive meshing environments. However, we believe that this work is the first step towards the development of the practical optical tomography system. This work is the first demonstration of the development and application of the adaptive tetrahedral dual-mesh FEM to the tomography, especially in terms of the fluorescence enhanced optical imaging in tissue. The fully adaptive tetrahedral dual-mesh based FEM developed by the authors has a wide range of applicability to the various tomography problems including the electrical impedance and microwave computed tomography, etc. Most importantly, it can be applied to the arbitrary curved geometries. Future works will include solution to the truncation issues, application to the experimental data, and hopefully the parallelization of the proposed fully adaptive tomography algorithm as an extension of this work.

## Acknowledgments

The authors acknowledge the support of the National Institute of Health (NIH Grant no. R01 CA112679).

## 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 **14**, 504–514 (1995). [CrossRef] [PubMed]

**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. **25**, 183–193 (1998) [CrossRef] [PubMed]

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

**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. **42**, 3117–3129 (2003). [CrossRef] [PubMed]

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

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

**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,” (*submitted*).

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

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

**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 *Biomedical Photonics Handbook* (CRC Press, 2003), Chap. 33.

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

**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 , **22**, 1215–1223 (2003). [CrossRef]

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