OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 15, Iss. 11 — May. 28, 2007
  • pp: 6955–6975
« Show journal navigation

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

Jae Hoon Lee, Amit Joshi, and Eva M. Sevick-Muraca  »View Author Affiliations


Optics Express, Vol. 15, Issue 11, pp. 6955-6975 (2007)
http://dx.doi.org/10.1364/OE.15.006955


View Full Text Article

Acrobat PDF (1364 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

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]

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

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]

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

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

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]

] 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

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

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

] 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

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

] 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 SUB8, SUB2, SUB4a, and SUB4b, respectively. The regular subdivision SUB8 creates eight subtetrahedrons and the resulting subtetrahedrons are called regular and labeled as S 8. SUB2, SUB4a, or SUB4b 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 SUB2, SUB4a, and SUB4b 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 S8; for details, see [9

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

].

Fig. 1. 8-subtetrahedron subdivision scheme for nested tetrahedral meshing [9]. The operations shown in (a) through (d) are termed SUB8, SUB2, SUB4a, and SUB4b, respectively. The subtetrahedrons resulting from the operations (a) to (d) are labeled as S 8, S 2, S 4a, and S 4b, respectively.

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

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

]. 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 2-, (b) S 4a-, and (c) S 4b-tetrahedron and the tetrahedrons given from the other mesh embedded therein.

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

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

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

minpI(p)=12i=1NSj=1ND[PjTΦmi(p)yji]2
(1)

subject to the simple lower bound constraint:

pk0,k=1,2,,Nμ
(2)

where NS is the number of sources; ND is the number of detectors; and Nμ is the number of nodes used for the discretization of μaxf. The term yji denotes the fluorescent emission field measured at jth detector for ith source illumination; Φ ix and Φ im 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 Φ im maps Φ im 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

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.

]:

{[Dx(r)Φxrω]+kxΦxrω=0[Dm(r)Φmrω]+kmΦmrω=βxmΦxrω
(3)

where Dx,m=1/[3(μaxi,ami+μaxf,amf+μ´sx,sm)], kx,m=iω/c+μaxi,ami(r)+μaxf,amf(r), and βxm=axf/[1- iωτ(r)]. An index x denotes the excitation field and m denotes the emission field; Dx,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:

{2DxnΦx+γΦx=Q(r)2DmnΦm+γΦm=0
(4)

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

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

].

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:

uα={φi(rα),φj(rα),φk(rα),φl(rα)},
vα={ψi'(rα),ψj'(rα),ψk'(rα),ψl'(rα)},
(5)

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

πα(rβ)=δαβ
(6)
Δ=πpaπqbπrcπsddv=6VΔa!b!c!d!(a+b+c+d+3)!,
(7)

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:

Unα=φn(rα),Vnα=ψn(rα).
(8)

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

{φi(r),ψi(r)}=α{Uiα,Viα}πα(r)
(9)

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 Dx,m is affected. Therefore we need to discretize both Dx,m and μaxf on the inverse mesh M´ giving rise to:

{Dx,m(r),μaxf(r)}=i{Dx,mi,μaxfi}ψi(r)
(10)

while the fields Φ x and Φ m are discretized on the forward mesh into:

Φx,m(r)=iΦx,miφi(r)
(11)

The standard Galerkin method using {φi,(r)} discretizes and arranges Eq. (3) into:

[KxOBxmKm][ΦxΦm]=[QO]
(12)

The entries of the matrices and the column vector become:

Kx,mij=kDx,mkZk,ij+kμaxf,amfkΘk,ij+kx,m0Eij+γ2Γij,
(13)
Bxmij=χkμaxfkΘk,ij,
(14)
Qi=12Ωφi(r)Q(r)da,
(15)

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

Zk,ij=Ωψk(r)φi(r).φj(r)dv,
(16)
Θk,ij=Ωψk(r)φi(r)φj(r)dv,
(17)
Eij=Ωφi(r)φj(r)dv,
(18)
Γij=Ωφi(r)φj(r)da,
(19)

with

kx,m0=iωc+μaxi,ami,χ=q(1iωτ).
(20)

Γ and Q can be assembled directly over the boundary triangle elements in the forward mesh M. On the other hand, the finite element assembly of the volume terms Z, Θ, and E are performed over the tetrahedron piece provided by FINT, that is:

{Z,Θ,E}=Δ{ZΔ,ΘΔ,EΔ}
(21)

where Z Δ, Θ Δ, and E Δ are given by:

Zk,ijΔ=αβγVkγUiαUjβΔπγπαπβdv
(22)
Θk,ijΔ=αβγVkγUiαUjβΔπγπαπβdv,
(23)
EijΔ=αβUiαUjβΔπαπβdv,
(24)

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:

Kx,mij=Δ(kDx,mkZk,ijΔ+kμaxf,amfkΘk,ijΔ)+kx,m0Δ(kEijΔ)+γ2Γij,
(25)
Bxmij=χΔ(kμaxfkΘk,ijΔ)
(26)

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 Dx,m and μaxf inside T´ by taking the average of the values at the four nodes of T´, that is:

Dx,m(r)14iDx,mi,μaxf(r)14iμaxfi.
(27)

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

ZijΔ=αβUiαUjβΔπα.πβdv,
(28)
Kx,mij=14Δ(kDx,mk)ZijΔ+14Δ(kμaxf,amfk)EijΔ+kx,m0Δ(kEijΔ)+γ2Γij,
(29)
Bxmij=14χΔ(kμaxfk)EijΔ,
(30)

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 pk at the kth node in the inverse mesh can be arranged into:

[KxOBxmKm][ΦxipkΦxipk]=[KxpkOBxmpkKmpk][ΦxiΦmi]
(31)

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

([KxBxmOKm][Ψ xΨ m])T[ΦxipkΦxipk]=[ΨxTΨmT][KxpkOBxmpkKmpk][ΦkiΦmi]
(32)

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 pk reads:

Jk,ij=PjT(Φmipk).
(33)

If Ψ x and Ψ m satisfy:

[KxBxmOKm][ΨxΨm]=[OPj],
(34)

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:

Jk,ij=ΨxjT(Kxpk)Φxi+ΨmjT(Bxmpk)ΦxiΨmjT(Kmpk)Φmi,
(35)

where Ψ ix and Ψ jm are the solutions of Eq. (34). For the point detection geometry, when the jth detector location is r j, the kth entry Pkj of the projection vector P j becomes:

Pjk=φk(rj)
(36)

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

Jk,ijΨxjT(Kxpk)Φxi+ΨmjT(Bxmpk)Φxi.
(37)

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

J=ΔJΔ
(38)

where J Δ is given by:

Jk,ijΔ=14(γDxγpkδkγ)αβ(Ψxj)αZαβΔ(Φxi)β14(γδkγ)αβ(Ψxj)αEαβΔ(Φxi)β+14χ(γδkγ)αβ(Ψmj)αEαβΔ(Φxi)β.
(39)

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

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]

] wherein the error of a physical quantity f for a tetrahedral element T is estimated by:

εT[f]=hTfn±2da, ,
(40)

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 eT Φ for the forward element T as:

eΦT=s=1NS(wx,sεT[Φxs]+wm,sεT[Φms])+d=1ND(wxa,dεT[Ψ xd]+wma,dεT[Ψ xd])
(41)

where wx,s, wm,s, wxa,d, and wma,d are normalizing factors chosen to be squares of the inverse of maximum nodal values in the amplitudes in Φ sx, Φ sm, Φ dx, and Φ dm, 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 eTμ as:

eμT=εT[μaxf].
(42)

Any tetrahedral element T is chosen for refinement in the corresponding meshes if:

eΦ,μT>ηΦ,μmaxT{eΦ,μT},0<ηΦ,μ<1.
(43)

κ=[pc+pd]2pmpmaxpmin
(44)

κ>θ,0<θ<12,
(45)

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

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

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.

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.

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

table-icon
View This 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 Dmax=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.

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 ∥d2<ε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.
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.
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.
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.

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.

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

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

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]

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

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]

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

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: March 8, 2007
Revised Manuscript: May 17, 2007
Manuscript Accepted: May 17, 2007
Published: May 22, 2007

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

Citation
Jae Hoon Lee, Amit Joshi, and Eva M. Sevick-Muraca, "Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue," Opt. Express 15, 6955-6975 (2007)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-11-6955


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. K. D. Paulsen, P. M. Meaney, M. J. Moskowitz, J. M. Sullivan, Jr., "A dual mesh scheme for finite element based reconstruction algorithms," IEEE Trans. Med. Imaging 14, 504-514 (1995). [CrossRef] [PubMed]
  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. Methods Eng. 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]

Cited By

Alert me when this paper is cited

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

Supplementary Material


» Media 1: AVI (161 KB)     
» Media 2: AVI (82 KB)     
» Media 3: AVI (68 KB)     
» Media 4: AVI (366 KB)     
» Media 5: AVI (681 KB)     
» Media 6: AVI (182 KB)     
» Media 7: AVI (90 KB)     
» Media 8: AVI (428 KB)     
» Media 9: AVI (477 KB)     
» Media 10: AVI (1399 KB)     
» Media 11: AVI (179 KB)     
» Media 12: AVI (88 KB)     
» Media 13: AVI (765 KB)     
» Media 14: AVI (570 KB)     
» Media 15: AVI (1697 KB)     
» Media 16: AVI (176 KB)     
» Media 17: AVI (89 KB)     
» Media 18: AVI (572 KB)     
» Media 19: AVI (422 KB)     
» Media 20: AVI (1091 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited