1. Introduction

Whole-body small animal fluorescence tomography (FT) is a pre-clinical imaging technique relying on the acquisition of diffuse near-infrared (NIR) light signals transmitted through the body which are used for reconstruction within the framework of an inverse problem [1]. Three-dimensional (3D) images are reconstructed by minimizing the difference between acquired measurements and theoretical measurements defined by analytical or numerical solutions derived from a light transport forward model using boundary conditions consistent with the geometry of the interrogated specimen. Five important factors affect image resolution and recovered contrast of FT images: noise, forward model approximations, detection geometry, anatomical priors, and data types. Understanding the interplay and tradeoffs between these factors within the confines of currently available photodetection instrumentation is central to system optimization. This study looks at an approach to compare the effects of these factors upon the tomography image quality.

The majority of research and commercial FT systems use either cooled charge-coupled device (CCD) cameras or photomultiplier tubes (PMT) for photo-detection, each associated with its own benefits and drawbacks. CCD cameras provide a very large number of optical projections; however, this wide-field detection suffers from intrinsic limitations such as the relatively small dynamical range and the photo-detection sensitivity, both related to the read-out noise. Conversely, PMT-based focused-detection systems are not as affected by read-out noise issues and therefore may be more sensitive to low photon fluxes and capable of a larger dynamical range. However, because detection is usually achieved one pixel per detector at a time, the optical projections attainable with this approach are typically much fewer than with CCDs. Because of this, acquiring large data sets requires either longer scan times or the inclusion of many detection channels, which can increase the price of a system beyond what is reasonable for biological sciences laboratories. Considering the tradeoffs of both PMT and CCD approaches, an in-depth study of the relative affects of noise and projection geometry on image resolution and recovered contrast was developed here in the context of potential anatomical prior and forward model approximation errors.

Previous studies have presented the issue of source-detector geometry optimization by examining the behavior of singular values associated with a specific forward model matrix. This paper re-examines this issue by introducing novel mathematical criteria that can be used to evaluate the merit of an imaging geometry based on singular vectors (both data and image singular vectors) as a way to complement the information provided by singular values. First, a metric is introduced that can be used to assess the spatial resolution of an imaging geometry based on the spatial frequency of the image singular vectors. Second, the use of data singular vectors is considered as a means to evaluate how stochastic noise and model mismatch errors contribute distinctly to the degradation of FT images. These novel singular vector metrics are then used to evaluate the spatial resolution and quantification potential of different imaging geometries, including focused and wide-field detection configurations. The simulation results presented in this work provide evidence that in the context of marginal noise, comparable levels of contrast recovery can be achieved for both focused and wide-field detection, even with sparsely sampled focused detection. Furthermore, higher density measurements do not necessarily translate into significant gains in terms of quantification. The motivation for the work is to provide a thorough theoretical method to analyze and optimize imaging geometry and system design for a focused detection tomography system to provide maximal sensitivity without limiting spatial resolution [2,3].

2. Background material

2.1 Imaging system and geometry

The mathematical analysis and simulations presented in this work are based on the detection geometry of a non-contact diffuse fluorescence tomography instrument, allows specimens to be interrogated with multiple laser projections [3]. In its current implementation, the system is co-registered with a microcomputed tomography (microCT) system allowing delineation of the specimen surface to be retrieved and used for light transport modeling using finite elements methods (FEM). Figure 1 presents a schematic of the generalized detection geometry considered in this work. The detection of filtered fluorescence light is achieved in transmission using a number N_d of photodetectors, separated by an angle θ, and opposing the point of illumination. The laser and photodetectors are fixed on a gantry that can rotate 360° around the specimen with an angular separation providing a number N_p of laser projections. In this theoretical study, fluorescence as well as transmission (excitation) measurements were simulated, allowing reconstructions to be performed based on fluorescence-to-transmission ratio data (Born normalization) [4] for a range of focused-detection data collection schemes (Table 1 ). A wide-field geometry (i.e., CCD detection) was emulated by minimizing θ and maximizing N_d.

 

Fig. 1 Optical detection geometries evaluated for diffuse fluorescence tomography simulations. In the limit that the number of detectors, N_d, is large and the angle between adjacent detection points, θ, is small, wide-field detection is achieved.

Download Full Size | PDF

Table 1. Technical specifications associated with the five representative imaging geometries for which the performance is compared. The last column in the table shows the numerical rank of the forward model fluorescence matrix associated with the geometries. The number of projections is in reference to the system described in Section 2.1, and corresponds to the number of laser positions used (manipulated by rotating the gantry). In other words, the number of measurements per 2D slice is equal to the number of projections times the number of detectors.

View Table | View all tables in this article

2.2 Forward modeling and image reconstruction

Diffuse fluorescence tomography relies on the resolution of an inverse problem of the form (see, e.g., Ref [5].),

where, for continuous-wave (CW) imaging, D_ρ is a column vector with N _m = N _d × N _p entries corresponding to the fluorescence-to-transmission ratio measurements obtained with the representative geometry shown in Fig. 1. The column vector C _a (a = 1,2,…,N _v) is composed of entries representing the concentration (e.g., in units of Molars or μg/ml) of fluorophores in each of the N _v voxels contained in the discretized volume (e.g., an FEM mesh) associated with the interrogated specimen. Finally, the N _m × N _v matrix A is derived from the light transport model of choice and relates the CW measurements vector to the fluorophore concentration vector as shown in Eq. (1). If the diffusion approximation to the radiative-transport equation (RTE) is used, the forward model matrix takes the form,

where Q _F is the quantum efficiency of the fluorophore, ε _F the extinction coefficient and τ the lifetime. The vectors r _s, r _d and r _a (a = 1,2,…,N _v) represent the location of a laser source point (labeled s) projected on the specimen, a detection point (labeled d) also projected on the specimen and the position of a voxel inside the interrogated medium, respectively. The functions Φ^x and Φ^e are fluence fields at the excitation and the fluorescence emission wavelengths, respectively. These fields are computed by solving the diffusion equation on a discretized mesh either numerically, or analytically in cases where simple imaging geometries are considered. Fluorescence tomography involves the retrieval of the vector C_a in Eq. (1) − a volumetric fluorescence image − through the resolution of an inverse problem. An analysis of the inverse problem for guiding source-detector geometry optimization is discussed in the following section.

2.3 Singular value decomposition analysis for imaging geometry optimization

Singular value analysis (SVA) is a tool that has been used in the past to derive mathematical criteria allowing relative a priori assessment of the performance of various imaging geometry setups in diffuse optical tomography (DOT) as well as in diffuse FT [6,7]. Thus far, the tool has predominantly been limited to the study of the singular values. Here, the evaluation is extended to a more detailed study of the singular mode vectors and how their behavior can translate into tomography features such as spatial resolution.

In singular value analysis, the inverse problem in Eq. (1) is solved by performing singular value decomposition (SVD) of the matrix A into two sets of singular vectors and a finite number of singular values (Eq. (2).1) in Ref [8].),

where, for N_v > N_m, U = (c₁, …,c_Nm) ∈ R^{Nm x Nv} and V = (d₁, …,d_Nv) ∈ R^{Nm x Nv} are orthonormal matrices, and where the diagonal matrix Σ = diag(s₁, …,s_Nm) has non-negative elements appearing in non-increasing order (s₁≥s₂≥…≥s_m≥0). A least-squares solution to the problem then takes the form,

where s_i are the singular values, d_i the singular data (SD) vectors and c_i the singular image (SI) vectors. Each singular value is associated with one SD vector and one SI vector, all of which are labeled with the index i running from 1 to N, where N is the maximal number of modes included in a reconstruction. Equation (4) represents a truncated singular value decomposition (TSVD), for which N ≤ N_rank. N_rank is the numerical rank of the matrix A, which is demarcated by the mode at which a significant and abrupt decrease in singular values is observed [8]. The singular value solution in Eq. (4) is roughly similar to the projection of a signal onto different spectral modes in Fourier analysis. According to this analogy, the vectors c_i are the Fourier modes, the coefficients F_i are regularization filter factors allowing for the relative contribution of the spectral components to be weighted, and the coefficients,

are factors weighting the contribution of each of the different spectral components according to how much overlap there is between the SD vectors and the tomography data vector D.

3. Methods

Evaluations of focused-detection and wide-field imaging geometries were performed on five potential imaging geometries described in Fig. 1 and Table 1. Geometries I to IV correspond to potential focused-detection applications, while geometry V corresponds to a wide-field application. As mentioned in Section 2.1, the focused-detection geometries were chosen as potential variants of an in-house developed non-contact diffuse FT instrument that allows specimens to be interrogated with multiple laser projections [2,3]. A rigorous analysis of all given imaging geometries was preceded by and completed in the context of an investigation into the relative contributions of forward model mismatch and stochastic noise to perturbations in the raw data set (Section 4.1). This was followed by an investigation into the characteristics of both singular value (Section 4.2) and novel singular vector analyses (Section 4.3) for evaluation of imaging geometries. Finally, both the model mismatch/noise analysis and the singular value/vector analyses were evaluated in the context of image quality, specifically, spatial resolution and contrast recovery (Section 4.4). All analyses were carried out assuming known anatomical priors for FEM boundary conditions. The benefits of Born normalization, a commonly used technique for improving fluorescence data reconstructions by normalizing to transmittance [4], were also investigated for completeness. In the past, fluorescence-to-transmission data sets have predominantly been shown to improve reconstructions by partially negating the impact of experimental factors, such as out-of-focus data acquisition for non-contact instruments. However, the emphasis here is on the further benefits this normalization scheme provides in terms of reducing the burden on precise light transport modeling in FT, a phenomenon referred to as model mismatch throughout this work.

3.1 The difference between model mismatch errors and stochastic noise

Analysis of the effects of model mismatch and stochastic noise on raw data vectors was carried out for geometry IV applied to a mouse head mesh created from a microCT data set acquired in vivo (Fig. 3c below). For fluorescence tomography simulations, two 4 mm-diameter inclusions were inserted with a 5 mm separation between their center-of-masses, as shown in Fig. 3c (right-most image). The fluorescence contrast between the inclusions and the homogeneously fluorescent background was set to 10:1. For all geometries, data vectors were simulated using two-dimensional diffusion theory using the simulation software NIRFAST [9], with an FEM mesh containing N_v = 1979 nodes. In this first study, five data vectors were simulated: three fluorescence data vectors, two of which were simulated assuming an optically homogeneous background with and without 10% stochastic (Gaussian) noise added and one assuming an optically heterogeneous background and 0% noise. The last two simulated data vectors were transmittance vectors, one assuming an optically heterogeneous medium and one assuming a homogeneous medium (noise-free). The left-most image in Fig. 3c illustrates the different anatomical regions imposed to provide heterogeneity (1: brain, 2: bone, 3: other soft tissue). Table 2 summarizes the optical properties (absorption and reduced scattering coefficients) that were assigned to each region in the scope of two different types of simulations: heterogeneous and homogeneous media.

 

Fig. 3 (a) Simulated data vectors associated with the distribution of inclusions shown in (c) (right-most image) for which a 10:1 fluorescent contrast with respect to the background has been assigned. Two of the curves in (a) represent raw fluorescence data sets, with one curve being associated with an homogenous medium (10% stochastic noise added) and the other one being associated with an heterogeneous medium (no stochastic noise added) where different optical properties (absorption and reduced scattering coefficients) were associated with the different regions shown in (c) (left-most image). The third curved in (a) (normalized by multiplying with the homogenous transmittance data set) represents the fluorescence-to-transmission ratio data set. (b) Graphs showing the percentage difference between a fluorescence data set generated for a homogenous medium (0% noise) and the data sets shown in (a).

Download Full Size | PDF

Table 2. Optical properties (absorption and reduced scattering coefficients) used for the light transport simulations presented in this work

View Table | View all tables in this article

3.2 Singular value decomposition for optimizing imaging system geometries

The theoretical formulation of singular value and singular vector analysis is introduced in Section 2.3. Typically, only the singular values are used as figures of merit for assessing imaging geometries [7,10]. In this study, characteristics of singular values and both SI and SD vectors were investigated using simulated data sets from the homogeneous mouse head mesh described in Section 3.1 and for the imaging geometries summarized in Table 1. For comparison with other studies investigating singular values, the magnitude of the singular values was plotted against mode order and the number of available modes in the presence of stochastic noise and model-mismatch was evaluated.

To provide a more rigorous analysis of the imaging geometries, the shapes of both the SI and the SD vectors were also investigated as a function of mode order for all geometries. The spatial frequency of the SI vector modes typically increases with the order, i, of the singular order (Fig. 2 , Fig. 5 ); whereas, the shape evolution of the SD vectors follows a more complex pattern with increasing mode order (Fig. 6 ). The resolution (inverse of spatial frequency) of each individual SI vector mode was defined as the average distance, in millimeters, between two consecutive extrema (maxima or minima) along one-dimensional cross-sections of SI vector images. An example for this is the dotted line shown in Fig. 2a. More precisely, a numerical value has been assigned to each two-dimensional SI mode (Fig. 5) corresponding to the number of extrema, N_ex^SI, along a one-dimensional cross-section divided by the length of the cross-section, L (in millimeters),

 

Fig. 2 (a) Singular image (SI) vectors associated with increasing singular order. Mode images are plotted on a log-scale to convey the idea that increasing the order of the mode typically leads to spatial frequency increases. (b) Cross-section of the modes shown in (a) along the dotted line.

Download Full Size | PDF

 

Fig. 5 Figure showing the metric labeled mode resolution as a function of singular mode order. The resolution of each SI mode is defined here as the average distance (in mm) between two consecutive maxima along the dotted line shown in Fig. 2. Higher spatial frequency is therefore equivalent to smaller mode resolution. The red dots represent (for each geometry) the maximum mode order (minimum mode resolution) that can be used to reconstruct an image when 1% stochastic noise is added to the simulated data set (see Table 3). The dotted line should be used as a visual guide indicating the approximate location of the noise floor threshold for all geometries.

Download Full Size | PDF

 

Fig. 6 (a) Representative singular data (SD) vectors associated with imaging geometry IV (see Table 1). Inspection of the modes of order i = 1, i = 10 and i = 20 show that the frequency (as a function of the source-detector pair index label measurement number) increases as a function of the singular order i. (b) Frequency of SD modes as a function of singular mode order for the five imaging geometries presented in Table 1. The curves demonstrate that the frequency, Eq. (6), tends to increase monotonically within a certain range, but that phases of sudden decrease always occur for large values of i. Red dots on the curves represent the order from which modes cannot be used for reconstruction because of noise propagation (see Table 3).

Download Full Size | PDF

The number of extrema was computed using an algorithm that automatically localizes and counts the total number of minima and maxima along cross-sections of the SI modes. As an example, Fig. 2b shows three such curves where the number of extrema are N_ex^SI = 3, N_ex^SI = 5 and N_ex^SI = 9 for SI modes of singular order i = 1, i = 5 and i = 30, respectively. Although Fig. 5 has been plotted using the one-dimensional cross-section (dotted line) shown in Fig. 2, numerical analysis (not presented) was performed demonstrating that the final conclusions derived from the mode resolution analysis in Section 4 are unchanged for different cross-sections. For this metric, favorable detection geometries were determined as those for which the modes available for reconstruction had the smallest possible spatial resolution.

Each SD vector (e.g., see Fig. 6a) was also characterized by a numerical criterion, in this case labeled mode frequency (Fig. 6b). The frequency for each SD mode was computed by evaluating the number of extrema N_ex^SD for each SD mode and dividing this number by the total number of measurements N_m,,

The number of extrema, N_ex^SD, was automatically computed using an algorithm similar to that described above for the SI vectors. While the SI mode resolution is an effective measure for the level of spatial resolution each image mode can contribute to a reconstructed FT image, the SD mode frequency does not lend itself to such a direct interpretation. However, as explained in Section 4, the numerical value of the SD mode frequency can be used as an indirect assessment of the propagation of stochastic noise and model mismatch error into FT images.

For the FT image reconstructions presented in this work, the number of usable modes per imaging geometry was determined by using the singular value decomposition to reconstruct images by sequentially including higher order modes in the presence of different levels of stochastic noise. This method, referred earlier, is defined as the TSVD [8]. The usable mode cutoff was empirically defined as the point when further inclusion of higher order modes did nothing to improve the image quality. The use of filter functions F corresponding to a Tikhonov regularization was also considered but did not provide significant image quality improvements.

Other aspects of the SVD analysis used in the a priori assessment of detector geometries included the steepness of the singular value vs. mode order and the number of singular modes above the noise floor (i.e., the modes available for image reconstruction).

3.3 Image quality: spatial resolution and recovered contrast

To evaluate the concepts put forth in the two previous subsections in the context of image quality – e.g., spatial resolution and recovered contrast – images were reconstructed from the simulated data sets of all detection geometries (Table 1) under various levels of noise and with respect to forward model mismatch. For each of these images, an assessment of quality was determined, favoring images with: (1) higher contrast recovery (with respect to the target simulated image), (2) increased ability to spatially resolve multiple targets, and (3) the absence of image artifacts. These a posteriori evaluations are then discussed in the context of the a priori evaluations introduced in Section 3.2.

4. Results and discussion

4.1 The difference between model mismatch errors and stochastic noise

Figure 3a displays the fluorescence data vectors described in Section 3.1 as a function of measurement number for the homogeneous medium with 10% noise and for the heterogeneous medium with no noise added. The third curve displayed is the fluorescence-to-transmittance (Born ratio) for the heterogeneous medium. The Born ratio curve was normalized by multiplication with the simulated transmittance data set associated with the homogeneous optical properties in Table 2. Based on the similarities between the normalized data set and the data set of the homogeneous medium data set, the results suggest that if unknown heterogeneous optical properties are suspected, a homogenous forward model may still be applicable as long as Born normalization is applied. This is further demonstrated in Fig. 3b.

Figure 3b displays the percentage difference between the three vectors plotted in Fig. 2a and the simulated fluorescence data vector for a homogeneous background with no added noise (not shown in Fig. 3a) as a function of measurement number. The differences between the heterogeneous and homogeneous data sets are seen to be up to 15-20% and the measurement-to-measurement variation is significant. In addition to corroborating the observations made about Fig. 3a, this result more precisely suggests that reconstructing a heterogeneous data set using a homogeneous forward model will lead to imprecise reconstructions in terms of center-of-mass localization, as well as potentially decreasing spatial resolution and contrast recovery. Again, befitting the preliminary conclusion drawn from Fig. 3a, the minimal difference between the Born curve and the homogeneous curve suggests that reconstructions based on normalized data sets are less reliant on the veracity of the forward model, which often inaccurately represents the imaging object: in vivo values for different organs are typically not available and in cases where they would be available, significant inter-specimen variations are expected. Additionally, different systems use different excitation and emission wavelengths, preventing use of generic values in the scope of routine experiments, and accurate co-registration with anatomical imaging modalities is difficult but required.

Figure 3b also displays the differences between the noise and noiseless homogeneous vectors (labeled homogenous medium - 10% noise), illustrating the characteristics of noise propagation into the data sets. As expected, the difference between the two vectors consists of stochastic noise normally distributed around zero, with a maximum amplitude of 10% of the fluorescence measurement amplitude. This curve demonstrates that sources of noise such as photon shot noise will always contribute high-frequency components to a data vector. This is in opposition to the predominantly low frequency contributions associated with model mismatch. Therefore, the effect of both stochastic noise and model mismatch on fluorescence tomography can be empirically written as,

where D_homo is the vector that would result for an homogenous medium while D_high-f and D_low-f are vectors corresponding high-frequency (stochastic noise) and low-frequency (predominantly model mismatch) contributions. This parametric classification of data vector contributions is clearly not exact; however, it will be conceptually useful in the upcoming discussion relating to singular value decomposition of the matrix A. It should also be noted here that model mismatch contributions associated with discretization in Eq. (2) will mainly contribute high-frequency components to the data vector in situations where the size of the individual mesh components is comparable or smaller than typical light mean free paths. On the other hand, the fact that most of the low frequency contributions is associated with model mismatch errors is a direct result of the diffusive nature of the light transport sensitivity functions, which leads to averaging of signals emerging from large regions within the interrogated medium.

4.2 Singular value analysis of image geometries

Figure 4 displays the singular value vs. mode order curves associated with the five imaging geometries described in Table 1 and the red dots indicate the maximum usable mode order as described in Section 3.2. The maximum number of modes used for different reconstruction scenarios are tabulated in Table 3 . The truncation value corresponds to the largest singular mode order where stochastic noise does not corrupt the reconstructed images. The values used in Fig. 4 (red dots) correspond to the three columns associated with a 1% stochastic noise level. Detection geometries are usually deemed favorable when the slope of the singular value curve is as flat as possible and when the number of singular values or modes above a certain noise threshold is maximized [7]. In other words, conventional wisdom in the fluorescence tomography community is that geometries providing a large value of N in Eq. (4) are superior because more modes can be used to reconstruct images. Based on those criteria, inspection of Fig. 4 indicates that the wide-field detection geometry (V) is superior compared to the other geometries. It can also be deduced that the ranking of the most favorable geometries goes like: V, IV, II, III, I.

 

Fig. 4 Singular values as a function of the singular mode order for the source-detection geometries presented in Table 1. The red dots represent (for each geometry) the minimum singular value (maximum mode order) that can be used to reconstruct an image when 1% stochastic noise is added to the simulated data set (see Table 3 for numerical values). The dotted line should be used as a visual guide indicating the approximate location of the noise floor threshold for all geometries.

Download Full Size | PDF

Table 3. Tabulation of the number of singular modes that can be used to reconstruct images in situations where different levels of noise and model mismatch are added to simulated tomography data sets. In cases where no noise is added, all singular modes (100% of rank value) can be used to reconstruct an image. In situations where noise is added there is a dramatic decrease in the number of useful singular modes.

View Table | View all tables in this article

The maximum number of singular modes that can in principle be used to reconstruct a fluorescence image is set by the rank of the matrix A, i.e., N = N _rank. However, in practice the number of modes that can actually be used to reconstruct an image is limited by, for example, the presence of high-frequency noise in the data vector D. The projection of stochastic noise contribution onto the data vector, D _high-f ⋅d_i^T, typically causes the coefficients K_i to contribute disproportionally to reconstructed images for some value i < N _rank, systematically leading to noisy reconstructed images. This effect might be perceived as a divergence in the singular coefficients caused by two competing effects: namely, monotonic decrease in singular values with increasing order and projection of high-frequency components of the data vector, e.g. stochastic noise or volume discretization noise, onto the SD vectors, which is disproportionately amplified for some values of i. For a better conditioned problem, sequentially smaller singular values (for i<N_rank) will not cause the K_i’s to diverge since the monotonic increase of 1/s_i is opposed by a continuously decreasing overlap of the SD and raw data vectors as i increases. However, in the presence of high-frequency noise, the large i behavior of the term D _high-f ⋅d_i^T introduces divergences, which hide the contribution of the lower frequency components, i.e., the modes containing the physically relevant biological information. In a nutshell, regularization in optical tomography amounts to minimizing the contribution of those singular vectors contributing disproportionally to reconstructed images due to the propagation of high-frequency noise. When using singular value decomposition to solve the inverse problem, regularization can be achieved by appropriately choosing the filter function F _i, or, by introducing a smoothing norm by using a generalized singular value decomposition (GSVD) [8].

Whereas singular value analysis can provide insight into the number of singular modes achievable for a given imaging geometry, it gives little insight into the shape of those modes, and therefore does not provide a full and fair analysis of the imaging geometry in question. In an attempt to utilize singular value decomposition to its full potential, the spatial frequency of the SI and SD vectors, novel metrics for analyzing imaging geometries, was investigated in the context of quantifying mode resolution characteristics for each of the described imaging geometries.

4.3 Singular vector analyses of image geometries

Figure 5 displays the mode resolution of the SI vectors (Eq. (6)) as a function of mode order for the five representative geometries described in Table 1. According to this figure of merit, favorable geometries will be those for which the modes used to reconstruct an image have the highest possible spatial frequency (smallest possible mode resolution). In other words, geometries with the steepest mode resolution versus mode order curves will be deemed favorable according to the SI mode resolution metric. Based on these criteria, the order of the most favorable metrics is (I, III, V, II, IV). Another criterion is the detection geometry-specific, minimum spatial resolution attainable (red dots in Fig. 5). The corresponding maximum number of modes used for different reconstruction scenarios are tabulated in Table 3. The values used in Fig. 5 correspond to the three columns associated with a 1% stochastic noise level. Using this last figure of merit and assuming equivalent levels of noise, geometry V (wide-field) seems to be favorable, although only by a small margin compared to the focused-based detection geometries. Perhaps a surprising observation derived from Fig. 5 is that the geometry with the smallest measurement density has resolution features comparable with geometries associated with a far larger number of measurements.

Figure 6b displays the frequency of SD vectors (Eq. (7)) as a function of mode order for individual imaging geometries. This metric was observed to generally increase monotonically with mode order, but was subject to interspersed dramatic decreases in resolution for some critical singular mode order values. This implies that low-frequency model mismatch errors, in addition to high-frequency stochastic noise errors, have to be taken into account at higher order modes. For a monotonically increasing SD mode frequency, it is expected that as the mode order increases, projection onto the high-frequency components in Eq. (8) will gradually make way to the projections onto the lower-frequency components. However, for geometries associated with interspersed, highly variable frequency-mode dependencies, this simple relationship between the contributions of high- and lower-frequency components is not expected to hold. An interesting observation here is that the truncation of useful modes for image reconstructions at 1% noise (red dots in Fig. 6b) always seems to occur shortly after (in terms of singular mode order) the critical points where low frequency noise starts to contribute to high order projection coefficients.

4.4 Image quality

Figure 7 displays reconstructed images for each geometry data set with increasing levels of stochastic noise. Inspection of this figure shows that even minute levels of noise degrade the image quality considerably in all cases. The imaging geometry emulating wide-field detection is able to maintain a level of spatial resolution that is clearly superior to the other four geometries for 0.1% levels of noise or less; however, for a more realistic level of noise (1%), geometry V remains only marginally superior to the focused-based geometries. One possible explanation for this phenomenon is the non-uniform loss of usable modes as stochastic noise increases. Employing the number of usable modes defined in previous sections, it can be demonstrated that for geometry I, where the rank of the forward matrix is the smallest, the percentage of modes (evaluated with respect to the rank of the matrix A) that can be retained for image reconstruction is 73% for 1% noise and 50% for 5% noise. In comparison, the wide-field detection geometry V allows only about 15% of its modes to be included for 1% and 5% noise. As for the other three geometries, they allow retention of between around 26-50% of the modes for 1% noise. In summary, more measurements do translate into more modes included in the images but in relative terms, high measurement density configurations may not be the most economical in terms of the percentage of useful modes.

 

Fig. 7 Images reconstructed from data simulated for different imaging geometries. Even minute levels of stochastic noise (0.001% - 1%) are shown to significantly degrade image quality for all detection geometries. The simulation target image used to generate the synthetic data is that shown in Fig. 3c. Images are scaled with respect to the maximum value of the target image.

Download Full Size | PDF

Figure 8 illustrates the affects of model mismatch and noise on image reconstruction for two different image geometries (II and V) using a forward model assuming homogeneous absorption and scattering properties when data was simulated for a heterogeneous environment (see Table 2). As shown in Fig. 3, this model mismatch is quite significant and as a result the quality of the reconstructed images (third column in Fig. 8) is significantly reduced for both imaging geometries. An important point to note here is that the optimal number of modes used to reconstruct these images is the same as that for the final column in Fig. 7 where the only source of noise was stochastic (see Table 3). This is a confirmation of the findings in earlier sections that model mismatch noise tends to contribute significantly only as low frequency components. The fourth column displays the corresponding images when Born normalization is employed. As illustrated in Fig. 3, using fluorescence-to-transmission data helps minimize the impact of forward model mismatch manifesting in a more balanced contribution of the two fluorescent inclusions in the presented reconstructed images.

 

Fig. 8 Images illustrating the impact of low and high frequency noise on fluorescence images reconstructed using a model assuming homogeneous absorption and scattering values for two different imaging geometries (Geometries II and V). The first column shows images where the forward model and reconstruction models were the same and no noise has been added to the synthetic data. The images found in the next three columns illustrate how different sources of noise degrade image quality: 1% stochastic noise added to the data vector (2nd column), 1% noise added plus low-frequency noise due to absorption heterogeneities in Table 2 (3rd column), same as column 3 but data vector used to reconstruct consists in fluorescence-to-transmission (F/T) ratio (column 4).

Download Full Size | PDF

Finally, Table 4 shows a summary of the image geometry performance predictions associated with both the singular value curve and the SI mode resolution curve. The ranking that is provided in the last column of the table orders the geometries starting with the most favorable geometry as assessed with either singular values or SI vectors. The last line of the table provides a ranking where image quality was assessed by evaluating the contrast recovery and the spatial resolution of the images presented in Fig. 7. The important thing to note here is that the ranking associated with either of the singular value analysis metrics (singular values or SI vectors) alone does not entirely agree with the a posteriori evaluation of the images based on image quality. For example, based on the singular value curve, geometry II is preferable to geometry III. However, based on the spatial resolution of the modes, geometry III is expected to perform better than geometry II. Similarly, based on singular values, the high density of measurements geometry IV is expected to be superior to geometry III, which is in disagreement with the prediction derived from the spatial resolution of the SI modes. Again, the prediction from the later criteria is more in tune with the ranking based on image quality.

Table 4. Figures of merit used to assess a priori and a posteriori performance of the imaging geometries presented in Table 1. For each set of figures of merit, a ranking is provided for geometries I-V starting with the one for which a figure of merit predicts the optimal performance.

View Table | View all tables in this article

4.5 Summary and conclusions

In this work, singular value decomposition analysis of image reconstruction matrices for fluorescence tomography was investigated in terms of its potential to predict the performance of various imaging geometries (corresponding to focused- and wide-field detection), quantified by reconstructed image quality. In the past, analysis of the number of usable singular values has been used to pre-evaluate the performance of proposed imaging geometries, generally predicting that, for larger fields of view, more source-detector projections are always preferable [7]. A singular value analysis was repeated in the present study for completeness, corroborating this finding; however, an analysis of the image quality improvements corresponding to increased source-detector projections demonstrated only very marginal improvements in the presence of realistic levels of noise and model-data mismatch. In fact, inspection of Fig. 7 provides support to the conservative claim that focused-based detection can achieve image quality close to what can be achieved with wide-field detection. To further support this claim, Fig. 9 presents simulated imaging scenarios illustrating how different imaging geometries perform when it comes to monitoring the evolution of a hypothetical fluorescent tumor. Three different evolution scenarios are presented where different tumor parameters are varied at four different stages: in Scenario A the contrast-to-background ratio remains constant (10:1) but the size of the tumor varies, in Scenario B the size of the inclusion is constant but the contrast varies, and, in Scenario C both the size and the contrast of the tumor vary. The tomography results for each scenario are presented in Fig. 9 where images were reconstructed using an homogenous matrix A but with fluorescence-to-transmission data generated with the heterogeneous optical properties found in Table 2. For all the reconstructions, 1% stochastic noise has been added to the data vector. The results presented in Fig. 9 show that under the same noise conditions, a wide-field detection geometry provides only minor imaging gains when it comes to spatial resolution and recovered contrast.

 

Fig. 9 Different simulation scenarios where FT is used to monitor the progress of a tumor with center-of-mass the same as top inclusion shown on the fluorescence image in Fig. 3c. Reconstruction results are presented as a cross-section of the corresponding tomography image through a horizontal line across the center-of-mass of the tumor. Scenario A shows a tumor of constant contrast-to-background ratio (10:1) but with decreasing size, Scenario B shows a tumors of decreasing contrast but with constant size, Scenario C shows tumor stages where both size and contrast vary. For each scenario, four tumor stages are shown and each image is reconstructed with three different imaging configurations, namely geometries I, II, V (Table 1).

Download Full Size | PDF

It should be mentioned here that the noise characteristics of wide-field detection, which are usually based on cooled CCD chips [11,12] are expected to be different than focused-detection geometries, typically based on PMTs, when it comes to the detection of low photon fluxes. For high fluorescence signals, the noise is expected in both cases to be shot-noise limited and to some extent the levels of noise could be equivalent. However, when it comes to the detection of very low photon fluxes, non-contact CCD-based systems will be limited by the read-out noise of the chip. This is particularly important given the limited dynamical range associated with non-contact CCD-based small animal imaging. In terms of noise characteristics for low photon fluence detection, the optimal systems are time-domain PMT systems using time correlated single photon counting [3]. In this case, detection is achieved in such a way that the signal-to-noise levels achievable for low photon fluxes encountered in fluorescence tomography applications can be orders of magnitude larger than what can be achieved with cooled CCD’s. In conclusion, when comparing wide-field detection with focused-based detection it might be more realistic to compare wide-field reconstructions at 1% noise with focused-based detection with 0.1% stochastic noise added.

Even assuming equal noise considerations for all imaging geometries, the geometry evaluation ranking of the singular values did not match the image quality rankings (Table 4). These findings suggest that the number of usable singular values alone cannot fully characterize the performance of a given imaging geometry for realistic data sets. In response, we developed new metrics, mode resolution and mode frequency, to quantify the characteristics of the singular image and singular data vectors, respectively, corresponding to each singular value. Like the singular value analysis, alone, these singular vector metrics did not adequately reflect the image quality ranking; however, their ranking could be used to account for the discrepancies between the singular value ranking and the image quality ranking, suggesting that together, singular value and singular vector metrics are complementary and could conceivably be used to fully characterize the efficacy of proposed imaging geometries. Future research efforts will be devoted to findings methodologies to better characterize the relative importance and interplay between the different singular value decomposition metrics.

4.6 Limitations and future work

Singular value analysis has traditionally been used in optical tomography as an analysis tool for imaging geometry optimization. This method is rarely used to actually reconstruct diffuse fluorescence images in the scope of biological experiments. One of the reasons for this is that singular value decomposition is a direct matrix inversion method that requires significant computer resources, which become particularly time intensive for reconstructing 3D volumes in whole body, small animal imaging or for larger volumes such as breast imaging. Though the current study focused on 2D reconstructions of a mouse head, the conclusion that increasing source-detector densities has diminishing gains in terms of quality of image reconstruction is expected to hold for 3D and larger volume reconstructions although further studies are warranted. Another reason why SVD is not commonly used for image reconstruction is that the approach is not flexible in terms of allowing for the inclusion of a generic penalty function when resolving the inverse problem. Although a conventional Tikhonov penalty term can be introduced using a specific filter function F _i in Eq. (4), an analysis of the problem based on a generalized SVD (GSVD) must be performed in order to allow penalty functions that would be consistent with, e.g., soft as well hard spatial priors. While the results presented in this paper provide useful information on the FT problem, they should be interpreted within the confines of those limitations. In fact, it is more commonplace in whole-body fluorescence tomography to use inverse problem resolution methods with regularization functions allowing certain aspects of the images to be favored as determined by prior knowledge (e.g., spectral and/or spatial) or pre-determined features that are expected for the images (e.g., edge preservation, absence of high spatial frequencies).

Finally, the analysis presented in this work was concerned with how different imaging configurations can affect FT results assuming that the data were from continuous-wave measurement systems. A conclusion that is reached here is that the relative gains, in terms of spatial resolution and recovered contrast, associated with using a high density of measurements, is not as impressive as one might expect based upon linear system performance. In fact, the results presented in Fig. 7 show only modest image quality gains associated with a system where the number of measurements was increased up to 10 - 50 times. Therefore, future developments in FT should be concerned with investigating imaging schemes that would allow more significant improvements to be achieved within the confine of the intrinsic limitations posed by light scattering. One way of achieving this would be to use the singular value analysis methods presented in this work to investigate how different non-redundant data types could decrease the impact of the ill-posedness of the problem. Approaches using tomography data sets derived from multi-spectral or time-domain systems constitute promising research directions.