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 (Fast Intersection on Nested Tetrahedron) 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₈, SUB₂, SUB_4a, and SUB_4b, respectively. The regular subdivision SUB₈ creates eight subtetrahedrons and the resulting subtetrahedrons are called regular and labeled as S ₈. SUB₂, 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₂, SUB_4a, and SUB_4b are labeled as S ₂, 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₈; for details, see [9].

 

Fig. 1. 8-subtetrahedron subdivision scheme for nested tetrahedral meshing [9]. The operations shown in (a) through (d) are termed SUB₈, SUB₂, SUB_4a, and SUB_4b, respectively. The subtetrahedrons resulting from the operations (a) to (d) are labeled as S ₈, S ₂, S _4a, and S _4b, respectively.

Download Full Size | PDF

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

 

Fig. 2. Nontrivial intersections resolved by FINT [8]. The figures show the disjoint tetrahedron pieces obtained by FINT for the intersection between the type (a) S ₂-, (b) S _4a-, and (c) S _4b-tetrahedron and the tetrahedrons given from the other mesh embedded therein.

Download Full Size | PDF

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 ₈-tetrahedrons that are regular because only SUB₈ provides the subdivision in the regular pattern. Whenever any tetrahedron on one of the two meshes is subdivided by SUB₈, 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 ₈-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 ₂, 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 μ_axf 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:

subject to the simple lower bound constraint:

where N_S is the number of sources; N_D is the number of detectors; and N_μ is the number of nodes used for the discretization of μ_axf. The term y_ji denotes the fluorescent emission field measured at jth detector for ith source illumination; Φ ⁱ_x and Φ ⁱ_m are the nodal solution vectors of the excitation and the emission fields, respectively, predicted by the model for the ith source illumination; P _j is the projection vector whose inner product with Φ ⁱ_m maps Φ ⁱ_m onto the measured fluorescence at they jth detector location for the ith 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_x,m=iω/c+μ_axi,ami(r)+μ_axf,amf(r), and β_xm=qμ_axf/[1- iωτ(r)]. An index x denotes the excitation field and m denotes the emission field; D_x,m are photon diffusion coefficients; μ_axi,ami is the absorption coefficient due to endogenous chromophores; μ_axf,amf is the absorption coefficient due to exogenous fluorophores; μ´_sx,sm is the reduced scattering coefficient; ω=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 μ_axf 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].

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:

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 μ_amf has a linear relationship with μ_axf given by μ_amf =ζ_axf, only D_x,m is affected. Therefore we need to discretize both D_x,m and μ_axf on the inverse mesh 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 μ_axf inside 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 ith source can be derived from Eq. (12) and the adjoint theorem. With p=μ_axf the derivative of Eq. (12) with respect to p_k at the kth 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 jth detector for the ith 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 Ψ ⁱ_x and Ψ ^j_m are the solutions of Eq. (34). For the point detection geometry, when the jth detector location is r _j, the kth entry P^k_j of the projection vector 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 μ_axf, respectively. For 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_xa,d, and w_ma,d are normalizing factors chosen to be squares of the inverse of maximum nodal values in the amplitudes in Φ ^s_x, Φ ^s_m, Φ ^d_x, and Φ ^d_m, respectively. For a tetrahedral element T in the inverse mesh (inverse element), we simply use Eq. (40) with f=μ_axf for the element error estimate e^T_μ as:

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_min} are the maximum and the minimum values of the parameter 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 ₈), T’s partner must exist in the forward mesh, and otherwise (ii) if T is nonregular (of the type S ₂, S _4a, or S _4b), the partner of T’s parent having 8 child tetrahedrons (of the type S ₈) 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∥₂ < ε_c 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 free and the bound parameters are updated; the maximum step length λ_max is initialized to 1. For a free index i satisfying x_i+d_i<ε, if x_i<ε and d_i<0 we set d_i=0 and delete i from the set of free indices and add i to the set of the bound indices and otherwise compute the maximum step length λ_i satisfying Eq. (2) and update λ_max =min{λ_i, λ_max}. (iv) Backtracking line search is performed for the interval (0, λ_max] and the ratio ρ 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_min} and 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{2D, D_max}. Else, x is updated and D remains unchanged.

 

Fig. 3. Flow diagram for adaptive reconstruction scheme. GN means Gauss-Newtonm, D is the radius of the trust-region and κ is the measure given by Eq. (44).

Download Full Size | PDF

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 ε^T_μ 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.

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

 

Fig. 4. Source/detector configurations for the breast simulating phantom where (a) shows the 27 source locations and (b) shows 128 detector locations.

Download Full Size | PDF

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, μ´_sx=10.8792 cm^-1, μ_ami=0.0322 cm^-1, μ´_sm=9.8241 cm^-1, and ζ= μ_amf/μ_axf =0.1692. The optical parameters specific to fluorescent targets are μ_axf =0.1 cm^-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 μ_axf. 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 Φ_x and Φ_m. To obtain the accurate fluence fields and the target locations with the specified diameter and gaps for μ_axf, the dual-mesh was refined to the fifth subdivision level using η _Φ,μ=0.1 and 7 refinement iterations.

 

Fig. 5. Configurations for spherical fluorescent targets in 5 mm diameter where (a) shows the location of single target and (b) shows the locations of two targets. The diameter is d=5 mm and g is the gap between the two targets.

Download Full Size | PDF

Table. 1. Locations of two spherical targets. All units are in cm.

View Table

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 Δ_max=0.1. The convergence tolerance for the 2-norm of the Gauss-Newton direction d is chosen as ε_c=10^-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.

 

Fig. 6. Initial tetrahedral dual-meshes: (a) and (b) are used for obtaining simulated boundary measurement data, while (c) and (d) are used for iterative reconstruction problem. In the figures, (a) and (c) are the initial meshes discretizing fluence fields, while (b) and (d) are the initial meshes discretizing the parameter μ_axf map. The separately discretized fluence fields and parameter μ_axf map are coupled by FINT.

Download Full Size | PDF

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∥₂<ε_c was encountered. 30 iterations suffice to locate the target. Quantitatively, 99.7 % of the target μ_axf 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.

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

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

 

Fig. 7. (1.39 MB) Movies of the adaptive refinements of the dual-mesh and iterative reconstructions for the single target. The movies show evolutions in the forward mesh in (a) [Media 1], the inverse mesh in (b) [Media 2], two isosurfaces of 80% and 50% of maximum of the reconstructed μ_axf value in (c) and (e) [Media 3] [Media 4], and target values above 0.1% of the maximum of the reconstructed μ_axf value for y=0 slice in (g) [Media 5]. Figures of ideal targets are also shown in (d), (f), and (h) for comparision.

Download Full Size | PDF

 

Fig. 8. (2.64 MB) Movies of the adaptive refinements of the dual-mesh and iterative reconstructions for the two targets with 1 cm gap. The movies show evolutions in the forward mesh in (a) [Media 6], the inverse mesh in (b) [Media 7], two isosurfaces of 80% and 50% of maximum of the reconstructed μ_axf value in (c) and (e) [Media 8] [Media 9], and target values above 0.1% of the maximum of the reconstructed μ_axf value for y=0 slice in (g) [Media 10]. Figures of ideal targets are also shown in (d), (f), and (h) for comparision.

Download Full Size | PDF

 

Fig. 9. (3.38 MB) Movies of the adaptive refinements of the dual-mesh and iterative reconstructions for the two targets with 5 mm gap. The movies show evolutions in the forward mesh in (a) [Media 11], the inverse mesh in (b) [Media 12], two isosurfaces of 80% and 50% of maximum of the reconstructed μ_axf value in (c) and (e) [Media 13] [Media 14], and target values above 0.1% of the maximum of the reconstructed μ_axf value for y=0 slice in (g) [Media 15]. Figures of ideal targets are also shown in (d), (f), and (h) for comparision.

Download Full Size | PDF

 

Fig. 10. (2.41 MB) Movies of the adaptive refinements of the dual-mesh and iterative reconstructions for the two targets with 2.5 mm gap. The movies show evolutions in the forward mesh in (a) [Media 16], the inverse mesh in (b) [Media 17], two isosurfaces of 80% and 50% of maximum of the reconstructed μ_axf value in (c) and (e) [Media 18] [Media 19], and target values above 0.1% of the maximum of the reconstructed μ_axf value for y=0 slice in (g) [Media 20]. Figures of ideal targets are also shown in (d), (f), and (h) for comparision.

Download Full Size | PDF

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

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

 

Fig. 11. Changes in the number of nodes in the forward/inverse meshes, total accumulated computing time, and maximum value of the reconstructed μ_axf distribution map, according to the number of iterations in reconstruction

Download Full Size | PDF

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

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.