|Home | About | Journals | Submit | Contact Us | Français|
Fluorescence molecular tomography is a powerful tool for 3D visualization of molecular targets and pathways in vivo in small animals. Owing to the high degrees of absorption and scattering of light through tissue, the fluorescence tomographic inverse problem is inherently ill-posed. In order to improve source localization and the conditioning of the light propagation model, multiple sets of data are acquired by illuminating the animal surface with different spatial patterns of near-infrared light. However, the choice of these patterns in most experimental setups is ad hoc and suboptimal. This paper presents a systematic approach for designing efficient illumination patterns for fluorescence tomography. Our objective here is to determine how to optimally illuminate the animal surface so as to maximize the information content in the acquired data. We achieve this by improving the conditioning of the Fisher information matrix. We parameterize the spatial illumination patterns and formulate our problem as a constrained optimization problem that, for a fixed number of illumination patterns, yields the optimal set of patterns. For geometric insight, we used our method to generate a set of three optimal patterns for an optically homogeneous, regular geometrical shape and observed expected symmetries in the result. We also generated a set of six optimal patterns for an optically homogeneous cuboidal phantom set up in the transillumination mode. Finally, we computed optimal illumination patterns for an optically inhomogeneous realistically shaped mouse atlas for different given numbers of patterns. The regularized pseudoinverse matrix, generated using the singular value decomposition, was employed to reconstruct the point spread function for each set of patterns in the presence of a sample fluorescent point source deep inside the mouse atlas. We have evaluated the performance of our method by examining the singular value spectra as well as plots of average spatial resolution versus estimator variance corresponding to different illumination schemes.
Optical imaging techniques have long been popular for generating contrast representing molecular processes in tissue. Over the past decade, their applications have extended from microscopy to macroscopic imaging of deep-tissue optical sources by exploiting the near-infrared (NIR) window where water, oxyhemoglobin, and deoxyhemoglobin, the primary absorbers in tissue, have relatively low absorption coefficients allowing light to penetrate several centimeters inside tissue (Weissleder and Ntziachristos 2003). The non-ionizing nature of NIR light, the availability of a variety of highly specific fluorescent probes, and, finally, low cost offer these methods some leverage over existing radiological techniques for molecular imaging (Hebden et al 1997). However, penetration depths of only a few centimeters are insufficient for most clinical deep-tissue imaging applications. Thus the success of 3D optical imaging methods in clinical diagnostics is restricted to only a few applications. These include NIR spectroscopy and diffuse optical tomography used for mapping brain structure and function in neonates (Koizumi et al 2003, Villringer and Chance 1997, Hebden et al 2002) and optical mammography for breast cancer screening (Ntziachristos and Chance 2001, Culver et al 2003). Optical tomographic methods show great promise in preclinical research, which is a valuable translational tool between in vitro studies and clinical applications. Fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT) have emerged as promising low-cost alternatives to PET and SPECT for functional imaging in small animals, thus greatly impacting diagnostics, drug discovery, and therapeutics (Ntziachristos et al 2005, Gibson et al 2005). Although plagued by tissue autofluorescence, FMT has many advantages over BLT. Compared to most bioluminescent probes, commonly used fluorophores emit light at longer wavelengths (where tissues are less absorbing) consequently offering higher detected signal strengths (Contag and Bachmann 2002, Hielscher 2005). Additionally multiple illumination patterns can be used to generate different mappings from the source space to the detector space making the FMT problem less ill-posed than the BLT problem (Chaudhari et al 2005). The success of FMT can be attributed to a number of breakthroughs.
Despite their tremendous potential and increasing popularity, fluorescence tomographic techniques are confounded by high degrees of absorption and scattering of photons propagating through tissue, making the FMT problem ill-posed. One approach for alleviating this problem and improving source localization is to harness spectral variations of tissue optical properties by using multispectral illumination and/or detection (Zacharakis et al 2005, Chaudhari et al 2009). Another approach is to exploit the degree of freedom offered by external illumination in FMT and design a set of spatial illumination patterns that improve the conditioning of the forward model matrix (Dutta et al 2009), and that is the focus of this paper.
FMT setups typically acquire data-sets corresponding to different surface illumination patterns. These patterns generate different excitation fields over the volume, which tune the system matrix. FMT setups available today employ illumination schemes chiefly guided by the availability and simplicity of the light source. Several of these use laser sources with focusing or diffuser lenses to generate point or distributed patch patterns (Graves et al 2003, Zavattini et al 2006, Li et al 2009). Other approaches include raster scanning (Patwardhan et al 2005, Joshi et al 2006) and structured light or spatially modulated illumination patterns (Lukic et al 2009). Most of these approaches are ad hoc. Although various performance metrics could be used to theoretically compare these standard approaches, it is impossible to make an exhaustive set of comparisons, since, for a given number of illumination patterns being used for an experiment, infinitely many designs exist. Therefore, the question we address in this paper is as follows: given a fixed number of illumination patterns, how do we design these patterns so as to maximize the information in the acquired data?
With the availability of Texas Instruments Digital Light Processor (DLP®) chips (Hornbeck 1996, Dudley et al 2003) which work in the near-infrared range and give us precise control over the spatial intensity distribution, it is feasible to generate any set of spatial illumination patterns with grayscale intensity variation (Gardner et al 2010, Bassi et al 2008, Konecky et al 2009, Bélanger et al 2010). The focus of this paper is to compute the set of illumination patterns for CW FMT that maximize the information content in the data by improving the conditioning of the Fisher information matrix. We formulate our problem as a constrained optimization problem that minimizes a cost function derived from the Fisher information matrix and computes the parameterized set of optimal spatial patterns.
Section 2 of this paper provides a description of the CW FMT problem. The formulation of the optimization problem that generates the optimal set of patterns is presented in section 3. In section 4, we describe the methods used to solve the forward and inverse problems, the optimization procedure, and the performance metrics used to evaluate different illumination schemes. Section 5 presents optimal patterns on a cylinder, a cuboidal tissue phantom, and a mouse atlas along with performance comparisons of illumination schemes for the atlas. Finally, a discussion of the results is presented in section 6.
In continuous wave fluorescence molecular tomography (CW FMT), the 3D biodistribution of the fluorophore inside an animal volume is computed from the steady-state surface fluorescence photon density measured using a CCD camera. With the diffusion approximation, the CW FMT problem can be described by a set of coupled partial differential equations (PDEs) (Hutchinson et al 1995, Milstein et al 2004, Pogue and Patterson 2006):
Here is the diffusion coefficient and μa (r, λ) and , respectively, are the absorption and reduced scattering coefficients at a position r and a wavelength λ. For a surface illumination pattern w(r, λex) at an excitation wavelength λex, the corresponding excitation field, Φ(r, λex), over an animal volume, Ω, can be computed by solving (1). For a known fluorophore distribution q(r) and emission spectral coefficient η(λem), the emission field, Φ(r, λem), can be computed by solving (2) for the previously computed excitation field. In FMT, we assume prior knowledge of the animal surface geometry and tissue optical properties, μa (r, λ) and . The diffusion approximation is valid under the condition . This assumption is valid for most tissue types at near-infrared wavelengths, except near boundaries and in non-scattering parts of some organs (e.g. ventricles in the brain and alveolar sacs in the lungs) (Dehghani et al 1999). Appropriate boundary conditions must be imposed while solving (1) and (2) (Haskell et al 1994, Aronson 1995). The Robin boundary condition, for example, requires that the total inwardly directed photon current at the boundary be zero:
Here denotes the outward unit normal at position r on the boundary, Ω, while G depends on the refractive index mismatch at the boundary (Schweiger et al 1995).
The PDEs (1) and (2) can be decoupled and solved subject to conditions (3) and (4), respectively, by replacing the source terms on the right-hand sides of (1) and (2) by point sources (delta functions) at different locations. In a discretized domain with nd detector nodes on the animal surface and ns point source locations distributed inside the volume, the excitation forward model at an excitation wavelength λex can be expressed as a matrix Aex ns×nd, with the j th column of this matrix representing the excitation field caused by illuminating the j th surface node. Similarly the emission forward model at an emission wavelength λem can be expressed as a matrix Aem nd×ns. The j th column of this matrix represents the surface fluorescence pattern due to the j th internal point source. For a set of p illumination patterns, we discretize the intensity distribution at the surface, w(r, λex), as a set of vectors wk nd, where k = 1: p. The solution to the forward problem is a system matrix dependent on Aex, Aem, and all the wk vectors.
The block of the system matrix corresponding to the kth illumination pattern can be obtained by diagonally scaling the emission forward model matrix by the excitation intensities at each internal point:
Here is a diagonal matrix representing the excitation field due to wk and is given by
where is the ith component of the vector dk computed from
The full system matrix, A pnd×ns, is obtained by vertically concatenating the individual forward model matrices for a set of p different illumination patterns:
Here the prime symbol (′) represents matrix transpose.
The measured steady-state surface fluorescence patterns corresponding to p different illumination patterns can be stacked up as a single data vector, b pnd. The unknown fluorophore distribution, q ns, can be computed by solving the following linear system of equations:
CW FMT setups commonly employ one or more mirrors that make the entire surface of the object (animal or phantom) visible in a single CCD camera image (Li et al 2009, Dutta et al 2008, Chaudhari et al 2005). Such setups require CCD cameras with a high dynamic range. Mapping of this data from the CCD image space to the object space requires knowledge of the object surface as well as calibration for various geometric and radiometric mapping parameters (Li et al 2009). The data vector, b, is obtained by interpolating CCD pixel-wise intensity values onto the surface mesh of the animal or phantom being imaged using the calibrated mapping scheme. Its size, therefore, depends on the mesh density.
We pose the image reconstruction problem as an optimization problem that seeks to minimize a regularized least-squares cost function:
Here is the reconstructed 3D fluorescent source distribution and α is the Tikhonov regularization parameter (Tikhonov and Arsenin 1977). The first term is a data-fitting term, while the second term is a regularization or a penalty term that smooths the reconstructed image. The problem in (10) can be solved analytically using singular value decomposition (SVD) or using a variety of iterative methods, the former being the approach adopted in this paper. Given the SVD, A = UΣV′, of the system matrix, this problem can be solved using the regularized pseudoinverse method:
The system matrix, A, is typically very large in size, with many more rows (pnd) than columns (ns). It is often infeasible to store this large matrix in computer memory. To obtain the solution in (11), it is sufficient to compute A′A and A′b as follows:
Here ‘○’ represents the Hadamard (entrywise) product and and dk can be computed from (6) and (7), respectively. The matrix A′A ns×ns is of the same size irrespective of the number of illumination patterns and/or the number of illumination/detection wavelengths. Its SVD, A′A = VΣ2V′, can be used in (11).
According to the Cramér Rao inequality, the inverse of the Fisher information sets a lower bound on the variance of an unbiased estimator (Kay 1993). This establishes a reciprocity between the variance and the information content of an estimator. For a scalar estimator, the Fisher information is a scalar, and the estimator variance can be minimized by maximizing the scalar Fisher information. However, when we are dealing with a vector estimator, the Fisher information is a matrix of elements, and there is no single scalar optimality criterion. In the field of optimal design of experiments, a variety of different functionals of the eigenvalues of the Fisher information matrix (FIM) are traditionally used as optimality criteria. For example, A-optimality minimizes the trace of the inverse of the FIM, D-optimality maximizes the determinant of the FIM, and E-optimality maximizes the smallest eigenvalue of the FIM (Pukelsheim 1993). The approach that we present here is inspired by E-optimality. Our goal is to tune the FIM so that its singular value spectrum approaches that for the perfectly conditioned case. Our method raises the entire spectrum in reference to the largest singular value. Then, for the same regularization parameter, a greater number of singular values (and the corresponding singular vectors) will affect the solution in the optimal case, leading to improved resolution characteristics.
In principle, an unlimited number of illumination patterns can be used. However, this would lead to lengthy data acquisition times (during which the animal needs to remain anesthetized), large data volumes, slow reconstruction, and possible redundancy in collected data, all of which are undesirable in practice. Therefore, we fix the number of patterns and seek to determine the optimal set for this given number.
For an unknown vector q and an observed data vector b with probability density p(b), the FIM is given by (Kay 1993)
We assume an additive white Gaussian noise model with the noise vector and the data vector b = Aq + n. Without loss of generality, we ignore the constant multiple, , in the covariance. Using (5) and (8), the Fisher information matrix for this system can be expressed as
Using (6), (7), and (15) the FIM can be expressed as a function of the set of illumination patterns, W = [w1w2 ··· wp], where W nd×p. We denote the ith row vector of Aex by and the ith column vector of Aem by . Then the ijth term of the FIM is given by
where we have used the notations . This implies that the FIM can be conveniently represented as a Hadamard matrix product of two matrices, one dependent on the emission model and the other on the excitation model and the illumination patterns:
Introducing the notation Fem = [(Aem)′Aem] and Gex = [AexWW′(Aex)′], we can rewrite (17) as
For the system matrix, A = UΣV′, to have a perfect condition number of 1, the diagonal matrix Σ of singular values should be σ I where σ is a singular value. Note that the singular value σ does not affect the condition number. Then the FIM for the ideal system is
Thus, for perfect conditioning, F should approach Fideal. We would like to find a set of illumination patterns, W, which minimizes the relative difference of the resultant FIM, F (W), and the ideal one, Fideal:
where F (W) is defined in (17) and ||·||F represents the Frobenius norm. The minimizer, (Ŵ, ), of the cost function, Θ(W, σ), with the nonnegativity constraint W ≥ 0 is not unique, since for any γ, Θ(γ Ŵ, γ) = Θ(Ŵ, ). We therefore first solve a simpler optimization problem:
Then we determine the scaling factor by imposing a power constraint, ||W||F = β, on the illumination patterns. The final solution is given by
In practice, the illumination is scaled to exploit the maximum available source power and thus ensure a high signal-to-noise ratio. Since the optimal solution, Ŵ, which minimizes Θ, is simply a scaled version of the solution W opt, we now focus on the simpler problem given in (21).
The formulated problem in (21) is fourth-order and non-convex with respect to the patterns W. W is of size nd ×p, which is typically large, and hence solving this problem is computationally intensive (Dutta et al 2008). Additionally, the solutions are not necessarily spatially smooth. Smooth patterns are convenient for experimental implementation since calibration errors in mapping intensities from the light source to the object surface are likely to affect the accuracy of the forward model less. They are also desirable because an optimal set of patterns computed for a mouse atlas can be meaningfully warped to any given mouse shape, as discussed later in section 6. To allow smooth solutions and to enable speedier computation, we parameterize the illumination patterns using a small number of spatially smooth basis functions. The transformation from the basis function domain to the spatial pattern domain is given by
Here yk m is a vector of the linear coefficients of the basis functions for the kth illumination pattern, L nd×m is a matrix whose columns are the basis vectors, and Y m×p with Y = [y1 y2 ··· yp], m being the number of basis vectors used to encode each spatial pattern. Then, the FIM can be represented in terms of Y as follows:
With the inclusion of basis functions for dimensionality reduction, the modified optimization problem is as follows:
A gradient-based approach can be used to solve the optimization problem in (25). Our cost function as a function of the argument yk, where k = 1: p, is as follows:
where δij is the Kronecker delta function and Ψij is given by
The gradient, k Φ, of the cost function with respect to the kth argument vector, yk, is given by
in (28) can be computed as follows:
Our formulation can be extended for application to multispectral detection. We assume that, for each illumination pattern, fluorescence data are collected over s wavelength bins. The fluorescent dye is assumed to have an emission spectral coefficient ηl at the lth wavelength bin. Then the block of the system matrix corresponding to the kth illumination pattern is given by
The Fisher information matrix for this case is
where . The FIM retains the Hadamard product form with the illumination-dependent matrix Gex unchanged. This implies that, once the multispectral emission model-dependent Fem−ms matrix is precomputed, the optimization problem can be solved for the multispectral emission case without any additional computational cost.
The optimization problem formulated in section 3 can be solved to obtain the optimal set of patterns for a given number of patterns, p. This problem is first solved for an optically homogeneous symmetric solid shape—a cylinder—to gather geometric intuition. Next, we solve this problem for an optically homogeneous cuboidal phantom set up in the transillumination mode with the top face illuminated and the bottom face used for detection. Finally, we solve the same problem for an optically inhomogeneous and realistically shaped mouse atlas. We perform source localization studies using the optimal patterns on the mouse atlas and evaluate the performance of the patterns based on appropriate metrics.
The finite element method (FEM) was used to solve the coupled PDEs (1) and (2) subject to boundary conditions (3) and (4), respectively, to precompute the excitation and emission forward model matrices Aex and Aem.
The cylinder is chosen since it is a regular geometrical shape without corners and elongated like a real mouse. The idea of this study is to ensure that the numerical scheme for computing the patterns preserves the natural symmetries expected in a regular geometrical shape when no spatial constraints have been imposed. For these simulations, we use a tessellated cylinder of height 40 mm and diameter 20 mm with 2373 tetrahedrons and 588 tessellation nodes. A set of 361 detector nodes on the surface are used for illumination and detection and 1564 point source locations on a uniform, volumetric grid with a 1.5 mm spacing are used for source localization. We use monochromatic excitation at λex = 650 nm and monochromatic detection at λem = 730 nm. Optical properties are assumed to be homogeneously distributed and are assumed to resemble those for muscle tissue (Alexandrakis et al 2005).
We simulate a cuboidal phantom of dimensions 40 mm × 30 mm × 20 mm. The phantom is set up in the transillumination mode, a planar imaging mode where the illumination source and the CCD camera are placed on opposite sides of the object (Ntziachristos 2006). The tessellated phantom contains 40 742 tetrahedrons and 7940 nodes. For generating the patterns, we illuminate over the 474 nodes on the top surface (z = 20 mm), detect over the 477 nodes on the bottom surface (z = 0 mm), and use 1155 internal grid points with a uniform 2.5 mm grid spacing for source localization. We use monochromatic excitation at λex = 650 nm and multispectral detection at 710 nm, 730 nm, and 750 nm. We assumed homogeneous optical properties similar to muscle tissue (Alexandrakis et al 2005). We use the emission spectrum of the Alexa Fluor 700 dye, with an emission peak at 719 nm, for computing the forward model.
For more realistic simulations, we use the Digimouse atlas (http://neuroimage.usc.edu/Digimouse.html) (Dogdas et al 2007, Stout et al 2002), a labeled atlas, based on co-registered CT and cryosection images of a 28 g normal male nude mouse. Tissue optical properties are assumed to be spatially inhomogeneous and assigned organ-wise, for 17 different organs, according to published values (Alexandrakis et al 2005). The tessellated atlas volume consists of 306 773 tetrahedrons and 58 244 nodes. For generating the patterns, 810 surface nodes and 1129 internal grid points with a uniform 2.4 mm grid spacing are used. We assume that all surface nodes except those lying on the limbs and the tip of the snout are used for illumination. The Alexa Fluor 700 dye (emission peak 719 nm) is used as the fluorophore. Accordingly, we use an excitation wavelength of λex = 650 nm and perform multispectral detection at 710 nm, 730 nm, and 750 nm.
Since our cost function in (26) is continuous and its gradient is computable using (30), we adopt a gradient-based approach for optimization. We use fmincon, an inbuilt function for constrained optimization in MATLAB® (The Mathworks Inc., Natick, MA, USA). This function chooses the active set method for optimization (Gill et al 1982). The non-convex nature of the cost function implies that several local minima could exist, and our gradient-based optimization approach could get stuck in one of these. Hence, we use a multi-start approach with random initializations. This is a heuristic but straightforward approach to the global minimization problem in which the local minimization algorithm is run from many different starting points, and then the best solution is picked (Liberti and Maculan 2009).
For dimensionality reduction, we define basis functions on the cylinder and atlas surfaces by treating them as smooth 2D manifolds. We use as our basis set eigenfunctions of the Laplace–Beltrami operator (Qiu et al 2006, Joshi 2008, Grenander and Miller 2007). The Laplace–Beltrami operator is a generalization of the Laplace operator on Riemannian and pseudo-Riemannian manifolds. These basis functions allow us to generate a spatially smooth solution. Their orthonormality allows efficient representation of the solution. Since the eigenvalues of a spatial differential operator are measures of the spatial rate of change of the corresponding eigenfunctions, the smaller the eigenvalue the smoother is the corresponding eigenfunction. Therefore, eigenfunctions corresponding to the m = 20 smallest eigenvalues of the operator were used as basis functions. Figure 1 shows front and back views of the tessellated cylinder, top views of the cuboidal phantom, and top and bottom views of the Digimouse atlas with the basis functions displayed on them.
The inverse problem is solved using the regularized pseudoinverse in (11). The illumination patterns on the Digimouse surface are interpolated onto a denser surface mesh with 3234 nodes. For source localization studies, a uniform grid of 9192 point sources with a 1.2 mm spacing is used. Our simulation setup uses surface data from all nodes on the atlas surface except those lying on the limbs and the tip of the snout.
In order to comparatively assess different illumination schemes, we need to establish appropriate performance metrics. We use two approaches for comparing the performance of different illumination patterns. The first approach is to examine the conditioning of the system matrix corresponding to a specific illumination pattern by looking at its singular value distribution. The second approach is to look at average resolution–variance curves to examine the behavior of the inverse solution.
The condition number of a matrix is given by the ratio of its largest to its smallest singular value. Although this is a good figure-of-merit for the robustness of the system, the point spread functions (and, hence, the reconstruction results for any source distribution) depend on the singular values as well as the singular vectors (Chaudhari et al 2005). So in addition to looking at the condition number, we analyze properties of the point spread functions corresponding to different source locations inside the animal volume.
The mean squared error of an estimator can be decomposed as a sum of the squared bias and the variance. As the regularization parameter is increased, the estimator variance decreases while the estimator bias increases (resolution worsens) and vice versa. Resolution–variance curves, which signify this inherent tradeoff, are commonly used to assess image reconstruction quality (Qi and Leahy 1999, Meng et al 2003, Chaudhari et al 2009). These curves are obtained by sweeping the regularization parameters over a range of values, computing the variance and resolution of the point spread functions, and plotting these two quantities against each other for the different values of the regularization parameter. A lower lying curve implies superior resolution for the same noise variance for the estimator corresponding to an illumination scheme. Since the regularized pseudoinverse operator in (11) is a linear operator, it can be completely characterized in terms of its point spread functions.
When the regularized pseudoinverse method is used for reconstruction, the spatial resolution and noise variance of the estimator can be computed in a closed form, thus eliminating the need for computationally expensive Monte Carlo simulations. For a true source distribution q, the mean value of the estimator for a regularization parameter α can be expressed in terms of the SVD, A = UσV′, of the system matrix as
We would like to compute the spatial resolution from the point spread function, PSFj, for the j th unit point source. For this source, q = ej, where ej is a unit vector with 1 as the j th element. Then PSFj can be computed by substituting q = ej in (33):
The spatial resolution is measured as the full width at half maximum (FWHM) computed from PSFj. Then the average spatial resolution over n different point source locations inside the object is a function of the regularization parameter, α:
For an additive white Gaussian noise model with the covariance , the estimator variance for the j th unit point source, q = ej, is given by
For fair comparison of illumination schemes involving different numbers of illumination patterns, we assume that the total acquisition time is always the same irrespective of the number of patterns used. This implies that, for a scheme with p patterns, the acquisition time per pattern is p times shorter than that for a scheme with 1 pattern. As a result, the individual data-sets are noisier by a factor of p. In other words, the data noise covariance is . With this compensation for acquisition time, the estimator variance averaged over n point source locations is given by
Using (35) and (37), the resolution and the variance for any illumination scheme (and hence any system matrix) averaged over different point source locations inside the volume can be computed for different values of α. The curves of resolution versus variance are obtained by plotting the two quantities against each other for different values of α. Comparison of these curves for different illumination patterns throws light on the relative magnitudes of mean squared error for estimators corresponding to different illumination patterns. It must be noted that (35) and (37) depend on only V and Σ and not on U. Therefore, as described in section 2.2, it is sufficient to compute the SVD A′A = VΣ2V′.
By solving the optimization problem described in section 3, one can determine an optimal set of illumination patterns for a given number of patterns, p. We generated optimal sets of patterns for the tessellated cylinder, the cuboidal phantom, and the Digimouse atlas described in section 4.1. Optimal sets of patterns for the atlas were generated for different values of p. The performance of these sets of patterns was then comparatively evaluated using the metrics described in section 4.3.
For the study using the tessellated cylinder, we set the number of patterns, p, to 3. The resulting set of three optimal patterns computed for the optically homogeneous tessellated cylinder is shown in figure 2. It is assumed that only the curved surface is used for illumination. The patterns are displayed on the curved surface by laying it out flat. We observe that the pattern (b) can be shifted by 120° (one-third of a full circle) right or left to obtain patterns (a) and (c), respectively. The sum of the three patterns in (d) looks quite uniform overall and additionally indicate rough bilateral symmetry about the horizontal plane cutting the cylinder lengthwise midway. Thus the optimal set of patterns exhibits radial and bilateral symmetries, characteristic of the cylindrical shape.
We generated a set of p = 6 patterns for the cuboidal phantom placed in the transillumination mode with illumination from the top and detection from the bottom. The resulting set of optimal illumination patterns on the top surface of the phantom is shown in figure 3. For reference, we generated another transillumination setup where 6 evenly spaced points on the top face (z = 20 mm) are illuminated. The point locations are (10, 10, 20), (20, 10, 20), (10, 20, 20), (20, 20, 20), (10, 30, 20), and (20, 30, 20), where all positions are in millimeters. The singular value spectra for the optimal scheme and the point illumination scheme compared in figure 4 show that the former is more well conditioned than the latter.
Optimal sets of patterns for the Digimouse atlas were generated for different numbers of patterns, p. We picked the values p = 1, 3, 6, and 12 and looked at the corresponding optimal sets of patterns. The average computation times for sets of 1, 3, 6, and 12 patterns were 5 min, 16 min, 2 h, and 20 h, respectively, on a 2.33 GHz quad-core Intel Xeon® processor. The optimal illumination schemes for the chosen numbers of patterns are shown in figure 5. For each set of patterns, the intensities are normalized with respect to the maximum intensity value for that set for display purposes. Unlike the cylinder, the mouse atlas is of an optically inhomogeneous, irregular shape. Accordingly, the patterns do not exhibit any clear symmetry.
We reconstructed a sample deep source inside the trunk of the mouse at a depth of 10.4 mm below the top surface and 7.6 mm from the bottom surface. The PSFs were generated using (34) for the optimal sets of patterns computed for p = 1, 3, 6, and 12. Figure 6 shows coronal slices of the Digimouse atlas displaying the point spread functions for the optimal illumination schemes for different values of p. For reference, we have studied a set of p = 3 patterns each of which uniformly illuminates one of three different longitudinal sections of the mouse. To generate these PSFs, we picked values of the regularization parameter that equalize the estimator variances as computed using (37) for the different illumination schemes. These values were 4.90×10−11, 3.26×10−12, 1.23×10−13, and 8.35×10−16, respectively, for optimal illumination schemes with p = 1, 3, 6, and 12, and 10−13 for the p = 3 uniform illumination scheme. We observe improvement in the measured FWHM from the cases p = 1 through 12. The FWHMs of the PSFs for the different cases are provided in table 1. Also, figure 6 and table 1 indicate that the reconstructed point source for the p = 3 optimal case appears to have better resolution than the p = 3 uniform case for the same estimator variance.
We used the performance metrics described in section 4.3 to investigate the benefit of having larger numbers of patterns in the optimal illumination scheme and to compare the optimal scheme for a set of three patterns with the uniform illumination scheme described in section 5.3.1.
Figure 7 shows the singular value distributions for the different illumination schemes. The singular values are plotted in descending order and normalized with respect to the largest singular value. The condition numbers of the system matrices are provided in table 2. The condition numbers clearly improve as p is increased from 1 through 12. The optimal illumination pattern for p = 3 is observed to generate a better conditioned system matrix than the p = 3 uniform illumination scheme.
The plots of average resolution versus average variance for the different illumination schemes are shown in figure 8. These are obtained by averaging resolution and variance for n = 250 randomly picked sources on the uniform source grid for different values of the regularization parameter. The regularization parameter was swept over values within the range 10−30–10−4. The best possible FWHM is 1.2 mm, which is the grid spacing. As explained in section 4.3.2, the idea underlying these curves is the trade-off between resolution and variance, which implies that, a decrease in estimator variance is accompanied by an increase in FWHM and vice versa. The intensities of the different sets were normalized to ensure that all sets of patterns have the same average intensity. As expected, as p goes from 1 to 12, the curves lie lower. In other words, the resolution–variance properties improve. Also, figure 8 indicates that, for the same average variance, the p = 3 optimal illumination pattern offers better average resolution than the p = 3 uniform illumination scheme.
We have developed an optimization framework for generating optimal spatial illumination patterns for CW FMT based on an approach that seeks to improve the condition number of the Fisher information matrix. We formulated our problem as a constrained optimization problem which, for a given number of patterns, can be solved to compute the set of optimal patterns which maximize the information content in the acquired data. Our formulation assumes an i.i.d. Gaussian noise model. The fourth-order, non-convex cost function was minimized in the presence of a non-negativity constraint using a gradient-based approach. A multi-start method with random initializations was used to ensure global convergence. The dimensionality of the optimization problem was reduced using a set of geometrical basis functions defined on the 2D manifold representing the surface of the object being imaged. Eigenfunctions of the Laplace–Beltrami operator were used as basis functions owing to their spatial smoothness and orthogonality. To ensure that this numerical approach is geometrically meaningful, we looked at optimal patterns generated for an optically homogeneous tessellated cylinder and observed radial and bilateral symmetries that are intrinsic to the cylindrical shape. We then applied this method to compute optimal patterns for an optically homogeneous cuboidal phantom in the transillumination mode. We compared this scheme with an evenly spaced point illumination scheme by examining the singular value spectra and showed that the optimal scheme is more well-conditioned than the point illumination scheme. Finally, we used our method to compute optimal patterns for the optically inhomogeneous and realistically shaped Digimouse atlas. The obtained patterns look smooth and physically realizable. For the same number of patterns, p = 3, the optimal set was shown to perform better than a uniform illumination scheme on the basis of condition number and average resolution versus variance curves. Since typical FMT experimental setups use a larger number of illumination patterns, we have generated optimal patterns for up to p = 12 and presented their singular value spectra and average resolution–variance characteristics. For fair comparison, the variance was computed based on the assumption of equality of total data acquisition times for different values of p.
In our simulations, we have assumed that the entire surface of the mouse except for the limbs and the tip of the snout is being illuminated. However, in an experimental setting, it may not be feasible to illuminate the entire mouse surface at one go. For such applications, our formulation can be easily modified to introduce spatial constraints on the illumination patterns. Another useful modification of our formulation would be its extension to include multispectral excitation. Owing to the larger variability of tissue optical properties at the excitation wavelengths of commonly used fluorescent dyes and proteins, for the same number of wavelength bins, multispectral excitation typically leads to a better conditioned system matrix than multispectral detection (Chaudhari et al 2009). We, therefore, believe that optimal illumination coupled with multispectral excitation will allow us to pack more information into the acquired data for a given number of illumination patterns and a given number of wavelengths.
Currently, the main limitation of our method is the large computation time required to generate large numbers of patterns. We, therefore, would like to explore more efficient numerical implementations that will enable us to compute larger sets of patterns. This will also allow us to investigate the benefits of using large (say > 20) sets of patterns and, in the process, determine the minimum number of patterns that can be used without making the acquired data overly redundant.
The framework for generating optimal patterns assumes prior knowledge of mouse surface topography and tissue optical properties and requires the forward problem to be solved prior to the optimization procedure. It might not be feasible to repeat the optimization procedure for each animal in between the surface profiling and fluorescence data acquisition steps of an experiment. An alternative approach would be to use the atlas as a surrogate and to warp its surface to match that of the target mouse being imaged (Joshi et al 2009), and, in doing so, we can warp the optimal patterns onto the surface of the target.
This work was supported by the National Cancer Institute under grants R01CA121783 and R44CA138243.