1. Introduction

Most developing optical tomography techniques, which employ time-dependent measurements and fluorescent contrast agents, depend upon (i) illuminating a tissue surface with incident points of light using either fiber optics or lens-focused light beams, and (ii) collecting the propagated light at various points on the tissue surface [1–16]. Besides creating practical problems for clinical implementation particularly when used with a fluorescent contrast agent, point illumination results in a limited tissue volume of interrogation, potentially missing the target tissue of interest. To avoid this problem, mechanical raster scanning has been proposed in animal studies [17–21]. In a simpler and more practical approach, investigators in our laboratories have experimentally demonstrated 3-D tomographic reconstruction from area measurement of a single projection of fluorescence resulting from planar illumination of modulated excitation light [13, 22–24]. Area-illumination and detection techniques are promising as the non-contact nature of data acquisition minimizes patient discomfort and reduces data acquisition time. While area-illumination provides enhanced volumes of excitation light interrogation, the tomography problem is hard to solve because of the increased ill-posedness introduced by the availability of fluorescence measurements on only a single projection plane. In potential clinical applications of fluorescence optical tomography involving sentinel lymph node mapping for imaging micrometastases in the lymph nodes of the axilla, internal mammary nodes, and supraclavicular nodes, one reflectance fluorescence measurement will yield insufficient information for accurate fluorescence yield map reconstructions. On the other hand, as will be shown in this contribution, the information content in the area measurements can be increased by employing spatially varying or patterned illumination, and taking multiple reflectance measurements corresponding to different illumination patterns. Spatially patterned illumination can be provided by scanning line sources, Gaussian spots, or patterns generated by diffractive optics components on different locations on the tissue phantom. Multiple frequency-domain fluorescence measurements are then acquired on the illumination surface corresponding to different illumination patterns. Fluorescence yield distributions in the imaged tissue volume can then be recovered by a model based tomography algorithm.

Model based tomography in an area illumination and detection framework requires accurate numerical solution of the photon fluence distribution in the tissue created by spatially patterned excitation illumination. In the past, we have proposed an efficient adaptive finite element tomography algorithmto solve the plane wave excitation based fluorescence tomography problem[23] and demonstrated its ability to invert experimental area fluorescence measurements [24]. In this contribution, we extend the adaptive finite element based tomography approach to handle multiple spatially patterned illumination sources. Adaptive finite element solutions are computed for multiple area sources independently of each other on separate finite element meshes, which are independently adapted to accurately resolve each illumination pattern. Information from all measurements is then combined to update the unknown interior fluorescence distribution in tissue on another adaptively refined finite element mesh to better resolve the embedded fluorescent targets. The proposed scheme is implemented in a parallelized framework, so that all computations related to different illumination patterns are solved simultaneously on different machines of a Linux Beowulf cluster computer, and sent back to a master node for joint reconstruction of the yield distribution.

The objective of this contribution is to introduce the approach and image reconstruction algorithm for non-contact spatially patterned excitation based fluorescence tomography. To this end, we perform a series of numerical experiments with different excitation patterns for reconstructing multiple fluorescence targets at varying depths from the illumination plane, and to quantify the advantage of using multiple measurement sets.

The outline of this article is as follows: In Section 2, we detail the formulation and implementation of the proposed fluorescence tomography scheme. Section 3 discusses the computational experiments conducted to demonstrate the parallel adaptive fluorescence tomography scheme for reconstructing multiple fluorescence targets embedded (at constant and varying depths from the illumination surface) in turbid media. We provide a quantitative assessment in Section 4. Finally, the results are summarized in Section 5. To the best of the authors’ knowledge, this contribution represents the first time that spatially patterned modulated area excitation sources have been demonstrated in optical tomography to improve fluorescence image reconstructions from reflectance measurements.

2. Methods and formulation

2.1. The photon transport model

Fluorescence optical tomography is typically performed in a model-based framework, wherein a photon transport model is used to generate predicted boundary fluorescence measurements for a given fluorescence absorption map in the tissue interior. The map of the absorption owing to fluorophore is then iteratively updated until the predicted boundary fluorescence measurements converge to the actual experimentally observed fluorescence measurements. For photon propagation in large tissue volumes, the following set of coupled photon diffusion equations is an accurate model:

Here,

$D_{x,m}=\frac{1}{3\left({\mu }_{\mathit{ax},\mathit{mi}}+{\mu }_{\mathit{ax},\mathit{mf}}+{\mu }_{\mathit{sx},m}^{\prime }\right)},k_{x,m}=\frac{i\omega }{c}+{\mu }_{\mathit{ax},\mathit{mi}}+{\mu }_{\mathit{ax},\mathit{mf}},{\beta }_{\mathit{xm}}=\frac{\varphi {\mu }_{\mathit{axf}}}{1-i\omega \tau },$

and subscripts x and m denote the excitation and the emission light fields, respectively. u,v are the complex-valued photon fluence fields at excitation and emission wavelengths, respectively; D _x,m are the photon diffusion coefficients; µ _ax,mi is the absorption coefficient due to endogenous chromophores; µ _ax,mf is the absorption coefficient due to exogenous fluorophore; ${\mathit{\mu }}_{\mathit{\text{sx}},m}^{\prime }$ is the reduced scattering coefficient; ω is the modulation frequency; ϕ is the quantum efficiency of the fluorophore, and finally, τ is the fluorophore lifetime associated with first order fluorescence decay kinetics. All these coefficients, and in particular the fluences u,v and the absorption/scattering coefficients are spatially variable.

Above equations are complemented by Robin-type boundary conditions on the boundary ∂Ω of the domain Ω modeling the NIR excitation source:

where n denotes the outward normal to the surface and γ is a constant depending on the optical refractive index mismatch at the boundary. The complex-valued function S=S(r) is the spatially variable excitation boundary source. There is no source term for the emission boundary condition. The goal of fluorescence tomography is to reconstruct the spatial map of coefficients µ _axf (r) and/or τ(r) from measurements of the complex emission fluence v(r) on the boundary. In this work we will focus on the recovery of only µ _axf (r), while all other coefficients are considered known a priori. For notational brevity, and to indicate the special role of this coefficient as the main unknown of the problem, we set q=µ _axf in the following paragraphs.

2.2. The inverse problem for multiple illumination patterns

We have previously proposed a novel fluorescence tomography algorithm utilizing adaptive finite element methods [23]. In the following, we briefly describe the formulation of the scheme and its extension to image reconstructions from multiple area illumination sources. To this end, assume that we employ W different area excitation light patterns S ⁱ (r), i=1,2, …, W, to excite the embedded fluorophore in the phantom. For each of these experiments, we can predict fluences u ⁱ ,v ⁱ satisfying (1)–(3) with S ⁱ (r) as source terms if we assume knowledge of the yield map q. In addition, we take fluorescence measurements on the measurement boundary G for each of these experiments that we will denote by z ⁱ . The fluorescence image reconstruction problem is then posed as a constrained optimization problem wherein an L ₂ norm based error functional of the distance between fluorescence measurements z={z ⁱ , i=1,2, …, W} and the diffusion model predictions v={v ⁱ , i=1,2, …, W} is minimized by variation of the parameter q; the diffusion model for each fluence prediction vⁱ is a constraint to this optimization problem. In a function space setting, the mathematical formulation of this minimization problem reads as follows:

Here, the error functional J(q,v) incorporates a least squares error term over the measurement part Γ of the boundary ∂Ω and a Tikhonov regularization term:

where β is the Tikhonov regularization parameter. The Tikhonov regularization term βr(q) is added to the minimization functional to control undesirable components in the map q(r) that result from a lack of resolvability. The constraint A ⁱ (q; [u ⁱ ,v ⁱ ])([ζ ⁱ ,ξ ⁱ ])=0 is the weak or variational form of the coupled photon diffusion Eq.s (1)–(3) with partial current boundary conditions for the i ^th excitation source, and with test functions [ζ,ξ] ∊H ¹(Ω):

The solution of minimization problem (4) is then a stationary point of the combined Lagrangian [25]

Here, ${\lambda }_i^{\mathit{\text{ex}}}$ ,${\lambda }_i^{\mathit{\text{em}}}$ are the Lagrange multipliers corresponding to the excitation and emission diffusion equation constraints for the i ^th source, respectively, and we have introduced the abbreviation x={u,v,λ ^ex ,λ ^em ,q} for simplicity; u={u ⁱ , i=1,2, …, W} and v={v ⁱ , i=1,2, …, W} are excitation and emission fluences for the W excitation sources; λ ^ex ={${\lambda }_i^{\mathit{\text{ex}}}$ , i=1,2, …, W} and λ ^em ={${\mathrm{\lambda }}_i^{\mathit{\text{em}}}$ , i=1,2, …, W} denote the Lagrange multipliers corresponding to the excitation and emission equation constraints for the W excitation sources.

In contrast with the output least squares based optical tomography techniques, we do not directly enforce the dependence of the state variables u and v on the parameter q by solving diffusion equations for a given parameter map and substiuting the solution in place of v in the error functional. Rather this dependence is implicitly enforced by treating u,v,q as independent variables and including diffusion equations as a constraint on the the error functional.

2.3. Solving the inverse problem

A stationary point of L(x), and therefore a solution of the constrained minimization problem (4), is found using the Gauss-Newton method wherein the update direction δx _k ={δ u _k ,δv _k ,${\delta \lambda }_k^{\mathit{\text{ex}}}$ ,${\delta \lambda }_k^{\mathit{\text{em}}}$ ,δ _qk } is determined by solving the linear system

where L _xx (x _k ) is the Gauss-Newton approximation to the Hessian matrix of second derivatives of L at point x _k , and y denotes possible test functions. These equations represent one condition for each variable in δx _k , i.e. they form a coupled system of equations for the 4W+1 variables involved: all W excitation and emission fluences, excitation and emission Lagrange multipliers, and the yield map q. Once the search direction is computed from Eq. (8), the actual update is determined by calculating a safeguarded step length α _k :

The step-length α _k can be computed from one of several methods, such as the Goldstein-Armijo backtracking line search [26, 27]. The Gauss-Newton equations are discretized by the finite element method. State and adjoint variables u,v,λ ^ex , and λ ^em are discretized and solved for on meshes with continuous finite elements, while the unknown parameter map q is discretized on a separate mesh with discontinuous finite elements. Whenever Gauss-Newton iterations on these meshes have reduced the error function by a factor of 10^-3 or the Gauss-Newton step length returned by the line search algorithm has fallen below 0.15, the meshes are refined using a posteriori refinement criteria. The advantage of function space Lagrangian framework lies in its ability to provide a mathematically cohesive view of the multiple experiment tomography problem. As the inverse iterations are derived in a discretization independent manner, it makes it easier to implement optimal finite elementmeshes corresponding to different excitation sources.

The discretized version of the Gauss-Newton system (8) leads to a matrix equation for the unknowns of our discrete finite element update δx _k . It has the following form:

where the updates for the primal and dual (Lagrange multiplier) variables are abbreviated as δp _k =[δ u _k ,δ v _k ] ^T , δd _k =[${\delta \lambda }_k^{\mathit{\text{ex}}}$ ,${\delta \lambda }_k^{\mathit{\text{em}}}$ ] ^T . Since each of these 4W+1 variables is discretized with several ten or hundred thousand unknowns on our finite element meshes, the matrix on the left hand side can easily have a dimension of several million to over 10 million. At first glance, it therefore seems infeasible or at least very expensive to solve such a system. In the past, this has led researchers to the following approach: use one experiment alone to solve for the fluorescence map, then use the result as the starting value for inverting the next data set and so on; one or several loops over all data sets may be performed. While this approach often works for problems that are only moderately ill-posed, it is inappropriate for problems with the severe ill-posedness of the one at hand. The reason is that if we scan the source over the surface, we will only be able to identify the yield map in the vicinity of the illuminated area. Far away from it, we have virtually no information on the map. We will therefore reconstruct invalid information away from the source, erasing all prior information we have already obtained there.

 

Fig. 1. Block structure for the (square) Gauss-Newton matrix (10) for multiple excitation sources. Indices denote the number of the experiment a matrix corresponds to. W is the total number of experiments (excitation illumination patterns or multiple imaging modalities) used. Blocks not shown correspond to matrices with all-zero entries.

Download Full Size | PDF

Consequently, it is mandatory that we use all available information at the same time, in a joint inversion scheme. Fortunately, we can make use of the structure of the problem: because experiments are independent of each other, the joint Gauss-Newton matrix is virtually empty, and in particular has no couplings between the entries corresponding to primal and dual variables of different illumination experiments. The only thing that keeps the matrix fully coupled is the presence of the yield map q on which all experiments depend.

This structure is manifested in the fact that M, the second derivative of the measurement error function with respect to state variables u ⁱ ,v ⁱ for all the excitation sources, is a block diagonal matrix {M ₁,M ₂,…,M _W }. Likewise P=blockdiag{P ₁,P ₂,…,P _W } is the representation of the discrete forward diffusion model for all the excitation sources. Finally, the matrix C=blockdiag{C ₁,C ₂,…,C _W } is obtained by differentiating the semi-linear form A ⁱ in Eq. (6) with respect to the parameter q. For reasons explained below, we choose different meshes for individual measurements. Consequently the individual blocks M _i ,P _i ,C _i all differ from each other. Finally, the right hand side F denotes the discretized form of -L _x (x _k )(y). The detailed formulation of the individual blocks M _i ,P _i ,C _i and the right hand side F is provided in Ref. [23] for the single excitation source based fluorescence measurements.

Given these considerations, the block structure of the Gauss-Newton KKT matrix (10) is shown in Fig. 1. Using this structure, we can form the Schur complement of this system with respect to the R block that reflects the independence of experiments:

Here, we first have to solve for the yield map updates δq _k , and then for updates of state and adjoint variables for all the experiments individually and independently, a task that is obviously simpler than solving for the one big and coupled matrix in (10). Note in particular that the Schur complement matrix

is symmetric and positive definite. We can therefore use the Conjugate Gradient (CG) method to invert it.

In the CG algorithm, solving for λq _k is done iteratively; in each iteration, one multiplication between S and a vector is required. Given the structure of the matrix, this can be implemented on separate computers or separate processors on a multiprocessor system, each of which is responsible for one or several of the experiments (and corresponding matrices C _i ,P _i ,M _i ). Since multiplication of a vector with the matrices $C_i^T$ $P_i^{-T}$ M _i $P_i^{-1}$ C _i is completely independent, a workstation cluster with W nodes is able to perform the image reconstruction task within approximately the same time as a single machine requires for inverting a single excitation source. For the examples shown in Section 3, reconstructions took between 10 and 20 minutes on a cluster of Opteron machines, independently of the number of experiments W.

Using this approach, Eq.s (11)–(13) are iteratively solved to convergence in order to obtain the unknown fluorescence map q. The results shown below are created using a program implementing this algorithm developed by Bangerth [27] with the help of the open source deal.II finite elements library [28]. Further implementation details including the line search procedure and the incorporation of bounds on the unknown parameters were described previously [23].

2.4. Choice of meshes

To date, finite element simulations for multiple illumination sources in a non-contact area illumination mode have not been reported in the literature. However, traditionally, optical tomography schemes for fiber illumination setups have utilized only one finite element mesh for the solution of diffusion equations corresponding to different illumination sources. If one does that, one ends up with a simpler matrix in (10) where all the block matrices M _i are equal, and similarly for the matrices and P _i . (The blocks C _i are still different since their elements are computed using the forward solutions u ⁱ ,v ⁱ .)

However, such a scheme is neither efficient nor accurate. Simulation of photon transport via area excitation illumination requires careful finite element mesh design to capture the photon fluence variation in the tissue media. A well-adapted, locally refined mesh tailored the illumination by a single source pattern would be fine only in the vicinity of illumination, and coarse far away from it. Such a mesh would therefore be unsuitable for illumination by a different source at a different position. Conversely, a uniform mesh would have to be prohibitively fine everywhere to accurately predict fluences for all possible illumination scenarios.

 

Fig. 2. Schematic of the multiple excitation illumination tomography algorithm. The implementation either runs on a single machine, or on a Linux Beowulf cluster where each node is responsible for one or several measurements.

Download Full Size | PDF

 

Fig. 3. Source patterns employed for multiple experiment fluorescence optical tomography: scanning lines (left column), scanning Gaussians (center column), diffractive optics patterns (right column)

Download Full Size | PDF

Consequently, efficient and accurate computations can only be performed if different meshes are chosen for each illumination pattern. Our code uses separate meshes for each experiment adapted to the illumination pattern employed (Fig. 4). Details of the refinement criteria are given in Ref. [29]. Such an algorithm runs efficiently on multiprocessor computers or workstation clusters wherein each source or a small number of sources can be simulated on a separate compute node. The flow of this algorithm for parallel computations is schematically shown in Fig. 2.

3. Image reconstruction simulations

In this section, non-contact fluorescence tomography with multiple excitation sources will be demonstrated. Synthetic measurements were generated on a 512 ml cubical tissue phantom with optical properties of 1% Liposyn by running a highly accurate finite element forward simulator with a yield map that is assumed to be known. For the simulated background absorption coefficient, we chose µ _axi =0.023cm ^-1 and µ _ami =0.0289cm ^-1 [30]. The absorption coefficient due to fluorophore at the excitation wavelength was set to µ _axf =0.23cm ^-1 in the target. The associated emission wavelength absorption coefficient was µ _amf =0.023cm ^-1 in the target. The lifetime of the fluorophore was taken to be τ=0.56ns and the quantum efficiency was ϕ=0.016 to match the corresponding properties of Indocyanine Green (ICG) dye used in experiments. The excitation wavelength for ICG is 785nm and the emission data is collected at 830nm. The reduced scattering coefficient was taken as ${\mathit{\mu }}_s^{\prime }$ =9.84cm ^-1 in both the target and the background, and for this study was taken to be the same at excitation and emission wavelengths.

The setup of these simulations is as follows: In a right handed coordinate system, our phantom occupies the volume [0,8cm]³. The top surface at z=8cm was used as the illumination and detection plane. We assumed a detector with an infinite number of pixels, by providing synthetic measurements at all the quadrature points when numerically evaluating integrals over the measurement surface. Excitation light modulated at 100MHz was delivered to the illumination plane via one of the following schemes: (i) four line sources, (ii) four Gaussian sources, or (iii) a combination of diffractive optics patterns. Figure 3 shows these patterns. With these assumed patterns used as sources S ⁱ (r), we computed fluorescence amplitude and phase across the domain containing fluorescent targets of 5mm diameter spheres filled with 1µM Indocyanine Green solution in 1% Liposyn. The phantom background was assumed to be devoid of fluorophore. As mentioned above, we use different meshes to compute fluences for different illumination patterns; Fig. 4 demonstrates how meshes are generated through adaptive mesh refinement for a single one of each of the three kinds of sources. It is clear that for resolving complex source patterns generated by diffractive optics, adaptive mesh refinement is a necessity since the incident excitation source is poorly resolved during the first couple of refinement stages, and at least 4 adaptive refinements are required. Resorting to a priori global mesh refinement would incur an exorbitant computational burden in the simulation of scanning and patterned area incident illumination.

All of the following computations were performed on a Beowulf cluster with 16 compute nodes. Each node consisted of dual Opteron 2.2 GHz 64 bit processors and 8GB of memory.

3.1. Single target reconstruction

For the first computational experiment, we position a single simulated fluorescent target at a depth of 1cm from the center of the illumination plane. We then run our image reconstruction algorithm using synthetic data generated for this situation.

Figure 5 depicts the reconstructed images obtained for source patterns (i)–(iii). All three illumination patterns are able to recover the embedded fluorescent target. However for the multiple Gaussian source patterns, the recovered target is shifted towards the illumination plane. The upwards shift of the reconstructed target may be attributed to the fact that the target is farthest from the excitation sources (none of the sources illuminates the area immediately above the target, see Fig. 3). Hence, less excitation light reaches the fluorophore, resulting in a lower signal for image reconstruction.

If the target is identified as the volume in which the yield map belongs to the top 10% of the recovered range, then the recovered target size varies with the type of the incident illumination patterns used (see Fig. 5).

 

Fig. 4. State mesh evolution: 5 automatic mesh refinements are depicted for the first illumination source employed in scanning lines (left), scanning Gaussians (center) and diffractive optics patterns (right). Colors indicate arbitrary units of excitation light fluence.

Download Full Size | PDF

 

Fig. 5. Single target reconstructions for scanning lines (left), scanning Gaussians (center), and diffractive optics patterns (right). Black wire-frames represent the true target location and size, while the colored blocks depict the reconstructed target. Top 10% of the reconstructed contour levels of fluorescence absorption are shown.

Download Full Size | PDF

 

Fig. 6. Three target reconstructions for targets at the same depth of 1cm from the illumination plane: (a) a single Gaussian source, (b) scanning line sources, (c) scanning Gaussian sources, (d) diffractive optics patterns. Black wire-frames represent the true target locations and sizes, while the colored blocks depict reconstructed targets. Top 50% of the contour levels of fluorescence absorption are shown.

Download Full Size | PDF

3.2. Multiple target reconstructions

Figure 6 depicts the image reconstruction of three 1cm deep fluorescent targets, performed with synthetic data generated from (a) a single Gaussian excitation source with half a width of 4cm centered on the illumination plane [23], (b) four line sources scanning the illumination plane of the phantom in equidistant steps, (c) four Gaussian excitation sources with half widths of 1cm focused on different parts of the surface, and (d) the four diffractive optics patterns.

The advantage gained by multiple area illumination patterns over a single Gaussian is obvious, as the latter lacks the information content to recover all three embedded fluorescent targets. This consequently shows that multiple area illumination are a necessity for unambiguous image reconstruction when the tissue contains either multiple fluorescent targets, or a distributed fluorescence source.

The performance of multiple source patterns can be further differentiated by placing the three fluorescent targets at different depths. We choose the same targets, but now place them at depths of 1cm, 1.5cm, and 2.0cm. Figure 7 shows the reconstructed images for the threemultiple source patterns. The reconstructed parameter contour level cutoff needed to be reduced to 80% of the maximum level for resolving the three targets. Since excitation light penetration descreases exponentially with depth, the deeper targets do not fluoresce as much as the shallower targets resulting in a nonuniformity in the reconstructed fluorescence absorption in the three targets. The scanning Gaussian source doesn’t resolve the 1.5cm and 2cm deep targets clearly, while the scanning line source illumination results in the detection of only the 1cm and 1.5cm deep targets. Diffractive optical pattern based illumination provides the best image reconstructions as all the targets can be identified and there is no overlap.

Figure 7 also shows the mesh used in the discretization of the unknown parameter q, for the three different illumination situations. Since the recovered parameter should be the same in all three situations, it is not surprising that the parameter meshes are similar although the illumination scenarios are very different. Note that the parameter mesh is refined only near the reconstructed fluorescent targets and remains coarse elsewhere, including close to the illumination and measurement surface, yielding an optimal distribution of unknowns for the reconstruction of the targets. Fluorescence optical tomography is an illposed inverse problem. In the numerical studies reported in the paper, the ill-posedness is further exacerbated by availability of only the reflectance plane measurements. An adaptive refinement based technique successfully recovers the location of embedded fluorescence targets but fluorescence absorption magnitudes are not uniquely quantified resulting in a variation in the reconstructed fluorescence absorption maps depending on the simulation configuration. Hence, for the ease of visualization of results different contour level cutoffs were used for the single target and three target image reconstructions.

The performance of patterned area illumination schemes can be contrasted with the traditional fiber optic based point illumination schemes. In order to assess the comparative performance of patterned illumination with point illumination, four point sources placed symmetrically over the illumination surface were used to generate synthetic measurements corresponding to the three fluorescent target simulations. Target positions and sizes were identical to that reported in the preceeding paragraphs. The software framework detailed in section-2 can be used to model point-illumination schemes by using a suitable mathematical function describing point illumination. Point sources were simulated by employing Gaussian source profiles with a narrow 1.5mm width. Figure 8(a) depicts the source positions. Figure 8(b) depicts the adaptively refined forward simulation mesh for the first source. The image reconstructions corresponding to the three fluorescent targets placed at the same and varying depths are reported in Fig. 8(c) and Fig. 8(d) respectively. The contour level cutoffs employed are identical to the patterned illumination based image reconstructions. The image reconstruction corresponding to the fluorescent targets at the same depth depicts significant image artifacts in the center of the imaged field of view, while the image reconstruction from the measurements generated for the three fluorescent targets at varying depths completely fails to accurately reconstruct the location and sizes of all the three fluorescent targets. The reason for the failure of point sources in locating three fluorescent targets is primarily because of insufficient excitation light penetration and consequently weaker fluorescence emission. This suggests the inadequacy of point illumination schemes for reconstructing fluorescence distributions in large tissue volumes with multiple targets, especially when only a limited number of sources can be employed to limit data acquisition time.

 

Fig. 7. Three target reconstructions for targets at varying depths from the illumination plane: (a) scanning line sources, (b) scanning Gaussian sources, (c) diffractive optics patterns. Left column: Meshes for the reconstruction of the parameter. Right column: Black wire-frames represent the true target locations and sizes, while the colored blocks depict reconstructed targets. Top 80% of the contour levels of fluorescence absorption are shown.

Download Full Size | PDF

4. Quantitative assessment

In the previous section, we have discussed qualitative measures for comparing different excitation patterns. In this section, we demonstrate quantitatively the advantages of using multiple area excitation sources by a series of computational experiments. Simulation studies have a long history in optical tomograpy (see, for example, Ref. [35]) and are ideally suited to determine optimal experimental setups before their implementation in hardware.

We focus on the scanning line setup: The phantom surface is scanned by up to W lines and we quantitate the advantages in image reconstruction with the increments in the number of line sources employed (W). The line sources are positioned symmetrically with respect to the center of the measurement surface (see Fig. 3). In the image reconstruction procedure detailed in the preceding sections, a linear problem defined by Eq. (11) is solved within each non-linear iteration to determine the update to the unknown parameter map q. We note that the once the solution has converged to within noise level, the right hand side of Eq. (11) only consists of contributions due to measurement noise, and further updates δq _k will only produce changes within the range of uncertainty. The matrix S defined in Eq. (14) determines this parameter update δ _q . Since it is symmetric and positive definite we can find a complete set of (normalized) eigenvectors v _ℓ and singular values σ _ℓ such that S=Σ _{ℓ σℓ} v _ℓ $v_{\mathit{\ell }}^T$. By convention the singular values σ _ℓ are arranged in decreasing order. We can then determine the update to be

$\delta q_k={\displaystyle \sum _{\ell }}\frac{1}{{\sigma }_{\ell }}{t}_{\ell }{v}_{\ell },$

where t _ℓ =$v_{\mathit{\ell }}^T$ [F ₂-${\mathrm{\Sigma }}_{i=1}^W$ $C_i^T$ $P_i^{-T}$ (F ₁-M _i $P_i^{\mathit{-}\mathit{1}}$ F ₃)]. Since we want small random updates δ _qk , it is important that the singular values σ _ℓ be as large as possible.

 

Fig. 8. Three target reconstructions with point illumination (a) point sources employed are numbered, (b) final adaptively refined forward mesh for source 1, (c) reconstructed fluorescence absorption map for the three targets placed at the depth of 1cm with a contour level cutoff at 50% of maximum, (d) reconstructed fluorescence absorption map for the three targets placed at depths of 1cm, 1.5cm, and 2cm with a contour level cutoff at 80% of maximum. Black wire-frames represent the true target locations and sizes, while the colored blocks depict reconstructed targets.

Download Full Size | PDF

Most large singular values of the resolution matrix S correspond to image modes v _ℓ that describe features in the solution close to the measurement surface, about which measurements provide enough information. On the other hand, small singular values correspond to high-frequency oscillations and deep objects that are badly constrained by the available measurements. Hence the decay profile of the singular value spectrum of the Gauss-Newton matrix

provides a measure of the ill-conditioning of the inverse image reconstruction problem. Similar analyses have been performed by other researchers to optimize source-detector geometries in both diffuse optical tomography [31–33] and fluorescence molecular tomography [34] for point source based illumination.

 

Fig. 9. Singular value spectra of the Gauss-Newton matrix S defined in Eq. (14) with increasing number of lines.

Download Full Size | PDF

In order to test our hypothesis that more measurements reduce ill-posedness, the matrices S for a varying number W of scanning line sources were computed for the case of a uniform fluorophore distribution (q=0.01cm ^-1). To remove the effect of adaptive mesh refinement on the singular value spectrum, only one iteration on uniform globally refined state and parameter meshes was computed. All the state/forward meshes consisted of 2.5mm cubical voxels, while the parameter mesh consisted of 1cm cubical voxels.

Figure 9 depicts the effect of increasing the number of measurements on the singular value spectrum of the matrix W. The larger the number of measurements the larger the number of singular values σℓ of S with significant magnitude, where the latter is determined by the fact that only singular values greater than the regularization parameter β(typically β=10^-10) impact image formation. Figure 10 shows this result in a different manner by plotting the number of singular values beyond the typical regularization parameter. From both these plots we conclude that by adding additional measurements, we can increase the number of resolvable modes by a factor of 4–5. If these modes were isotropically distributed in the domain, this result would imply that 16 measurements will yield a resolution that is a factor of $\sqrt[3]{4}-\sqrt[3]{5}$ better than when a single measurement is used. In addition, we see that the absolute value of the largest eigenvalue also increases by a factor of at least 10 with the number of measurements, meaning that the uncertainty on the resolvable image modes v _ℓ is reduced by a factor of √10.

Figure 10 indicates that even a further increase in the number of measurements will not increase the number of significant singular values appreciably beyond around 70. This should not surprise: the amount of information that can be extracted from a system is bounded by physical considerations (for example, no features smaller than the scattering length scale can be recovered assuming the validity of the diffusion equation). Even optimally designed experiments will therefore not be able to increase the number of significant singular values beyond limit. On the other hand, it is conceivable that a set of experiments different from the ones considered here may extract more information from the system than shown in Figs. 9 and 10.

 

Fig. 10. Increase of useful singular values σ _ℓ >10^-10 for resolving the reconstructed fluorophore image with increasing number W of lines.

Download Full Size | PDF

5. Conclusion and future implications

In this contribution, we have demonstrated a novel fluorescence optical tomography framework for recovering the location of single and multiple fluorescent targets embedded in tissue, from time dependent, non-contact boundary fluorescence measurements generated by multiple area excitation illumination sources. As with fiber optic based point illumination and detection schemes, multiple area illumination patterns enable wide spatial sampling of the tissue surface, but unlike fiber optic sources and detectors, multiple area illumination techniques generate dense measurement datasets and enhance excitation light penetration in the tissue. As a consequence, we obtain an increase in sensitivity from the patterned modulated illumination.

Multiple area illumination sources can be generated by scanning a simple source geometry, such as a line or a Gaussian, across the illumination plane or by utilizing complex patterns generated by diffractive optical components. The proposed multiple experiment fluorescence tomography scheme implements an algorithm in which each of these measurements is represented by an adaptively refined mesh. Contributions from each experiment are assembled in the nonlinear scheme using a parallel framework that results in image reconstruction times that are equivalent to single excitation source computations on a single machine.

As we have shown both qualitatively and quantitatively, multiple area illumination has distinct advantages over single illumination schemes in the image reconstruction of multiple fluorescent targets.

Adaptive mesh refinement tomography schemes are essential to generate optimal and efficient finite element meshes in clinically relevant situations, wherein large (>100cm ²) tissue surfaces may need to be sampled, and no a priori information about the fluorescent target locations is available. Different illumination patterns perform differently in the challenging image reconstruction problems involving multiple fluorescent targets placed at varying depths. Hence, excitation source patterns need to be optimized for covering large tissue surfaces and simultaneously handling fluorescent target distributions varying both in lateral and vertical directions. Since the numerical computations corresponding to multiple sources are distributed to separate compute nodes of a cluster, the proposed tomography scheme enables us to optimize excitation sources as the number of patterns which can be employed is only limited by the amount of computational resources available. The optimal source patterns for fluorescence optically tomography can be determined by an exhaustive search procedure over the illumination surface pixels by maximizing a suitable measure of the conditioning of the Gauss-Newton hessian matrix S with respect to the variation of illumination patterns. Finally, numerical studies such as those presented in this contribution are significantly faster than the comparable tests using actual measurements in the laboratory, and may therefore accelerate the development and deployment of optimal measurement schemes for fluorescence optical tomography in clinic.